diff --git a/src/DFTsolver.cc b/src/DFTsolver.cc
index 034366fd..5e512090 100644
--- a/src/DFTsolver.cc
+++ b/src/DFTsolver.cc
@@ -316,7 +316,7 @@ int DFTsolver<OrbitalsType>::solve(OrbitalsType& orbitals,
         const bool ortho
             = (ct.getOrthoType() == OrthoType::Eigenfunctions || orthof);
 
-        if (!ortho)
+        if (!ortho || !ct.fullyOccupied())
         {
             // strip dm from the overlap contribution
             // dm <- Ls**T * dm * Ls
@@ -337,7 +337,14 @@ int DFTsolver<OrbitalsType>::solve(OrbitalsType& orbitals,
             }
             else
             {
-                orbitals.orthonormalizeLoewdin();
+                bool updateDM = false;
+                if (!ct.fullyOccupied())
+                {
+                    orbitals.computeGramAndInvS();
+                    dm_strategy_->dressDM();
+                    updateDM = true;
+                }
+                orbitals.orthonormalizeLoewdin(true, nullptr, updateDM);
 
                 orbitals_stepper_->restartMixing();
             }
diff --git a/src/ProjectedMatrices.cc b/src/ProjectedMatrices.cc
index d3707657..6d6b6339 100644
--- a/src/ProjectedMatrices.cc
+++ b/src/ProjectedMatrices.cc
@@ -448,12 +448,14 @@ void ProjectedMatrices<MatrixType>::updateDM(const int iterative_index)
 
 #ifndef NDEBUG
     double nel = getNel();
-    std::cout << "ProjectedMatrices<MatrixType>::updateDM(), nel = " << nel
-              << std::endl;
+    if (mmpi.instancePE0())
+        std::cout << "ProjectedMatrices<MatrixType>::updateDM(), nel = " << nel
+                  << std::endl;
     assert(std::isfinite(nel));
     double energy = getExpectationH();
-    std::cout << "ProjectedMatrices<MatrixType>::updateDM(), energy = "
-              << energy << std::endl;
+    if (mmpi.instancePE0())
+        std::cout << "ProjectedMatrices<MatrixType>::updateDM(), energy = "
+                  << energy << std::endl;
 #endif
 }
 
@@ -935,6 +937,11 @@ void ProjectedMatrices<MatrixType>::computeLoewdinTransform(
 
     if (transform_matrices)
     {
+        Control& ct     = *(Control::instance());
+        MGmol_MPI& mmpi = *(MGmol_MPI::instance());
+        if (mmpi.instancePE0() && ct.verbose > 1)
+            std::cout << "Transform DM to reflect Loewdin orthonormalization"
+                      << std::endl;
         assert(sqrtMat);
 
         // transform DM to reflect Loewdin orthonormalization
diff --git a/tests/MVP/mvp.cfg b/tests/MVP/mvp.cfg
index e34d3513..25d51ed2 100644
--- a/tests/MVP/mvp.cfg
+++ b/tests/MVP/mvp.cfg
@@ -20,6 +20,7 @@ type=QUENCH
 solver=PSD
 max_steps=300
 atol=1.e-7
+ortho_freq=10
 [Orbitals]
 nempty=10
 initial_type=random