llnl · siuwuncheung · May 15, 2025 · Sep 26, 2024 · Sep 30, 2024 · Sep 30, 2024
diff --git a/examples/PinnedH2O/job.basis_1_50 b/examples/PinnedH2O/job.basis_1_50
@@ -10,7 +10,7 @@ setenv OMP_NUM_THREADS 1
 
 set ncpus = 8
 
-set maindir = /p/lustre2/cheung26/mgmol-20250219
+set maindir = /p/lustre2/cheung26/mgmol
 
 setenv LD_LIBRARY_PATH ${maindir}/build_quartz/libROM/build/lib:$LD_LIBRARY_PATH
 

diff --git a/examples/PinnedH2O/job.offline b/examples/PinnedH2O/job.offline
@@ -10,7 +10,7 @@ setenv OMP_NUM_THREADS 1
 
 set ncpus = 8
 
-set maindir = /p/lustre2/cheung26/mgmol-20250219
+set maindir = /p/lustre2/cheung26/mgmol
 
 setenv LD_LIBRARY_PATH ${maindir}/build_quartz/libROM/build/lib:$LD_LIBRARY_PATH
 

diff --git a/examples/PinnedH2O/job.ref b/examples/PinnedH2O/job.ref
@@ -11,7 +11,7 @@ setenv OMP_NUM_THREADS 1
 set ncpus = 8
 set case = 2
 
-set maindir = /p/lustre2/cheung26/mgmol-20250219
+set maindir = /p/lustre2/cheung26/mgmol
 
 setenv LD_LIBRARY_PATH ${maindir}/build_quartz/libROM/build/lib:$LD_LIBRARY_PATH
 

diff --git a/examples/PinnedH2O/job.rom b/examples/PinnedH2O/job.rom
@@ -10,7 +10,7 @@ setenv OMP_NUM_THREADS 1
 
 set ncpus = 8
 
-set maindir = /p/lustre2/cheung26/mgmol-20250219
+set maindir = /p/lustre2/cheung26/mgmol
 
 setenv LD_LIBRARY_PATH ${maindir}/build_quartz/libROM/build/lib:$LD_LIBRARY_PATH
 

diff --git a/examples/PinnedH2O/job.rom_3DOF b/examples/PinnedH2O/job.rom_3DOF
@@ -11,7 +11,7 @@ setenv OMP_NUM_THREADS 1
 set ncpus = 8
 set case = 2
 
-set maindir = /p/lustre2/cheung26/mgmol-20250219
+set maindir = /p/lustre2/cheung26/mgmol
 
 setenv LD_LIBRARY_PATH ${maindir}/build_quartz/libROM/build/lib:$LD_LIBRARY_PATH
 

diff --git a/examples/PinnedH2O/mgmol_rom_3DOF_test1.cfg b/examples/PinnedH2O/mgmol_rom_3DOF_test1.cfg
@@ -1,4 +1,4 @@
-verbosity=1
+verbosity=2
 xcFunctional=PBE
 FDtype=4th
 [Mesh]
@@ -34,7 +34,9 @@ initial_width=2.
 nempty=30
 [DensityMatrix]
 solver=MVP
-nb_inner_it=100
+nb_inner_it=20
+mixing=0.8
+tol=1.e-8
 [Restart]
 output_level=4
 [ROM]

diff --git a/examples/PinnedH2O/mgmol_rom_3DOF_test2.cfg b/examples/PinnedH2O/mgmol_rom_3DOF_test2.cfg
@@ -1,4 +1,4 @@
-verbosity=1
+verbosity=2
 xcFunctional=PBE
 FDtype=4th
 [Mesh]
@@ -30,7 +30,9 @@ initial_width=2.
 nempty=30
 [DensityMatrix]
 solver=MVP
-nb_inner_it=100
+nb_inner_it=20
+mixing=0.8
+tol=1.e-8
 [Restart]
 output_level=4
 [ROM]

diff --git a/src/Control.cc b/src/Control.cc
@@ -224,6 +224,7 @@ void Control::print(std::ostream& os)
        << conv_tol << std::endl;
     os << std::fixed;
     os << " Density matrix mixing = " << dm_mix << std::endl;
+    os << std::setprecision(4) << std::scientific << " Density matrix tol = " << dm_tol << std::endl;
     if (DMEigensolver() == DMEigensolverType::Eigensolver)
     {
         os << " Density matrix computation algorithm = "
@@ -439,7 +440,7 @@ void Control::sync(void)
         memset(&int_buffer[0], 0, size_int_buffer * sizeof(int));
     }
 
