fix: fix a bug in Nose-Hoover Chain (#1470)

* refactor: rename md_mnhc to md_tchain * refactor: update the default value of md_tchain * refactor: rename NVT_NHC to Nose_Hoover * fix: fix a bug in Nose-Hoover Chain * refactor: rename some para and functions in md * test: update result.ref due to bug fix * test: update MD_func_test due to refactor * test: update result.ref of 183_PW_MD_SDFT_10S
deepmodeling · Nov 4, 2022 · 1905280 · 1905280
1 parent a326912
commit 1905280
Show file tree

Hide file tree

Showing 42 changed files with 817 additions and 469 deletions.
diff --git a/docs/advanced/input_files/input-main.md b/docs/advanced/input_files/input-main.md
@@ -33,7 +33,7 @@
 
 - [Molecular dynamics](#molecular-dynamics)
 
-  [md_type](#md_type) | [md_thermostat](#md_thermostat) | [md_nstep](#md_nstep) | [md_ensolver](#md_ensolver) | [md_restart](#md_restart) | [md_dt](#md_dt) | [md_tfirst, md_tlast](#md_tfirst-md_tlast) | [md_dumpfreq](#md_dumpfreq) | [md_restartfreq](#md_restartfreq) | [md_seed](#md_seed) | [md_tfreq](#md_tfreq) | [md_mnhc](#md_mnhc) | [lj_rcut](#lj_rcut) | [lj_epsilon](#lj_epsilon) | [lj_sigma](#lj_sigma) | [pot_file](#pot_file) | [msst_direction](#msst_direction) | [msst_vel](#msst_vel) | [msst_vis](#msst_vis) | [msst_tscale](#msst_tscale) | [msst_qmass](#msst_qmass) | [md_damp](#md_damp) | [md_tolerance](#md_tolerance) | [md_nraise](#md_nraise)
+  [md_type](#md_type) | [md_thermostat](#md_thermostat) | [md_nstep](#md_nstep) | [md_ensolver](#md_ensolver) | [md_restart](#md_restart) | [md_dt](#md_dt) | [md_tfirst, md_tlast](#md_tfirst-md_tlast) | [md_dumpfreq](#md_dumpfreq) | [md_restartfreq](#md_restartfreq) | [md_seed](#md_seed) | [md_tfreq](#md_tfreq) | [md_tchain](#md_tchain) | [lj_rcut](#lj_rcut) | [lj_epsilon](#lj_epsilon) | [lj_sigma](#lj_sigma) | [pot_file](#pot_file) | [msst_direction](#msst_direction) | [msst_vel](#msst_vel) | [msst_vis](#msst_vis) | [msst_tscale](#msst_tscale) | [msst_qmass](#msst_qmass) | [md_damp](#md_damp) | [md_tolerance](#md_tolerance) | [md_nraise](#md_nraise)
 
 - [vdW correction](#vdw-correction)
 
@@ -1472,11 +1472,11 @@ This part of variables are used to control the molecular dynamics calculations.
 - **Description**: control the frequency of the temperature oscillations during the simulation. If it is too large, the temperature will fluctuate violently; if it is too small, the temperature will take a very long time to equilibrate with the atomic system.
 - **Default**: 1/40/md_dt
 
-### md_mnhc
+### md_tchain
 
 - **Type**: Integer
 - **Description**: Number of Nose-Hoover chains.
-- **Default**: 4
+- **Default**: 1
 
 ### lj_rcut
 

diff --git a/examples/md/lcao_gammaonly_Sn64/INPUT_1 b/examples/md/lcao_gammaonly_Sn64/INPUT_1
@@ -33,4 +33,4 @@ md_nstep               10
 md_dt                  1
 md_tfirst              300
 md_tfreq               1
-md_mnhc                1
+md_tchain                1
diff --git a/source/Makefile.Objects b/source/Makefile.Objects
@@ -240,7 +240,7 @@ OBJS_MD=MD_func.o\
 	mdrun.o\
 	verlet.o\
     MSST.o\
-    NVT_NHC.o\
+    Nose_Hoover.o\
     FIRE.o\
     Langevin.o\
 

diff --git a/source/input.cpp b/source/input.cpp
@@ -1212,9 +1212,9 @@ bool Input::Read(const std::string &fn)
         {
             read_value(ifs, mdp.md_dt);
         }
-        else if (strcmp("md_mnhc", word) == 0)
+        else if (strcmp("md_tchain", word) == 0)
         {
-            read_value(ifs, mdp.md_mnhc);
+            read_value(ifs, mdp.md_tchain);
         }
         else if (strcmp("md_tfirst", word) == 0)
         {
@@ -2289,7 +2289,7 @@ void Input::Bcast()
     Parallel_Common::bcast_string(mdp.md_thermostat);
     Parallel_Common::bcast_int(mdp.md_nstep);
     Parallel_Common::bcast_double(mdp.md_dt);
-    Parallel_Common::bcast_int(mdp.md_mnhc);
+    Parallel_Common::bcast_int(mdp.md_tchain);
     Parallel_Common::bcast_double(mdp.msst_qmass);
     Parallel_Common::bcast_double(mdp.md_tfirst);
     Parallel_Common::bcast_double(mdp.md_tlast);
@@ -2645,17 +2645,24 @@ void Input::Check(void)
             out_level = "m"; // zhengdy add 2019-04-07
 
         // deal with input parameters , 2019-04-30
-        // if(basis_type == "pw" ) ModuleBase::WARNING_QUIT("Input::Check","calculate = MD is only availble for LCAO.");
         if (mdp.md_dt < 0)
             ModuleBase::WARNING_QUIT("Input::Check", "time interval of MD calculation should be set!");
         if (mdp.md_tfirst < 0 && tddft==0)
             ModuleBase::WARNING_QUIT("Input::Check", "temperature of MD calculation should be set!");
         if (mdp.md_tlast < 0.0)
             mdp.md_tlast = mdp.md_tfirst;
+        if(mdp.md_pmode != "none" && mdp.md_pfirst < 0)
+            ModuleBase::WARNING_QUIT("Input::Check", "pressure of MD calculation should be set!");
+        if (mdp.md_plast < 0.0)
+            mdp.md_plast = mdp.md_pfirst;
 
         if(mdp.md_tfreq == 0)
         {
-            mdp.md_tfreq = 1.0/40.0/mdp.md_dt;
+            mdp.md_tfreq = 1.0/40/mdp.md_dt;
+        }
+        if(mdp.md_pfreq == 0)
+        {
+            mdp.md_pfreq = 1.0/400/mdp.md_dt;
         }
         if(mdp.md_restart) 
         {
@@ -2679,23 +2686,6 @@ void Input::Check(void)
                 ModuleBase::WARNING_QUIT("Input::Check", "Can not find DP model !");
             }
         }
-        // if(mdp.md_tfirst!=mdp.md_tlast)
-        // {
-        //     std::ifstream file1;
-        //     file1.open("ChangeTemp.dat");
-        //     if(!file1)                      // Peize Lin fix bug 2016-08-06
-        //    {
-        //         std::ofstream file;
-        //         file.open("ChangeTemp.dat");
-        //         for(int ii=0;ii<30;ii++)
-        //         {
-        //             file<<mdp.md_tfirst+(mdp.md_tlast-mdp.md_tfirst)/double(30)*double(ii+1)<<" ";
-        //         }
-        //         file.close();
-        //     }
-        //     else
-        //         file1.close();
-        // }
     }
     else if (calculation == "cell-relax") // mohan add 2011-11-04
     {

diff --git a/source/module_md/CMakeLists.txt b/source/module_md/CMakeLists.txt
@@ -6,7 +6,7 @@ add_library(
 	mdrun.cpp
 	verlet.cpp
     MSST.cpp
-    NVT_NHC.cpp
+    Nose_Hoover.cpp
     FIRE.cpp
     Langevin.cpp
 )

diff --git a/source/module_md/MD_func.cpp b/source/module_md/MD_func.cpp
@@ -51,48 +51,29 @@ double MD_func::GetAtomKE(
 	return ke;
 }
 
-void MD_func::kinetic_stress(
+void MD_func::compute_stress(
 		const UnitCell_pseudo &unit_in,
 		const ModuleBase::Vector3<double> *vel, 
 		const double *allmass, 
-		double &kinetic,
+        const ModuleBase::matrix &virial,
 		ModuleBase::matrix &stress)
 {
-//---------------------------------------------------------------------------
+//--------------------------------------------------------------------------------------------
 // DESCRIPTION:
-//   This function calculates the classical kinetic energy of atoms
-//   and its contribution to stress.
-//----------------------------------------------------------------------------
-    kinetic = MD_func::GetAtomKE(unit_in.nat, vel, allmass);
+//   This function calculates the contribution of classical kinetic energy of atoms to stress.
+//--------------------------------------------------------------------------------------------
 
     if(GlobalV::CAL_STRESS)
     {
-        ModuleBase::matrix temp;
-        temp.create(3,3);    // initialize
+        ModuleBase::matrix t_vector;
 
-        for(int ion=0; ion<unit_in.nat; ++ion)
-        {
-            for(int i=0; i<3; ++i)
-            {
-                for(int j=i; j<3; ++j)
-                {
-                    temp(i, j) += allmass[ion] * vel[ion][i] * vel[ion][j];
-                }
-            }
-        }
+        temp_vector(unit_in.nat, vel, allmass, t_vector);
 
         for(int i=0; i<3; ++i)
         {
             for(int j=0; j<3; ++j)
             {
-                if(j<i) 
-                {
-                    stress(i, j) = stress(j, i);
-                }
-                else
-                {
-                    stress(i, j) = temp(i, j)/unit_in.omega;
-                }
+                stress(i, j) = virial(i, j) + t_vector(i, j) / unit_in.omega;
             }
         }
     }
@@ -227,7 +208,7 @@ void MD_func::force_virial(
 		UnitCell_pseudo &unit_in,
 		double &potential,
 		ModuleBase::Vector3<double> *force,
-		ModuleBase::matrix &stress)
+		ModuleBase::matrix &virial)
 {
 	ModuleBase::TITLE("MD_func", "force_stress");
     ModuleBase::timer::tick("MD_func", "force_stress");
@@ -241,14 +222,14 @@ void MD_func::force_virial(
 
     if(GlobalV::CAL_STRESS)
     {
-        p_esolver->cal_Stress(stress);
+        p_esolver->cal_Stress(virial);
     }
 
     if(mdp.md_ensolver == "FP")
     {
         potential *= 0.5;
         force_temp *= 0.5;
-        stress *= 0.5;
+        virial *= 0.5;
     }
 
     for(int i=0; i<unit_in.nat; ++i)
@@ -384,12 +365,32 @@ double MD_func::target_temp(const int &istep, const double &tfirst, const double
     return tfirst + delta * (tlast - tfirst);
 }
 
-double MD_func::current_temp(const int &natom, 
+double MD_func::current_temp(double &kinetic,
+            const int &natom, 
             const int &frozen_freedom, 
             const double *allmass,
             const ModuleBase::Vector3<double> *vel)
 {
-    double ke = GetAtomKE(natom, vel, allmass);
+    kinetic = GetAtomKE(natom, vel, allmass);
 
-    return 2 * ke / (3 * natom - frozen_freedom) * ModuleBase::Hartree_to_K;
+    return 2 * kinetic / (3 * natom - frozen_freedom) * ModuleBase::Hartree_to_K;
+}
+
+void MD_func::temp_vector(const int &natom, 
+            const ModuleBase::Vector3<double> *vel, 
+            const double *allmass, 
+            ModuleBase::matrix &t_vector)
+{
+    t_vector.create(3, 3);
+
+    for(int ion=0; ion<natom; ++ion)
+    {
+        for(int i=0; i<3; ++i)
+        {
+            for(int j=0; j<3; ++j)
+            {
+                t_vector(i, j) += allmass[ion] * vel[ion][i] * vel[ion][j];
+            }
+        }
+    }
 }
diff --git a/source/module_md/MD_func.h b/source/module_md/MD_func.h
@@ -43,18 +43,18 @@ class MD_func
 		UnitCell_pseudo &unit_in,
 		double &potential,
 		ModuleBase::Vector3<double> *force,
-		ModuleBase::matrix &stress);
+		ModuleBase::matrix &virial);
 
 	static double GetAtomKE(
 		const int &numIon,
 		const ModuleBase::Vector3<double> *vel, 
 		const double *allmass);
 
-	static void kinetic_stress(
+	static void compute_stress(
 		const UnitCell_pseudo &unit_in,
 		const ModuleBase::Vector3<double> *vel, 
 		const double *allmass, 
-		double &kinetic,
+        const ModuleBase::matrix &virial,
 		ModuleBase::matrix &stress);
 
 	static void outStress(const ModuleBase::matrix &virial, const ModuleBase::matrix &stress);
@@ -72,9 +72,15 @@ class MD_func
 
     static double target_temp(const int &istep, const double &tfirst, const double &tlast);
 
-    static double current_temp(const int &natom, 
+    static double current_temp(double &kinetic,
+            const int &natom, 
             const int &frozen_freedom, 
             const double *allmass,
             const ModuleBase::Vector3<double> *vel);
+
+    static void temp_vector(const int &natom, 
+            const ModuleBase::Vector3<double> *vel, 
+            const double *allmass, 
+            ModuleBase::matrix &t_vector);
 };
 #endif
diff --git a/source/module_md/MD_parameters.h b/source/module_md/MD_parameters.h
@@ -34,8 +34,14 @@ class MD_parameters
 		msst_tscale = 0.01;
 
 		// NHC
-		md_tfreq = 0.0;
-		md_mnhc = 4;
+        md_pmode = "none";
+        md_pcouple = "none";
+        md_tfreq = 0.0;
+        md_pfirst = -1.0;
+        md_plast = -1.0;
+        md_pfreq = 0.0;
+        md_tchain = 1;
+        md_pchain = 1;
 
 		// Langevin
 		md_damp = 1.0;
@@ -72,8 +78,14 @@ class MD_parameters
 	double msst_tscale;              // reduction in initial temperature (0~1)
 
 	// NHC
-	double md_tfreq;               // Oscillation frequency, used to determine msst_qmass of NHC
-	int md_mnhc;                   // num of NHC
+    std::string md_pmode;          // NPT ensemble mode: none, iso, aniso, tri
+    std::string md_pcouple;        // whether couple different components: xyz, xy, yz, xz, none
+    double md_tfreq;               // Oscillation frequency, used to determine qmass of thermostats coupled with particles
+    double md_pfirst;             // Initial pressure
+    double md_plast;              // Final pressure
+    double md_pfreq;               // Oscillation frequency, used to determine qmass of thermostats coupled with barostat
+    int md_tchain;                   // num of thermostats coupled with particles
+    int md_pchain;                   // num of thermostats coupled with barostat
 
 	// Langevin
 	double md_damp;                // damping parameter (time units)

diff --git a/source/module_md/MSST.cpp b/source/module_md/MSST.cpp
@@ -63,8 +63,8 @@ void MSST::setup(ModuleESolver::ESolver *p_esolver)
             }
         }
 
-        MD_func::kinetic_stress(ucell, vel, allmass, kinetic, stress);
-        stress += virial;
+        MD_func::compute_stress(ucell, vel, allmass, virial, stress);
+        t_current = MD_func::current_temp(kinetic, ucell.nat, frozen_freedom_, allmass, vel);
     }
 
     ModuleBase::timer::tick("MSST", "setup");
@@ -146,8 +146,8 @@ void MSST::second_half()
     propagate_vel();
 
     vsum = vel_sum();
-    MD_func::kinetic_stress(ucell, vel, allmass, kinetic, stress);
-    stress += virial;
+    MD_func::compute_stress(ucell, vel, allmass, virial, stress);
+    t_current = MD_func::current_temp(kinetic, ucell.nat, frozen_freedom_, allmass, vel);
 
     // propagate the time derivative of volume 1/2 step
     propagate_voldot();

diff --git a/source/module_md/Makefile.Objects b/source/module_md/Makefile.Objects
@@ -12,7 +12,7 @@ run_md_classic.o\
 mdrun.o\
 verlet.o\
 MSST.o\
-NVT_NHC.o\
+Nose_Hoover.o\
 FIRE.o\
 Langevin.o\