read velocity from SRTU controlled by parameter set_vel

doc/examples/md.md

@@ -34,7 +34,6 @@ mixing_beta         0.4
 charge_extrap       second-order
 
 md_mdtype           1    //choose ensemble
-md_potential        0    //choose potential  0:Pseudopotential  1:LJ potential
 md_dt               1    //time step
 md_tfirst           700  //the first target temperature
 md_rstmd            0    //whether restart md
@@ -47,7 +46,6 @@ These MD parameters means that ABACUS will use NVT ensemble with Nosé-hoover th
 Note: *Please turn off symmetry when do MD simulation.*
 
 - md_mdtype : 0, NVE; 1, NVT; 2, velocity scaling
-- md_potential : 0, Pseudopotential  1, LJ potential
 - md_dt : time step in md simulation (fs)
 - md_tfirst : target temperature in md simulation(K), you should set parameter md_tlast and md_fixtemperature when you want to change temperature during md simulation.
 - md_rstmd : 0, no need of restart ; 1, restart with restart file, you must repalce STRU file with STRU_MD before you run the restart task.

source/driver.cpp

@@ -81,17 +81,8 @@ void Driver::atomic_world(void)
 	// lcao_in_pw: LCAO expaned by plane wave basis set
 	// lcao: linear combination of atomic orbitals
 	//--------------------------------------------------
-#ifndef __NOMD //qianrui do not want to add md to pw module temporarily before md is compeletly built.
-	if(CALCULATION=="md" || CALCULATION=="md-sto") // Yu Liu 2021-07-12
-	{
-		Run_md::md_line();
-	}
-	else if(BASIS_TYPE=="pw" || BASIS_TYPE=="lcao_in_pw")
-#else
 	if(BASIS_TYPE=="pw" || BASIS_TYPE=="lcao_in_pw")
-#endif
 	{
-
 		Run_pw::plane_wave_line();
 	}
 #ifdef __LCAO

source/input.cpp

@@ -165,6 +165,7 @@ void Input::Default(void)
     search_radius=-1.0; // unit: a.u. -1.0 has no meaning.
     search_pbc=true;
     symmetry=false;
+	set_vel=false;
 	mlwf_flag=false;
     force=0;
     force_set=false;