-    const short size_float_buffer = 43;
+    const short size_float_buffer = 44;
     float* float_buffer           = new float[size_float_buffer];
     if (mype_ == 0)
     {
@@ -485,6 +486,7 @@ void Control::sync(void)
         float_buffer[40] = threshold_eigenvalue_gram_quench_;
         float_buffer[41] = pair_mlwf_distance_threshold_;
         float_buffer[42] = e0_;
+        float_buffer[43] = dm_tol;
     }
     else
     {
@@ -680,6 +682,7 @@ void Control::sync(void)
     threshold_eigenvalue_gram_quench_ = float_buffer[40];
     pair_mlwf_distance_threshold_     = float_buffer[41];
     e0_                               = float_buffer[42];
+    dm_tol                            = float_buffer[43];
     max_electronic_steps_loose_       = max_electronic_steps;
 
     delete[] short_buffer;
@@ -699,8 +702,8 @@ void Control::setDefaultValues()
 
 void Control::adjust()
 {
-    // change dm_mix default to 1. if not using Davidson
-    if (it_algo_type_ != 2 && dm_mix < 0.) dm_mix = 1.;
+    // change dm_mix default to 1. if not using Davidson or MVP
+    if ((it_algo_type_ != 2 && DM_solver_ != 1) && dm_mix < 0.) dm_mix = 1.;
 
     if (nel_ - 2 * numst == 0)
     {
@@ -1720,6 +1723,7 @@ void Control::setOptions(const boost::program_options::variables_map& vm)
             lrs_extrapolation = 10;
 
         dm_mix         = vm["DensityMatrix.mixing"].as<float>();
+        dm_tol         = vm["DensityMatrix.tol"].as<float>();
         dm_inner_steps = vm["DensityMatrix.nb_inner_it"].as<short>();
         dm_use_old_    = vm["DensityMatrix.use_old"].as<bool>() ? 1 : 0;
         str            = vm["DensityMatrix.algo"].as<std::string>();
@@ -1739,8 +1743,6 @@ void Control::setOptions(const boost::program_options::variables_map& vm)
         else
             dm_algo_ = 2;
 
-        dm_tol = vm["DensityMatrix.tol"].as<float>();
-
         str = vm["DensityMatrix.solver"].as<std::string>();
         if (str.compare("Mixing") == 0) DM_solver_ = 0;
         if (str.compare("MVP") == 0) DM_solver_ = 1;

diff --git a/src/Ions.cc b/src/Ions.cc
@@ -1315,22 +1315,6 @@ void Ions::setLocalForces(
         {
             if (ion->compareName(*s))
             {
-                //std::cout << "Ion found: " << ion->name() << std::endl;
-                //std::cout << "Ion force: (" << *it << ", " << *(it + 1) << ", " << *(it + 2) << ")" << std::endl;
-                //std::cout << "names: ";
-                //for (int i = 0; i < names.size(); i++)
-                //{
-                //    std::cout << names[i];
-                //    if (i == forces.size() - 1) std::cout << std::endl;
-                //    else std::cout << ", ";
-                //}
-                //std::cout << "forces: ";
-                //for (int i = 0; i < forces.size(); i++)
-                //{
-                //    std::cout << forces[i];
-                //    if (i == forces.size() - 1) std::cout << ")" << std::endl;
-                //    else std::cout << ", ";
-                //}
                 ion->set_force(0, *it);
                 ion->set_force(1, *(it + 1));
                 ion->set_force(2, *(it + 2));
@@ -1364,9 +1348,6 @@ void Ions::setLocalForces(
             double d = std::sqrt(d2);
             if (d < tol)
             {
-                std::cout << "Ion found: " << ion->name() << std::endl;
-                std::cout << "Ion position:( " << *cit << ", " << *(cit + 1) << ", " << *(cit + 2) << ")" << std::endl;
-                std::cout << "Ion force: (" << *fit << ", " << *(fit + 1) << ", " << *(fit + 2) << ")" << std::endl;
                 ion->set_force(0, *fit);
                 ion->set_force(1, *(fit + 1));
                 ion->set_force(2, *(fit + 2));

diff --git a/src/MGmol.cc b/src/MGmol.cc
@@ -102,6 +102,7 @@ extern Timer md_moveVnuc_tm;
 extern Timer md_updateMasks_tm;
 extern Timer md_extrapolateOrbitals_tm;
 extern Timer md_updateRhoAndPot_tm;
+extern Timer md_updateDMandEnergy_tm;
 extern Timer quench_tm;
 extern Timer ions_setupInteractingIons_tm;
 extern Timer ions_setup_tm;
@@ -928,6 +929,7 @@ void MGmol<OrbitalsType>::printTimers()
     init_nuc_tm_.print(os_);
     md_updateMasks_tm.print(os_);
     md_extrapolateOrbitals_tm.print(os_);
+    md_updateDMandEnergy_tm.print(os_);
     quench_tm.print(os_);
     evnl_tm_.print(os_);
     ions_setupInteractingIons_tm.print(os_);

diff --git a/src/MVPSolver.cc b/src/MVPSolver.cc
@@ -47,13 +47,16 @@ MVPSolver<OrbitalsType, MatrixType>::MVPSolver(MPI_Comm comm, std::ostream& os,
     Electrostatic* electrostat, MGmol<OrbitalsType>* mgmol_strategy,
     const int numst, const double kbT,
     const std::vector<std::vector<int>>& global_indexes,
-    const short n_inner_steps, const bool use_old_dm)
+    const short n_inner_steps, const double mixing, const double tol_de0,
+    const bool use_old_dm)
     : comm_(comm),
       os_(os),
       n_inner_steps_(n_inner_steps),
       use_old_dm_(use_old_dm),
       ions_(ions),
-      numst_(numst)
+      numst_(numst),
+      mixing_(mixing),
+      tol_de0_(tol_de0)
 {
     Control& ct = *(Control::instance());
     if (onpe0 && ct.verbose > 0)
@@ -207,7 +210,6 @@ int MVPSolver<OrbitalsType, MatrixType>::solve(OrbitalsType& orbitals)
 
         kbpsi.computeHvnlMatrix(&kbpsi, ions_, h11_nl);
 
-        const double tol_de0 = 1.e-12;
         for (int inner_it = 0; inner_it < n_inner_steps_; inner_it++)
         {
             if (onpe0 && ct.verbose > 1)
@@ -270,62 +272,79 @@ int MVPSolver<OrbitalsType, MatrixType>::solve(OrbitalsType& orbitals)
 
                 double de0 = evaluateDerivative(dmInit, delta_dm, ts0);
 
-                if (std::abs(de0) < tol_de0 && inner_it > 0)
+                // check for convergence
+                if (std::abs(de0) < tol_de0_ && inner_it > 0)
                 {
                     if (onpe0 && ct.verbose > 0)
                         std::cout << "MVP: de0 = " << de0
                                   << ", convergence achieved" << std::endl;
                     break;
                 }
 
-                //
-                // evaluate free energy at beta=1
-                //
-                if (onpe0 && ct.verbose > 2)
-                    std::cout << "MVP --- Target energy..." << std::endl;
-                proj_mat_work_->setDM(target, orbitals.getIterativeIndex());
-                proj_mat_work_->computeOccupationsFromDM();
-                if (ct.verbose > 2) proj_mat_work_->printOccupations(os_);
-                const double nel = proj_mat_work_->getNel();
-                if (onpe0 && ct.verbose > 1)
-                    os_ << "MVP --- Number of electrons at beta=1 : " << nel
-                        << std::endl;
-
-                rho_->computeRho(orbitals, target);
-
-                mgmol_strategy_->update_pot(vh_init, ions_);
-
-                energy_->saveVofRho();
-
-                // update h11
+                double beta = 0.;
+                if (mixing_ > 0.)
                 {
-                    h11 = h11_nl;
-                    mgmol_strategy_->addHlocal2matrix(orbitals, orbitals, h11);
+                    beta = mixing_;
+                    if (onpe0 && ct.verbose > 0)
+                    {
+                        if (ct.verbose > 1)
+                            os_ << "MVP with beta = " << beta << std::endl;
+                        os_ << std::setprecision(12);
+                        os_ << std::fixed << "MVP inner iteration " << inner_it
+                            << ", E0=" << e0 << std::endl;
+                    }
                 }
-
-                proj_mat_work_->assignH(h11);
-                proj_mat_work_->setHB2H();
-
-                const double ts1
-                    = evalEntropyMVP(proj_mat_work_, (ct.verbose > 2), os_);
-                const double e1 = energy_->evaluateTotal(
-                    ts1, proj_mat_work_, ions_, orbitals, ct.verbose - 1, os_);
-
-                // line minimization
-                const double beta
-                    = minQuadPolynomial(e0, e1, de0, (ct.verbose > 2), os_);
-                assert(!std::isnan(beta));
-
-                if (onpe0 && ct.verbose > 0)
+                else
                 {
-                    os_ << std::setprecision(12);
-                    os_ << std::fixed << "MVP inner iteration " << inner_it
-                        << ", E0=" << e0 << ", E1=" << e1;
-                    os_ << std::scientific << ", E0'=" << de0
-                        << " -> beta=" << beta;
-                    os_ << std::endl;
-                }
+                    //
+                    // evaluate free energy at beta=1
+                    //
+                    if (onpe0 && ct.verbose > 2)
+                        std::cout << "MVP --- Target energy..." << std::endl;
+                    proj_mat_work_->setDM(target, orbitals.getIterativeIndex());
+                    proj_mat_work_->computeOccupationsFromDM();
+                    if (ct.verbose > 2) proj_mat_work_->printOccupations(os_);
+                    const double nel = proj_mat_work_->getNel();
+                    if (onpe0 && ct.verbose > 1)
+                        os_ << "MVP --- Number of electrons at beta=1 : " << nel
+                            << std::endl;
+
+                    rho_->computeRho(orbitals, target);
+
+                    mgmol_strategy_->update_pot(vh_init, ions_);
+
+                    energy_->saveVofRho();
+
+                    // update h11
+                    {
+                        h11 = h11_nl;
+                        mgmol_strategy_->addHlocal2matrix(
+                            orbitals, orbitals, h11);
+                    }
+
+                    proj_mat_work_->assignH(h11);
+                    proj_mat_work_->setHB2H();
+
+                    const double ts1
+                        = evalEntropyMVP(proj_mat_work_, (ct.verbose > 2), os_);
+                    const double e1 = energy_->evaluateTotal(ts1,
+                        proj_mat_work_, ions_, orbitals, ct.verbose - 1, os_);
+
+                    // line minimization
+                    beta
+                        = minQuadPolynomial(e0, e1, de0, (ct.verbose > 2), os_);
+                    assert(!std::isnan(beta));
 
+                    if (onpe0 && ct.verbose > 0)
+                    {
+                        os_ << std::setprecision(12);
+                        os_ << std::fixed << "MVP inner iteration " << inner_it
+                            << ", E0=" << e0 << ", E1=" << e1;
+                        os_ << std::scientific << ", E0'=" << de0
+                            << " -> beta=" << beta;
+                        os_ << std::endl;
+                    }
+                }
                 // update DM
                 *work_ = dmInit;
                 work_->axpy(beta, delta_dm);

diff --git a/src/MVPSolver.h b/src/MVPSolver.h
@@ -33,6 +33,12 @@ class MVPSolver
     Ions& ions_;
 
     int numst_;
+    double mixing_;
+
+    /*!
+     * tolerance on energy slope in inner iterations
+     */
+    double tol_de0_;
 
     Rho<OrbitalsType>* rho_;
     Energy<OrbitalsType>* energy_;
@@ -56,7 +62,8 @@ class MVPSolver
         Electrostatic* electrostat, MGmol<OrbitalsType>* mgmol_strategy,
         const int numst, const double kbT,
         const std::vector<std::vector<int>>& global_indexes,
-        const short n_inner_steps, const bool use_old_dm);
+        const short n_inner_steps, const double mixing, const double tol_de0,
+        const bool use_old_dm);
     ~MVPSolver();
 
     int solve(OrbitalsType& orbitals);

diff --git a/src/MVP_DMStrategy.cc b/src/MVP_DMStrategy.cc
@@ -54,7 +54,7 @@ int MVP_DMStrategy<OrbitalsType, MatrixType>::update(OrbitalsType& orbitals)
 
     MVPSolver<OrbitalsType, MatrixType> solver(comm_, os_, ions_, rho_, energy_,
         electrostat_, mgmol_strategy_, ct.numst, ct.occ_width, global_indexes_,
-        ct.dm_inner_steps, use_old_dm_);
+        ct.dm_inner_steps, ct.dm_mix, ct.dm_tol, use_old_dm_);
 
     return solver.solve(orbitals);
 }

diff --git a/src/md.cc b/src/md.cc
@@ -50,6 +50,7 @@ Timer md_tau_tm("md_tau");
 Timer md_moveVnuc_tm("md_moveVnuc");
 Timer md_updateMasks_tm("md_updateMasks");
 Timer md_extrapolateOrbitals_tm("md_extrapolateOrbitals");
+Timer md_updateDMandEnergy_tm("md_updateDMandEnergy");
 
 #define DUMP_MAX_NUM_TRY 5
 
@@ -485,7 +486,9 @@ void MGmol<OrbitalsType>::md(OrbitalsType** orbitals, Ions& ions)
 
         if (ROM_MVP)
         {
+            md_updateDMandEnergy_tm.start();
             updateDMandEnergy(**orbitals, *ROM_ions, eks);
+            md_updateDMandEnergy_tm.stop();
         }
         else
         {