Feat: NPT ensemble (#1491)

* Feat: implement NPT ensemble * test: add test for iso NPT * Feat: restart npt
deepmodeling · Nov 14, 2022 · 3cce826 · 3cce826
1 parent 61843c9
commit 3cce826
Show file tree

Hide file tree

Showing 18 changed files with 755 additions and 54 deletions.
diff --git a/source/input.cpp b/source/input.cpp
@@ -1200,8 +1200,10 @@ bool Input::Read(const std::string &fn)
         {
             read_value(ifs, lcao_rmax);
         }
-        // about molecular dynamics
-        // added begin by zheng daye
+        //----------------------------------------------------------
+        // Molecule Dynamics
+        // Yu Liu add 2021-07-30
+        //----------------------------------------------------------
         else if (strcmp("md_type", word) == 0)
         {
             read_value(ifs, mdp.md_type);
@@ -1254,11 +1256,30 @@ bool Input::Read(const std::string &fn)
         {
             read_value(ifs, mdp.md_restart);
         }
-        // added by zheng daye
-        //----------------------------------------------------------
-        // Classic MD
-        // Yu Liu add 2021-07-30
-        //----------------------------------------------------------
+        else if (strcmp("md_pmode", word) == 0)
+        {
+            read_value(ifs, mdp.md_pmode);
+        }
+        else if (strcmp("md_pcouple", word) == 0)
+        {
+            read_value(ifs, mdp.md_pcouple);
+        }
+        else if (strcmp("md_pchain", word) == 0)
+        {
+            read_value(ifs, mdp.md_pchain);
+        }
+        else if (strcmp("md_pfirst", word) == 0)
+        {
+            read_value(ifs, mdp.md_pfirst);
+        }
+        else if (strcmp("md_plast", word) == 0)
+        {
+            read_value(ifs, mdp.md_plast);
+        }
+        else if (strcmp("md_pfreq", word) == 0)
+        {
+            read_value(ifs, mdp.md_pfreq);
+        }
         else if (strcmp("lj_rcut", word) == 0)
         {
             read_value(ifs, mdp.lj_rcut);
@@ -2320,6 +2341,12 @@ void Input::Bcast()
     Parallel_Common::bcast_string(mdp.pot_file);
     Parallel_Common::bcast_int(mdp.md_nraise);
     Parallel_Common::bcast_string(mdp.md_thermostat);
+    Parallel_Common::bcast_string(mdp.md_pmode);
+    Parallel_Common::bcast_string(mdp.md_pcouple);
+    Parallel_Common::bcast_int(mdp.md_pchain);
+    Parallel_Common::bcast_double(mdp.md_pfirst);
+    Parallel_Common::bcast_double(mdp.md_plast);
+    Parallel_Common::bcast_double(mdp.md_pfreq);
     // Yu Liu add 2022-05-18
     Parallel_Common::bcast_bool(efield_flag);
     Parallel_Common::bcast_bool(dip_cor_flag);
@@ -2651,7 +2678,7 @@ void Input::Check(void)
             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)
+        if(mdp.md_type == 1 && 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;
@@ -2668,7 +2695,7 @@ void Input::Check(void)
         {
             init_vel = 1;
         }
-        if(esolver_type == "lj" || esolver_type == "dp" || mdp.md_type == 4)
+        if(esolver_type == "lj" || esolver_type == "dp" || mdp.md_type == 4 || (mdp.md_type == 1 && mdp.md_pmode != "none"))
         {
             cal_stress = 1;
         }

diff --git a/source/module_md/MD_func.cpp b/source/module_md/MD_func.cpp
@@ -367,7 +367,7 @@ double MD_func::current_temp(double &kinetic,
 {
     kinetic = GetAtomKE(natom, vel, allmass);
 
-    return 2 * kinetic / (3 * natom - frozen_freedom) * ModuleBase::Hartree_to_K;
+    return 2 * kinetic / (3 * natom - frozen_freedom);
 }
 
 void MD_func::temp_vector(const int &natom, 

diff --git a/source/module_md/MSST.cpp b/source/module_md/MSST.cpp
@@ -251,9 +251,6 @@ void MSST::rescale(double volume)
     ucell.latvec.e11 *= dilation[0];
     ucell.latvec.e22 *= dilation[1];
     ucell.latvec.e33 *= dilation[2];
-    ucell.a1 *= dilation[0];
-    ucell.a2 *= dilation[1];
-    ucell.a3 *= dilation[2];
 
     ucell.setup_cell_after_vc(GlobalV::ofs_running);
     MD_func::InitPos(ucell, pos);