Reconstruction MD:

ABACUS.develop/source/Makefile.Objects

@@ -53,9 +53,10 @@ chi0_hilbert.o\
 chi0_standard.o\
 epsilon0_pwscf.o\
 epsilon0_vasp.o\
-md.o\
-mdNVE.o\
-mdNVT.o\
+MD_basic.o\
+MD_thermo.o\
+MD_fire.o\
+MD_func.o\
 exx_lip.o\
 soc.o\
 to_wannier90.o \

ABACUS.develop/source/input.cpp

@@ -292,27 +292,6 @@ void Input::Default(void)
 	md_tstep=1; //reduec md_delt every md_tstep step.
 	md_delt=1.0;
 */
-	//md and related parameters(added by zheng da ye)
-	md_mdtype=1;
-	md_tauthermo=0;
-	md_taubaro=0;
-	md_dt=-1;
-	md_nresn=3;
-	md_nyosh=3;
-	md_qmass=1;
-	md_tfirst=-1;         //kelvin
-	md_tlast=md_tfirst;
-	md_dumpmdfred=1;
-	md_mdoutpath="mdoutput";
-	md_domsd=0;
-	md_domsdatom=0;
-	md_rstmd=0;
-	md_outputstressperiod=1;
-	md_fixtemperature=1;
-	md_ediff=1e-4;
-	md_ediffg=1e-3;
-	md_msdstartTime=1;
-	//end of zhengdaye's add.
 
 /* //----------------------------------------------------------
 // vdwD2									//Peize Lin add 2014-03-31, update 2015-09-30
@@ -1128,79 +1107,59 @@ bool Input::Read(const string &fn)
 //added begin by zheng daye
 		else if (strcmp("md_mdtype",word) == 0)
 		{
-			read_value(ifs, md_mdtype);
+			read_value(ifs, mdp.mdtype);
 		}
-		else if (strcmp("md_tauthermo",word) == 0)
+		else if (strcmp("NVT_tau",word) == 0)
 		{
-			read_value(ifs, md_tauthermo);
+			read_value(ifs, mdp.NVT_tau);
 		}
-		else if (strcmp("md_taubaro",word) == 0)
+		else if (strcmp("NVT_control",word) == 0)
 		{
-			read_value(ifs,md_taubaro );
+			read_value(ifs,mdp.NVT_control );
 		}
 		else if (strcmp("md_dt",word) == 0)
 		{
-			read_value(ifs, md_dt);
+			read_value(ifs, mdp.dt);
 		}
-		else if (strcmp("md_nresn",word) == 0)
+		else if (strcmp("mnhc",word) == 0)
 		{
-			read_value(ifs,md_nresn );
-		}
-		else if (strcmp("md_nyosh",word) == 0)
-		{
-			read_value(ifs, md_nyosh);
+			read_value(ifs,mdp.MNHC );
 		}
 		else if (strcmp("md_qmass",word) == 0)
 		{
-			read_value(ifs,md_qmass );
+			read_value(ifs,mdp.Qmass );
 		}
 		else if (strcmp("md_tfirst",word) == 0)
 		{
-			read_value(ifs, md_tfirst);
+			read_value(ifs, mdp.tfirst);
 		}
 		else if (strcmp("md_tlast",word) == 0)
 		{
-			read_value(ifs,md_tlast );
+			read_value(ifs,mdp.tlast );
 		}
 		else if (strcmp("md_dumpmdfred",word) == 0)
 		{
-			read_value(ifs, md_dumpmdfred);
+			read_value(ifs, mdp.recordFreq);
 		}
 		else if (strcmp("md_mdoutpath",word) == 0)
 		{
-			read_value(ifs,md_mdoutpath );
-		}
-		else if (strcmp("md_domsd",word) == 0)
-		{
-			read_value(ifs, md_domsd);
-		}
-		else if (strcmp("md_domsdatom",word) == 0)
-		{
-			read_value(ifs, md_domsdatom);
+			read_value(ifs,mdp.mdoutputpath );
 		}
 		else if (strcmp("md_rstmd",word) == 0)
 		{
-			read_value(ifs,md_rstmd );
-		}
-		else if (strcmp("md_outputstressperiod",word) == 0)
-		{
-			read_value(ifs,md_outputstressperiod );
+			read_value(ifs,mdp.rstMD );
 		}
 		else if (strcmp("md_fixtemperature",word) == 0)
 		{
-			read_value(ifs,md_fixtemperature );
+			read_value(ifs,mdp.fixTemperature );
 		}
 		else if (strcmp("md_ediff",word) == 0)
 		{
-			read_value(ifs,md_ediff );
+			read_value(ifs,mdp.ediff );
 		}
 		else if (strcmp("md_ediffg",word) == 0)
 		{
-			read_value(ifs,md_ediffg );
-		}
-		else if (strcmp("md_msdstarttime",word) == 0)
-		{
-			read_value(ifs,md_msdstartTime );
+			read_value(ifs,mdp.ediffg );
 		}
 //added by zheng daye
 //----------------------------------------------------------
@@ -2087,7 +2046,7 @@ void Input::Bcast()
 	Parallel_Common::bcast_int( selinv_niter);
 /*
 	// mohan add 2011-11-07
-	Parallel_Common::bcast_double( md_dt );
+	Parallel_Common::bcast_double( mdp.dt );
 	Parallel_Common::bcast_int( md_restart );
 	Parallel_Common::bcast_double( md_tolv );
 	Parallel_Common::bcast_string( md_thermostat );
@@ -2096,25 +2055,20 @@ void Input::Bcast()
 	Parallel_Common::bcast_double( md_delt );
 */
 	//zheng daye add 2014/5/5
