Test the correctness of the modified MD module

source/module_md/md_base.cpp

@@ -29,7 +29,7 @@ MD_base::MD_base(MD_para& MD_para_in, UnitCell& unit_in)
     assert(ModuleBase::AU_to_FS!=0.0);
     assert(ModuleBase::Hartree_to_K!=0.0);
 
-    /// convert to a.u. unit
+    /// convert time and temperatrue to a.u. unit
     mdp.md_dt /= ModuleBase::AU_to_FS;
     mdp.md_tfirst /= ModuleBase::Hartree_to_K;
     mdp.md_tlast /= ModuleBase::Hartree_to_K;

source/module_md/md_base.h

@@ -8,7 +8,7 @@
  * @brief base class of md
  *
  * This class implements the velocity-Verlet method.
- * The system is assumed to be isolated in the sense that it cannot exchange energy or particles with its environment,
+ * The corresponding system(nve) is assumed to be isolated in the sense that it cannot exchange energy or particles with its environment,
  * so that the energy of the system does not change with time.
  */
 class MD_base
@@ -88,4 +88,4 @@ class MD_base
     double energy_;  ///< total energy of the system
 };
 
-#endif // MD_BASE_H
+#endif // MD_BASE_H

source/module_md/msst.cpp

@@ -269,38 +269,82 @@ void MSST::rescale(std::ofstream& ofs, const double& volume)
 
 void MSST::propagate_vel(void)
 {
-    if (mdp.my_rank == 0)
+#ifndef __MPI
+    const int sd = mdp.msst_direction;
+    const double dthalf = 0.5 * mdp.md_dt;
+    const double fac = mdp.msst_vis * pow(omega[sd], 2) / (vsum * ucell.omega);
+
+    for (int i = 0; i < ucell.nat; ++i)
     {
-        const int sd = mdp.msst_direction;
-        const double dthalf = 0.5 * mdp.md_dt;
-        const double fac = mdp.msst_vis * pow(omega[sd], 2) / (vsum * ucell.omega);
+        ModuleBase::Vector3<double> const_C = force[i] / allmass[i];
+        ModuleBase::Vector3<double> const_D;
+        const_D.set(fac / allmass[i], fac / allmass[i], fac / allmass[i]);
+        const_D[sd] -= 2 * omega[sd] / ucell.omega;
 
-        for (int i = 0; i < ucell.nat; ++i)
+        for (int k = 0; k < 3; ++k)
         {
-            ModuleBase::Vector3<double> const_C = force[i] / allmass[i];
-            ModuleBase::Vector3<double> const_D;
-            const_D.set(fac / allmass[i], fac / allmass[i], fac / allmass[i]);
-            const_D[sd] -= 2 * omega[sd] / ucell.omega;
-
-            for (int k = 0; k < 3; ++k)
+            if (fabs(dthalf * const_D[k]) > 1e-6)
+            {
+                double expd = exp(dthalf * const_D[k]);
+                vel[i][k] = expd * (const_C[k] + const_D[k] * vel[i][k] - const_C[k] / expd) / const_D[k];
+            }
+            else
             {
-                if (fabs(dthalf * const_D[k]) > 1e-6)
-                {
-                    double expd = exp(dthalf * const_D[k]);
-                    vel[i][k] = expd * (const_C[k] + const_D[k] * vel[i][k] - const_C[k] / expd) / const_D[k];
-                }
-                else
-                {
-                    vel[i][k]
-                        += (const_C[k] + const_D[k] * vel[i][k]) * dthalf
-                           + 0.5 * (const_D[k] * const_D[k] * vel[i][k] + const_C[k] * const_D[k]) * dthalf * dthalf;
-                }
+                vel[i][k]
+                    += (const_C[k] + const_D[k] * vel[i][k]) * dthalf
+                        + 0.5 * (const_D[k] * const_D[k] * vel[i][k] + const_C[k] * const_D[k]) * dthalf * dthalf;
             }
         }
     }
-
+#endif
+
 #ifdef __MPI
-    MPI_Bcast(vel, ucell.nat * 3, MPI_DOUBLE, 0, MPI_COMM_WORLD);
+    int size;
+    MPI_Comm_size(MPI_COMM_WORLD, &size);
+    const int sd = mdp.msst_direction;
+    const double dthalf = 0.5 * mdp.md_dt;
+    const double fac = mdp.msst_vis * pow(omega[sd], 2) / (vsum * ucell.omega);
+    int each_ucell_nat = ucell.nat / size;
+    int ucell_nat_begin = (each_ucell_nat) * mdp.my_rank;
+    int ucell_nat_end = ucell_nat_begin + each_ucell_nat;
+    if(mdp.my_rank == size - 1)
+    {
+        ucell_nat_end = ucell.nat;
+    }
+    for(int i = ucell_nat_begin; i < ucell_nat_end; i++)
+    { 
+        ModuleBase::Vector3<double> const_C = force[i] / allmass[i];
+        ModuleBase::Vector3<double> const_D;
+        const_D.set(fac / allmass[i], fac / allmass[i], fac / allmass[i]);
+        const_D[sd] -= 2 * omega[sd] / ucell.omega;
+
+        for (int k = 0; k < 3; ++k)
+        {
+            if (fabs(dthalf * const_D[k]) > 1e-6)
+            {
+                double expd = exp(dthalf * const_D[k]);
+                vel[i][k] = expd * (const_C[k] + const_D[k] * vel[i][k] - const_C[k] / expd) / const_D[k];
+            }
+            else
+            {
+                vel[i][k]
+                    += (const_C[k] + const_D[k] * vel[i][k]) * dthalf
+                        + 0.5 * (const_D[k] * const_D[k] * vel[i][k] + const_C[k] * const_D[k]) * dthalf * dthalf;
+            }
+        }
+    }
+
+    for(int i = 0;i < size; i++){
+        int each_ucell_nat = ucell.nat / size;
+        int ucell_nat_begin = (each_ucell_nat) * i;
+        int ucell_nat_end = ucell_nat_begin + each_ucell_nat;
+        if(i == size - 1)
+        {
+            ucell_nat_end = ucell.nat;
+        }
+        MPI_Bcast(vel + ucell_nat_begin, (ucell_nat_end - ucell_nat_begin) * 3, MPI_DOUBLE, i, MPI_COMM_WORLD);
+    }
+
 #endif
 
     return;

source/module_md/run_md.h

@@ -7,6 +7,7 @@
 /**
  * @brief the md loop line
  *
+ * Determine different classes and methods based on md_type.
  */
 namespace Run_MD
 {
@@ -20,4 +21,4 @@ namespace Run_MD
 void md_line(UnitCell& unit_in, ModuleESolver::ESolver* p_esolver, MD_para& md_para);
 } // namespace Run_MD
 
-#endif
+#endif

source/module_md/verlet.cpp

@@ -75,21 +75,7 @@ void Verlet::apply_thermostat(void)
     {
         if (mdp.my_rank == 0)
         {
-            double deviation;
-            for (int i = 0; i < ucell.nat; ++i)
-            {
-                if (static_cast<double>(std::rand()) / RAND_MAX <= 1.0 / mdp.md_nraise)
-                {
-                    deviation = sqrt(mdp.md_tlast / allmass[i]);
-                    for (int k = 0; k < 3; ++k)
-                    {
-                        if (ionmbl[i][k])
-                        {
-                            vel[i][k] = deviation * MD_func::gaussrand();
-                        }
-                    }
-                }
-            }
+            anderson_set_vel();
         }
 #ifdef __MPI
         MPI_Bcast(vel, ucell.nat * 3, MPI_DOUBLE, 0, MPI_COMM_WORLD);
@@ -145,3 +131,23 @@ void Verlet::restart(const std::string& global_readin_dir)
     MD_base::restart(global_readin_dir);
     return;
 }
+
+void Verlet::anderson_set_vel(void)
+{
+    double deviation;
+    for (int i = 0; i < ucell.nat; ++i)
+    {
+        if (static_cast<double>(std::rand()) / RAND_MAX <= 1.0 / mdp.md_nraise)
+        {
+            deviation = sqrt(mdp.md_tlast / allmass[i]);
+            for (int k = 0; k < 3; ++k)
+            {
+                if (ionmbl[i][k])
+                {
+                    vel[i][k] = deviation * MD_func::gaussrand();
+                }
+            }
+        }
+    }
+    return;
+}

source/module_md/verlet.h

@@ -6,6 +6,7 @@
 /**
  * @brief the md methods based on the velocity-Verlet equation
  *
+ * This method corresponds to the NVE and the NVT ensemble。
  */
 class Verlet : public MD_base
 {
@@ -35,6 +36,12 @@ class Verlet : public MD_base
      * @param target_temp the target temperature
      */
     void thermalize(const int& nraise, const double& current_temp, const double& target_temp);
-};
 
-#endif
+    /**
+    * @brief set verlet when themostat is anderson
+    *
+    */
+    void anderson_set_vel();
+
+};
+#endif