diff --git a/candidateset.cpp b/candidateset.cpp
index f0319738..589e97b7 100644
--- a/candidateset.cpp
+++ b/candidateset.cpp
@@ -150,7 +150,7 @@ bool CandidateSet::update(string tree, double score) {
 }
 
 vector<double> CandidateSet::getBestScores(int numBestScore) {
-	if (numBestScore == 0)
+	if (numBestScore == 0 || numBestScore > size())
 		numBestScore = size();
 	vector<double> res;
 	for (reverse_iterator rit = rbegin(); rit != rend() && numBestScore > 0; rit++, numBestScore--) {
diff --git a/iqtree.cpp b/iqtree.cpp
index 405411e6..1a32ef90 100644
--- a/iqtree.cpp
+++ b/iqtree.cpp
@@ -1539,7 +1539,7 @@ string IQTree::optimizeModelParameters(bool printInfo) {
 }
 
 void IQTree::printBestScores(int numBestScore) {
-	vector<double> bestScores = candidateTrees.getBestScores(candidateTrees.popSize);
+	vector<double> bestScores = candidateTrees.getBestScores(numBestScore);
 	for (vector<double>::iterator it = bestScores.begin(); it != bestScores.end(); it++)
 		cout << (params->maximum_parsimony ? -(*it) : (*it)) << " ";
 	cout << endl;
diff --git a/phyloanalysis.cpp b/phyloanalysis.cpp
index d96dfc0f..f7dd7bf2 100644
--- a/phyloanalysis.cpp
+++ b/phyloanalysis.cpp
@@ -1851,8 +1851,13 @@ void runTreeReconstruction(Params &params, string &original_model, IQTree &iqtre
 	if (iqtree.isSuperTree())
 			((PhyloSuperTree*) &iqtree)->mapTrees();
 	if (params.snni && params.min_iterations) {
-		cout << (params.maximum_parsimony ? "Scores" : "Log-likelihoods") << " of best " << params.popSize << " trees: " << endl;
-		iqtree.printBestScores(iqtree.candidateTrees.popSize);
+		if (params.savek == 0) {
+			cout << (params.maximum_parsimony ? "Scores" : "Log-likelihoods") << " of best " << params.popSize << " trees: " << endl;
+			iqtree.printBestScores(iqtree.candidateTrees.popSize);
+		} else {
+			cout << (params.maximum_parsimony ? "Scores" : "Log-likelihoods") << " of best " << params.savek << " trees: " << endl;
+			iqtree.printBestScores(params.savek);
+		}
 	}
 
 	/******** Performs final model parameters optimization ******************/
diff --git a/tools.cpp b/tools.cpp
index d58577a6..859d49c3 100644
--- a/tools.cpp
+++ b/tools.cpp
@@ -805,6 +805,7 @@ void parseArg(int argc, char *argv[], Params &params) {
     params.do_first_rell = false;
     params.remove_dup_seq = false;
     params.test_mode = false;
+	params.savek = 0;
 
 #ifdef _OPENMP
     params.num_threads = 0;
@@ -2263,6 +2264,15 @@ void parseArg(int argc, char *argv[], Params &params) {
 				assert(params.popSize < params.numParsTrees);
 				continue;
 			}
+			if (strcmp(argv[cnt], "-savek") == 0) {
+				cnt++;
+				if (cnt >= argc)
+					throw "Use -savek <number_of_saving_trees>";
+				params.savek = convert_int(argv[cnt]);
+				params.numParsTrees = max(params.savek, params.numParsTrees);
+				assert(params.savek > 0);
+				continue;
+			}
 			if (strcmp(argv[cnt], "-beststart") == 0) {
 				params.bestStart = true;
 				cnt++;
diff --git a/tools.h b/tools.h
index 315db2dd..92875d03 100644
--- a/tools.h
+++ b/tools.h
@@ -1625,6 +1625,11 @@ struct Params {
 	 */
 	bool test_mode;
 
+        /**
+         * number of search trees saved
+         */
+        int savek;
+
 	/*
 	 * Diep:
 	 * Use with -test_mode