-        Parallel_Common::bcast_int(md_mdtype);
-        Parallel_Common::bcast_double(md_tauthermo);
-        Parallel_Common::bcast_double(md_taubaro);
-        Parallel_Common::bcast_double(md_dt);
-        Parallel_Common::bcast_int(md_nresn);
-        Parallel_Common::bcast_int(md_nyosh);
-        Parallel_Common::bcast_double(md_qmass);
-        Parallel_Common::bcast_double(md_tfirst);
-        Parallel_Common::bcast_double(md_tlast);
-        Parallel_Common::bcast_int(md_dumpmdfred);
-        Parallel_Common::bcast_string(md_mdoutpath);
-        Parallel_Common::bcast_bool(md_domsd);
-        Parallel_Common::bcast_bool(md_domsdatom);
-        Parallel_Common::bcast_int(md_rstmd);
-        Parallel_Common::bcast_int(md_outputstressperiod);
-        Parallel_Common::bcast_int(md_fixtemperature);
-        Parallel_Common::bcast_double(md_ediff);
-        Parallel_Common::bcast_double(md_ediffg);
-        Parallel_Common::bcast_int(md_msdstartTime);
+        Parallel_Common::bcast_int(mdp.mdtype);
+        Parallel_Common::bcast_double(mdp.NVT_tau);
+        Parallel_Common::bcast_int(mdp.NVT_control);
+        Parallel_Common::bcast_double(mdp.dt);
+        Parallel_Common::bcast_int(mdp.MNHC);
+        Parallel_Common::bcast_double(mdp.Qmass);
+        Parallel_Common::bcast_double(mdp.tfirst);
+        Parallel_Common::bcast_double(mdp.tlast);
+        Parallel_Common::bcast_int(mdp.recordFreq);
+        Parallel_Common::bcast_string(mdp.mdoutputpath);
+        Parallel_Common::bcast_int(mdp.rstMD);
+        Parallel_Common::bcast_int(mdp.fixTemperature);
+        Parallel_Common::bcast_double(mdp.ediff);
+        Parallel_Common::bcast_double(mdp.ediffg);
 /* 	// Peize Lin add 2014-04-07
 	Parallel_Common::bcast_bool( vdwD2 );
 	Parallel_Common::bcast_double( vdwD2_scaling );
@@ -2483,10 +2437,10 @@ void Input::Check(void)
 
         //deal with input parameters , 2019-04-30
         if(basis_type == "pw" ) WARNING_QUIT("Input::Check","calculate = MD is only availble for LCAO.");
-        if(md_dt == -1) WARNING_QUIT("Input::Check","time interval of MD calculation should be set!");
-        if(md_tfirst == -1) WARNING_QUIT("Input::Check","temperature of MD calculation should be set!");
-        if(md_tlast  == -1) md_tlast = md_tfirst;
-        if(md_tfirst!=md_tlast)
+        if(mdp.dt < 0) WARNING_QUIT("Input::Check","time interval of MD calculation should be set!");
+        if(mdp.tfirst < 0) WARNING_QUIT("Input::Check","temperature of MD calculation should be set!");
+        if(mdp.tlast  < 0.0) mdp.tlast = mdp.tfirst;
+        if(mdp.tfirst!=mdp.tlast)
         {
             ifstream file1;
             file1.open("ChangeTemp.dat");
@@ -2496,7 +2450,7 @@ void Input::Check(void)
                 file.open("ChangeTemp.dat");
                 for(int ii=0;ii<30;ii++)
                 {
-                    file<<md_tfirst+(md_tlast-md_tfirst)/double(30)*double(ii+1)<<" ";
+                    file<<mdp.tfirst+(mdp.tlast-mdp.tfirst)/double(30)*double(ii+1)<<" ";
                 }
                 file.close();
             }

ABACUS.develop/source/input.h

@@ -7,6 +7,7 @@
 
 #include "src_global/vector3.h"
 #include <string>
+#include "src_pw/MD_parameters.h"
 using namespace std;
 
 class Input
@@ -244,7 +245,7 @@ class Input
 // molecular dynamics 
 // added by Daye Zheng
 //==========================================================
-    int md_mdtype;                   //choose ensemble
+/*    int md_mdtype;                   //choose ensemble
 	double md_tauthermo;
 	double md_taubaro;
 	double md_dt;                    //time step
@@ -262,7 +263,8 @@ class Input
 	int md_fixtemperature;          //period to change temperature
 	double md_ediff;             //parameter for constraining total energy change
 	double md_ediffg;             //parameter for constraining max force change
-	int md_msdstartTime;            //choose which step that msd be calculated 
+	int md_msdstartTime;            //choose which step that msd be calculated */
+	MD_parameters mdp;
 
 //==========================================================
 // vdw

