fix: a bug that MD stuck

It is because different cores get different different postions of atoms.
fix issue #932

Author: Qianruipku

source/module_md/FIRE.cpp

@@ -36,7 +36,8 @@ void FIRE::first_half()
 {
     ModuleBase::TITLE("FIRE", "first_half");
     ModuleBase::timer::tick("FIRE", "first_half");
-
+    if(GlobalV::MY_RANK == 0)
+    {
     for(int i=0; i<ucell.nat; ++i)
     {
         for(int k=0; k<3; ++k)
@@ -60,6 +61,10 @@ void FIRE::first_half()
             }
         }
     }
+    }
+#ifdef __MPI
+    MPI_Bcast(pos , ucell.nat*3,MPI_DOUBLE,0,MPI_COMM_WORLD);
+#endif
 
     ucell.update_pos_tau(pos);
     ucell.periodic_boundary_adjustment();

source/module_md/MD_func.cpp

@@ -75,7 +75,8 @@ void MD_func::kinetic_stress(
 //   This function calculates the classical kinetic energy of atoms
 //   and its contribution to stress.
 //----------------------------------------------------------------------------
-
+	if(GlobalV::MY_RANK==0) //only first rank do md
+	{
 	kinetic = MD_func::GetAtomKE(unit_in.nat, vel, allmass);
 
 	ModuleBase::matrix temp;
@@ -106,6 +107,7 @@ void MD_func::kinetic_stress(
 			}
 		}
 	}
+	}
 }
 
 // Read Velocity from STRU liuyu 2021-09-24

source/module_md/MSST.cpp

@@ -83,7 +83,9 @@ void MSST::first_half()
 
     const int sd = mdp.msst_direction;
     const double dthalf = 0.5 * mdp.md_dt;
-
+    double vol;
+    if( GlobalV::MY_RANK == 0 )
+    {
     energy_ = potential + kinetic;
 
     // propagate the time derivative of volume 1/2 step
@@ -112,7 +114,7 @@ void MSST::first_half()
     propagate_vel();
 
     // propagate volume 1/2 step
-    double vol = ucell.omega + omega[sd] * dthalf;
+    vol = ucell.omega + omega[sd] * dthalf;
 
     // rescale positions and change box size
     rescale(vol);
@@ -122,6 +124,11 @@ void MSST::first_half()
     {
         pos[i] += vel[i] * mdp.md_dt;
     }
+    }
+#ifdef __MPI
+    MPI_Bcast(pos , ucell.nat*3,MPI_DOUBLE,0,MPI_COMM_WORLD);
+#endif
+
     ucell.update_pos_tau(pos);
     ucell.periodic_boundary_adjustment();
 
@@ -141,7 +148,8 @@ void MSST::second_half()
 
     const int sd = mdp.msst_direction;
     const double dthalf = 0.5 * mdp.md_dt;
-
+    if( GlobalV::MY_RANK == 0 )
+    {
     energy_ = potential + kinetic;
 
     // propagate velocities 1/2 step
@@ -156,6 +164,7 @@ void MSST::second_half()
 
     // calculate Lagrangian position
     lag_pos -= mdp.msst_vel * ucell.omega / v0 * mdp.md_dt;
+    }
 
     ModuleBase::timer::tick("MSST", "second_half");
 }

source/module_md/verlet.cpp

@@ -72,6 +72,7 @@ void Verlet::setup(ModuleESolver::ESolver *p_esolver)
 
 void Verlet::first_half()
 {
+    if(GlobalV::MY_RANK==0) //only first rank do md
     for(int i=0; i<ucell.nat; ++i)
     {
         for(int k=0; k<3; ++k)
@@ -83,6 +84,9 @@ void Verlet::first_half()
             }
         }
     }
+#ifdef __MPI
+    MPI_Bcast(pos , ucell.nat*3,MPI_DOUBLE,0,MPI_COMM_WORLD);
+#endif
 
     ucell.update_pos_tau(pos);
     ucell.periodic_boundary_adjustment();
@@ -91,6 +95,7 @@ void Verlet::first_half()
 
 void Verlet::second_half()
 {
+    if(GlobalV::MY_RANK==0) //only first rank do md
     for(int i=0; i<ucell.nat; ++i)
     {
         for(int k=0; k<3; ++k)