@@ -710,6 +711,10 @@ bool Input::Read(const string &fn)
         else if (strcmp("symmetry", word) == 0)
         {
             read_value(ifs, symmetry);
+        }
+		else if (strcmp("set_vel", word) == 0)
+        {
+            read_value(ifs, set_vel);
         }
         else if (strcmp("mlwf_flag", word) == 0)
         {
@@ -1148,10 +1153,10 @@ bool Input::Read(const string &fn)
 		{
 			read_value(ifs, mdp.mdtype);
 		}
-		else if (strcmp("md_potential",word) == 0)
-		{
-			read_value(ifs, mdp.md_potential);
-		}
+		//else if (strcmp("md_potential",word) == 0)
+		//{
+		//	read_value(ifs, mdp.md_potential);
+		//}
 		else if (strcmp("NVT_tau",word) == 0)
 		{
 			read_value(ifs, mdp.NVT_tau);
@@ -2010,6 +2015,7 @@ void Input::Bcast()
 	Parallel_Common::bcast_bool( search_pbc );
     Parallel_Common::bcast_double( search_radius );
     Parallel_Common::bcast_bool( symmetry );
+	Parallel_Common::bcast_bool( set_vel );  //liuyu 2021-07-14
     Parallel_Common::bcast_bool( mlwf_flag );
     Parallel_Common::bcast_int( force );
     Parallel_Common::bcast_bool( force_set );
@@ -2127,7 +2133,7 @@ void Input::Bcast()
 */
 	//zheng daye add 2014/5/5
         Parallel_Common::bcast_int(mdp.mdtype);
-		Parallel_Common::bcast_int(mdp.md_potential);
+		//Parallel_Common::bcast_int(mdp.md_potential);
         Parallel_Common::bcast_double(mdp.NVT_tau);
         Parallel_Common::bcast_int(mdp.NVT_control);
         Parallel_Common::bcast_double(mdp.dt);

source/input.h

@@ -51,6 +51,8 @@ class Input
 	double emin_sto;
 	string stotype;
 
+	bool set_vel;           // read velocity from STRU or not  liuyu 2021-07-14
+
     bool symmetry;			// turn on symmetry or not
 	int npool; 				// ecch pool is for one k point
 

source/module_cell/atom_spec.cpp

@@ -13,6 +13,7 @@ Atom::Atom()
     stapos_wf = 0;
     tau = new Vector3<double>[1];
     taud = new Vector3<double>[1];
+    vel = new Vector3<double>[1];
     mag = new double[1];
     l_nchi = new int[1];
     iw2l = new int[1];
@@ -26,6 +27,7 @@ Atom::~Atom()
 {
     delete[] tau;
     delete[] taud;
+    delete[] vel;
     delete[] mag;
     delete[] l_nchi;
     delete[] iw2l;

source/module_cell/atom_spec.h

@@ -33,6 +33,7 @@ class Atom: public Atom_pseudo
     string label; // atomic symbol
     Vector3<double> *tau;// Cartesian coordinates of each atom in this type.
 	Vector3<double> *taud;// Direct coordinates of each atom in this type.
+    Vector3<double> *vel;// velocities of each atom in this type.
 
 	double* mag;
 	double angle1;//spin angle, added by zhengdy-soc

source/module_cell/read_atoms.cpp

@@ -398,10 +398,12 @@ bool UnitCell_pseudo::read_atom_positions(ifstream &ifpos)
 			{
        			delete[] atoms[it].tau;
 				delete[] atoms[it].taud;
+				delete[] atoms[it].vel;
        			delete[] atoms[it].mbl;
 				delete[] atoms[it].mag;
        			atoms[it].tau = new Vector3<double>[na];
        			atoms[it].taud = new Vector3<double>[na];
+				atoms[it].vel = new Vector3<double>[na];
        			atoms[it].mbl = new Vector3<int>[na];
 				atoms[it].mag = new double[na];
 				atoms[it].mass = this->atom_mass[it]; //mohan add 2011-11-07 
@@ -436,6 +438,14 @@ bool UnitCell_pseudo::read_atom_positions(ifstream &ifpos)
 					{
 						ifpos >> v.x >> v.y >> v.z
 							>> mv.x >> mv.y >> mv.z;
+						if(INPUT.set_vel)
+						{
+							ifpos >> atoms[it].vel[ia].x >> atoms[it].vel[ia].y >> atoms[it].vel[ia].z;
+						}
+						else
+						{
+							atoms[it].vel[ia].set(0,0,0);
+						}
 						string mags;
 						std::getline( ifpos, mags );
 						// change string to double.
@@ -675,7 +685,8 @@ void UnitCell_pseudo::print_stru_file(const string &fn, const int &type)const
 			for(int ia=0; ia<atoms[it].na; ia++)
 			{
 				ofs << atoms[it].tau[ia].x << " " << atoms[it].tau[ia].y << " " << atoms[it].tau[ia].z << " " 
-					<< atoms[it].mbl[ia].x << " " << atoms[it].mbl[ia].y << " " << atoms[it].mbl[ia].z << endl;
+					<< atoms[it].mbl[ia].x << " " << atoms[it].mbl[ia].y << " " << atoms[it].mbl[ia].z << " "
+					<< atoms[it].vel[ia].x << " " << atoms[it].vel[ia].y << " " << atoms[it].vel[ia].z << endl;
 			}
 		}
 	}
@@ -697,7 +708,8 @@ void UnitCell_pseudo::print_stru_file(const string &fn, const int &type)const
 			for(int ia=0; ia<atoms[it].na; ia++)
 			{
 				ofs << atoms[it].taud[ia].x << " " << atoms[it].taud[ia].y << " " << atoms[it].taud[ia].z << " " 
-					<< atoms[it].mbl[ia].x << " " << atoms[it].mbl[ia].y << " " << atoms[it].mbl[ia].z << endl;
+					<< atoms[it].mbl[ia].x << " " << atoms[it].mbl[ia].y << " " << atoms[it].mbl[ia].z << " "
+					<< atoms[it].vel[ia].x << " " << atoms[it].vel[ia].y << " " << atoms[it].vel[ia].z << endl;
 			}
 		}
 	}
@@ -723,6 +735,9 @@ void UnitCell_pseudo::print_tau(void)const
         << setw(20) << "y"
         << setw(20) << "z"
         << setw(20) << "mag"
+		<< setw(20) << "vx"
+        << setw(20) << "vy"
+        << setw(20) << "vz"
         << endl;
         ofs_running << setprecision(12);
 
@@ -744,6 +759,9 @@ void UnitCell_pseudo::print_tau(void)const
                 << setw(20) << atoms[it].tau[ia].z
                 //<< setw(20) << atoms[it].mag[ia]
                 << setw(20) << magnet.start_magnetization[it]
+				<< setw(20) << atoms[it].vel[ia].x
+                << setw(20) << atoms[it].vel[ia].y
+                << setw(20) << atoms[it].vel[ia].z
                 << endl;
 
                 ++iat;
@@ -763,6 +781,9 @@ void UnitCell_pseudo::print_tau(void)const
         << setw(20) << "y"
         << setw(20) << "z"
         << setw(20) << "mag"
+		<< setw(20) << "vx"
+        << setw(20) << "vy"
+        << setw(20) << "vz"
         << endl;
 
         int iat=0;
@@ -783,6 +804,9 @@ void UnitCell_pseudo::print_tau(void)const
                 << setw(20) << atoms[it].taud[ia].z
                 //<< setw(20) << atoms[it].mag[ia]
                 << setw(20) << magnet.start_magnetization[it]
+				<< setw(20) << atoms[it].vel[ia].x
+                << setw(20) << atoms[it].vel[ia].y
+                << setw(20) << atoms[it].vel[ia].z
                 << endl;
 
                 ++iat;

source/module_cell/unitcell_pseudo.cpp

@@ -397,7 +397,6 @@ void UnitCell_pseudo::setup_cell_classic(const string &fn, output &outp, ofstrea
 
 	ofs_running << endl;
 	outp.printM3(ofs_running,"Lattice vectors: (Cartesian coordinate: in unit of a_0)",latvec); 
-
 }
 
 

source/module_md/MD_basic.cpp

@@ -22,7 +22,18 @@ MD_basic::MD_basic(MD_parameters& MD_para_in, UnitCell_pseudo &unit_in):
     step_=0;
 //	ucell.latvec=ucell.latvec;
 
-    vel=new Vector3<double>[ucell.nat];
+    vel=new Vector3<double>[ucell.nat]; //initialize velocity of atoms from STRU  Yu Liu 2021-07-14
+    int iat=0;
+    for(int it=0; it<ucell.ntype; ++it)
+    {
+        for(int ia=0; ia<ucell.atoms[it].na; ++ia)
+        {
+            vel[iat] = ucell.atoms[it].vel[ia];
+            ++iat;
+        }
+    }
+    assert(iat==ucell.nat);
+
 	cart_change=new Vector3<double>[ucell.nat];
 	tauDirectChange=new Vector3<double>[ucell.nat];
 	allmass=new double[ucell.nat];

source/module_md/MD_parameters.h

@@ -7,7 +7,7 @@ class MD_parameters
     MD_parameters()
 	{
 		mdtype=1;
-		md_potential=0;
+		//md_potential=0;
 		NVT_tau=0;
 		NVT_control=1;
 		dt=1;
@@ -27,7 +27,7 @@ class MD_parameters
     int rstMD;          // 1 : restart MD, vel.restart and ion.restart files will be read
                         // 0 : not restart from ion and vel files
 	int mdtype;         // 1: NVT:Nose Hoover, 2:NVT:velocity scaling 3: NPT, 0: NVE 
-	int md_potential;   // 0: Pseudopotential, 1: LJ potential   Yu Liu 2021-07-12
+	//int md_potential;   // 0: Pseudopotential, 1: LJ potential   Yu Liu 2021-07-12
     double Qmass;       // Inertia of extended system variable
     double dt ;         // Time increment (hbar/E_hartree)
     double tfirst; //Temperature (in Hartree, 1 Hartree ~ 3E5 K)