ABACUS.develop/source/src_io/print_info.cpp

@@ -37,16 +37,17 @@ void Print_Info::setup_parameters(void)
 
 			cout << " ---------------------------------------------------------" << endl;
 
-			if(INPUT.md_mdtype ==1 || INPUT.md_mdtype==2)
+			if(INPUT.mdp.mdtype ==1 || INPUT.mdp.mdtype==2)
 			{
 				cout << " ENSEMBLE                 : " << "NVT" << endl;
+				cout << " Qmass for NVT(a.u.)      : " << INPUT.mdp.Qmass/6.02/9.109*1e5 << endl;
 			}
-			else if(INPUT.md_mdtype==0)
+			else if(INPUT.mdp.mdtype==0)
 			{
 				cout << " ENSEMBLE                 : " << "NVE" << endl;
 			}
-			cout << " Qmass for NVT(a.u.)      : " << INPUT.md_qmass/6.02/9.109*1e5 << endl;
-			cout << " Time interval(fs)        : " << INPUT.md_dt << endl;
+
+			cout << " Time interval(fs)        : " << INPUT.mdp.dt << endl;
 		}
 		cout << " ---------------------------------------------------------" << endl;
 

ABACUS.develop/source/src_io/write_input.cpp

@@ -141,22 +141,20 @@ void Input::Print(const string &fn)const
 	OUTP(ofs,"selinv_niter",selinv_niter,"max number of steps to update mu");
 
 	ofs << "\n#Parameters (10.Molecular dynamics)" << endl;
-	OUTP(ofs,"md_mdtype",md_mdtype,"choose ensemble");
-	OUTP(ofs,"md_dt",md_dt,"time step");
-	OUTP(ofs,"md_nresn",md_nresn,"parameter during integrater");
-	OUTP(ofs,"md_nyosh",md_nyosh,"parameter during integrater");
-	OUTP(ofs,"md_qmass",md_qmass,"mass of thermostat");
-	OUTP(ofs,"md_tfirst",md_tfirst,"temperature first");
-	OUTP(ofs,"md_tlast",md_tlast,"temperature last");
-	OUTP(ofs,"md_dumpmdfred",md_dumpmdfred,"The period to dump MD information for monitoring and restarting MD");
-	OUTP(ofs,"md_mdoutpath",md_mdoutpath,"output path of md");
-	OUTP(ofs,"md_domsd",md_domsd,"whether compute <r(t)-r(0)>");
-	OUTP(ofs,"md_domsdatom",md_domsdatom,"whether compute msd for each atom");
-	OUTP(ofs,"md_rstmd",md_rstmd,"whether restart");
-	OUTP(ofs,"md_fixtemperature",md_fixtemperature,"period to change temperature");
-	OUTP(ofs,"md_ediff",md_ediff,"parameter for constraining total energy change");
-	OUTP(ofs,"md_ediffg",md_ediffg,"parameter for constraining max force change");
-	OUTP(ofs,"md_msdstarttime",md_msdstartTime,"choose which step that msd be calculated");
+	OUTP(ofs,"md_mdtype",mdp.mdtype,"choose ensemble");
+	OUTP(ofs,"md_dt",mdp.dt,"time step");
+	OUTP(ofs,"mnhc",mdp.MNHC,"number of Nose-Hoover chains");
+	OUTP(ofs,"md_qmass",mdp.Qmass,"mass of thermostat");
+	OUTP(ofs,"md_tfirst",mdp.tfirst,"temperature first");
+	OUTP(ofs,"md_tlast",mdp.tlast,"temperature last");
+	OUTP(ofs,"md_dumpmdfred",mdp.recordFreq,"The period to dump MD information for monitoring and restarting MD");
+	OUTP(ofs,"md_mdoutpath",mdp.mdoutputpath,"output path of md");
+	OUTP(ofs,"md_rstmd",mdp.rstMD,"whether restart");
+	OUTP(ofs,"md_fixtemperature",mdp.fixTemperature,"period to change temperature");
+	OUTP(ofs,"md_ediff",mdp.ediff,"parameter for constraining total energy change");
+	OUTP(ofs,"md_ediffg",mdp.ediffg,"parameter for constraining max force change");
+	OUTP(ofs,"NVT_tau",mdp.NVT_tau,"parameter for adjust effect of thermostat");
+	OUTP(ofs,"NVT_control",mdp.NVT_control,"choose which thermostat used in NVT ensemble");
 
 	ofs << "\n#Parameters (11.Efield)" << endl;
 	OUTP(ofs,"efield",efield,"add electric field");

ABACUS.develop/source/src_ions/variable_cell.cpp

@@ -58,7 +58,7 @@ void Variable_Cell::final_calculation_after_vc(void)
     cout<<" -----------------------------------------------------------------"<<endl;
 
     cout<<"\n -----------------------------------------------------------------"<<endl;
-    cout<<" The structure has been fully relaxed, and the following is a scf"<<endl;
+    cout<<" The structure has been fully relaxed or MD finished, and the following is a scf"<<endl;
     cout<<" calculation at the final structure. The fft grids and G-vectors "<<endl;
     cout<<" are recalculated for the final relaxed unit cell."<<endl;
     cout<<" -----------------------------------------------------------------"<<endl;

ABACUS.develop/source/src_lcao/ELEC_evolve.cpp

@@ -217,8 +217,8 @@ void ELEC_evolve::using_LAPACK_complex(const int &ik, complex<double>** c, compl
         double delta_t;
 //      delta_t = 0.2;	//identity: fs;
         ComplexMatrix Numerator(NLOCAL,NLOCAL);
-        Numerator = Idmat - 0.5*INPUT.md_dt*41.34*Denominator;
-        Denominator = Idmat + 0.5*INPUT.md_dt*41.34*Denominator;
+        Numerator = Idmat - 0.5*INPUT.mdp.dt*41.34*Denominator;
+        Denominator = Idmat + 0.5*INPUT.mdp.dt*41.34*Denominator;
 
 	int info;
         int lwork=3*NLOCAL-1; //tmp
@@ -368,8 +368,8 @@ void ELEC_evolve::using_LAPACK_complex_2(const int &ik, complex<double>** c, com
         }
         double delta_t;
         ComplexMatrix Numerator(NLOCAL,NLOCAL);
-        Numerator = Idmat - 0.5*INPUT.md_dt*41.34*Denominator;
-        Denominator = Idmat + 0.5*INPUT.md_dt*41.34*Denominator;
+        Numerator = Idmat - 0.5*INPUT.mdp.dt*41.34*Denominator;
+        Denominator = Idmat + 0.5*INPUT.mdp.dt*41.34*Denominator;
 
 	int info;
         int lwork=3*NLOCAL-1; //tmp

ABACUS.develop/source/src_lcao/FORCE_STRESS.h

@@ -16,6 +16,7 @@ class Force_Stress_LCAO
 	friend void Input_Conv::Convert();
 	friend class Update_input;
 	friend class LOOP_ions;
+	friend class MD_func;
 
 	public :
 

ABACUS.develop/source/src_lcao/LCAO_evolve.cpp

@@ -150,8 +150,8 @@ void Evolve_LCAO_Matrix::using_LAPACK_complex(const int &ik, complex<double>** c
         double delta_t;
 //      delta_t = 0.2;	//identity: fs;
         ComplexMatrix Numerator(NLOCAL,NLOCAL);
-        Numerator = Idmat - 0.5*INPUT.md_dt*41.34*Denominator;
-        Denominator = Idmat + 0.5*INPUT.md_dt*41.34*Denominator;
+        Numerator = Idmat - 0.5*INPUT.mdp.dt*41.34*Denominator;
+        Denominator = Idmat + 0.5*INPUT.mdp.dt*41.34*Denominator;
 
 	int info;
         int lwork=3*NLOCAL-1; //tmp