diff --git a/ABACUS.develop/source/Makefile.Objects b/ABACUS.develop/source/Makefile.Objects
index 27e60059ce..d241353e70 100644
--- a/ABACUS.develop/source/Makefile.Objects
+++ b/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 \
diff --git a/ABACUS.develop/source/input.cpp b/ABACUS.develop/source/input.cpp
index 25fa136250..a56a80ca24 100644
--- a/ABACUS.develop/source/input.cpp
+++ b/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();
             }
diff --git a/ABACUS.develop/source/input.h b/ABACUS.develop/source/input.h
index 3385c146ac..0250d12f94 100644
--- a/ABACUS.develop/source/input.h
+++ b/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
diff --git a/ABACUS.develop/source/src_io/print_info.cpp b/ABACUS.develop/source/src_io/print_info.cpp
index 515cca97db..993816b4e3 100644
--- a/ABACUS.develop/source/src_io/print_info.cpp
+++ b/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;
 
diff --git a/ABACUS.develop/source/src_io/write_input.cpp b/ABACUS.develop/source/src_io/write_input.cpp
index 5db6d49b03..877a6cd7de 100644
--- a/ABACUS.develop/source/src_io/write_input.cpp
+++ b/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");
diff --git a/ABACUS.develop/source/src_ions/variable_cell.cpp b/ABACUS.develop/source/src_ions/variable_cell.cpp
index 009e700837..7ba80704cc 100644
--- a/ABACUS.develop/source/src_ions/variable_cell.cpp
+++ b/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;
diff --git a/ABACUS.develop/source/src_lcao/ELEC_evolve.cpp b/ABACUS.develop/source/src_lcao/ELEC_evolve.cpp
index af59b42099..0db138bfca 100644
--- a/ABACUS.develop/source/src_lcao/ELEC_evolve.cpp
+++ b/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
diff --git a/ABACUS.develop/source/src_lcao/FORCE_STRESS.h b/ABACUS.develop/source/src_lcao/FORCE_STRESS.h
index d27cb86bcb..acfe68c784 100644
--- a/ABACUS.develop/source/src_lcao/FORCE_STRESS.h
+++ b/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 :
 	
diff --git a/ABACUS.develop/source/src_lcao/LCAO_evolve.cpp b/ABACUS.develop/source/src_lcao/LCAO_evolve.cpp
index 21b0971822..286095aaa8 100644
--- a/ABACUS.develop/source/src_lcao/LCAO_evolve.cpp
+++ b/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
diff --git a/ABACUS.develop/source/src_lcao/run_md.cpp b/ABACUS.develop/source/src_lcao/run_md.cpp
index 4fee78666d..28bcae3e9c 100644
--- a/ABACUS.develop/source/src_lcao/run_md.cpp
+++ b/ABACUS.develop/source/src_lcao/run_md.cpp
@@ -13,7 +13,6 @@
 #include "ELEC_scf.h"
 #include "src_global/sltk_atom_arrange.h"
 #include "src_pw/vdwd2.h"
-#include "src_pw/vdwd3.h"
 
 Run_MD::Run_MD()
 {}
@@ -88,20 +87,8 @@ void Run_MD::opt_ions(void)
         LCM.allocate();
     }
 
-	mdtype= INPUT.md_mdtype;//control md
-
-	// NVT ensemble
-	if(mdtype==1||mdtype==2) 
-	{
-		MDNVT.md_allocate();
-		MDNVT.initMD();
-	}
-	// NVE ensemble
-	else if(mdtype==0)
-	{
-		MDNVE.md_allocate();
-		MDNVE.initMD();
-	}
+    MD_basic mdb(INPUT.mdp, ucell);
+    int mdtype = INPUT.mdp.mdtype;
 
     this->istep = 1;
     int force_step = 1;
@@ -115,10 +102,7 @@ void Run_MD::opt_ions(void)
         if(OUT_LEVEL=="ie" || OUT_LEVEL=="m")
         {
 			cout << " ---------------------------------------------------------" << endl;
-			if(mdtype==1||mdtype==2)
-			{
-				cout<<" Molecular Dynamics (NVT) STEP "<< MDNVT.step_rst + istep<<endl;
-			}
+			cout<<" Molecular Dynamics STEP "<< mdb.getRealStep()<<endl;
 			cout << " ---------------------------------------------------------" << endl;
         }
 
@@ -133,20 +117,20 @@ void Run_MD::opt_ions(void)
 
 		if(mdtype==1||mdtype==2)   
 		{
-			MDNVT.runnvt(istep);
+			mdb.runNVT(istep);
 		}
 		else if(mdtype==0)  
 		{
-			MDNVE.runNVE(istep);
-		}
-		if((mdtype==1||mdtype==2)&&(istep==NSTEP||stop))
-		{
-			MDNVT.md_release();
-		}
-		if(mdtype==0&&(istep==NSTEP||stop))
-		{
-			MDNVE.md_release();
+			mdb.runNVE(istep);
 		}
+        else if(mdtype==-1)
+        {
+            stop = mdb.runFIRE(istep);
+        }
+        else
+        {
+            WARNING_QUIT("opt_ions", "mdtype should be -1~2!");
+        }
 
         if(pot.out_potential == 2)
         {
@@ -203,6 +187,8 @@ void Run_MD::opt_ions(void)
         ++istep;
     }
 
+    if(istep>1) final_scf();
+
 	// mohan update 2021-02-10
     hm.orb_con.clear_after_ions(UOT, ORB, INPUT.out_descriptor);
 
@@ -210,190 +196,6 @@ void Run_MD::opt_ions(void)
     return;
 }
 
-//bool Run_MD::force_stress(void)
-bool Run_MD::force_stress(const int &istep, int &force_step, int &stress_step)
-{
-    TITLE("Run_MD","force_stress");
-    if(!FORCE && !STRESS)
-    {
-        return 1;
-    }
-    timer::tick("Run_MD","force_stress",'D');
-	matrix fcs;
-	matrix scs;
-	Force_Stress_LCAO FSL;
-	FSL.allocate (); 
-	FSL.getForceStress(FORCE, STRESS, TEST_FORCE, TEST_STRESS, fcs, scs);
-	//--------------------------------------------------
-	// only forces are needed, no stresses are needed
-	//--------------------------------------------------
-    if(FORCE && !STRESS)
-    {
-
-#ifdef __MPI //2015-10-01, xiaohui
-        atom_arrange::delete_vector( SEARCH_RADIUS );
-#endif //2015-10-01, xiaohui
-
-        if(CALCULATION=="relax") 
-        {
-            IMM.cal_movement(istep, istep, fcs, en.etot);
-
-            if(IMM.get_converged() || (istep==NSTEP))
-            {
-                return 1; // 1 means converged
-            }
-            else // ions are not converged
-            {
-                CE.update_istep(istep);
-                CE.extrapolate_charge();
-
-                if(pot.extra_pot=="dm")
-                {
-                }
-                else
-                {
-                    pot.init_pot( istep, pw.strucFac );
-                }
-            }
-            return 0;
-        }
-        else
-        {
-            return 1;
-        }
-
-        // mohan update 2013-04-11
-        // setup the structure factor
-        // and do the density extraploation.
-        // for both ionic iteratoin and
-        // force calculations.
-
-        //xiaohui modify 2014-08-09
-        //pw.setup_structure_factor();
-
-        // charge extrapolation if istep>0.
-        //xiaohui modify 2014-08-09
-        //CE.extrapolate_charge();
-
-/*xiaohui modify 2014-08-09
-        if(pot.extra_pot==4)
-        {
-            // done after grid technique.
-        }
-        else
-        {
-            pot.init_pot( istep );
-        }
-xiaohui modify 2014-08-09*/
-    }
-
-    static bool converged_force = false;
-    static bool converged_stress = false;
-
-    if(!FORCE&&STRESS)
-    {
-
-#ifdef __MPI
-		atom_arrange::delete_vector( SEARCH_RADIUS );
-#endif
-		if(CALCULATION=="cell-relax")
-		{
-           	LCM.cal_lattice_change(stress_step, scs, en.etot);
-           	converged_stress = LCM.get_converged();
-           	if(converged_stress)
-           	{
-               	return 1;
-           	}
-           	else
-           	{
-               	Variable_Cell::init_after_vc();
-               	pot.init_pot(stress_step, pw.strucFac);
-
-               	++stress_step;
-               	return 0;
-           	}
-		}
-        else
-        {
-            return 1;
-        }
-	}
-
-    if(FORCE&&STRESS)
-    {
-
-//#ifdef __MPI
-        atom_arrange::delete_vector( SEARCH_RADIUS );
-//#endif
-        
-        //if(CALCULATION=="relax") IMM.cal_movement(istep, FL.fcs, en.etot);
-        if(CALCULATION=="relax" || CALCULATION=="cell-relax")
-        {
-            IMM.cal_movement(istep, force_step, fcs, en.etot);
-
-            if(IMM.get_converged())
-            {
-                force_step = 1;
-
-
-			    if(CALCULATION=="cell-relax")
-			    {
-            	    LCM.cal_lattice_change(stress_step, scs, en.etot);
-            	    converged_stress = LCM.get_converged();
-            	    if(converged_stress)
-            	    {
-                	    return 1;
-            	    }
-            	    else
-            	    {
-                	    Variable_Cell::init_after_vc();
-                	    pot.init_pot(stress_step, pw.strucFac);
-
-                	    ++stress_step;
-                	    return 0;
-                    }
-                }
-                else
-                {
-                    return 1;
-                }
-
-#ifdef __MPI
-            //atom_arrange::delete_vector( SEARCH_RADIUS );
-#endif
-            }
-            else
-            {
-#ifdef __MPI
-            //atom_arrange::delete_vector( SEARCH_RADIUS );
-#endif
-                //CE.istep = istep;
-                CE.update_istep(force_step);
-                CE.extrapolate_charge();
-
-                if(pot.extra_pot=="dm")//xiaohui modify 2015-02-01
-                {
-                    // done after grid technique.
-                }
-                else
-                {
-                    pot.init_pot( istep, pw.strucFac );
-                }
-                ++force_step;
-                return 0;
-            }
-        }
-        else
-        {
-            return 1;
-        }
-    }
-
-    return 0;
-
-    timer::tick("Run_MD","force_stress",'D');
-}
-
 void Run_MD::final_scf(void)
 {
     TITLE("Run_MD","final_scf"); 
@@ -453,18 +255,15 @@ void Run_MD::final_scf(void)
         Vdwd3 vdwd3(ucell,vdwd3_para);
         vdwd3.cal_energy();
         en.evdw = vdwd3.get_energy();
-    }											  
+    }												  
     
 	ELEC_scf es;
 	es.scf(0);
 
-    if(CALCULATION=="scf" || CALCULATION=="relax")
-    {
-        ofs_running << "\n\n --------------------------------------------" << endl;
-        ofs_running << setprecision(16);
-        ofs_running << " !FINAL_ETOT_IS " << en.etot * Ry_to_eV << " eV" << endl; 
-        ofs_running << " --------------------------------------------\n\n" << endl;
-    }
+    ofs_running << "\n\n --------------------------------------------" << endl;
+    ofs_running << setprecision(16);
+    ofs_running << " !FINAL_ETOT_IS " << en.etot * Ry_to_eV << " eV" << endl; 
+    ofs_running << " --------------------------------------------\n\n" << endl;
 
     return;
 }
diff --git a/ABACUS.develop/source/src_lcao/run_md.h b/ABACUS.develop/source/src_lcao/run_md.h
index 96e7b42dc8..93c17b5c39 100644
--- a/ABACUS.develop/source/src_lcao/run_md.h
+++ b/ABACUS.develop/source/src_lcao/run_md.h
@@ -4,9 +4,7 @@
 #include "LOOP_elec.h"
 #include "../src_ions/ions_move_methods.h"
 #include "../src_pw/charge_extra.h"
-#include "../src_pw/md.h"
-#include "../src_pw/mdNVT.h"
-#include "../src_pw/mdNVE.h"
+#include "../src_pw/MD_basic.h"
 #include "../src_ions/lattice_change_methods.h"
 
 class Run_MD 
@@ -22,10 +20,6 @@ class Run_MD
 	void opt_cell(void);
 	void opt_ions(void);
 
-	//2014-06-06, xiaohui
-	mdnvt MDNVT ;
-	mdNVE MDNVE ;
-
 	private:
 
 	Ions_Move_Methods IMM;
@@ -33,16 +27,11 @@ class Run_MD
 	//bool force_stress(void);
 	Lattice_Change_Methods LCM;
 
-	bool force_stress(const int &istep, int &force_step, int &stress_step);
-
 	int istep;
 
 	// electron charge density extropolation method
 	Charge_Extra CE;
 
-	//choose md ensemble, zheng daye
-	int mdtype;
-
 	void final_scf(void);
 
 };
diff --git a/ABACUS.develop/source/src_pw/MD_basic.cpp b/ABACUS.develop/source/src_pw/MD_basic.cpp
new file mode 100644
index 0000000000..01a3ee009d
--- /dev/null
+++ b/ABACUS.develop/source/src_pw/MD_basic.cpp
@@ -0,0 +1,606 @@
+#include "src_pw/MD_basic.h"
+#include "src_pw/global.h"
+
+//define in MD_basic.h
+//class MD_basic
+
+const double MD_basic::fundamentalTime = 2.418884326505e-17;
+
+MD_basic::MD_basic(MD_parameters& MD_para_in, UnitCell_pseudo &unit_in):
+    mdp(MD_para_in),
+    ucell(unit_in)
+{	
+	mdp.dt=mdp.dt/fundamentalTime/1e15;
+	temperature_=mdp.tfirst/3.1577464e5;
+
+
+	//other parameter default
+	outputstressperiod_ = 1 ;// The period to output stress
+//	ucell.nat=ucell.nat;
+//	ntype=ucell.ntype;
+	step_rst_=0;
+    step_=0;
+//	ucell.latvec=ucell.latvec;
+
+    vel=new Vector3<double>[ucell.nat];
+	cart_change=new Vector3<double>[ucell.nat];
+	tauDirectChange=new Vector3<double>[ucell.nat];
+	allmass=new double[ucell.nat];
+	ionmbl=new Vector3<int>[ucell.nat];
+	force=new Vector3<double>[ucell.nat];
+
+	energy_=en.etot/2;
+
+	//MD starting setup
+	if(!mdp.rstMD){
+		nfrozen_ = mdf.getMassMbl(ucell, allmass, ionmbl);
+		//mdf.gettauDirectChange(ucell, tauDirectChange);
+		//A fresh new MD: Do not restart MD
+		mdf.InitVelocity(ucell.nat, temperature_, fundamentalTime, allmass, vel) ;
+		// Initialize thermostat, and barostat
+		
+	}
+	else{
+		nfrozen_ = mdf.getMassMbl(ucell, allmass, ionmbl);
+		//mdf.gettauDirectChange(ucell, tauDirectChange);
+		mdf.InitVelocity(ucell.nat, temperature_, fundamentalTime, allmass, vel) ;
+		if(!mdf.RestartMD(ucell.nat, vel, step_rst_)){
+			cout<<"error in restart MD!"<<endl;
+			return;
+		}
+	}
+	if(mdp.NVT_tau<1e-10)
+	{
+		if(mdp.NVT_control == 1) mdp.NVT_tau = 20 * mdp.dt; //NHC
+		else if(mdp.NVT_control == 2) mdp.NVT_tau = 20 * mdp.dt; //LGV
+		else if(mdp.NVT_control == 3) mdp.NVT_tau = 20 * mdp.dt; //ADS
+		else WARNING_QUIT("md:InitMD", "please choose available reservoir!!!");
+	}
+	else
+	{
+		mdp.NVT_tau /= fundamentalTime*1e15;
+	}
+	//tau = 1.0/(NVT_tau*fundamentalTime*1e15);
+	// Q=KbT*tau**2
+    mdp.Qmass=mdp.Qmass/6.02/9.109*1e5;
+	if(mdp.Qmass<1e-10)//Qmass=dt*fundamentalTime*1e15/6.02/9.109*1e5;
+	mdp.Qmass = pow(mdp.NVT_tau,2)*temperature_;///beta;
+	//NHC setup
+	if(mdp.mdtype==1 && mdp.NVT_control == 1)
+	{
+		mdt.init_NHC(
+        mdp.MNHC, 
+        mdp.Qmass, 
+        mdp.NVT_tau, 
+        mdp.dt,
+        mdp.NVT_control, 
+        ofs_running, 
+        ucell.nat,
+        temperature_,
+        vel,
+        allmass
+        );
+        if(mdp.rstMD) mdt.NHC_restart();
+	}
+}
+
+MD_basic::~MD_basic()
+{
+    delete []force;
+    delete []ionmbl;
+    delete []cart_change;
+    delete []vel;
+    delete []tauDirectChange;
+    delete []allmass;
+}
+
+void MD_basic::runNVT(int step1){
+//------------------------------------------------------------------------------
+// DESCRIPTION:
+// Molecular dynamics calculation with fixed Volume and slight fluctuated temperature
+// Using thermostat : 1, Nose-Hoover Chains; 2, Langevin; 3, Anderson
+// Normal Nose-Hoover thermostat method is retained for test.
+//------------------------------------------------------------------------------
+
+	TITLE("MD_basic","runnvt");
+	timer::tick("MD_basic","runnvt",'C');
+	step_=step1+step_rst_;
+	//the real MD step
+	
+	//change target temperature
+//	if (step_!=1)mdf.ReadNewTemp( step_ );
+	
+	cout << " ---------------------------------------------------------" << endl;
+	
+	cout<< " Target temperature  : "<< temperature_/K_BOLTZMAN_AU<< " (K)"<<endl;
+	
+	if(step_==1||step_%mdp.fixTemperature==1)
+    {
+		for(int k=0;k<ucell.nat;k++)
+        {
+			if(ionmbl[k].x==0)vel[k].x=0;
+			if(ionmbl[k].y==0)vel[k].y=0;
+			if(ionmbl[k].z==0)vel[k].z=0;
+		}
+		mdf.scalevel(ucell.nat, nfrozen_, temperature_, vel, allmass);
+	}
+
+	// get the kinetic energy
+	double twiceKE = mdf.GetAtomKE(ucell.nat, vel, allmass);
+	twiceKE = twiceKE * 2;
+	
+	//Set up forces and stress
+    if(INPUT.basis_type == "lcao")
+    {
+	    mdf.callInteraction_LCAO(ucell.nat, force, stress);
+    }
+    else if(INPUT.basis_type == "pw")
+    {
+        mdf.callInteraction_PW(ucell.nat, force, stress);
+    }
+	//print total stress + stress_MD
+	if(STRESS)
+	{
+		outStressMD(stress, twiceKE);
+	}
+
+	for(int k=0;k<ucell.nat;k++){
+		if(ionmbl[k].x==0)force[k].x=0;
+		if(ionmbl[k].y==0)force[k].y=0;
+		if(ionmbl[k].z==0)force[k].z=0;
+	}
+	double maxForce = mdf.MAXVALF(ucell.nat, force);
+	energy_=en.etot/2;
+    double hamiltonian;
+	//----------------------------------------------
+	// big loop
+	//-----------------------------------------------
+	cout<<" "<<std::left<<setw(12)<<"MD_STEP"<<std::left<<setw(12)<< "SystemE"<<std::left<<setw(12)<< "Conserved"<<std::left<<setw(12)<< "DeltaE"<<std::left<<setw(12)<< "Temperature"<<endl;
+	
+	cout<<" "<<std::left<<setw(12)<<step_<<std::left<<setw(12)<< energy_<<std::left<<setw(12)<< hamiltonian<<std::left<<setw(12)<< energy_-oldEtot_<<std::left<<setw(12)<<twiceKE/(3*double(ucell.nat-nfrozen_))/K_BOLTZMAN_AU<<endl;
+	
+	oldEtot_=energy_;
+
+	if (!MY_RANK)
+    {
+        ofs_running<<" --------------------------------------------------"<<endl;
+        ofs_running<<" Molecular Dynamics (NVT) STEP "<< step_<<endl;
+        ofs_running<< " --------------------------------------------------"<<endl;
+	}
+	
+	
+	// Calculate the Mean-Square-Displacement.
+	if(step_==1&&mdp.rstMD==0)
+    { 
+        //Note: side scheme position before
+        //Now turn to middle scheme.
+	
+	    // New vel obtained (Verlet-like step)
+        if(mdp.NVT_control == 1)
+        {
+            hamiltonian = mdt.NHChamiltonian(
+                twiceKE/2, 
+                energy_,
+                temperature_,
+                nfrozen_);
+        }
+        else hamiltonian = mdf.Conserved(twiceKE/2, energy_, ucell.nat-nfrozen_);
+	    this->update_half_velocity();
+	}
+	else 
+    {
+		this->update_half_velocity();
+		
+		twiceKE=mdf.GetAtomKE(ucell.nat, vel, allmass);
+		twiceKE = 2 * twiceKE;
+		if(mdp.NVT_control==1)
+        {
+            hamiltonian = mdt.NHChamiltonian(
+                twiceKE/2, 
+                energy_,
+                temperature_,
+                nfrozen_);
+        } 
+		else hamiltonian = mdf.Conserved(twiceKE/2, energy_, ucell.nat-nfrozen_);
+        //Note: side scheme position before
+        //Now turn to middle scheme. 
+		this->update_half_velocity();
+	}
+
+	// Update the Non-Wrapped cartesion coordinates
+    if(mdp.mdtype==2) mdf.scalevel(ucell.nat, nfrozen_, temperature_, vel, allmass);//choose velocity scaling method
+
+    update_half_direct(1);
+	
+	mdt.Integrator(mdp.NVT_control, temperature_, vel, allmass);//thermostat interact with velocity
+    twiceKE=mdf.GetAtomKE(ucell.nat, vel, allmass);
+    twiceKE = 2 * twiceKE;
+	
+	update_half_direct(0);
+
+	// Calculate the maximal velocities.
+	// The step in fractional coordinates
+	double maxStep = 0;
+    for(int  ii=0;ii<ucell.nat;ii++)
+    { 
+		if((pow(vel[ii].x,2)+pow(vel[ii].y,2)+pow(vel[ii].z,2))>maxStep)
+        {
+		    maxStep = pow(vel[ii].x,2)+pow(vel[ii].y,2)+pow(vel[ii].z,2);
+        }
+    }
+	getTaudUpdate();
+
+	//save the atom position change to DFT module
+	save_output_position();
+	maxStep = sqrt(maxStep)*mdp.dt;
+
+    if (!MY_RANK){
+        ofs_running<<setw(15)<<"maxForce=     "<<setw(15)<<"maxstep=      "<<setw(15)<<"step=     "<<endl;
+        ofs_running<<setw(15)<<maxForce<<setw(15)<<maxStep<<setw(15)<<step_<<endl;
+        ofs_running<<step_<<" "<< energy_<<" "<< hamiltonian<<" "<< energy_-oldEtot_<<" "<< twiceKE/(3*double(ucell.nat-nfrozen_))/K_BOLTZMAN_AU<<endl;
+    }
+    oldEtot_=energy_;
+    //output basic restart info
+    mdf.mdRestartOut(step_, mdp.recordFreq, ucell.nat, vel);
+    //output NHC thermo restart info
+    if(mdp.NVT_control==1) mdt.NHC_info_out(step_, mdp.recordFreq, ucell.nat);
+
+//2015-09-25, xiaohui
+#ifdef __MPI
+    if(step_%100==0)
+    {
+        MPI_Bcast(vel,ucell.nat*3,MPI_DOUBLE,0,MPI_COMM_WORLD);
+    //    MPI_Bcast(&xLogS,1,MPI_DOUBLE,0,MPI_COMM_WORLD);
+    //    MPI_Bcast(&vLogS,1,MPI_DOUBLE,0,MPI_COMM_WORLD);
+    }
+#endif
+    timer::tick("MD_basic","runnvt",'C');
+    return;
+}
+
+void MD_basic::runNVE(int step1){
+//-------------------------------------------------------------------------------
+// Adiabatic ensemble 
+// Molecular dynamics calculation with Verlet algorithm
+//-------------------------------------------------------------------------------
+
+    TITLE("MD_basic","runNVE");
+    timer::tick("MD_basic","runNVE",'C');
+    step_=step1+step_rst_;
+    if(step_==1)
+    {
+        cout<<" (NVE) Start."<<endl;
+    }
+    else 
+    {
+        cout<<" (NVE) step: "<<step_<<endl;
+    }
+
+  // Set up extra output to ion optimizer / MD header
+  cout<<"Time interval       : "<< mdp.dt*fundamentalTime*1e15<< " (fs)"<<endl;
+  cout<<"Target temperature  : "<< temperature_/K_BOLTZMAN_AU<< " (K)"<<endl;
+
+    if(step_==1){
+        for(int k=0;k<ucell.nat;k++)
+        {
+            if(ionmbl[k].x==0)vel[k].x=0;
+            if(ionmbl[k].y==0)vel[k].y=0;
+            if(ionmbl[k].z==0)vel[k].z=0;
+        }
+        mdf.scalevel(ucell.nat, nfrozen_, temperature_, vel, allmass);
+    }
+    double twiceKE=mdf.GetAtomKE(ucell.nat, vel, allmass);
+    twiceKE = twiceKE * 2;
+
+    double tempNow = twiceKE/(3*double(ucell.nat-nfrozen_))/K_BOLTZMAN_AU;
+    cout<<" start temperature = "<< tempNow/K_BOLTZMAN_AU<< " (k)"<<endl;
+
+    // Set up forces 
+    mdf.callInteraction_LCAO(ucell.nat, force, stress);
+	if(STRESS)
+	{
+		outStressMD(stress, twiceKE);
+	}
+  
+    for(int k=0;k<ucell.nat;k++)
+    {
+        if(ionmbl[k].x==0)force[k].x=0;
+        if(ionmbl[k].y==0)force[k].y=0;
+        if(ionmbl[k].z==0)force[k].z=0;
+    }
+
+    //cout<<"begin maxForce"<<endl;
+    double maxForce = mdf.MAXVALF(ucell.nat, force);
+    cout<<"maxForce: "<<sqrt(maxForce)<<endl; 
+    energy_=en.etot/2;
+    double conservedE = mdf.Conserved(twiceKE/2, energy_, ucell.nat-nfrozen_);
+
+    cout<< "NVE_STEP "<<" "<<"SystemEnergy"<<" "<< "Conserved"<<" "<< "DeltaE"<<" "<< "Temperature"<<endl;
+
+    // Output the message to the screen.
+    cout<<step_<<" "<< energy_<<" "<< conservedE<<" "
+        << energy_-oldEtot_<<" "<< twiceKE/(3*double(ucell.nat-nfrozen_))/K_BOLTZMAN_AU<<endl;
+
+    cout<<"------------------------------------------------------------------------------"<<endl;
+    cout<< "MD(NVE) STEP "<< step_<<endl;
+    cout<<"------------------------------------------------------------------------------"<<endl;
+    if (!MY_RANK){
+        ofs_running<<"------------------------------------------------------------------------------"<<endl;
+        ofs_running<< "MD(NVE) STEP "<< step_<<endl;
+        ofs_running<<"------------------------------------------------------------------------------"<<endl;
+    }
+  
+    // (1) 1st step of Verlet-Velocity
+    // New velocity are obtained
+	this->update_half_velocity();
+    // (2) 2nd step of Verlet-Velocity 
+    // Update the Non-Wrapped cartesion coordinate
+    twiceKE = mdf.GetAtomKE(ucell.nat, vel, allmass);
+    twiceKE = 2 * twiceKE;
+    if(step_!=1||mdp.rstMD==1)this->update_half_velocity();
+    update_half_direct(1);
+    update_half_direct(0);
+    // Calculate the maximal velocities.
+    // The step in fractional coordinates
+    double maxStep = 0;
+    for(int  i = 0;i< ucell.nat;i++)
+    {
+        if((pow(vel[i].x,2.0)+pow(vel[i].y,2.0)+pow(vel[i].z,2.0))>maxStep)
+        {
+            maxStep = pow(vel[i].x,2.0)+pow(vel[i].y,2.0)+pow(vel[i].z,2.0);
+        }
+    }
+    getTaudUpdate();
+    save_output_position();
+    maxStep = sqrt(maxStep)*mdp.dt;
+
+
+    // calculate the conserved quantity during MD 
+    double hamiltonian = mdf.Conserved(twiceKE/2, energy_, ucell.nat-nfrozen_);
+
+    cout<< setprecision (9)<<hamiltonian<<" "<< setprecision (9)<<twiceKE/2<<endl;
+
+    
+    // Output the message to the screen.
+    if (!MY_RANK)
+    { 
+        ofs_running<<step_<<" "<< energy_<<" "<< hamiltonian<<" "<< energy_-oldEtot_<<" "<< twiceKE/(3*double(ucell.nat-nfrozen_))/K_BOLTZMAN_AU<<endl;
+        ofs_running<<step_<<" "<< energy_<<" "<< conservedE<<" "<< energy_-oldEtot_<<" "<< twiceKE/(3*double(ucell.nat-nfrozen_))/K_BOLTZMAN_AU<<endl;
+    }
+    oldEtot_=energy_;
+    mdf.mdRestartOut(step_, mdp.recordFreq, ucell.nat, vel);
+
+    cout<<"(NVE): this step finished."<<endl;
+    if (!MY_RANK)
+    {
+        ofs_running<<"(NVE): this step finished."<<endl;
+    }
+//2015-09-25, xiaohui
+#ifdef __MPI
+    if(step_%100==0)
+    {
+        MPI_Bcast(vel,ucell.nat*3,MPI_DOUBLE,0,MPI_COMM_WORLD);
+ //       MPI_Bcast(&xLogS,1,MPI_DOUBLE,0,MPI_COMM_WORLD);
+ //       MPI_Bcast(&vLogS,1,MPI_DOUBLE,0,MPI_COMM_WORLD);
+    }
+#endif
+
+    timer::tick("MD_basic","runNVE",'C');
+    return;
+}
+
+bool MD_basic::runFIRE(int step1){
+//-------------------------------------------------------------------------------
+// REFERENCES:
+//   
+    TITLE("MD_basic","runFIRE");
+    timer::tick("MD_basic","runFIRE",'C');
+    step_ = step1;
+    if(step_==1)
+    {
+        cout<<" FIRE Start."<<endl;
+    }
+    else 
+    {
+        cout<<" FIRE step: "<<step_<<endl;
+    }
+  
+    cout<<" FIRE Fronzen Atom Number : "<<nfrozen_<<endl;
+    
+
+    // Set up extra output to ion optimizer / MD header
+    cout<<"Time interval       : "<< mdp.dt*fundamentalTime*1e15<< " (fs)"<<endl;
+
+    if(step_==1)
+    {
+        for(int k=0;k<ucell.nat;k++)
+        {
+            if(ionmbl[k].x==0)vel[k].x=0;
+            if(ionmbl[k].y==0)vel[k].y=0;
+            if(ionmbl[k].z==0)vel[k].z=0;
+        }
+    //      if(nfrozen_==0)RemoveMovementOfCenterOfMass();
+        mdf.scalevel(ucell.nat, nfrozen_, temperature_, vel, allmass);
+    }
+    double twiceKE=mdf.GetAtomKE(ucell.nat, vel, allmass);
+    twiceKE = twiceKE * 2;
+
+    double tempNow = twiceKE/(3*double(ucell.nat-nfrozen_))/K_BOLTZMAN_AU;
+    cout<<" start temperature = "<< tempNow/K_BOLTZMAN_AU<< " (k)"<<endl;
+
+
+
+
+
+    // Set up forces 
+    mdf.callInteraction_LCAO(ucell.nat, force, stress);
+    
+    for(int k=0;k<ucell.nat;k++)
+    {
+        if(ionmbl[k].x==0)force[k].x=0;
+        if(ionmbl[k].y==0)force[k].y=0;
+        if(ionmbl[k].z==0)force[k].z=0;
+    }
+
+    //cout<<"begin maxForce"<<endl;
+    double maxForce = mdf.MAXVALF(ucell.nat, force);
+    cout<<"maxForce: "<<sqrt(maxForce)<<endl; 
+    energy_ = en.etot/2;
+    double conservedE = mdf.Conserved(twiceKE/2, energy_, ucell.nat-nfrozen_);
+
+    cout<<"------------------------------------------------------------------------------"<<endl;
+    cout<< "FIRE STEP "<< step_<<endl;
+    cout<<"------------------------------------------------------------------------------"<<endl;
+    if (!MY_RANK){
+         ofs_running<<"------------------------------------------------------------------------------"<<endl;
+         ofs_running<< "FIRE STEP "<< step_<<endl;
+         ofs_running<<"------------------------------------------------------------------------------"<<endl;
+    }
+  
+    // (1) 1st step of Verlet-Velocity
+    // New velocity are obtained
+	this->update_half_velocity();
+
+    // (2) 2nd step of Verlet-Velocity 
+    // Update the Non-Wrapped cartesion coordinate
+    twiceKE = mdf.GetAtomKE(ucell.nat, vel, allmass);
+    twiceKE = 2 * twiceKE;
+    if(step_!=1)this->update_half_velocity();
+
+	double largest_grad_FIRE = 0.0;
+	for(int i=0;i<ucell.nat;i++)
+    {
+        if(largest_grad_FIRE < abs(force[i].x))
+        {
+            largest_grad_FIRE= abs(force[i].x);
+        } 
+        if(largest_grad_FIRE < abs(force[i].y))
+        {
+            largest_grad_FIRE= abs(force[i].y);
+        } 
+        if(largest_grad_FIRE < abs(force[i].z))
+        {
+            largest_grad_FIRE= abs(force[i].z);
+        } 
+	}
+	
+	ofs_running << " LARGEST GRAD (eV/A)  : " << largest_grad_FIRE * Ry_to_eV / 0.529177 << endl;
+
+	if(largest_grad_FIRE*Ry_to_eV/0.529177 < 0.01)
+    {
+	//"convergency reach"
+		//cout <<"CONVERGENCY REACH of FIRE in the "<<step <<" steps " <<endl;
+		timer::tick("MD_basic","runFIRE",'C');
+		return true;
+	}
+
+	fire.check_FIRE(ucell.nat, force, mdp.dt, vel);
+    
+    update_half_direct(1);
+    update_half_direct(0);
+    // Calculate the maximal velocities.
+    // The step in fractional coordinates
+    double maxStep = 0;
+    for(int  i = 0;i< ucell.nat;i++)
+    {
+        if((pow(vel[i].x,2.0)+pow(vel[i].y,2.0)+pow(vel[i].z,2.0))>maxStep)
+        {
+            maxStep = pow(vel[i].x,2.0)+pow(vel[i].y,2.0)+pow(vel[i].z,2.0);
+        }
+    }
+    getTaudUpdate();
+
+    save_output_position();
+    maxStep = sqrt(maxStep)*mdp.dt;
+
+    double hamiltonian = mdf.Conserved(twiceKE/2, energy_, ucell.nat-nfrozen_);
+
+
+
+    
+    // Output the message to the screen.
+    oldEtot_=energy_;
+
+    cout<<"(NVE): this step finished."<<endl;
+    if (!MY_RANK){
+        ofs_running<<"(NVE): this step finished."<<endl;
+    }
+//2015-09-25, xiaohui
+#ifdef __MPI
+    if(step1%100==0){
+        MPI_Bcast(vel,ucell.nat*3,MPI_DOUBLE,0,MPI_COMM_WORLD);
+    }
+#endif
+
+    timer::tick("MD_basic","runFIRE",'C');
+    return false;
+}
+
+//update velocities of ions for half MD step
+void MD_basic::update_half_velocity()
+{
+    for(int  ii=0;ii<ucell.nat;ii++){ 
+        vel[ii] = vel[ii] + force[ii]/allmass[ii]*mdp.dt/2.0;
+    }
+}
+
+//update cartesian coordinate changes for half MD step, would turn to direct coordinate changes later
+void MD_basic::update_half_direct(const bool is_restart)
+{
+    if(is_restart)
+    {//first half
+        for(int  ii=0;ii<ucell.nat;ii++)
+        { 
+            cart_change[ii] = vel[ii]*mdp.dt/2.0/ucell.lat0;
+        }
+    }
+    else
+    {//second half
+        for(int  ii=0;ii<ucell.nat;ii++)
+        { 
+            cart_change[ii] += vel[ii]*mdp.dt/2.0/ucell.lat0;
+        }
+    }
+}
+
+//output Structure information for each MD step
+void MD_basic::save_output_position()
+{
+    string posOutName("md_pos_");
+	posOutName=mdf.intTurnTostring(step_,posOutName);
+	ucell.update_pos_taud(tauDirectChange);
+	mdf.printpos(posOutName, step_, mdp.recordFreq, ucell);
+}
+
+int MD_basic::getRealStep()
+{
+    return this->step_;
+}
+
+//output pressure of total MD system, P = tr(stress) + P_kin
+void MD_basic::outStressMD(const matrix& stress, const double& twiceKE)
+{
+    ofs_running<<"output Pressure for check!"<<endl;
+    double press;
+    for(int i=0;i<3;i++)
+        press += stress(i,i)/3;
+    press += twiceKE/3/ucell.omega; //output virtual press = 2/3 *Ek/V + sum(sigma[i][i])/3
+    double unit_transform = RYDBERG_SI / pow(BOHR_RADIUS_SI,3) * 1.0e-8 ;
+    ofs_running<<"Virtual Pressure is "<<press*unit_transform<<" Kbar "<<endl;
+}
+
+//turn cartesian coordinate changes to direct changes
+void MD_basic::getTaudUpdate()
+{
+    Vector3<double> fracStep;
+	for(int  ii=0;ii<ucell.nat;ii++)
+    { 
+		Mathzone::Cartesian_to_Direct(cart_change[ii].x,cart_change[ii].y,cart_change[ii].z,
+					ucell.latvec.e11,ucell.latvec.e12,ucell.latvec.e13,
+					ucell.latvec.e21,ucell.latvec.e22,ucell.latvec.e23,
+					ucell.latvec.e31,ucell.latvec.e32,ucell.latvec.e33,
+					fracStep.x,fracStep.y,fracStep.z);
+
+		tauDirectChange[ii] = fracStep;
+	}
+}
\ No newline at end of file
diff --git a/ABACUS.develop/source/src_pw/MD_basic.h b/ABACUS.develop/source/src_pw/MD_basic.h
new file mode 100644
index 0000000000..3646126b85
--- /dev/null
+++ b/ABACUS.develop/source/src_pw/MD_basic.h
@@ -0,0 +1,56 @@
+#ifndef MD_BASIC_H
+#define MD_BASIC_H
+
+#include "src_pw/MD_func.h"
+#include "src_pw/MD_thermo.h"
+#include "src_pw/MD_parameters.h"
+#include "src_pw/unitcell_pseudo.h"
+#include "src_pw/MD_fire.h"
+
+//using namespace std;
+class MD_basic
+{
+	public:	
+	MD_basic(MD_parameters& MD_para_in, UnitCell_pseudo &unit_in);
+	~MD_basic();
+	void runNVT(int step1);//NVT ensemble MD
+	void runNVE(int step1); //NVE ensemble MD
+	bool runFIRE(int step1); //relax method FIRE
+	int getRealStep();
+
+	private:
+    MD_func mdf;
+    MD_thermo mdt;
+	MD_parameters &mdp;
+	UnitCell_pseudo &ucell;
+	MD_fire fire;
+
+	double temperature_;
+	const static double fundamentalTime;
+	int outputstressperiod_;
+	int step_rst_;
+	int step_;
+	double energy_;
+	int nfrozen_;
+	double oldEtot_;
+
+	//allocable variaints
+	Vector3<double> *vel;//   velocity of each atom, unit is a.u.
+	double *allmass;     //mass of each atom
+	Vector3<double> *force;  //force of each atom
+	matrix stress;           //stress for this lattice
+	Vector3<int> *ionmbl;    //if frozen of each atom
+	Vector3<double> *tauDirectChange;           //change of dirac coord of atoms
+	Vector3<double> *cart_change;        //cartensian coord of atoms, *not* wrapped
+
+	//repete code
+	void update_half_velocity();
+	void update_half_direct(const bool is_restart);
+	void save_output_position();
+	void outStressMD(const matrix& stress, const double& twiceKE);
+	void getTaudUpdate();
+	
+
+};
+
+#endif
\ No newline at end of file
diff --git a/ABACUS.develop/source/src_pw/MD_fire.cpp b/ABACUS.develop/source/src_pw/MD_fire.cpp
new file mode 100644
index 0000000000..91a1c9f10c
--- /dev/null
+++ b/ABACUS.develop/source/src_pw/MD_fire.cpp
@@ -0,0 +1,68 @@
+#include "src_pw/MD_fire.h"
+
+MD_fire::MD_fire()
+{
+    dt_max=-1.0;
+    alpha_start=0.10;
+    alpha = alpha_start;
+
+    finc=1.1;
+    fdec=0.5;
+    f_alpha=0.99;
+    N_min=4;
+    negative_count=0;
+}
+
+void MD_fire::check_FIRE(
+    const int& numIon, 
+    const Vector3<double>const* force, 
+    double& deltaT, 
+    Vector3<double>* vel)
+{
+    if(dt_max<0) dt_max = deltaT * 2.5; //initial dt_max
+	double P(0.0);
+	double sumforce(0.0);
+	double normvel(0.0);
+	
+	for(int iatom =0;iatom<numIon;iatom++)
+    {
+		P += vel[iatom].x*force[iatom].x+vel[iatom].y*force[iatom].y+vel[iatom].z*force[iatom].z;
+		
+		sumforce += force[iatom].x*force[iatom].x+force[iatom].y*force[iatom].y+force[iatom].z*force[iatom].z;
+		
+		normvel += vel[iatom].x*vel[iatom].x +vel[iatom].y*vel[iatom].y +vel[iatom].z*vel[iatom].z ;
+	}
+	
+	sumforce = sqrt(sumforce);
+	normvel = sqrt(normvel);
+	
+	for(int iatom=0; iatom<numIon; iatom ++)
+    {
+		vel[iatom].x = (1.0-alpha)*vel[iatom].x+alpha*force[iatom].x/sumforce*normvel;
+		vel[iatom].y = (1.0-alpha)*vel[iatom].y+alpha*force[iatom].y/sumforce*normvel;
+		vel[iatom].z = (1.0-alpha)*vel[iatom].z+alpha*force[iatom].z/sumforce*normvel;
+	}
+	
+	if(P > 0 ){
+	    negative_count +=1 ;
+		if(negative_count >=N_min)
+        {
+			deltaT=min(deltaT*finc,dt_max);
+			alpha= alpha *f_alpha;
+		}
+	}
+	else if( P<=0)
+	{
+		deltaT=deltaT*fdec;
+		negative_count = 0;
+		
+		for(int iatom=0; iatom<numIon ; iatom++)
+        {
+			vel[iatom].x = 0.0;
+			vel[iatom].y = 0.0;
+			vel[iatom].z = 0.0;
+		}
+		
+		alpha=alpha_start;
+	}
+}
\ No newline at end of file
diff --git a/ABACUS.develop/source/src_pw/MD_fire.h b/ABACUS.develop/source/src_pw/MD_fire.h
new file mode 100644
index 0000000000..88e07dcc15
--- /dev/null
+++ b/ABACUS.develop/source/src_pw/MD_fire.h
@@ -0,0 +1,27 @@
+#ifndef MD_FIRE_H
+#define MD_FIRE_H
+
+#include "src_pw/tools.h"
+
+class MD_fire
+{
+    public:
+        MD_fire();
+        ~MD_fire(){};
+        void check_FIRE(
+            const int& numIon, 
+            const Vector3<double>const* force, 
+            double& deltaT, 
+            Vector3<double>* vel);
+    private:
+        double alpha_start ;//alpha_start begin
+        double alpha;//alpha begin
+        double finc;//finc begin
+        double fdec;//fdec begin
+        double f_alpha;
+        int N_min;
+        double dt_max;//dt_max
+        int negative_count;//Negative count
+};
+
+#endif
\ No newline at end of file
diff --git a/ABACUS.develop/source/src_pw/MD_func.cpp b/ABACUS.develop/source/src_pw/MD_func.cpp
index 9905f14eed..9756d7b4a7 100644
--- a/ABACUS.develop/source/src_pw/MD_func.cpp
+++ b/ABACUS.develop/source/src_pw/MD_func.cpp
@@ -1,225 +1,85 @@
-#include "MD_func.h"
-
-
-
-void MD_func::initMD(){
-	//parameter from input file
-	mdtype=INPUT.md_mdtype;
-	Qmass=INPUT.md_qmass/6.02/9.109*1e5;
-	dt=INPUT.md_dt/fundamentalTime/1e15;
-	temperature=INPUT.md_tfirst/3.1577464e5;
-	rstMD=INPUT.md_rstmd;
-	dumpmdfreq=INPUT.md_dumpmdfred;
-	fixTemperature=INPUT.md_fixtemperature;
-	ediff=INPUT.md_ediff;
-	ediffg=INPUT.md_ediffg;
-	mdoutputpath=INPUT.md_mdoutpath;
-	NVT_control = INPUT.md_nvtcontrol;
-   	NVT_tau = INPUT.md_nvttau;
-	MNHC = INPUT.md_mnhc;
-
-
-	//other parameter default
-	outputstressperiod = 1 ;// The period to output stress
-	numIon=ucell.nat;
-	ntype=ucell.ntype;
-	step_rst=0;
-	ionlatvec=ucell.latvec;
-
-	mdenergy=en.etot/2;
-
-	//allocate for MD
-	this->md_allocate();
-
-	//MD starting setup
-	if(rstMD==0){
-		connection0();
-		connection1();
-		//A fresh new MD: Do not restart MD
-		InitVelocity() ;
-		// Initialize thermostat, and barostat
-		
-	}
-	else if ( rstMD>0 ){
-		connection0();
-		connection1();
-		InitVelocity() ;
-		if(!RestartMD()){
-			cout<<"error in restart MD!"<<endl;
-			return;
-		}
-		connection2();
-	}
-	if(NVT_tau<1e-10)
-	{
-		if(NVT_control == 1) NVT_tau = 20 * dt; //NHC
-		else if(NVT_control == 2) NVT_tau = 20 * dt; //LGV
-		else if(NVT_control == 3) NVT_tau = 20 * dt; //ADS
-		else WARNING_QUIT("md:InitMD", "please choose available reservoir!!!");
-	}
-	else
-	{
-		NVT_tau *= fundamentalTime*1e15;
-	}
-	//tau = 1.0/(NVT_tau*fundamentalTime*1e15);
-	// Q=KbT*tau**2
-	if(Qmass<1e-10)//Qmass=dt*fundamentalTime*1e15/6.02/9.109*1e5;
-	Qmass = pow(1.0/NVT_tau,2)*temperature;///beta;
-	//NHC setup
-	if(mdtype==1 && NVT_control == 1)
-	{
-		unsigned long init[4]={0x123, 0x234, 0x345, 0x456}, length=4;
-		init_by_array(init, length);
-		ofs_running<<" ...............Nose-Hoover Chain parameter initialization...............  " << endl;
-		ofs_running<<" Temperature =    "<< temperature << endl;
-		ofs_running<<" Temperature2 =    "<< temperature/K_BOLTZMAN_AU << endl;
-		ofs_running<<" NHC frequency =    "<< 1.0/NVT_tau << endl;
-		ofs_running<<" NHC chain =    "<< MNHC << endl;
-		ofs_running<<" Qmass  =    "<< Qmass << endl;
-		ofs_running<<" ...............................................................  " << endl;
-		w[0]=0.7845136105;
-		w[6]=0.7845136105;
-		w[1]=0.2355732134;
-		w[5]=0.2355732134;
-		w[2]=-1.177679984;
-		w[4]=-1.177679984;
-		w[3]=1-w[0]-w[1]-w[2]-w[4]-w[5]-w[6];
-		//
-		for(int i=0;i<nsy;i++){
-			delta[i]=w[i]*dt;
-		}
-		
-		
-		for(int j=0;j<MNHC;j++){
-			for(int k=0;k<numIon;k++){
-				NHCeta[j*numIon+k].x=genrand_real2()-0.5; 
-				NHCeta[j*numIon+k].y=genrand_real2()-0.5; 
-				NHCeta[j*numIon+k].z=genrand_real2()-0.5;
-				NHCpeta[j*numIon+k].x=gaussrand()*sqrt(Qmass * temperature);
-				NHCpeta[j*numIon+k].y=gaussrand()*sqrt(Qmass * temperature);
-				NHCpeta[j*numIon+k].z=gaussrand()*sqrt(Qmass * temperature);
-			}
-		}
-		
-		for(int j=MNHC-1;j>=1;j--){
-			for(int k=0;k<numIon;k++){
-				G[j*numIon+k].x=pow(NHCpeta[(j-1)*numIon+k].x,2)/Qmass - 1.0 * temperature;
-				G[j*numIon+k].y=pow(NHCpeta[(j-1)*numIon+k].y,2)/Qmass - 1.0 * temperature;
-				G[j*numIon+k].z=pow(NHCpeta[(j-1)*numIon+k].z,2)/Qmass - 1.0 * temperature;
-			}
-		}
-		
-		for(int k=0;k<numIon;k++){
-			G[0*numIon+k].x=pow(vel[k].x,2)*allmass[k]-1.0 * temperature;
-			G[0*numIon+k].y=pow(vel[k].y,2)*allmass[k]-1.0 * temperature;
-			G[0*numIon+k].z=pow(vel[k].z,2)*allmass[k]-1.0 * temperature;
-		}
-		ofs_running << "finish NHC thermostat define" << endl; //zifei
-	}
-
-	return;
-}
-	
-/*void md::runMD(int step1){
-       if(rstMD==1)MDNVT.runnvt();
-       else if(mdtype==3)MDNVE.runNVE(step1);
-       return;
-}*///we need this form
-
-bool  MD_func::RestartMD(){
-	int i;
+#include "src_pw/MD_func.h"
+#include "src_pw/tools.h"
+#include "src_lcao/FORCE_STRESS.h"
+#include "src_pw/forces.h"
+#include "src_pw/stress_pw.h"
+#include "src_global/sltk_atom_arrange.h"
+
+bool MD_func::RestartMD(const int& numIon, Vector3<double>* vel, int& step_rst)
+{
+	int error(0);
 	double *vell=new double[numIon*3];
-	double *cart=new double[numIon*3];
-	double *dira=new double[numIon*3];
-	double *nhcg=new double[numIon*3*MNHC];
-	double *nhcp=new double[numIon*3*MNHC];
-	double *nhce=new double[numIon*3*MNHC];
-	if (!MY_RANK){
+	if (!MY_RANK)
+	{
 		stringstream ssc;
-		ssc << global_out_dir << "Restart_md.dat";
+		ssc << global_readin_dir << "Restart_md.dat";
 		ifstream file(ssc.str().c_str());
 
-		if(!file){
-			cout<<"please ensure whether 'Restart_md.dat' exists"<<endl;
-			exit(0);
+		if(!file)
+		{
+			cout<<"please ensure whether 'Restart_md.dat' exists!"<<endl;
+			error = 1;
 		}
-		if (file.good())
+		if (!error)
 		{
 			file.ignore(11, '\n');
-		}
 
-//----------------------------------------------------------
-// main parameters
-//----------------------------------------------------------
-		file.ignore(7, '\n');
-		//                file>>xLogS;
-		file.get();	
-		file.ignore(7, '\n');
-		//                file>>vLogS;
-		file.get();
-		file.ignore(23, '\n');
-		for(i = 0;i<numIon*3;i++){
-			file>>vell[i];
+        //----------------------------------------------------------
+        // main parameters
+        //----------------------------------------------------------
+			file.ignore(13, '\n');//read and assert number of atoms
+			int num;
+			file>>num;
+			if(num != numIon)
+			{
+				cout<<"please ensure whether 'Restart_md.dat' right!"<<endl;
+				error = 1;
+			}
 		}
-		/*file.get();	file.ignore(17, '\n');
-		for(i=0;i<numIon*3;i++){
-			file>>cart[i];
+		if (!error)
+		{
+			file.get();
+			file.ignore(23, '\n');//read velocities
+			for(int i = 0;i<numIon*3;i++)
+			{
+				file>>vell[i];
+			}
+			file.get();
+			file.ignore(6, '\n');//read start step of MD
+			file>>step_rst;
+			file.close();
 		}
-		file.get();	file.ignore(13, '\n');
-		for(i=0;i<numIon*3;i++){
-			file>>dira[i];
-		}*/
-		file.get();     
-		file.ignore(6, '\n');
-		file>>step_rst;
-		file.close();
 	}
-//2015-09-05, xiaohui
+#ifdef __MPI
+	MPI_Bcast(&error,1,MPI_INT,0,MPI_COMM_WORLD);
+#endif
+	if(error)
+	{
+		delete[] vell;
+		exit(0);
+	}
 #ifdef __MPI
 	MPI_Bcast(&step_rst,1,MPI_INT,0,MPI_COMM_WORLD);
 	MPI_Bcast(vell,numIon*3,MPI_DOUBLE,0,MPI_COMM_WORLD);
 #endif
-	for(i=0;i<numIon;i++){
+	for(int i=0;i<numIon;i++)
+	{
 		vel[i].x=vell[i*3];
 		vel[i].y=vell[i*3+1];
 		vel[i].z=vell[i*3+2];
 	}
-	// DEBUG
-	//int debug = 0;
-	//if (debug) {
-	//  cout<<" dump cell + ion infomation file to -> echo_md_restart.dat, for debug purpose. "<<endl;
-	//	fstream debugfile;
-	//	debugfile.open("echo_md_restart.dat");
-	//  debugfile<<"  "<<endl;
-	//debugfile<<"cellReal[0] = "<<a.latvec.e11<<" "<<a.latvec.e12<<" "<<a.latvec.e13<<endl;
-	//debugfile<<"cellReal[1] = "<< a.latvec.e21<<" "<<a.latvec.e22<<" "<<a.latvec.e23<<endl;
-	//debugfile<<"cellReal[2] = "<< a.latvec.e31<<" "<<a.latvec.e32<<" "<<a.latvec.e33<<endl;
-	//debugfile<<"ION FRACTIONAL COORDINATES:"<<endl;
-	//for(ii=0;ii<numIon;ii++){
-	//  debugfile<<taudirac[ii].x<<" "<<taudirac[ii].y<<" "<<taudirac[ii].z<<endl;
-	//}
-	//debugfile<<"ION VELOCITIES:"<<endl;
-	//for( ii=0;ii<ionIon;ii++){
-	// debugfile<<vel[ii].x<<" "<<vel[ii].y<<" "<<vel[ii].z<<endl;
-	//	}
-	//  debugfile.close();
-	//}
-	// END DEBUG
-	if(vell!=NULL)delete []vell;
-	if(cart!=NULL)delete []cart;
-	if(dira!=NULL)delete []dira;
-	if(nhcg!=NULL)delete []nhcg;
-	if(nhcp!=NULL)delete []nhcp;
-	if(nhce!=NULL)delete []nhce;
+
+	delete []vell;
 	return true;
 }
 
-void MD_func::mstout(int step){
+void MD_func::mdRestartOut(const int& step, const int& recordFreq, const int& numIon, Vector3<double>* vel)
+{
 //this function used for outputting the information of restart file
 	bool pass;
 	pass = 0;
 
-	if (dumpmdfreq==1||step==1||( dumpmdfreq > 1&& step%dumpmdfreq==0 ) )
+	if (recordFreq==1||step==1||( recordFreq > 1&& step%recordFreq==0 ) )
 		pass =1;
 	if (!pass) return;
 
@@ -228,25 +88,18 @@ void MD_func::mstout(int step){
 		ssc << global_out_dir << "Restart_md.dat";
 		ofstream file(ssc.str().c_str());
 		file<<"MD_RESTART"<<endl;
+		file<<"ATOM_NUMBERS: "<<numIon<<endl;
 		file<<"ION_VELOCITIES_(a.u.): "<<endl;
 		for(int i=0;i<numIon;i++){
 			file<<setprecision (12)<<vel[i].x<<" "<<setprecision (12)<<vel[i].y<<" "<<setprecision (12)<<vel[i].z<<endl;
 		}
-/*	file<<"CARTESIAN_COORD: "<<endl;
-	for(i=0;i<numIon;i++){
-		file<<cartNoWrap[i].x<<" "<<cartNoWrap[i].y<<" "<<cartNoWrap[i].z<<endl;
-	}
-	file<<"DIRAC_COORD: "<<endl;
-	for(i=0;i<numIon;i++){
-		file<<taudirac[i].x<<" "<<taudirac[i].y<<" "<<taudirac[i].z<<endl;
-	}*/
 		file<<"step: "<<step<<endl;                
 		file.close();
 	}
 	return;
 }
 
-double MD_func::GetAtomKE(){
+double MD_func::GetAtomKE(const int& numIon, const Vector3<double>const* vel, const double const* allmass){
 //---------------------------------------------------------------------------
 // DESCRIPTION:
 //   This function calculates the classical kinetic energy of a group of atoms.
@@ -264,7 +117,12 @@ double MD_func::GetAtomKE(){
 	return ke;
 }
 
-void MD_func::InitVelocity()
+void MD_func::InitVelocity(
+	const int& numIon, 
+	const double& temperature, 
+	const double& fundamentalTime, 
+	const double const* allmass,
+	Vector3<double>* vel)
 {
 	if(!MY_RANK){ //xiaohui add 2015-09-25
 		srand(time(0));
@@ -323,82 +181,7 @@ void MD_func::InitVelocity()
 	return;
 }
 
-
-void MD_func::OutputMDHeader()
-{
-//---------------------------------------------------
-// Output the Header for MD simulation
-//---------------------------------------------------
-
-	// Set up extra output to ion optimizer / MD header
-	if (!MY_RANK)
-		out<<"Q ="<<Qmass<<"a.u.\n"<<" dt = "<<dt*fundamentalTime*1e15<<"fs\n"<< "T ="<<temperature/K_BOLTZMAN_AU<<endl;
-		ofs_running<<"Q ="<<Qmass<<"a.u.\n"<<" dt = "<<dt*fundamentalTime*1e15<<"fs\n"<< "T ="<<temperature/K_BOLTZMAN_AU<<endl;
-		return;
-} 
-
-void MD_func::OutputMDGeom(int iter)
-{
-//---------------------------------------------------------------------------
-// DESCRIPTION:
-//  Output the ion's position and cell's structure for monitoring and MD
-//  restart.
-//---------------------------------------------------------------------------
-  
-	bool pass;
-	int i,it;
-	ofstream uf;
-	string iteration_number;
-	string eleName;
-	   
-	//! >> FUNCTION << !
-	pass = 0;
-	// if(rankGlobal!=0) return;
-	if(dumpmdfreq==1||iter==1||( dumpmdfreq > 1&&(iter%dumpmdfreq)==0 ) ) 
-		pass = 1;
-	if(!pass)return;
-	// cout<<"test the function OutputMDGeom"<<endl;
-
-	//the final file name is the md_geometry plus the iteraction number
-	iteration_number=intTurnTostring(iter,iteration_number);
-	string  message=mdoutputpath+"_ion"+iteration_number+".txt";
-	//  cout<<"file name is:"<<message<<endl;
-	// open the file: 'ion.1.dat'
-	uf.open(message.c_str());
-	//  cout<<"already open file"<<endl;
-	uf<< "BLOCK LATTICE_CART"<<endl;
-	uf<<"    "<<ionlatvec.e11*ucell.lat0<<"  "<<ionlatvec.e12*ucell.lat0<<"  "<<ionlatvec.e13*ucell.lat0<<endl;
-	  
-	uf<<"    "<<ionlatvec.e21*ucell.lat0<<"  "<<ionlatvec.e22*ucell.lat0<<"  "<<ionlatvec.e23*ucell.lat0<<endl;
-
-	uf<<"    "<<ionlatvec.e31*ucell.lat0<<"  "<<ionlatvec.e32*ucell.lat0<<"  "<<ionlatvec.e33*ucell.lat0<<endl;
-
-	uf<< "END BLOCK LATTICE_CART"<<endl;
-	//  cout<<"already out put lattice cart"<<endl;
-  
-	//output the atom positions.
-	uf<<"BLOCK POSITIONS_CART"<<endl;
-	int ion=0;
-	for(it=0;it<ntype;it++){
-		for(i = 0;i< na[it];i++){
-			eleName=ucell.atoms[it].label; 
-			uf<<"  "<<eleName 
-				<<"  "<<cartNoWrap[ion].x  // unit of cartNorap is 'Bohr' 
-				<<"  "<<cartNoWrap[ion].y   
-				<<"  "<<cartNoWrap[ion++].z<<endl;
-		}
-	}
-	uf<< "END BLOCK POSITIONS_CART"<<endl;
-	uf<<"already output positions cart"<<endl;
-	double distance;
-	distance=pow(cartNoWrap[1].x-cartNoWrap[0].x,2.0)+pow(cartNoWrap[1].y-cartNoWrap[0].y,2.0)+pow(cartNoWrap[1].z-cartNoWrap[0].z,2);
-	distance=sqrt(distance);
-	//  cout<<"distance is "<<distance<<" ";
-	uf.close();
-	return;
-}
-
-void MD_func::ReadNewTemp(int step )
+/*void MD_func::ReadNewTemp(int step )
 {
 //-----------------------------------------------------------------
 //  If fixTemperature == 0, then this subroutine will be skipped
@@ -446,7 +229,7 @@ void MD_func::ReadNewTemp(int step )
 	}
 
 	return;
-}
+}*/
 
 //int to string and add to output path
 string MD_func::intTurnTostring(long int iter, string path)
@@ -469,117 +252,68 @@ string MD_func::intTurnTostring(long int iter, string path)
 	return path;
 }
 
-void MD_func::connection0(){
+int MD_func::getMassMbl(const UnitCell_pseudo &unit_in, double* allmass, Vector3<int>* ionmbl)
+{
 //some prepared information
 //mass and degree of freedom
 	int ion=0;
-	nfrozen=0;
-	for(int it=0;it<ntype;it++){
-		for(int i=0;i<na[it];i++){
-			allmass[ion]=ucell.atoms[it].mass/6.02/9.109*1e5;
-			ionmbl[ion]=ucell.atoms[it].mbl[i];
+	int nfrozen=0;
+	for(int it=0;it<unit_in.ntype;it++){
+		for(int i=0;i<unit_in.atoms[it].na;i++){
+			allmass[ion]=unit_in.atoms[it].mass/6.02/9.109*1e5;
+			ionmbl[ion]=unit_in.atoms[it].mbl[i];
 			if (ionmbl[ion].x==0||ionmbl[ion].y==0||ionmbl[ion].z==0) nfrozen++;
 
 			ion++;
 		}
 	}
-	return;
-}
-	
-void MD_func::connection1()
-{
-//call each atom's position
-	int i,it,ion;
-	ion=0;
-	for(it=0;it<ntype;it++){
-		for(i=0;i<na[it];i++){
-			cartNoWrap[ion]=ucell.atoms[it].tau[i]*ucell.lat0;
-			taudirac[ion]=ucell.atoms[it].taud[i];
-			ion++;
-		}
-	}
-	return;
+	return nfrozen;
 }
 
-void MD_func::connection2()
-{
-//to refresh the atoms' positions
-	int i,it,ion;
-	ion=0;
-	for(it=0;it<ntype;it++){
-		for(i=0;i<na[it];i++){
-			ucell.atoms[it].tau[i]=cartNoWrap[ion]/ucell.lat0;
-			ucell.atoms[it].taud[i]=taudirac[ion];
-			ion++;
-		}
-	}
-	return;
-}
-
-void MD_func::callforce()
+void MD_func::callInteraction_LCAO(const int& numIon, Vector3<double>* force, matrix& stress_lcao)
 {
 //to call the force of each atom
-	int ion;
-	Force_LCAO FL;
-	FL.allocate (); 
-	FL.start_force();
-	for(ion=0;ion<numIon;ion++){
-		force[ion].x=FL.fcs(ion,0)/2.0;
-		force[ion].y=FL.fcs(ion,1)/2.0;
-		force[ion].z=FL.fcs(ion,2)/2.0;
-		// cout<<force[ion].x/0.0195<<" "<<force[ion].y/0.0195<<" "<<force[ion].z/0.0195<<endl;
+	matrix fcs;//temp force matrix
+	Force_Stress_LCAO FSL;
+	FSL.allocate (); 
+	FSL.getForceStress(FORCE, STRESS, TEST_FORCE, TEST_STRESS, fcs, stress_lcao);
+	for(int ion=0;ion<numIon;ion++){
+		force[ion].x =fcs(ion, 0)/2.0;
+		force[ion].y =fcs(ion, 1)/2.0;
+		force[ion].z =fcs(ion, 2)/2.0;
 	}
 #ifdef __MPI //2015-10-01, xiaohui
 	atom_arrange::delete_vector( SEARCH_RADIUS );
 #endif //2015-10-01, xiaohui
-	//to call the stress
-	if(STRESS)
-	{
-	//	matrix stress_lcao;//this is the stress matrix same as src_pw/ion.cpp
-		stress_lcao.create(3,3);
-		FL.cal_stress(stress_lcao);
-	}
+
 	return;
 }
-
-void MD_func::moveatoms(int step){
-//intent to ensure that no atom is out of the cell
-
-	for(int ia=0;ia<numIon;ia++)
+void MD_func::callInteraction_PW(const int& numIon, Vector3<double>* force, matrix& stress_pw)
+{
+//to call the force of each atom
+	matrix fcs;//temp force matrix
+	Forces ff;
+	ff.init(fcs);
+	for(int ion=0;ion<numIon;ion++){
+		force[ion].x =fcs(ion, 0)/2.0;
+		force[ion].y =fcs(ion, 1)/2.0;
+		force[ion].z =fcs(ion, 2)/2.0;
+	}
+	if(STRESS)
 	{
-//changed by a constrain box
-                 
-		if(taudirac[ia].x<0) taudirac[ia].x += 1.0;
-		if(taudirac[ia].y<0) taudirac[ia].y += 1.0;
-		if(taudirac[ia].z<0) taudirac[ia].z += 1.0;
-		if(taudirac[ia].x>=1.0) taudirac[ia].x -= 1.0;
-		if(taudirac[ia].y>=1.0) taudirac[ia].y -= 1.0;
-		if(taudirac[ia].z>=1.0) taudirac[ia].z -= 1.0;
-
-		if(taudirac[ia].x<0 || taudirac[ia].y<0
-			|| taudirac[ia].z<0 ||
-			taudirac[ia].x>=1.0 ||
-			taudirac[ia].y>=1.0 ||
-			taudirac[ia].z>=1.0)
-		{
-			cout << " ion: "<< ia  << endl;
-			cout << "d=" << taudirac[ia].x << " " <<
-				taudirac[ia].y << " " << taudirac[ia].z << endl;
-			WARNING_QUIT("md::moveatoms","the movement of atom is larger than the length of cell.");
-		}
-
-		cartNoWrap[ia] = taudirac[ia] * ionlatvec*ucell.lat0;
+		Stress_PW ss;
+		ss.cal_stress(stress_pw);
 	}
 	return;
 }
 
-void MD_func::printpos(string file,int iter)
+void MD_func::printpos(const string& file, const int& iter, const int& recordFreq, const UnitCell_pseudo& unit_in)
 {
 //intend to output the positions of atoms to ordered file
 	bool pass;
 	pass = 0;
 
-	if (dumpmdfreq==1||iter==1||( dumpmdfreq > 1&& iter%dumpmdfreq==0 ) )
+	if (recordFreq==1||iter==1||( recordFreq > 1&& iter%recordFreq==0 ) )
 		pass =1;
 	if (!pass) return;
 
@@ -587,26 +321,68 @@ void MD_func::printpos(string file,int iter)
 	string file2=file+".cif";
 
 	//xiaohui add 'OUT_LEVEL', 2015-09-16
-	if(OUT_LEVEL == "i"||OUT_LEVEL == "ie") ucell.print_tau();
-	if(OUT_LEVEL == "i"||OUT_LEVEL == "ie") ucell.print_cell_xyz(file1);
-	ucell.print_cell_cif(file2);
+	if(OUT_LEVEL == "i"||OUT_LEVEL == "ie") unit_in.print_tau();
+	if(OUT_LEVEL == "i"||OUT_LEVEL == "ie") unit_in.print_cell_xyz(file1);
+	unit_in.print_cell_cif(file2);
 	stringstream ss;
 
 	ss << global_out_dir << "STRU_MD";
 
 	//zhengdy modify 2015-05-06, outputfile "STRU_Restart"
-	ucell.print_stru_file(ss.str(),2);
+	unit_in.print_stru_file(ss.str(),2);
 
 	return;
 }
 
 //rescale velocities to target temperature.
-void MD_func::scalevel()
+void MD_func::scalevel(
+	const int& numIon,
+	const int& nfrozen,
+	const double& temperature,
+	Vector3<double>* vel,
+	const double* allmass
+)
 {
-	double ke=GetAtomKE();
+	double ke=GetAtomKE(numIon, vel, allmass);
 	if(ke>1e-9)
-		for(int i=0;i<numIon;i++){
+	{
+		for(int i=0;i<numIon;i++)
+		{
 			vel[i]*=sqrt(3*(numIon-nfrozen)*temperature/ke/2);
 		}
+	}
 	return;
+}
+
+double MD_func::Conserved(const double KE, const double PE, const int number){
+//---------------------------------------------------------------------------
+//   This function calculates the conserved quantity for the NVE system. 
+//----------------------------------------------------------------------------
+
+   	// KE   // Kinetic energy of particles (due to their velocity)
+   	// PE   // Potential energy (DFT total energy)
+	// number //number of atoms which have full freedoms
+
+  	double Conserved; // The conserved quantity
+
+   	Conserved = KE + PE ;
+   
+	if (!MY_RANK)
+	{                 
+		ofs_running<< "NVE Conservation     : "<< Conserved<<" (Hartree)"<<endl;
+		ofs_running<< "NVE Temperature      : "<< 2*KE/(3*number)/K_BOLTZMAN_AU<<" (K)"<<endl;
+		ofs_running<< "NVE Kinetic energy   : "<< KE<<" (Hartree)"<<endl;
+		ofs_running<< "NVE Potential energy : "<< PE<<" (Hartree)"<<endl;
+	}
+   	return Conserved;
+}
+
+double MD_func::MAXVALF(const int numIon, const Vector3<double>const* force){
+	//cout<<"enter in MAXVALF"<<endl;
+	double max=0;
+	for(int i=0;i<numIon;i++){
+		double force0 = pow(force[i].x,2)+pow(force[i].y,2)+pow(force[i].z,2);
+		if(max<force0) max = force0;
+	}
+	return max;
 }
\ No newline at end of file
diff --git a/ABACUS.develop/source/src_pw/MD_func.h b/ABACUS.develop/source/src_pw/MD_func.h
index d0abfd3cf6..f1b45dedaa 100644
--- a/ABACUS.develop/source/src_pw/MD_func.h
+++ b/ABACUS.develop/source/src_pw/MD_func.h
@@ -1,33 +1,36 @@
 #ifndef MD_FUNC_H
 #define MD_FUNC_H
 
+#include "src_pw/unitcell_pseudo.h"
+
+
 class MD_func
 {
     public:
     MD_func(){};
     ~MD_func(){};
-    void initMD();
-	bool RestartMD();
-    void mstout(int step);
-//	void RemoveMovementOfCenterOfMass();
-	double GetAtomKE();
-//	void MakeIntCoeff();
-	void InitVelocity();
-//	void MonitorMeanSquareDisplacement(int step);
-//	void OutputIonVelocityAndStats(int iter);
-	void OutputMDHeader();
-	void OutputMDGeom(int iter);
-	void ReadNewTemp(int step);
-	string intTurnTostring(int iter,string path);
-	void connection0();
-	void connection1();
-	void connection2();
-	void callforce();
-    void moveatoms(int step);
-    void printpos(string file,int iter);
-//        void md_release();
-    void scalevel();
-	double MAXVALF();
-	double Conserved(double KE, double PE);
+	bool RestartMD(const int& numIon, Vector3<double>* vel, int& step_rst);
+    void mdRestartOut(const int& step, const int& recordFreq, const int& numIon, Vector3<double>* vel);
+	double GetAtomKE(const int& numIon, const Vector3<double>const* vel, const double const* allmass);
+	void InitVelocity(
+		const int& numIon, 
+		const double& temperature, 
+		const double& fundamentalTime, 
+		const double const* allmass,
+		Vector3<double>* vel);
+//	void ReadNewTemp(int step);
+	string intTurnTostring(long int iter,string path);
+	int getMassMbl(const UnitCell_pseudo &unit_in, double* allmass, Vector3<int>* ionmbl);
+	void callInteraction_LCAO(const int& numIon, Vector3<double>* force, matrix& stress_lcao);
+	void callInteraction_PW(const int& numIon, Vector3<double>* force, matrix& stress_pw);
+    void printpos(const string& file, const int& iter, const int& recordFreq, const UnitCell_pseudo& unit_in);
+    void scalevel(
+		const int& numIon,
+		const int& nfrozen,
+		const double& temperature,
+		Vector3<double>* vel,
+		const double* allmass);
+	double MAXVALF(const int numIon, const Vector3<double>const* force);
+	double Conserved(const double KE, const double PE, const int number);
 };
 #endif
\ No newline at end of file
diff --git a/ABACUS.develop/source/src_pw/MD_parameters.h b/ABACUS.develop/source/src_pw/MD_parameters.h
index 73fe75ff50..d985329b45 100644
--- a/ABACUS.develop/source/src_pw/MD_parameters.h
+++ b/ABACUS.develop/source/src_pw/MD_parameters.h
@@ -3,25 +3,43 @@
 
 class MD_parameters
 {
-    MD_parameters(){};
+public:
+    MD_parameters()
+	{
+		mdtype=1;
+		NVT_tau=0;
+		NVT_control=1;
+		dt=1;
+		Qmass=1;
+		tfirst=0;         //kelvin
+		tlast=-1;
+		recordFreq=1;
+		mdoutputpath="mdoutput";
+		rstMD=0;
+		fixTemperature=1;
+		ediff=1e-4;
+		ediffg=1e-3;
+		MNHC=4;
+	};
     ~MD_parameters(){};
 
-    private:
     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 
     double Qmass;       // Inertia of extended system variable
     double dt ;         // Time increment (hbar/E_hartree)
-    double temperature; //Temperature (in Hartree, 1 Hartree ~ 3E5 K)
-    int ntype;          //
-    int dumpmdfreq;     // The period to dump MD information for monitoring and restarting MD
+    double tfirst; //Temperature (in Hartree, 1 Hartree ~ 3E5 K)
+    double tlast;
+    int recordFreq;     // The period to dump MD information for monitoring and restarting MD
     int NVT_control;
     double NVT_tau;
     int MNHC;
     string mdoutputpath;// output directory of md files: .ion .vel
 	double ediff;       //parameter for constrain 
 	double ediffg;      //parameter for constrain
+    int fixTemperature;
 
 };
 
+
 #endif
\ No newline at end of file
diff --git a/ABACUS.develop/source/src_pw/MD_run.cpp b/ABACUS.develop/source/src_pw/MD_run.cpp
deleted file mode 100644
index e40142cd43..0000000000
--- a/ABACUS.develop/source/src_pw/MD_run.cpp
+++ /dev/null
@@ -1,540 +0,0 @@
-#include "MD_run.h"
-
-//define in MD_run.h
-//class MD_run
-
-
-void MD_run::runnvt(int step1){
-//------------------------------------------------------------------------------
-// DESCRIPTION:
-// Molecular dynamics calculation with fixed Volume and slight fluctuated temperature
-// Using thermostat : 1, Nose-Hoover Chains; 2, Langevin; 3, Anderson
-// Normal Nose-Hoover thermostat method is retained for test.
-//------------------------------------------------------------------------------
-
-	TITLE("MD_run","runnvt");
-	timer::tick("MD_run","runnvt",'C');
-	step=step1+step_rst;
-	//the real MD step
-	if (!MY_RANK)out<<"step: " << step <<endl;
-	
-	//change target temperature
-	if (step!=1)ReadNewTemp( step );
-	
-	cout << " ---------------------------------------------------------" << endl;
-	
-	cout<< " Target temperature  : "<< temperature/K_BOLTZMAN_AU<< " (K)"<<endl;
-	
-	if(step==1||step%fixTemperature==1)
-    {
-		for(int k=0;k<numIon;k++)
-        {
-			if(ionmbl[k].x==0)vel[k].x=0;
-			if(ionmbl[k].y==0)vel[k].y=0;
-			if(ionmbl[k].z==0)vel[k].z=0;
-		}
-		scalevel();
-	}
-
-	// get the kinetic energy
-	twiceKE = GetAtomKE();
-	twiceKE = twiceKE * 2;
-	
-	//Set up forces 
-	callforce();
-
-	//print total stress + stress_MD
-	if(STRESS)
-	{
-		cout<<"output Pressure for check!"<<endl;
-		double press;
-		for(int i=0;i<3;i++)
-			press += stress_lcao(i,i)/3;
-		press += twiceKE/3/ucell.omega; //output virtual press = 2/3 *Ek/V + sum(sigma[i][i])/3
-		double unit_transform = RYDBERG_SI / pow(BOHR_RADIUS_SI,3) * eps8 ;
-		cout<<"Virtual Pressure is "<<press*unit_transform<<" Kbar "<<endl;
-	}
-
-	for(int k=0;k<numIon;k++){
-		if(ionmbl[k].x==0)force[k].x=0;
-		if(ionmbl[k].y==0)force[k].y=0;
-		if(ionmbl[k].z==0)force[k].z=0;
-	}
-	maxForce = MAXVALF();
-	mdenergy=en.etot/2;
-	if(NVT_control == 1) hamiltonian = NHChamiltonian(twiceKE/2, mdenergy);
-	else hamiltonian = Conserved(twiceKE/2, mdenergy);
-	//----------------------------------------------
-	// big loop
-	//-----------------------------------------------
-	cout<<" "<<std::left<<setw(12)<<"MD_STEP"<<std::left<<setw(12)<< "SystemE"<<std::left<<setw(12)<< "Conserved"<<std::left<<setw(12)<< "DeltaE"<<std::left<<setw(12)<< "Temperature"<<endl;
-	
-	cout<<" "<<std::left<<setw(12)<<step<<std::left<<setw(12)<< mdenergy<<std::left<<setw(12)<< hamiltonian<<std::left<<setw(12)<< mdenergy-oldEtot<<std::left<<setw(12)<<twiceKE/(3*double(numIon-nfrozen))/K_BOLTZMAN_AU<<endl;
-	
-	oldEtot=mdenergy;
-
-	if (!MY_RANK)
-    {
-	      out<<" --------------------------------------------------"<<endl;
-	      out<<" Molecular Dynamics (NVT) STEP "<< step<<endl;
-	      out<< " --------------------------------------------------"<<endl;
-	      ofs_running<<" --------------------------------------------------"<<endl;
-	      ofs_running<<" Molecular Dynamics (NVT) STEP "<< step<<endl;
-	      ofs_running<< " --------------------------------------------------"<<endl;
-	}
-	
-	
-	// Calculate the Mean-Square-Displacement.
-	if(step==1&&rstMD==0)
-    { 
-        //Note: side scheme position before
-        //Now turn to middle scheme.
-	
-	    // New vel obtained (Verlet-like step)
-	    for(int  ii=0;ii<numIon;ii++){ 
-			mass = allmass[ii];
-			vel[ii] = vel[ii] + force[ii]/mass*dt/2.0;
-		}
-	}
-	else 
-    {
-		for(int  ii=0;ii<numIon;ii++)
-        { 
-			mass = allmass[ii];
-			vel[ii] = vel[ii] + force[ii]/mass*dt/2.0;
-		}
-		
-		twiceKE=GetAtomKE();
-		twiceKE = 2 * twiceKE;
-		if(NVT_control==1) hamiltonian = NHChamiltonian(twiceKE/2, mdenergy);
-		else hamiltonian = Conserved(twiceKE/2, mdenergy);
-        //Note: side scheme position before
-        //Now turn to middle scheme. 
-		for(int  ii=0;ii<numIon;ii++){ 
-			mass = allmass[ii];
-			vel[ii] = vel[ii] + force[ii]/mass*dt/2.0;
-		}
-	}
-
-	// Update the Non-Wrapped cartesion coordinates
-    if(mdtype==2) scalevel();//choose velocity scaling method
-
-	for(int  ii=0;ii<numIon;ii++)
-    { 
-        cartNoWrap[ii] = vel[ii]*dt/2.0 +cartNoWrap[ii];
-	}
-	
-	//by zifei
-	Integrator(NVT_control);
-    twiceKE=GetAtomKE();
-    twiceKE = 2 * twiceKE;
-	
-	for(int  ii=0;ii<numIon;ii++){ 
-		cartNoWrap[ii] = vel[ii]*dt/2.0 +cartNoWrap[ii];
-	}
-
-	// Calculate the maximal velocities.
-	// The step in fractional coordinates
-	maxStep = 0;
-	Vector3<double> fracStep;
-	for(int  ii=0;ii<numIon;ii++)
-    { 
-		if((pow(vel[ii].x,2)+pow(vel[ii].y,2)+pow(vel[ii].z,2))>maxStep)
-        {
-		    maxStep = pow(vel[ii].x,2)+pow(vel[ii].y,2)+pow(vel[ii].z,2);
-        }
-		Mathzone::Cartesian_to_Direct(vel[ii].x*dt/ucell.lat0,vel[ii].y*dt/ucell.lat0,vel[ii].z*dt/ucell.lat0,
-					ionlatvec.e11,ionlatvec.e12,ionlatvec.e13,
-					ionlatvec.e21,ionlatvec.e22,ionlatvec.e23,
-					ionlatvec.e31,ionlatvec.e32,ionlatvec.e33,
-					fracStep.x,fracStep.y,fracStep.z);
-
-		taudirac[ii] = taudirac[ii] + fracStep;
-	}
-
-	//save the atom position change to DFT module
-	moveatoms(step);
-	string t("md_pos_");
-	t=intTurnTostring(step,t);
-	connection2();
-	printpos(t,step);
-	maxStep = sqrt(maxStep)*dt;
-
-    if (!MY_RANK){
-         out<<setw(15)<<"maxForce=     "<<setw(15)<<"maxstep=      "<<setw(15)<<"step=     "<<endl;
-         out<<setw(15)<<maxForce<<setw(15)<<maxStep<<setw(15)<<step<<endl;
-         out<<step<<" "<< mdenergy<<" "<< hamiltonian<<" "<< mdenergy-oldEtot<<" "<< twiceKE/(3*double(numIon-nfrozen))/K_BOLTZMAN_AU<<endl;
-    }
-    oldEtot=mdenergy;
-    mstout(step);
-
-//2015-09-25, xiaohui
-#ifdef __MPI
-    if(step%100==0)
-    {
-        MPI_Bcast(vel,numIon*3,MPI_DOUBLE,0,MPI_COMM_WORLD);
-        MPI_Bcast(cartNoWrap,numIon*3,MPI_DOUBLE,0,MPI_COMM_WORLD);
-        MPI_Bcast(taudirac,numIon*3,MPI_DOUBLE,0,MPI_COMM_WORLD);
-    //    MPI_Bcast(&xLogS,1,MPI_DOUBLE,0,MPI_COMM_WORLD);
-    //    MPI_Bcast(&vLogS,1,MPI_DOUBLE,0,MPI_COMM_WORLD);
-    }
-#endif
-    timer::tick("MD_run","runnvt",'C');
-    return;
-}
-
-void MD_run::runNVE(int step1){
-//-------------------------------------------------------------------------------
-// Adiabatic ensemble 
-// Molecular dynamics calculation with Verlet algorithm
-//-------------------------------------------------------------------------------
-
-    TITLE("MD_run","runNVE");
-    timer::tick("MD_run","runNVE",'C');
-    step=step1+step_rst;
-    if(step==1)
-    {
-        cout<<" (NVE) Start."<<endl;
-    }
-    else 
-    {
-        cout<<" (NVE) step: "<<step<<endl;
-    }
-  
-  cout<<" (NVE) Fronzen Atom Number : "<<nfrozen<<endl;
-
-  // Set up extra output to ion optimizer / MD header
-  cout<<"Time interval       : "<< dt*fundamentalTime*1e15<< " (fs)"<<endl;
-  cout<<"Target temperature  : "<< temperature/K_BOLTZMAN_AU<< " (K)"<<endl;
-
-    if(step==1){
-        for(int k=0;k<numIon;k++)
-        {
-            if(ionmbl[k].x==0)vel[k].x=0;
-            if(ionmbl[k].y==0)vel[k].y=0;
-            if(ionmbl[k].z==0)vel[k].z=0;
-        }
-        scalevel();
-    }
-    moveatoms(step);
-    twiceKE=GetAtomKE();
-    twiceKE = twiceKE * 2;
-
-    tempNow = twiceKE/(3*double(numIon-nfrozen))/K_BOLTZMAN_AU;
-    cout<<" start temperature = "<< tempNow/K_BOLTZMAN_AU<< " (k)"<<endl;
-
-    // Set up forces 
-    callforce();
-	if(STRESS)
-	{
-		cout<<"output Pressure for check!"<<endl;
-		double press;
-		for(int i=0;i<3;i++)
-			press += stress_lcao(i,i)/3;
-		press += twiceKE/3/ucell.omega; //output virtual press = 2/3*Ek/V + sum(sigma[i][i])/3
-		double unit_transform = RYDBERG_SI / pow(BOHR_RADIUS_SI,3) * eps8 ;
-		cout<<"Virtual Pressure is "<<press*unit_transform<<" Kbar "<<endl;
-	}
-  
-    for(int k=0;k<numIon;k++)
-    {
-        if(ionmbl[k].x==0)force[k].x=0;
-        if(ionmbl[k].y==0)force[k].y=0;
-        if(ionmbl[k].z==0)force[k].z=0;
-    }
-
-    //cout<<"begin maxForce"<<endl;
-    maxForce = MAXVALF();
-    cout<<"maxForce: "<<sqrt(maxForce)<<endl; 
-    mdenergy=en.etot/2;
-    conservedE = Conserved(twiceKE/2, mdenergy);
-
-    cout<< "NVE_STEP "<<" "<<"SystemEnergy"<<" "<< "Conserved"<<" "<< "DeltaE"<<" "<< "Temperature"<<endl;
-
-    // Output the message to the screen.
-    cout<<step<<" "<< mdenergy<<" "<< conservedE<<" "
-        << mdenergy-oldEtot<<" "<< twiceKE/(3*double(numIon-nfrozen))/K_BOLTZMAN_AU<<endl;
-
-    cout<<"------------------------------------------------------------------------------"<<endl;
-    cout<< "MD(NVE) STEP "<< step<<endl;
-    cout<<"------------------------------------------------------------------------------"<<endl;
-    if (!MY_RANK){
-         out<<"------------------------------------------------------------------------------"<<endl;
-         out<< "MD(NVE) STEP "<< step<<endl;
-         out<<"------------------------------------------------------------------------------"<<endl;
-         ofs_running<<"------------------------------------------------------------------------------"<<endl;
-         ofs_running<< "MD(NVE) STEP "<< step<<endl;
-         ofs_running<<"------------------------------------------------------------------------------"<<endl;
-    }
-  
-    // (1) 1st step of Verlet-Velocity
-    // New velocity are obtained
-	for(int k=0;k<numIon;k++)
-    {       
-        mass = allmass[k];
-        vel[k] = vel[k] + force[k]/mass*dt/2.0;
-    } 
-    // (2) 2nd step of Verlet-Velocity 
-    // Update the Non-Wrapped cartesion coordinate
-    twiceKE = GetAtomKE();
-    twiceKE = 2 * twiceKE;
-    if(step!=1||rstMD==1)for(int k=0;k<numIon;k++)
-    {       
-        mass = allmass[k];
-        vel[k] = vel[k] + force[k]/mass*dt/2.0;
-    } 
-    for(int i=0;i<numIon;i++)
-    {
-        cartNoWrap[i] = vel[i]*dt + cartNoWrap[i];	
-	}
-    // Calculate the maximal velocities.
-    // The step in fractional coordinates
-    maxStep = 0;
-    for(int  i = 0;i< numIon;i++)
-    {
-        if((pow(vel[i].x,2.0)+pow(vel[i].y,2.0)+pow(vel[i].z,2.0))>maxStep)
-        {
-            maxStep = pow(vel[i].x,2.0)+pow(vel[i].y,2.0)+pow(vel[i].z,2.0);
-        }
-        Mathzone::Cartesian_to_Direct(vel[i].x*dt/ucell.lat0,vel[i].y*dt/ucell.lat0,vel[i].z*dt/ucell.lat0,
-                                        ionlatvec.e11,ionlatvec.e12,ionlatvec.e13,
-                                        ionlatvec.e21,ionlatvec.e22,ionlatvec.e23,
-                                        ionlatvec.e31,ionlatvec.e32,ionlatvec.e33,
-                                        fracStep.x,fracStep.y,fracStep.z);
-        taudirac[i] = taudirac[i] + fracStep;
-    }
-    moveatoms(step);
-    string t("md_pos_");
-    t=intTurnTostring(step,t);
-	connection2();
-    printpos(t,step);
-    maxStep = sqrt(maxStep)*dt;
-
-
-    // calculate the conserved quantity during MD 
-    hamiltonian = Conserved(twiceKE/2, mdenergy);
-
-    cout<< setprecision (9)<<hamiltonian<<" "<< setprecision (9)<<twiceKE/2<<endl;
-
-    
-    // Output the message to the screen.
-    if (!MY_RANK)
-    { 
-         out<<step<<" "<< mdenergy<<" "<< hamiltonian<<" "<< mdenergy-oldEtot<<" "<< twiceKE/(3*double(numIon-nfrozen))/K_BOLTZMAN_AU<<endl;
-         out<<step<<" "<< mdenergy<<" "<< conservedE<<" "<< mdenergy-oldEtot<<" "<< twiceKE/(3*double(numIon-nfrozen))/K_BOLTZMAN_AU<<endl;
-         ofs_running<<step<<" "<< mdenergy<<" "<< hamiltonian<<" "<< mdenergy-oldEtot<<" "<< twiceKE/(3*double(numIon-nfrozen))/K_BOLTZMAN_AU<<endl;
-         ofs_running<<step<<" "<< mdenergy<<" "<< conservedE<<" "<< mdenergy-oldEtot<<" "<< twiceKE/(3*double(numIon-nfrozen))/K_BOLTZMAN_AU<<endl;
-    }
-    oldEtot=mdenergy;
-    mstout(step);
-
-    cout<<"(NVE): this step finished."<<endl;
-    if (!MY_RANK)
-    {
-        out<<"(NVE): this step finished."<<endl;
-        ofs_running<<"(NVE): this step finished."<<endl;
-    }
-//2015-09-25, xiaohui
-#ifdef __MPI
-    if(step%100==0)
-    {
-        MPI_Bcast(vel,numIon*3,MPI_DOUBLE,0,MPI_COMM_WORLD);
-        MPI_Bcast(cartNoWrap,numIon*3,MPI_DOUBLE,0,MPI_COMM_WORLD);
-        MPI_Bcast(taudirac,numIon*3,MPI_DOUBLE,0,MPI_COMM_WORLD);
- //       MPI_Bcast(&xLogS,1,MPI_DOUBLE,0,MPI_COMM_WORLD);
- //       MPI_Bcast(&vLogS,1,MPI_DOUBLE,0,MPI_COMM_WORLD);
-    }
-#endif
-
-    timer::tick("MD_run","runNVE",'C');
-    return;
-}
-
-bool MD_run::runFIRE(int step1){
-//-------------------------------------------------------------------------------
-// REFERENCES:
-//   
-    TITLE("MD_run","runFIRE");
-    timer::tick("MD_run","runFIRE",'C');
-    if(step==1)
-    {
-        cout<<" FIRE Start."<<endl;
-    }
-    else 
-    {
-        cout<<" FIRE step: "<<step<<endl;
-    }
-  
-    cout<<" FIRE Fronzen Atom Number : "<<nfrozen<<endl;
-    
-
-    // Set up extra output to ion optimizer / MD header
-    cout<<"Time interval       : "<< dt*fundamentalTime*1e15<< " (fs)"<<endl;
-
-    if(step==1)
-    {
-        for(int k=0;k<numIon;k++)
-        {
-            if(ionmbl[k].x==0)vel[k].x=0;
-            if(ionmbl[k].y==0)vel[k].y=0;
-            if(ionmbl[k].z==0)vel[k].z=0;
-        }
-    //      if(nfrozen==0)RemoveMovementOfCenterOfMass();
-        scalevel();
-    }
-    moveatoms(step);
-    twiceKE=GetAtomKE();
-    twiceKE = twiceKE * 2;
-
-    tempNow = twiceKE/(3*double(numIon-nfrozen))/K_BOLTZMAN_AU;
-    cout<<" start temperature = "<< tempNow/K_BOLTZMAN_AU<< " (k)"<<endl;
-
-
-
-
-
-    // Set up forces 
-    callforce();
-    
-    for(int k=0;k<numIon;k++)
-    {
-        if(ionmbl[k].x==0)force[k].x=0;
-        if(ionmbl[k].y==0)force[k].y=0;
-        if(ionmbl[k].z==0)force[k].z=0;
-    }
-
-    //cout<<"begin maxForce"<<endl;
-    maxForce = MAXVALF();
-    cout<<"maxForce: "<<sqrt(maxForce)<<endl; 
-    mdenergy = en.etot/2;
-    conservedE = Conserved(twiceKE/2, mdenergy);
-
-    cout<<"------------------------------------------------------------------------------"<<endl;
-    cout<< "FIRE STEP "<< step<<endl;
-    cout<<"------------------------------------------------------------------------------"<<endl;
-    if (!MY_RANK){
-         ofs_running<<"------------------------------------------------------------------------------"<<endl;
-         ofs_running<< "FIRE STEP "<< step<<endl;
-         ofs_running<<"------------------------------------------------------------------------------"<<endl;
-    }
-  
-    // (1) 1st step of Verlet-Velocity
-    // New velocity are obtained
-	for(int k=0;k<numIon;k++)
-    {       
-        mass = allmass[k];
-        vel[k] = vel[k] + force[k]/mass*dt/2.0;
-    } 
-
-	//addedy by zifei
-	//initialise FIRE parameters
-	if(step1 == 1 ){
-		dt_max=2.5*dt;
-		alpha_start=0.10;
-		alpha = alpha_start;
-
-		finc=1.1;
-		fdec=0.5;
-		f_alpha=0.99;
-		N_min=4;
-		negative_count=0;
-
-		ofs_running<< "dt_max = " <<dt_max <<endl;
-		ofs_running<< "alpha_start begin = " <<alpha_start <<endl;
-		ofs_running<< "alpha begin = " <<alpha <<endl;
-		ofs_running<< "finc begin = " <<finc <<endl;
-		ofs_running<< "fdec begin = " <<fdec <<endl;
-		ofs_running<< "N_min = " <<N_min <<endl;
-		ofs_running<< "Negative count = " <<negative_count <<endl;
-	}
-
-    // (2) 2nd step of Verlet-Velocity 
-    // Update the Non-Wrapped cartesion coordinate
-    twiceKE = GetAtomKE();
-    twiceKE = 2 * twiceKE;
-    if(step!=1)for(int k=0;k<numIon;k++)
-    {       
-        mass = allmass[k];
-        vel[k] = vel[k] + force[k]/mass*dt/2.0;
-    } 
-
-	double largest_grad_FIRE = 0.0;
-	for(int i=0;i<numIon;i++)
-    {
-        if(largest_grad_FIRE < abs(force[i].x))
-        {
-            largest_grad_FIRE= abs(force[i].x);
-        } 
-        if(largest_grad_FIRE < abs(force[i].y))
-        {
-            largest_grad_FIRE= abs(force[i].y);
-        } 
-        if(largest_grad_FIRE < abs(force[i].z))
-        {
-            largest_grad_FIRE= abs(force[i].z);
-        } 
-	}
-	
-	ofs_running << " LARGEST GRAD (eV/A)  : " << largest_grad_FIRE * Ry_to_eV / 0.529177 << endl;
-
-	if(largest_grad_FIRE*Ry_to_eV/0.529177 < 0.01)
-    {
-	//"convergency reach"
-		//cout <<"CONVERGENCY REACH of FIRE in the "<<step <<" steps " <<endl;
-		timer::tick("MD_run","runFIRE",'C');
-		return true;
-	}
-
-	check_FIRE();
-    for(int i=0;i<numIon;i++)
-    {
-        cartNoWrap[i] = vel[i]*dt + cartNoWrap[i];	
-	}
-    // Calculate the maximal velocities.
-    // The step in fractional coordinates
-    maxStep = 0;
-    for(int  i = 0;i< numIon;i++)
-    {
-        if((pow(vel[i].x,2.0)+pow(vel[i].y,2.0)+pow(vel[i].z,2.0))>maxStep)
-        {
-            maxStep = pow(vel[i].x,2.0)+pow(vel[i].y,2.0)+pow(vel[i].z,2.0);
-        }
-        Mathzone::Cartesian_to_Direct(vel[i].x*dt/ucell.lat0,vel[i].y*dt/ucell.lat0,vel[i].z*dt/ucell.lat0,
-                                    ionlatvec.e11,ionlatvec.e12,ionlatvec.e13,
-                                    ionlatvec.e21,ionlatvec.e22,ionlatvec.e23,
-                                    ionlatvec.e31,ionlatvec.e32,ionlatvec.e33,
-                                    fracStep.x,fracStep.y,fracStep.z);
-        taudirac[i] = taudirac[i] + fracStep;
-    }
-
-    moveatoms(step);
-    string t("md_pos_");
-    t=intTurnTostring(step,t);
-	connection2();
-    printpos(t,step);
-    maxStep = sqrt(maxStep)*dt;
-
-    hamiltonian = Conserved(twiceKE/2, mdenergy);
-
-
-
-    
-    // Output the message to the screen.
-    oldEtot=mdenergy;
-
-    cout<<"(NVE): this step finished."<<endl;
-    if (!MY_RANK){
-        ofs_running<<"(NVE): this step finished."<<endl;
-    }
-//2015-09-25, xiaohui
-#ifdef __MPI
-    if(step%100==0){
-        MPI_Bcast(vel,numIon*3,MPI_DOUBLE,0,MPI_COMM_WORLD);
-        MPI_Bcast(cartNoWrap,numIon*3,MPI_DOUBLE,0,MPI_COMM_WORLD);
-        MPI_Bcast(taudirac,numIon*3,MPI_DOUBLE,0,MPI_COMM_WORLD);
-    }
-#endif
-
-    timer::tick("MD_run","runFIRE",'C');
-    return false;
-}
diff --git a/ABACUS.develop/source/src_pw/MD_run.h b/ABACUS.develop/source/src_pw/MD_run.h
deleted file mode 100644
index cebf606e35..0000000000
--- a/ABACUS.develop/source/src_pw/MD_run.h
+++ /dev/null
@@ -1,24 +0,0 @@
-#ifndef MD_RUN_H
-#define MD_RUN_H
-
-#include "MD_func.h"
-#include "MD_thermo.h"
-#include "MD_parameters.h"
-
-//using namespace std;
-class MD_run
-{
-	public:	
-	MD_run(){};
-	~MD_run(){};
-	void runnvt(int step1);
-	void runNVE(int step1); //NVE ensemble MD
-	bool runFIRE(int step1); //relax method FIRE
-
-	private:
-    MD_func mdf;
-    MD_thermo mdt;
-	MD_parameters &mdp;
-};
-
-#endif
\ No newline at end of file
diff --git a/ABACUS.develop/source/src_pw/MD_thermo.cpp b/ABACUS.develop/source/src_pw/MD_thermo.cpp
index fe55fcf2ff..cf1164090b 100644
--- a/ABACUS.develop/source/src_pw/MD_thermo.cpp
+++ b/ABACUS.develop/source/src_pw/MD_thermo.cpp
@@ -1,21 +1,45 @@
 #include"MD_thermo.h"
 
-MD_thermo::init_NHC(
+MD_thermo::MD_thermo()
+{
+    G = new Vector3<double>[1];
+    NHCeta = new Vector3<double>[1];
+    NHCpeta = new Vector3<double>[1];
+}
+
+MD_thermo::~MD_thermo()
+{
+    delete[] G;
+    delete[] NHCeta;
+    delete[] NHCpeta;
+}
+
+void MD_thermo::init_NHC(
     const int &MNHC_in, 
+    const double &Qmass_in, 
+    const double &NVT_tau_in, 
+    const double &dt_in,
     const int &NVT_control, 
     ofstream &ofs, 
-    const int &numIon
+    const int &numIon,
+    const double &temperature,
+    const Vector3<double>* vel,
+    const double* allmass
     )
 {
-    this->MNHC = MNHC_in;
+    this->MNHC_ = MNHC_in;
+    this->Qmass_ = Qmass_in;
+    this->NVT_tau_ = NVT_tau_in;
+    this->dt_ = dt_in;
+    this->numIon_ = numIon;
     unsigned long init[4]={0x123, 0x234, 0x345, 0x456}, length=4;
 	init_by_array(init, length);
 	ofs<<" ...............Nose-Hoover Chain parameter initialization...............  " << endl;
 	ofs<<" Temperature =    "<< temperature << endl;
 	ofs<<" Temperature2 =    "<< temperature/K_BOLTZMAN_AU << endl;
-	ofs<<" NHC frequency =    "<< 1.0/NVT_tau << endl;
-	ofs<<" NHC chain =    "<< MNHC << endl;
-	ofs<<" Qmass  =    "<< Qmass << endl;
+	ofs<<" NHC frequency =    "<< 1.0/NVT_tau_ << endl;
+	ofs<<" NHC chain =    "<< MNHC_ << endl;
+	ofs<<" Qmass  =    "<< Qmass_ << endl;
 	ofs<<" ...............................................................  " << endl;
 	w[0]=0.7845136105;
 	w[6]=0.7845136105;
@@ -27,44 +51,44 @@ MD_thermo::init_NHC(
 	//
 	for(int i=0;i<nsy;i++)
     {
-		delta[i]=w[i]*dt;
+		delta[i]=w[i]*dt_;
 	}
 		
     delete[] G;	
-	G=new Vector3<double>[MNHC*numIon];
+	G=new Vector3<double>[MNHC_*numIon_];
 	delete[] NHCeta;	
-	NHCeta=new Vector3<double>[MNHC*numIon];
+	NHCeta=new Vector3<double>[MNHC_*numIon_];
 	delete[] NHCpeta;	
-	NHCpeta=new Vector3<double>[MNHC*numIon];
+	NHCpeta=new Vector3<double>[MNHC_*numIon_];
 
-	for(int j=0;j<MNHC;j++)
+	for(int j=0;j<MNHC_;j++)
     {
-		for(int k=0;k<numIon;k++)
+		for(int k=0;k<numIon_;k++)
         {
-			NHCeta[j*numIon+k].x=genrand_real2()-0.5; 
-			NHCeta[j*numIon+k].y=genrand_real2()-0.5; 
-			NHCeta[j*numIon+k].z=genrand_real2()-0.5;
-			NHCpeta[j*numIon+k].x=gaussrand()*sqrt(Qmass * temperature);
-			NHCpeta[j*numIon+k].y=gaussrand()*sqrt(Qmass * temperature);
-			NHCpeta[j*numIon+k].z=gaussrand()*sqrt(Qmass * temperature);
+			NHCeta[j*numIon_+k].x=genrand_real2()-0.5; 
+			NHCeta[j*numIon_+k].y=genrand_real2()-0.5; 
+			NHCeta[j*numIon_+k].z=genrand_real2()-0.5;
+			NHCpeta[j*numIon_+k].x=gaussrand()*sqrt(Qmass_ * temperature);
+			NHCpeta[j*numIon_+k].y=gaussrand()*sqrt(Qmass_ * temperature);
+			NHCpeta[j*numIon_+k].z=gaussrand()*sqrt(Qmass_ * temperature);
 		}
 	}
 		
-	for(int j=MNHC-1;j>=1;j--)
+	for(int j=MNHC_-1;j>=1;j--)
     {
-		for(int k=0;k<numIon;k++)
+		for(int k=0;k<numIon_;k++)
         {
-			G[j*numIon+k].x=pow(NHCpeta[(j-1)*numIon+k].x,2)/Qmass - 1.0 * temperature;
-			G[j*numIon+k].y=pow(NHCpeta[(j-1)*numIon+k].y,2)/Qmass - 1.0 * temperature;
-			G[j*numIon+k].z=pow(NHCpeta[(j-1)*numIon+k].z,2)/Qmass - 1.0 * temperature;
+			G[j*numIon_+k].x=pow(NHCpeta[(j-1)*numIon_+k].x,2)/Qmass_ - 1.0 * temperature;
+			G[j*numIon_+k].y=pow(NHCpeta[(j-1)*numIon_+k].y,2)/Qmass_ - 1.0 * temperature;
+			G[j*numIon_+k].z=pow(NHCpeta[(j-1)*numIon_+k].z,2)/Qmass_ - 1.0 * temperature;
 		}
 	}
 		
-	for(int k=0;k<numIon;k++)
+	for(int k=0;k<numIon_;k++)
     {
-		G[0*numIon+k].x=pow(vel[k].x,2)*allmass[k]-1.0 * temperature;
-		G[0*numIon+k].y=pow(vel[k].y,2)*allmass[k]-1.0 * temperature;
-		G[0*numIon+k].z=pow(vel[k].z,2)*allmass[k]-1.0 * temperature;
+		G[0*numIon_+k].x=pow(vel[k].x,2)*allmass[k]-1.0 * temperature;
+		G[0*numIon_+k].y=pow(vel[k].y,2)*allmass[k]-1.0 * temperature;
+		G[0*numIon_+k].z=pow(vel[k].z,2)*allmass[k]-1.0 * temperature;
 	}
 	ofs << "finish NHC thermostat define" << endl; //zifei
 }
@@ -72,7 +96,6 @@ MD_thermo::init_NHC(
 double MD_thermo::NHChamiltonian(
     const double &KE,
     const double &PE,
-    const int &numIon,
     const double &temperature,
     const int &nfrozen )
 {
@@ -81,32 +104,23 @@ double MD_thermo::NHChamiltonian(
 	NHChamiltonian0 = KE + PE ;
 //	cout << "hamiltonian0 =    "<< NHChamiltonian0 <<endl;
 
-	for(int i=0;i<numIon;i++)
+	for(int i=0;i<numIon_;i++)
     {
-		for(int j=0;j<MNHC;j++)
+		for(int j=0;j<MNHC_;j++)
         {
-			NHChamiltonian0 += pow(NHCpeta[i*MNHC+j].x,2)/2.0/Qmass+temperature*NHCeta[i*MNHC+j].x;
-			NHChamiltonian0 += pow(NHCpeta[i*MNHC+j].y,2)/2.0/Qmass+temperature*NHCeta[i*MNHC+j].y;
-			NHChamiltonian0 += pow(NHCpeta[i*MNHC+j].z,2)/2.0/Qmass+temperature*NHCeta[i*MNHC+j].z;
+			NHChamiltonian0 += pow(NHCpeta[i*MNHC_+j].x,2)/2.0/Qmass_+temperature*NHCeta[i*MNHC_+j].x;
+			NHChamiltonian0 += pow(NHCpeta[i*MNHC_+j].y,2)/2.0/Qmass_+temperature*NHCeta[i*MNHC_+j].y;
+			NHChamiltonian0 += pow(NHCpeta[i*MNHC_+j].z,2)/2.0/Qmass_+temperature*NHCeta[i*MNHC_+j].z;
 		}
 	}
 
 	if (!MY_RANK)
     {
-		out<< " "<<endl;
-        out<< " --------------------------------------------------"<<endl;
-        out<< " SUMMARY OF NVT CALCULATION"<<endl;
-        out<<" --------------------------------------------------"<<endl;
-        out<<" NVT Conservation     : "<<setw(10)<< NHChamiltonian0*2<<" (Rydberg)"<<endl;
-        out<<" NVT Temperature      : "<<setw(10)<< KE*2/(3*double(numIon-nfrozen))/K_BOLTZMAN_AU<<" (K)"<<endl;
-        out<<" NVT Kinetic energy   : "<<setw(10)<< KE*2<<" (Rydberg)"<<endl;
-        out<<" NVT Potential energy : "<<setw(10)<< PE*2<<" (Rydberg)"<<endl;
-        ofs_running<< " "<<endl;
         ofs_running<< " --------------------------------------------------"<<endl;
         ofs_running<< " SUMMARY OF NVT CALCULATION"<<endl;
         ofs_running<<" --------------------------------------------------"<<endl;
         ofs_running<<" NVT Conservation     : "<<setw(10)<< NHChamiltonian0*2<<" (Rydberg)"<<endl;
-        ofs_running<<" NVT Temperature      : "<<setw(10)<< KE*2/(3*double(numIon-nfrozen))/K_BOLTZMAN_AU<<" (K)"<<endl;
+        ofs_running<<" NVT Temperature      : "<<setw(10)<< KE*2/(3*double(numIon_-nfrozen))/K_BOLTZMAN_AU<<" (K)"<<endl;
         ofs_running<<" NVT Kinetic energy   : "<<setw(10)<< KE*2<<" (Rydberg)"<<endl;
         ofs_running<<" NVT Potential energy : "<<setw(10)<< PE*2<<" (Rydberg)"<<endl;
     }
@@ -272,16 +286,24 @@ double MD_thermo::genrand_res53(void)
 } 
 
 
-void MD_thermo::Integrator(const int control)
+void MD_thermo::Integrator(
+    const int control,
+    const double &temperature,
+    Vector3<double>* vel,
+    const double* allmass)
 {
-	if(control == 1) NHCIntegrator();
-	else if(control == 2) LGVIntegrator();
-	else if(control == 3) ADSIntegrator();
+	if(control == 1) NHCIntegrator(temperature, vel, allmass);
+	else if(control == 2) LGVIntegrator(temperature, vel, allmass);
+	else if(control == 3) ADSIntegrator(temperature, vel, allmass);
 	else WARNING_QUIT("MD_thermo:Integrator", "please choose available reservoir!!!");
 	return;
 }
 
-void MD_thermo::LGVIntegrator()
+void MD_thermo::LGVIntegrator(
+    const double &temperature,
+    Vector3<double>* vel,
+    const double* allmass
+)
 {
 //---------------------------------------------------------------------------
 // DESCRIPTION:
@@ -289,11 +311,11 @@ void MD_thermo::LGVIntegrator()
 //--------------------------------------------------------------------
 	double randomx,randomy,randomz;
 	double tempV;
-	double c1k = exp(-1/NVT_tau * dt);
+	double c1k = exp(-1/NVT_tau_ * dt_);
 	//c1k=e^(-gamma*dt)
 	double c2k=sqrt(1.0-c1k*c1k);
 
-    for(int iatom=0;iatom<numIon;iatom++)
+    for(int iatom=0;iatom<numIon_;iatom++)
     {
         randomx = gaussrand();
         tempV = vel[iatom].x;
@@ -314,7 +336,11 @@ void MD_thermo::LGVIntegrator()
 
 
 //added by zifei
-void MD_thermo::ADSIntegrator()
+void MD_thermo::ADSIntegrator(
+    const double &temperature,
+    Vector3<double>* vel,
+    const double* allmass
+)
 {
 //---------------------------------------------------------------------------
 // DESCRIPTION:
@@ -326,13 +352,13 @@ void MD_thermo::ADSIntegrator()
     double rany;
     double ranz;
 
-	double nu = 1/ NVT_tau;
+	double nu = 1/ NVT_tau_;
 
-    for(int iatom=0;iatom<numIon;iatom++)
+    for(int iatom=0;iatom<numIon_;iatom++)
     {
         uniform_random = genrand_real2() ;
 
-        if(uniform_random < 1.0-exp(-nu*dt) ) 
+        if(uniform_random < 1.0-exp(-nu*dt_) ) 
         {
         //if(uniform_random < 0.3 ) {
             ranx=sqrt(allmass[iatom] * temperature);
@@ -359,7 +385,11 @@ void MD_thermo::ADSIntegrator()
 
 
 //zifei
-void MD_thermo::NHCIntegrator(){
+void MD_thermo::NHCIntegrator(
+    const double &temperature,
+    Vector3<double>* vel,
+    const double* allmass
+){
 //---------------------------------------------------------------------------
 // DESCRIPTION:
 //   This function propagates the Nose-Hoover Chain
@@ -368,85 +398,212 @@ void MD_thermo::NHCIntegrator(){
 //   Q : NHC mass
 //--------------------------------------------------------------------
 
-    for(int i=0;i<numIon;i++)
+    for(int i=0;i<numIon_;i++)
     {
-        G[0*numIon+i].x=pow(vel[i].x,2)*allmass[i]-1.0 * temperature;
-        G[0*numIon+i].y=pow(vel[i].y,2)*allmass[i]-1.0 * temperature;
-        G[0*numIon+i].z=pow(vel[i].z,2)*allmass[i]-1.0 * temperature;
+        G[0*numIon_+i].x=pow(vel[i].x,2)*allmass[i]-1.0 * temperature;
+        G[0*numIon_+i].y=pow(vel[i].y,2)*allmass[i]-1.0 * temperature;
+        G[0*numIon_+i].z=pow(vel[i].z,2)*allmass[i]-1.0 * temperature;
     }
   
     for(int i=0;i<nsy;i++)
     {
-        for(int j=0;j<numIon;j++)
+        for(int j=0;j<numIon_;j++)
         {
-            NHCpeta[(MNHC-1)*numIon+j] += G[(MNHC-1)*numIon+j]*delta[i]/2.0;
+            NHCpeta[(MNHC_-1)*numIon_+j] += G[(MNHC_-1)*numIon_+j]*delta[i]/2.0;
         }
 
-        for(int j=MNHC-2;j>=0;j--)
+        for(int j=MNHC_-2;j>=0;j--)
         {
-            for(int k=0;k<numIon;k++)
+            for(int k=0;k<numIon_;k++)
             {
-                NHCpeta[j*numIon+k].x *=exp(-NHCpeta[(j+1)*numIon+k].x/Qmass*delta[i]/4.0);
-                NHCpeta[j*numIon+k].y *=exp(-NHCpeta[(j+1)*numIon+k].y/Qmass*delta[i]/4.0);
-                NHCpeta[j*numIon+k].z *=exp(-NHCpeta[(j+1)*numIon+k].z/Qmass*delta[i]/4.0);
+                NHCpeta[j*numIon_+k].x *=exp(-NHCpeta[(j+1)*numIon_+k].x/Qmass_*delta[i]/4.0);
+                NHCpeta[j*numIon_+k].y *=exp(-NHCpeta[(j+1)*numIon_+k].y/Qmass_*delta[i]/4.0);
+                NHCpeta[j*numIon_+k].z *=exp(-NHCpeta[(j+1)*numIon_+k].z/Qmass_*delta[i]/4.0);
 
-                NHCpeta[j*numIon+k]+=G[j*numIon+k]*delta[i]/2.0;
+                NHCpeta[j*numIon_+k]+=G[j*numIon_+k]*delta[i]/2.0;
 
-                NHCpeta[j*numIon+k].x *=exp(-NHCpeta[(j+1)*numIon+k].x/Qmass*delta[i]/4.0);
-                NHCpeta[j*numIon+k].y *=exp(-NHCpeta[(j+1)*numIon+k].y/Qmass*delta[i]/4.0);
-                NHCpeta[j*numIon+k].z *=exp(-NHCpeta[(j+1)*numIon+k].z/Qmass*delta[i]/4.0);
+                NHCpeta[j*numIon_+k].x *=exp(-NHCpeta[(j+1)*numIon_+k].x/Qmass_*delta[i]/4.0);
+                NHCpeta[j*numIon_+k].y *=exp(-NHCpeta[(j+1)*numIon_+k].y/Qmass_*delta[i]/4.0);
+                NHCpeta[j*numIon_+k].z *=exp(-NHCpeta[(j+1)*numIon_+k].z/Qmass_*delta[i]/4.0);
             }
         }//MNHC-2
   	
-        for(int j=MNHC-1;j>=0;j--)
+        for(int j=MNHC_-1;j>=0;j--)
         {
-            for(int k=0;k<numIon;k++)
+            for(int k=0;k<numIon_;k++)
             {
-                NHCeta[j*numIon+k].x += NHCpeta[j*numIon+k].x/Qmass*delta[i];
-                NHCeta[j*numIon+k].y += NHCpeta[j*numIon+k].y/Qmass*delta[i];
-                NHCeta[j*numIon+k].z += NHCpeta[j*numIon+k].z/Qmass*delta[i];
+                NHCeta[j*numIon_+k].x += NHCpeta[j*numIon_+k].x/Qmass_*delta[i];
+                NHCeta[j*numIon_+k].y += NHCpeta[j*numIon_+k].y/Qmass_*delta[i];
+                NHCeta[j*numIon_+k].z += NHCpeta[j*numIon_+k].z/Qmass_*delta[i];
             }
         }
         //p=p*dexp(-peta(1,:)/Q(1,:)*delta(i))
-        for(int k=0;k<numIon;k++)
+        for(int k=0;k<numIon_;k++)
         {
-            vel[k].x *= exp(-NHCpeta[0*numIon+k].x/Qmass*delta[i]);
-            vel[k].y *= exp(-NHCpeta[0*numIon+k].y/Qmass*delta[i]);
-            vel[k].z *= exp(-NHCpeta[0*numIon+k].z/Qmass*delta[i]);
+            vel[k].x *= exp(-NHCpeta[0*numIon_+k].x/Qmass_*delta[i]);
+            vel[k].y *= exp(-NHCpeta[0*numIon_+k].y/Qmass_*delta[i]);
+            vel[k].z *= exp(-NHCpeta[0*numIon_+k].z/Qmass_*delta[i]);
         }
 
   	//G(1,:)=p**2/mass-1/beta
-  	    for(int k=0;k<numIon;k++){
-            G[0*numIon+k].x=pow(vel[0*numIon+k].x,2)*allmass[k]-1.0 * temperature;
-            G[0*numIon+k].y=pow(vel[0*numIon+k].y,2)*allmass[k]-1.0 * temperature;
-            G[0*numIon+k].z=pow(vel[0*numIon+k].z,2)*allmass[k]-1.0 * temperature;
+  	    for(int k=0;k<numIon_;k++){
+            G[0*numIon_+k].x=pow(vel[0*numIon_+k].x,2)*allmass[k]-1.0 * temperature;
+            G[0*numIon_+k].y=pow(vel[0*numIon_+k].y,2)*allmass[k]-1.0 * temperature;
+            G[0*numIon_+k].z=pow(vel[0*numIon_+k].z,2)*allmass[k]-1.0 * temperature;
   	    }
   		
-        for(int j=0;j<MNHC-1;j++)
+        for(int j=0;j<MNHC_-1;j++)
         {
-            for(int k=0;k<numIon;k++)
+            for(int k=0;k<numIon_;k++)
             {
-                NHCpeta[j*numIon+k].x *=exp(-NHCpeta[(j+1)*numIon+k].x/Qmass*delta[i]/4.0);
-                NHCpeta[j*numIon+k].y *=exp(-NHCpeta[(j+1)*numIon+k].y/Qmass*delta[i]/4.0);
-                NHCpeta[j*numIon+k].z *=exp(-NHCpeta[(j+1)*numIon+k].z/Qmass*delta[i]/4.0);
+                NHCpeta[j*numIon_+k].x *=exp(-NHCpeta[(j+1)*numIon_+k].x/Qmass_*delta[i]/4.0);
+                NHCpeta[j*numIon_+k].y *=exp(-NHCpeta[(j+1)*numIon_+k].y/Qmass_*delta[i]/4.0);
+                NHCpeta[j*numIon_+k].z *=exp(-NHCpeta[(j+1)*numIon_+k].z/Qmass_*delta[i]/4.0);
 
-                NHCpeta[j*numIon+k]+=G[j*numIon+k]*delta[i]/2.0;
+                NHCpeta[j*numIon_+k]+=G[j*numIon_+k]*delta[i]/2.0;
 
-                NHCpeta[j*numIon+k].x *=exp(-NHCpeta[(j+1)*numIon+k].x/Qmass*delta[i]/4.0);
-                NHCpeta[j*numIon+k].y *=exp(-NHCpeta[(j+1)*numIon+k].y/Qmass*delta[i]/4.0);
-                NHCpeta[j*numIon+k].z *=exp(-NHCpeta[(j+1)*numIon+k].z/Qmass*delta[i]/4.0);
+                NHCpeta[j*numIon_+k].x *=exp(-NHCpeta[(j+1)*numIon_+k].x/Qmass_*delta[i]/4.0);
+                NHCpeta[j*numIon_+k].y *=exp(-NHCpeta[(j+1)*numIon_+k].y/Qmass_*delta[i]/4.0);
+                NHCpeta[j*numIon_+k].z *=exp(-NHCpeta[(j+1)*numIon_+k].z/Qmass_*delta[i]/4.0);
 
-                G[(j+1)*numIon+k].x = pow(NHCpeta[j*numIon+k].x,2)/Qmass-1.0 * temperature;
-                G[(j+1)*numIon+k].y = pow(NHCpeta[j*numIon+k].y,2)/Qmass-1.0 * temperature;
-                G[(j+1)*numIon+k].z = pow(NHCpeta[j*numIon+k].z,2)/Qmass-1.0 * temperature;
+                G[(j+1)*numIon_+k].x = pow(NHCpeta[j*numIon_+k].x,2)/Qmass_-1.0 * temperature;
+                G[(j+1)*numIon_+k].y = pow(NHCpeta[j*numIon_+k].y,2)/Qmass_-1.0 * temperature;
+                G[(j+1)*numIon_+k].z = pow(NHCpeta[j*numIon_+k].z,2)/Qmass_-1.0 * temperature;
             }
         }
 
-        for(int k=0;k<numIon;k++)
+        for(int k=0;k<numIon_;k++)
         {
-            NHCpeta[(MNHC-1)*numIon+k]+=G[(MNHC-1)*numIon+k]*delta[i]/2.0;
+            NHCpeta[(MNHC_-1)*numIon_+k]+=G[(MNHC_-1)*numIon_+k]*delta[i]/2.0;
         }			
     }//nsy
 
+    return;
+}
+
+void MD_thermo::NHC_info_out(const int& step, const int& recordFreq, const int& numIon)
+{
+    bool pass;
+	pass = 0;
+
+	if (recordFreq==1||step==1||( recordFreq > 1&& step%recordFreq==0 ) )
+		pass =1;
+	if (!pass) return;
+
+	if(!MY_RANK){
+		stringstream ssc;
+		ssc << global_out_dir << "Restart_md.dat";
+		ofstream file(ssc.str().c_str(), ios::app);
+        file<<'\n';
+		file<<"MD_THERMOSTAT"<<endl;
+		file<<"MNHC: "<<MNHC_<<endl;
+		file<<"G: "<<endl;
+		for(int i=0;i<numIon_*MNHC_;i++){
+			file<<setprecision (12)<<G[i].x<<" "<<setprecision (12)<<G[i].y<<" "<<setprecision (12)<<G[i].z<<endl;
+		}
+        file<<"NHCeta: "<<endl;
+		for(int i=0;i<numIon_*MNHC_;i++){
+			file<<setprecision (12)<<NHCeta[i].x<<" "<<setprecision (12)<<NHCeta[i].y<<" "<<setprecision (12)<<NHCeta[i].z<<endl;
+		}
+        file<<"NHCpeta: "<<endl;
+		for(int i=0;i<numIon_*MNHC_;i++){
+			file<<setprecision (12)<<NHCpeta[i].x<<" "<<setprecision (12)<<NHCpeta[i].y<<" "<<setprecision (12)<<NHCpeta[i].z<<endl;
+		}                
+		file.close();
+	}
+	return;
+    return;
+}
+
+#include<cstring>
+void MD_thermo::NHC_restart()
+{
+    int error=0;
+    double *nhcg=new double[numIon_*3*MNHC_];
+    double *nhcp=new double[numIon_*3*MNHC_];
+    double *nhce=new double[numIon_*3*MNHC_];
+    if (!MY_RANK)
+	{
+		stringstream ssc;
+		ssc << global_readin_dir << "Restart_md.dat";
+		ifstream file(ssc.str().c_str());
+
+        char word[80];
+        int mnhc;
+        while(file.good())//search and read MNHC
+        {
+            file>>word;
+            if(std::strcmp("MNHC:",word)==0)
+            {
+                file>>mnhc;
+                break;
+            }
+            else
+            {
+                file.ignore(150, '\n');
+            }
+        }
+
+        if(mnhc!=MNHC_)
+        {
+            cout<<"please ensure whether 'Restart_md.dat' right!"<<endl;
+            error = 1;
+        }
+
+        if(!error)
+        {
+            file.get();
+            file.ignore(3, '\n');//read G
+            for(int i = 0;i<numIon_*3*MNHC_;i++)
+            {
+                file>>nhcg[i];
+            }
+            file.get();
+            file.ignore(8, '\n');//read eta
+            for(int i = 0;i<numIon_*3*MNHC_;i++)
+            {
+                file>>nhce[i];
+            }
+            file.get();
+            file.ignore(9, '\n');//read peta
+            for(int i = 0;i<numIon_*3*MNHC_;i++)
+            {
+                file>>nhcp[i];
+            }
+            file.close();
+        }
+
+    }
+#ifdef __MPI
+    MPI_Bcast(&error,1,MPI_INT,0,MPI_COMM_WORLD);
+#endif
+    if(error)
+    {
+        delete[] nhcg;
+        delete[] nhcp;
+        delete[] nhce;
+        exit(0);
+    }
+#ifdef __MPI
+	MPI_Bcast(nhcg,numIon_*3*MNHC_,MPI_DOUBLE,0,MPI_COMM_WORLD);
+	MPI_Bcast(nhcp,numIon_*3*MNHC_,MPI_DOUBLE,0,MPI_COMM_WORLD);
+    MPI_Bcast(nhce,numIon_*3*MNHC_,MPI_DOUBLE,0,MPI_COMM_WORLD);
+#endif
+
+    for(int i=0;i<numIon_*MNHC_;i++)
+    {
+        G[i].x = nhcg[i*3];
+        G[i].y = nhcg[i*3+1];
+        G[i].z = nhcg[i*3+2];
+        NHCpeta[i].x = nhcp[i*3];
+        NHCpeta[i].y = nhcp[i*3+1];
+        NHCpeta[i].z = nhcp[i*3+2];
+        NHCeta[i].x = nhce[i*3];
+        NHCeta[i].y = nhce[i*3+1];
+        NHCeta[i].z = nhce[i*3+2];
+    }
+    delete[] nhcg;
+    delete[] nhcp;
+    delete[] nhce;
     return;
 }
\ No newline at end of file
diff --git a/ABACUS.develop/source/src_pw/MD_thermo.h b/ABACUS.develop/source/src_pw/MD_thermo.h
index 18329f4f04..6a3cb286b7 100644
--- a/ABACUS.develop/source/src_pw/MD_thermo.h
+++ b/ABACUS.develop/source/src_pw/MD_thermo.h
@@ -1,30 +1,58 @@
 #ifndef MD_THERMO_H
 #define MD_THERMO_H
+#include "tools.h"
 
 class MD_thermo 
 {
     public:
     MD_thermo();
-    ~MD_thrmo();
+    ~MD_thermo();
 
-    void Integrator(const int control);
+    void Integrator(
+        const int control,
+        const double &temperature,
+        Vector3<double>* vel,
+        const double* allmass
+        );
     void init_NHC(
         const int &MNHC_in, 
+        const double &Qmass_in, 
+        const double &NVT_tau_in, 
+        const double &dt_in,
         const int &NVT_control, 
         ofstream &ofs, 
-        const int &numIon
+        const int &numIon,
+        const double &temperature,
+        const Vector3<double>* vel,
+        const double* allmass
         );
     double NHChamiltonian(
         const double &KE,
         const double &PE,
-        const int &numIon,
         const double &temperature,
-        const int &nfrozen );
+        const int &nfrozen
+        );
+
+    void NHC_info_out(const int& step, const int& recordFreq, const int& numIon);
+    void NHC_restart();
 
     private:    
-    void ADSIntegrator();
-	void LGVIntegrator();
-	void NHCIntegrator();
+    void ADSIntegrator(
+        const double &temperature,
+        Vector3<double>* vel,
+        const double* allmass
+    );
+	void LGVIntegrator(
+        const double &temperature,
+        Vector3<double>* vel,
+        const double* allmass
+    );
+	void NHCIntegrator(
+        const double &temperature,
+        Vector3<double>* vel,
+        const double* allmass
+    );
+
 
 	double gaussrand();
 	
@@ -35,7 +63,7 @@ class MD_thermo
 	double genrand_real1(void);
 	double genrand_real2(void);
 	double genrand_real3(void);
-	double genrand_res53(void) ;
+	double genrand_res53(void);
 
     //NVT thermostat parameters
 	const static int N_mt19937=624;
@@ -47,19 +75,24 @@ class MD_thermo
 	int mti=625; /* mti==N+1 means mt[N] is not initialized */
 	long double gamma; //langevin friction coefficient
 	long double nu; //Andersen collision frequency
-	long long double c1k; //parameter in Langevin
+	long double c1k; //parameter in Langevin
 	double c2k; //parameter in Langevin
     const static int nsy=7; //parameter in NHC, constant, no need to modification
-    long double w[nsy]; //parameters in NHC
-    long double delta[nsy]; //parameters in NHC
+    double w[nsy]; //parameters in NHC
+    double delta[nsy]; //parameters in NHC
 
     //input parameters
-    int MNHC;
-    int Qmass;
+    int MNHC_;
+    double Qmass_;
+    double NVT_tau_;
+    double dt_;
+    int numIon_;
 
     //need to be allocated
     Vector3<double> *G; //parameter in NHC
     Vector3<double> *NHCeta; //NHC position
     Vector3<double> *NHCpeta; //NHC momentum
 
-};
\ No newline at end of file
+};
+
+#endif
\ No newline at end of file
diff --git a/ABACUS.develop/source/src_pw/md.cpp b/ABACUS.develop/source/src_pw/md.cpp
deleted file mode 100644
index be99fd8763..0000000000
--- a/ABACUS.develop/source/src_pw/md.cpp
+++ /dev/null
@@ -1,1211 +0,0 @@
-//Description:this module is the base of MD simulation
-#include<iostream>
-#include<string>
-#include<fstream>
-#include<stdlib.h>
-#include<time.h>
-#include<cmath>
-#include"md.h"
-#include"../input.h"
-#include"../src_lcao/LOOP_ions.h"
-#include "../src_global/sltk_atom_arrange.h" //2015-10-01, xiaohui
-
-md::md(int n)
-{
-	xLogS=0; // thermostat's position
-	vLogS=0 ;// thermostat's speed
-	tauthermo = 0;  // The thermostat fluctuation period, for Al fcc at 300K ~ 1000fs
-	taubaro   = 0 ; // The barostat   fluctuation period, for Al fcc at 300K ~ 500 fs
-	Qmass      = -1;  // Inertia of extended system variable
-	dt         = -1; // Time increment (hbar/E_hartree)
-	temperature= -1;  //Temperature (in Hartree, 1 Hartree ~ 3E5 K)
-	doMSD = 1;   // compute the mean square displacement.
-	doMSDAtom = 1;        // compute the mean square displacement for each atom.
-	msdstartTime = 1;   //beyond this time, msd is going to be computed
-	diffuCoeff    = 0;
-	nResn=3;   // Coarseness of Trotter factorization of Nose-Hoover thermostat
-	nYosh=3;  // Order of thermostat scheme corresponding with "intCoeff"
-	outputstressperiod = 1 ;// The period to output stress
-	dumpmdfreq = 1  ;       // The period to dump MD information for monitoring and restarting MD
-	fixTemperature = 1;       // The period to read in new temperature during MD run.
-	rstMD = 0  ;// 1 : restart MD, vel.restart and ion.restart files will be read
-	// 0 : not restart from ion and vel files
-	// -n: restart from vel.n and ion.n in directory
-	// "md_output_path"
-
-	startStep=0 ;// the first step number of MD
-	velRescale = 0 ;// whether to rescale velocities at the first step
-	mdoutputpath="mdDataOut";
-	step_rst=0;
-	numIon=ucell.nat;
-	ntype=ucell.ntype;
-};
-
-void md::md_allocate()
-{
-	na=new int[ntype];
-	if (!MY_RANK)
-	{
-		stringstream ssc;
-		ssc << global_out_dir << "MD_OUT";
-		out.open(ssc.str().c_str());
-	}
-	for(int i=0;i<ntype;i++)
-	{
-		na[i]=ucell.atoms[i].na;
-	} 
-	ionlatvec=ucell.latvec;
-	msd=new double[ntype];
-	vel=new Vector3<double>[numIon];
-	cartNoWrap=new Vector3<double>[numIon];
-	taudirac=new Vector3<double>[numIon];
-	allmass=new double[numIon];
-	ionmbl=new Vector3<int>[numIon];
-	force=new Vector3<double>[numIon];
-	return;	
-}
-
-const double md::fundamentalTime = 2.418884326505e-17;
-
-void md::initMD(void)
-{
-
-	mdtype=INPUT.md_mdtype;
-
-	mdenergy=en.etot/2;
-
-	Qmass=INPUT.md_qmass/6.02/9.109*1e5;
-	dt=INPUT.md_dt/fundamentalTime/1e15;
-	temperature=INPUT.md_tfirst/3.1577464e5;
-	rstMD=INPUT.md_rstmd;
-
-	nResn=INPUT.md_nresn;
-	nYosh=INPUT.md_nyosh;
-	dumpmdfreq=INPUT.md_dumpmdfred;
-	fixTemperature=INPUT.md_fixtemperature;
-	ediff=INPUT.md_ediff;
-	ediffg=INPUT.md_ediffg;
-	doMSD=INPUT.md_domsd;
-	doMSDAtom=INPUT.md_domsdatom;
-	mdoutputpath=INPUT.md_mdoutpath;
-	msdstartTime=INPUT.md_msdstartTime;
-
-	if(rstMD==0)
-	{
-		connection0();
-		connection1();
-		//A fresh new MD: Do not restart MD
-		InitVelocity() ;
-		// Initialize thermostat, and barostat
-		xLogS = 0.0 ;      // position of thermostat
-		vLogS = 0.0;        // velocity of thermostat
-	}
-	else if( rstMD>0 )
-	{
-		connection0();
-		connection1();
-		InitVelocity() ;
-		if(!RestartMD()){
-			cout<<"error in restart MD!"<<endl;
-			return;
-		}
-		connection2();
-	}
-	if(Qmass==0)
-	{
-		Qmass=dt*fundamentalTime*1e15/6.02/9.109*1e5;
-	}
-
-	return;
-}
-	
-/*void md::runMD(int step1){
-       if(rstMD==1)MDNVT.runnvt();
-       else if(mdtype==3)MDNVE.runNVE(step1);
-       return;
-}*///we need this form
-
-bool md::RestartMD(void)
-{
-	int i;
-	double *vell=new double[numIon*3];
-	double *cart=new double[numIon*3];
-	double *dira=new double[numIon*3];
-
-	if (!MY_RANK)
-	{
-		stringstream ssc;
-		ssc << global_out_dir << "Restart_md.dat";
-		ifstream file(ssc.str().c_str());
-
-		if(!file)
-		{
-			cout<<"please ensure whether 'Restart_md.dat' exists"<<endl;
-			exit(0);
-		}
-	    if (file.good())
-	    {
-	        file.ignore(11, '\n');
-	    }
-//----------------------------------------------------------
-// main parameters
-//----------------------------------------------------------
-		file.ignore(7, '\n');
-		file>>xLogS;
-		file.get();	file.ignore(7, '\n');
-		file>>vLogS;
-		file.get();	file.ignore(23, '\n');
-		for(i = 0;i<numIon*3;i++)
-		{
-			file>>vell[i];
-		}
-		/*file.get();	file.ignore(17, '\n');
-		  for(i=0;i<numIon*3;i++){
-		  file>>cart[i];
-		  }
-		  file.get();	file.ignore(13, '\n');
-		  for(i=0;i<numIon*3;i++){
-		  file>>dira[i];
-		  }*/
-		file.get();     
-		file.ignore(6, '\n');
-		file>>step_rst;
-		file.close();
-	}
-
-//2015-09-05, xiaohui
-#ifdef __MPI
-	MPI_Bcast(&xLogS,1,MPI_DOUBLE,0,MPI_COMM_WORLD);
-	MPI_Bcast(&vLogS,1,MPI_DOUBLE,0,MPI_COMM_WORLD);
-	MPI_Bcast(&step_rst,1,MPI_INT,0,MPI_COMM_WORLD);
-	MPI_Bcast(vell,numIon*3,MPI_DOUBLE,0,MPI_COMM_WORLD);
-	//      MPI_Bcast(cart,numIon*3,MPI_DOUBLE,0,MPI_COMM_WORLD);
-	//      MPI_Bcast(dira,numIon*3,MPI_DOUBLE,0,MPI_COMM_WORLD);
-#endif
-
-	for(i=0;i<numIon;i++)
-	{
-		vel[i].x=vell[i*3];
-		vel[i].y=vell[i*3+1];
-		vel[i].z=vell[i*3+2];
-		//     cartNoWrap[i].x=cart[i*3];
-		//     cartNoWrap[i].y=cart[i*3+1];
-		//     cartNoWrap[i].z=cart[i*3+2];
-		//     taudirac[i].x=dira[i*3];
-		//     taudirac[i].y=dira[i*3+1];
-		//     taudirac[i].z=dira[i*3+2];
-	}
-
-	// DEBUG
-	//int debug = 0;
-	//if (debug) {
-	//  cout<<" dump cell + ion infomation file to -> echo_md_restart.dat, for debug purpose. "<<endl;
-	//	fstream debugfile;
-	//	debugfile.open("echo_md_restart.dat");
-	//  debugfile<<"  "<<endl;
-	//debugfile<<"cellReal[0] = "<<a.latvec.e11<<" "<<a.latvec.e12<<" "<<a.latvec.e13<<endl;
-	//debugfile<<"cellReal[1] = "<< a.latvec.e21<<" "<<a.latvec.e22<<" "<<a.latvec.e23<<endl;
-	//debugfile<<"cellReal[2] = "<< a.latvec.e31<<" "<<a.latvec.e32<<" "<<a.latvec.e33<<endl;
-	//debugfile<<"ION FRACTIONAL COORDINATES:"<<endl;
-	//for(ii=0;ii<numIon;ii++){
-	//  debugfile<<taudirac[ii].x<<" "<<taudirac[ii].y<<" "<<taudirac[ii].z<<endl;
-	//}
-	//debugfile<<"ION VELOCITIES:"<<endl;
-	//for( ii=0;ii<ionIon;ii++){
-	// debugfile<<vel[ii].x<<" "<<vel[ii].y<<" "<<vel[ii].z<<endl;
-	//	}
-	//  debugfile.close();
-	//}
-  	// END DEBUG
-
-	if(vell!=NULL)delete []vell;
-	if(cart!=NULL)delete []cart;
-	if(dira!=NULL)delete []dira;
-
-	return true;
-}
-
-
-void md::mstout(int step)
-{
-	//this function used for outputting the information of restart file
-	bool pass;
-	pass = 0;
-
-
-	if (dumpmdfreq==1||step==1||( dumpmdfreq > 1&& step%dumpmdfreq==0 ) )
-		pass =1;
-	if (!pass) return;
-
-	int i;
-	if(!MY_RANK)
-	{
-		stringstream ssc;
-		ssc << global_out_dir << "Restart_md.dat";
-		ofstream file(ssc.str().c_str());
-		file<<"MD_RESTART"<<endl;
-		file<<"xLogS: "<<setprecision (12)<<xLogS<<endl;
-		file<<"vLogS: "<<setprecision (12)<<vLogS<<endl;
-		file<<"ION_VELOCITIES_(a.u.): "<<endl;
-		for( i=0;i<numIon;i++)
-		{
-			file<<setprecision (12)<<vel[i].x
-			<<" "<<setprecision (12)
-			<<vel[i].y<<" "<<setprecision (12)<<vel[i].z<<endl;
-		}
-
-		/*	file<<"CARTESIAN_COORD: "<<endl;
-			for(i=0;i<numIon;i++){
-			file<<cartNoWrap[i].x<<" "<<cartNoWrap[i].y<<" "<<cartNoWrap[i].z<<endl;
-			}
-			file<<"DIRAC_COORD: "<<endl;
-			for(i=0;i<numIon;i++){
-			file<<taudirac[i].x<<" "<<taudirac[i].y<<" "<<taudirac[i].z<<endl;
-			}*/
-		file<<"step: "<<step<<endl;                
-		file.close();
-	}
-
-	return;
-}
-
-double md::GetAtomKE()
-{
-//---------------------------------------------------------------------------
-// DESCRIPTION:
-//   This function calculates the classical kinetic energy of a group of atoms.
-//----------------------------------------------------------------------------
-//!>> EXTERNAL VARIABLES <<!
-//vel[3][size(cell%ionTable,1)]  ! Velocities of atoms
-// ke ! kinetic energy we want.
-
-//!>> INTERIOR VARIABLES <<!
-	double mass;
-	int it,ii,ion;
-
-	//!>> INITIALIZATION <<!
-	//!>> FUNCTION BODY <<!
-	double ke = 0.0;
-
-	//  WRITE(outputUnit,*) "velocities"
-	// kinetic energy = \sum_{i} 1/2*m_{i}*(vx^2+vy^2+vz^2)
-	for(ion=0;ion<numIon;ion++){
-		mass = allmass[ion];
-		ke +=   0.5 * mass * (vel[ion].x*vel[ion].x+vel[ion].y*vel[ion].y+vel[ion].z*vel[ion].z);
-		//    WRITE(outputUnit,*) vel(1,ion), vel(2,ion), vel(3,ion)
-	}
-	//cout<<"in GetAtomKE KE="<< ke<<endl;
-	return ke;
-}
-
-
-void md::MakeIntCoeff()
-{
-//---------------------------------------------------------------------------
-// DESCRIPTION:
-//  Initialize the Yoshida/Suzuki coefficient
-//----------------------------------------------------------------------------
-	// nYosh ! approximation order
-	int i;
-	switch (nYosh)
-	{
-		case 1:
-			intCoeff[0] = 1.0;
-			break;
-		case 3:
-			for(i=0;i<3;i++)intCoeff[i] = 1.0/(2.0-pow(2,1/3.0));
-			intCoeff[1] = 1.0 - 2.0*intCoeff[1];
-			break;
-		case 5:
-			for(i=0;i<5;i++)intCoeff[i] = 1.0/(4.0-pow(4,1/3.0));
-			intCoeff[2] = 1.0 - 4.0*intCoeff[2];
-			break;
-		default:
-			cout<<"ERROR: NYOSH MUST BE 1, 3, or 5.  STOPPING."<<endl;
-			exit(0);
-	}
-
-	return;
-}
-
-
-//  double md::InitRandomSeed(){
-//---------------------------------------------------------------------------
-// DESCRIPTION:
-//  Initialize the random seed used for random number generator.
-//----------------------------------------------------------------------------
-
-  //}//don't know whether it is necessary,so I neglect it.
-void md::InitVelocity(void)
-{
-	if(!MY_RANK)
-	{ //xiaohui add 2015-09-25
-		srand(time(0));
-		int iy=rand()%21,im=rand()%11,id=rand()%21,ih=rand()%31,in=rand()%61,is=rand()%41;
-		int jy=rand()%21,jm=rand()%11,jd=rand()%21,jh=rand()%31,jn=rand()%61,js=rand()%41;
-		//    int iy=13,im=10,id=1,ih=21,in=54,is=23;
-		//  int jy=14,jm=6,jd=19,jh=4,jn=29,js=32;
-		long int i,j;
-		long int i0,i1,i2,j0,j1,j2;
-		long int a=16807,m=2147483647,q=127773,r=2836;
-		long int M;  
-		double x,y,u,v;
-		i=iy+70*(im+12*(id+31*(ih+23*(in+59*is))));				//initial the seeds
-		j=jy+70*(jm+12*(jd+31*(jh+23*(jn+59*js))));	
-		i0=a*(i%q)-r*((i-i%q)/q);
-		j0=a*(j%q)-r*((j-j%q)/q);
-		if(i0>=0)	i1=i0;  	else i1=i0+m;	
-		i0=a*(j%q)-r*((j-j%q)/q);
-		if(j0>=0)	j1=j0;  	else j1=j0+m;	
-		for(M=0;M<3*numIon;)
-		{	
-			i0=a*(i1%q)-r*((i1-i1%q)/q);                        
-			if(i0>=0)	i2=i0;  
-			else i2=i0+m;									
-			u=((double)i1)/((double)m);		       //first ramdon number        
-			i1=i2;
-			j0=a*(j1%q)-r*((j1-j1%q)/q);
-			if(j0>=0)	j2=j0;  
-			else j2=j0+m;										
-			v=((double)j1)/((double)m);		        //second ramdon number
-			j1=j2;
-			x=tan(PI*(u-0.5));
-			y=1.6*v/(PI*(x*x+1.0));
-			if(y<=(1.0/sqrt(2.0*PI)*exp(-0.5*x*x)))
-			{
-				if(M<numIon) vel[M].x=x;
-				else if(M<2*numIon) vel[M-numIon].y=x;
-				else vel[M-2*numIon].z=x;
-				// cout<<"this vel is:"<<x<<endl;
-				M++;
-			}
-		}
-
-		int ion;
-		for(ion=0;ion<numIon;ion++)
-		{
-			vel[ion]/=sqrt(allmass[ion]*9.01938e-31);
-		}
-
-		for(ion=0;ion<numIon;ion++)
-		{
-			vel[ion]*=sqrt(temperature/K_BOLTZMAN_AU * K_BOLTZMAN_SI);
-			vel[ion]*= fundamentalTime/BOHR_RADIUS_SI;
-			//rescale the velocities to a.u.
-		}
-	} //xiaohui add 2 lines 2015-09-25, bcast vel
-
-	//2015-09-25, xiaohui
-#ifdef __MPI
-	MPI_Bcast(vel,numIon*3,MPI_DOUBLE,0,MPI_COMM_WORLD);
-#endif
-
-    return;
-}
-
-
-/*void md::InitVelocity(){
-//----------------------------------------------------------------------------
-// DESCRIPTION:
-//   This function initializes the velocities of atom before 
-//   doing MD simulation.
-//----------------------------------------------------------------------------
-  
-  //! >> EXTERNAL VARIABLES <<!
-  //temperature ! temperature to be simulated
-  //double v[3][numIon] // velocity of each atom, in a.u.
-
-  //! >> LOCAL VARIABLES << !!
-  int i, ii, it,ion,nfrozen;
-  double  vel1, vel2, mass,  KE  ;              //KE: to store kinetic energy
-
-  TITLE("md","InitVelocity");
-
-  // note from mohan 2013-03
-  // get the random wave functions from processor 1!!
-  // otherwise the velocity on each node will be different!
-  // then the geometry will be different at each MD step
-  // the results are totally wrong.
-      srand(   (unsigned)time(   NULL   )   );
-    // Assign random numbers to velocities
-	  for(i = 1; i<numIon*3/2 + 2;i++){
-      
-      vel1 = sqrt(-2*log((rand()%10000)/10000.0))*cos(2*PI*(rand()%10000)/10000.0);
-      vel2 = sqrt(-2*log((rand()%10000)/10000.0))*sin(2*PI*(rand()%10000)/10000.0);
-      if ((2*i-2)/3 < numIon){
-		  if(2*i-2%3==0)vel[(2*i-2)/3].x = vel1;
-		  else if((2*i-2)%3==1)vel[(2*i-2)/3].y = vel1;
-		  else if((2*i-2)%3==2)vel[(2*i-2)/3].z = vel1;
-	  }
-	  if ((2*i-1)/3 < numIon){ 
-		  if(2*i-1%3==0)vel[(2*i-1)/3].x = vel2;
-		  else if((2*i-1)%3==1)vel[(2*i-1)/3].y = vel2;
-		  else if((2*i-1)%3==2)vel[(2*i-1)/3].z = vel2;
-	  }
-
-    // Scale velocity according to mass of atoms
-	  for(ion = 0;ion<numIon;ion++){
-              mass = allmass[ion];
-	          vel[ion] = vel[ion]/mass;			  
-	  }
-	  
-  
-    RemoveMovementOfCenterOfMass();
-  
-    // Do not move atoms which are frozen
-    nfrozen = 0;
-	for(ion=0;ion<numIon;ion++){
-         if(ionmbl[ion].x==0)vel[ion].x=0;
-         if(ionmbl[ion].y==0)vel[ion].y=0;
-         if(ionmbl[ion].z==0)vel[ion].z=0;
-		 nfrozen++;
-	}
-
-  // Rescale all velocities to start to exact temperature desired
-    KE=GetAtomKE();
-    cout<<"Kinetic energy from random distribution of velocities = "<< KE;
-    for(ion=0;ion<numIon;ion++){
-       vel[ion] = vel[ion] * sqrt(temperature*1.5 * double( numIon-nfrozen) / KE);
-	}
-  //mohan add 2013-04-17
-    KE= GetAtomKE();
-    out<<"In order to start from the input temperature, kinetic energy should be = "<<KE<<endl;
-
-  // mohan fix bug 2013-03
-  // distribute the random velocity from first processor to other processors
-
-//  WRITE(outputUnit,*) " init velocity:"
-// DO i=1, numIon
-//   WRITE(outputUnit,*) v(1,i), v(2,i), v(3,i)
-// ENDDO
-    
-  return ;
-  }
-}
-*/
-
-
-void md::RemoveMovementOfCenterOfMass(void)
-{
-//---------------------------------------------------------------------------
-// DESCRIPTION:
-//    Compute the velocity for the center of mass, and remove it
-//----------------------------------------------------------------------------
-  //double vel[3][numIon] // velocties of each atom.
-
-  //!! >> LOCAL VARIABLES << !!
-  int ii,it, id,j,ion;
-  Vector3<double> avgvel;
-  double mass,totMass;
-
-  //!! >> FUNCTION << !!
-  
-  avgvel.x = 0;
-  avgvel.y=0;
-  avgvel.z=0;
-  totMass = 0;
-
-  // sum up the total momentum of all the atoms
-  // in the system. (Px,Py,Pz) =\sum m_{i}*(vx,vy,vz)_{i}
-  for(ion=0;ion<numIon;ion++){
-        if(ionmbl[ion].x==0||ionmbl[ion].y==0||ionmbl[ion].z==0)
-        {
-        continue;
-        }
-
-        mass = allmass[ion];
-        avgvel += mass * vel[ion];
-        totMass = totMass + mass;
-  }
-	//cout<<"avgvel: "<<avgvel<<endl;
-    //    cout<<"numIon:"<<numIon<<endl;
-    //    cout<<"last mass:"<<mass<<endl;
-    //    cout<<"last vel:"<<vel[ion].x<<" "<<vel[ion].y<<endl;
-//	cout<<"totMass: "<<totMass<<endl;
-    avgvel = avgvel / totMass;
-//  cout<<setw(15)<<" Average velocity : "<<avgvel.x<<" "<< avgvel.y<<" "<< avgvel.z<<endl;
-  for(ion=0;ion<numIon;ion++){
-    if(ionmbl[ion].x==0||ionmbl[ion].y==0||ionmbl[ion].z==0) 
-    {
-        continue;
-    }
-    vel[ion] = vel[ion] - avgvel;
-  }
-  return ;
-}
-
-
-  void md::OutputMDHeader(){
-//---------------------------------------------------
-// Output the Header for MD simulation
-//---------------------------------------------------
-
-  
-  
-  //!!>>> FUNCTION BODY <<<!!
-
-  // Set up extra output to ion optimizer / MD header
-  if (!MY_RANK)
-  out<<"Q ="<<Qmass<<"a.u.\n"<<" dt = "<<dt*fundamentalTime*1e15<<"fs\n"<< "T ="<<temperature/K_BOLTZMAN_AU<<endl;
-  ofs_running<<"Q ="<<Qmass<<"a.u.\n"<<" dt = "<<dt*fundamentalTime*1e15<<"fs\n"<< "T ="<<temperature/K_BOLTZMAN_AU<<endl;
-  return;
- } 
-
-
-
-  void md::MonitorMeanSquareDisplacement(int step){
-//--------------------------------------------------------------------------
-// DESCRIPTION:
-//  Calculate the Mean Square Displcacement of particles
-//  Compute the mean square displacement of atoms according to formula:
-//    MSD = <|r(t)-r(0)|^2>
-//  <...> denotes the averaging over all the atoms
-//
-//  For solids, MSD saturates to a finte value, for liquids, 
-//  MSD grows linearly with time.
-// 
-//  Also compute the diffusion coefficient D:  
-//   D = MSD/(6*t)
-//----------------------------------------------------------------------------
-
-  //! >> LOCAL VARIABLES << !
-        Matrix3 box;
-        int it,ii,jj,ion,i; 
-        double timePassed;
-        double *tmp;
-        tmp =new double[ntype];
-        Vector3<double> diff; 
-  // INITIALIZATION
-   //       box.e11 =ionlatvec.e11;
-  //	  box.e12 =ionlatvec.e12;
-//	  box.e13=ionlatvec.e13;
-//	  box.e21 =ionlatvec.e21;
-//	  box.e22 =ionlatvec.e22;
-//	  box.e23 =ionlatvec.e23;
-//	  box.e31 =ionlatvec.e31;
-//	  box.e32 =ionlatvec.e32;
-//	  box.e33 =ionlatvec.e33;
-        box=ionlatvec;
-
-
-
-  //! >> FUNCTION << !
-        if(msdstartTime>step-step_rst)return;
-        if (doMSD ==1) {
-    
-               if (step-step_rst == msdstartTime) {
-                    msdcoords0=new Vector3<double>[numIon];
-                    msdcoordst=new Vector3<double>[numIon];
-                    if(rstMD==1){
-                         double *msd0=new double[numIon*3];
-                         double *msdt=new double[numIon*3];
-                         if (!MY_RANK){
-                              stringstream ssc;
-                              ssc << global_out_dir << "Restart_msd.dat";
-                              ifstream file(ssc.str().c_str());
-                              if(!file){
-                                   cout<<"please ensure whether 'Restart_msd.dat' exists"<<endl;
-                                   rstMD=0;
-                              }
-                              if(rstMD==1){
-                                   for(i=0;i<numIon*3;i++){
-                                        file>>msd0[i];
-                                   }
-                                   for(i=0;i<numIon*3;i++){
-                                        file>>msdt[i];
-                                   }
-                              }
-                              file.close();
-                         }
-//2015-09-05, xiaohui
-#ifdef __MPI
-                         MPI_Bcast(msd0,numIon*3,MPI_DOUBLE,0,MPI_COMM_WORLD);
-                         MPI_Bcast(msdt,numIon*3,MPI_DOUBLE,0,MPI_COMM_WORLD);
-#endif
-                         for(i=0;i<numIon;i++){
-                              msdcoords0[i].x=msd0[i*3];
-                              msdcoords0[i].y=msd0[i*3+1];
-                              msdcoords0[i].z=msd0[i*3+2];
-                              msdcoordst[i].x=msdt[i*3];
-                              msdcoordst[i].y=msdt[i*3+1];
-                              msdcoordst[i].z=msdt[i*3+2];
-                         } 
-                   }
-                   if(rstMD==0){for(ion=0;ion<numIon;ion++){
-                         msdcoords0[ion] = cartNoWrap[ion];
-                         msdcoordst[ion] = cartNoWrap[ion];
-		   }
-                   }
-	      }
-
-	      if (step-step_rst > msdstartTime){
-		   for(ion=0;ion<numIon;ion++){
-                         msdcoordst[ion] += vel[ion]*dt;
-		   }
-	      }
-      //-----------------------------------
-      // compute msd and diffuCoeff 
-              for(ii=0;ii<ntype;ii++){ 
-                   msd[ii] = 0;
-                   tmp[ii]=0;
-              }
-
-              ion=0;
-              for( ii=0;ii<ntype;ii++){
-                   for(jj=0;jj<na[ii];jj++){// loop over ions
-                       diff = msdcoords0[ion] - msdcoordst[ion];
-                       // the mean square displacement has unit of Bohr^2,
-                       // so we should not use 'SQRT' here.
-                       // mohan note 2013-05-01
-                       //tmp = SQRT( SUM(diff**2) )
-                       tmp[ii] += diff*diff;
-                       ion++;
-                   }
-                   msd[ii] = msd[ii] + tmp[ii];
-              }
-              if (!MY_RANK){
-                   for(ii=0;ii<ntype;ii++){ 
-                        out<<"msd_single_atom:"<< ii<<" "<< tmp[ii]<<endl;
-			//xiaohui remove, 2015-09-16
-                        //ofs_running<<"msd_single_atom:"<< ii<<" "<< tmp[ii]<<endl;
-
-                   }
-                   if(rstMD==0)
-                        timePassed = dt * ( step - step_rst - msdstartTime );
-                   else timePassed = dt * ( step - msdstartTime );
-                   for(ii=0;ii<ntype;ii++){
-                        msd[ii] = msd[ii] / double(na[ii]);
-                        diffuCoeff = msd[ii] / (6.0*timePassed);
-
-			//xiaohui move 2 lines, add 2 another lines 2015-09-25, setprecision
-                        //out<<ii<< "msd:STEP="<<step<< " msd= "<<msd[ii]<<" (Bohr^2), D= "<<diffuCoeff<< " a.u."<<endl;
-                        //ofs_running<<" "<<ii<< "msd:STEP="<<step<< " msd= "<<msd[ii]<<" (Bohr^2), D= "<<diffuCoeff<< " a.u."<<endl;
-                        out<<ii<< "msd:STEP="<<step<< " msd= "<<setprecision (12)<<msd[ii]<<" (Bohr^2), D= "<<setprecision (12)<<diffuCoeff<< " a.u."<<endl;
-                        ofs_running<<ii<< "msd:STEP="<<step<< " msd= "<<setprecision (12)<<msd[ii]<<" (Bohr^2), D= "<<setprecision (12)<<diffuCoeff<< " a.u."<<endl;
-                   }
-              }
-
-// NOTE !
-// to convert the diffuCoeff to A^2/PIcosecond, 
-// D * bohr * bohr / ( fundamentalTime / 1.0e-12 )
-// where fundamentalTime is 2.418884326505E-17_DP
-//!DEBUG          
-//          DO ii = 1,numIon
-//            WRITE(message,'(a,I6,a,ES12.4)') & 
-//              ' ion:',ii,' ion_sd:', & 
-//              SQRT(SUM((msd_coordst(:,ii)-msd_coords0(:,ii))**2))
-//            CALL WRTOUT(6,Message)
-//            WRITE(message,'(a,I6,a,3ES12.4)') & 
-//              ' ion:',ii,' ion_fracCoord:', cell%ionTable(ii)%coord(:)
-//              CALL WRTOUT(6,Message)
-//          ENDDO ! loop over ions
-//!~DEBUG
-
-	 // if (step > msd_startTime)
-               if(step==NSTEP){
-                     if(msdcoords0!=NULL)delete []msdcoords0;
-                     if(msdcoordst!=NULL)delete []msdcoordst;
-               }
-	  } // doMSD
-     return;
-  }
-
-
-  void md::OutputIonVelocityAndStats(int iter){
-//---------------------------------------------------------------------------
-// DESCRIPTION:
-//  Output the 
-//   (1) ion's velocity 
-//   (2) thermostat's xLogS, and vLogs (for NVT, NPT)
-//  (3) barostat's vbox (box's velocity)  (NPT only)
-// for restart purpose
-//---------------------------------------------------------------------------
-// REVISION LOG:
-//    Dec/29/2008: Created by Chen Huang
-//----------------------------------------------------------------------------
-  
-  // iter, mdtype //number of ions, iteration number, NVT/NPT
-
-//2015/4/30:change to output msd data for restart.
-  int i,j;
-
-  bool pass;
-  pass = 0;
-
-
-  if (dumpmdfreq==1||iter==1||( dumpmdfreq > 1&& iter%dumpmdfreq==0 ) )
-    pass =1;
-  if (!pass) return;
-  
-  if (msdstartTime>iter)return;
-  if (!MY_RANK){
-        stringstream ssc;
-        ssc << global_out_dir << "Restart_msd.dat";
-        ofstream vfile(ssc.str().c_str());
-
-/*  vfile<<"MD STEP= "<<iter<<endl;
-  //----- thermo-stat's velocity and positions
-  vfile<<" xLogS: "<< xLogS<< " <= thermostat position :(a.u.)"<<endl;
-  vfile<<" vLogS: "<< vLogS<< " <= thermostat velocity :(a.u.)"<<endl;
-	  switch(mdtype){
-  case 1: // NVT
-    vfile<<" 0.0 0.0 0.0 <= barostat velocity 1"<<endl;
-    vfile<< " 0.0 0.0 0.0 <= barostat velocity 2"<<endl;
-    vfile<< " 0.0 0.0 0.0 <= barostat velocity 3"<<endl;    // add addiontional outputs if needed
-    break;  
-  case 2: // NPT
-    vfile<<vBoxG.e11<<" "<<vBoxG.e12<<" "<<vBoxG.e13<< " <= barostat velocity : vBoxG[0][]"<<endl;
-    vfile<<vBoxG.e21<<" "<<vBoxG.e22<<" "<<vBoxG.e23<< " <= barostat velocity : vBoxG[1][]"<<endl; 
-    vfile<<vBoxG.e31<<" "<<vBoxG.e32<<" "<<vBoxG.e33<< " <= barostat velocity : vBoxG[2][]"<<endl; 
-    break; 
-  case 3: // NVE
-    vfile<< " 0.0 0.0 0.0 <= barostat velocity 1"<<endl;
-    vfile<< " 0.0 0.0 0.0 <= barostat velocity 2"<<endl;
-    vfile<< " 0.0 0.0 0.0 <= barostat velocity 3"<<endl;
-	break;
-  default: break;
-  }
-  // output velocities
-  vfile<<"ION VELOCITIES (a.u.):"<<endl;
-	  for( i = 0;i<numIon;i++){
-          vfile<<"  "<< vel[i].x<<"  "<<vel[i].y<<"  "<<vel[i].z;
-	  }
-*/
-   //vfile<<"msd0: " <<endl;
-        for(i=0;i<numIon;i++){
-                vfile<<setprecision (12)<<msdcoords0[i].x<<" "<<setprecision (12)<<msdcoords0[i].y<<" "<<setprecision (12)<<msdcoords0[i].z<<endl;
-        }
-     //   vfile<<"msdt: " <<endl;
-        for(i=0;i<numIon;i++){
-                vfile<<setprecision (12)<<msdcoordst[i].x<<" "<<setprecision (12)<<msdcoordst[i].y<<" "<<setprecision (12)<<msdcoordst[i].z<<endl;
-        }
-
-  vfile.close();
-  }
-   return;
-}
-
-
-
-void md::OutputMDGeom(int iter){
-//---------------------------------------------------------------------------
-// DESCRIPTION:
-//  Output the ion's position and cell's structure for monitoring and MD
-//  restart.
-//---------------------------------------------------------------------------
-// REVISION LOG:
-//    Dec/29/2008: Created by Chen Huang
-//----------------------------------------------------------------------------
-  
-  bool pass;
-  int i,it;
-  ofstream uf;
-  string iteration_number;
-  string eleName;
-   
-  //! >> FUNCTION << !
-  pass = 0;
- // if(rankGlobal!=0) return;
-  if(dumpmdfreq==1||iter==1||( dumpmdfreq > 1&&(iter%dumpmdfreq)==0 ) ) 
-    pass = 1;
-  if(!pass)return;
- // cout<<"test the function OutputMDGeom"<<endl;
-
-  //the final file name is the md_geometry plus the iteraction number
-  iteration_number=intTurnTostring(iter,iteration_number);
-  string  message=mdoutputpath+"_ion"+iteration_number+".txt";
-//  cout<<"file name is:"<<message<<endl;
-  // open the file: 'ion.1.dat'
-  uf.open(message.c_str());
-//  cout<<"already open file"<<endl;
-  uf<< "BLOCK LATTICE_CART"<<endl;
-  uf<<"    "<<ionlatvec.e11*ucell.lat0<<"  "<<ionlatvec.e12*ucell.lat0<<"  "<<ionlatvec.e13*ucell.lat0<<endl;
-  
-  uf<<"    "<<ionlatvec.e21*ucell.lat0<<"  "<<ionlatvec.e22*ucell.lat0<<"  "<<ionlatvec.e23*ucell.lat0<<endl;
-
-  uf<<"    "<<ionlatvec.e31*ucell.lat0<<"  "<<ionlatvec.e32*ucell.lat0<<"  "<<ionlatvec.e33*ucell.lat0<<endl;
-
-  uf<< "END BLOCK LATTICE_CART"<<endl;
-//  cout<<"already out put lattice cart"<<endl;
-  
-  //output the atom positions.
-  uf<<"BLOCK POSITIONS_CART"<<endl;
-  int ion=0;
-  for(it=0;it<ntype;it++){
-      for(i = 0;i< na[it];i++){
-          eleName=ucell.atoms[it].label; 
-          uf<<"  "<<eleName 
-            <<"  "<<cartNoWrap[ion].x  // unit of cartNorap is 'Bohr' 
-            <<"  "<<cartNoWrap[ion].y   
-            <<"  "<<cartNoWrap[ion++].z<<endl;
-	  }
-  }
-  uf<< "END BLOCK POSITIONS_CART"<<endl;
-  uf<<"already output positions cart"<<endl;
-  double distance;
-  distance=pow(cartNoWrap[1].x-cartNoWrap[0].x,2.0)+pow(cartNoWrap[1].y-cartNoWrap[0].y,2.0)+pow(cartNoWrap[1].z-cartNoWrap[0].z,2);
-  distance=sqrt(distance);
-//  cout<<"distance is "<<distance<<" ";
-  uf.close();
-    return;
-}
-
-
-
-
-void md::ReadNewTemp(int step ){
-//-----------------------------------------------------------------
-//  If fixTemperature == 0, then this subroutine will be skipped
-//  otherwise, we are going to read a new temperature from the disk.
-//  You only need to create a file named temperature.dat, and then 
-//  change the temperature in it. 
-//-----------------------------------------------------------------
- 
-
-    double intemp;
-    int sstep=0;
-
-    if ( fixTemperature > 0 ){
-
-    // Change the temperature every 'fixTemperature' steps.
-	  if( (step!=1)&&(step%fixTemperature == 1) ){
-    
-      // Read in new temperature from file.
-		  ifstream file;
-		  file.open("ChangeTemp.dat");
-		  if (!file){
-                         cout<<"ERROR IN OPENING ChangeTemp.dat, CODE STOP!"<<endl;
-                         exit(0);
-		  }
-                  while((sstep+fixTemperature)<step){
-                  file>> intemp;
-                  sstep+=fixTemperature;
-                  }
-                  file.close();
-
-      // Renew information.
-                  intemp =  intemp * K_BOLTZMAN_AU;
-                  if ( fabs(intemp-temperature) >1e-6 ) {
-                       cout <<"(ReadNewTemp): Read in new temp:"<< intemp/K_BOLTZMAN_AU 
-                            <<" previous temp:"<< temperature/K_BOLTZMAN_AU<<endl;
-                       temperature = intemp;
-                  }
-    	          else{
-                       cout<<"(ReadNewTemp): new temp:"<< intemp/K_BOLTZMAN_AU
-     		      	   <<" previous temp:"<<temperature/K_BOLTZMAN_AU
-    			   << ". No change of temp."<<endl;
-                  }
-	  }
-    }
-
-    return;
-}
-string md::intTurnTostring(int iter,string path){
-//int number turn to string term
-	int i[6],k;
-	if(iter>99999){
-        k=6;
-        i[5]=iter%10;
-		iter=iter/10;
-		i[4]=iter%10;
-		iter=iter/10;
-		i[3]=iter%10;
-		iter=iter/10;
-		i[2]=iter%10;
-		iter=iter/10;
-		i[1]=iter%10;
-		iter=iter/10;
-		i[0]=iter%10;
-	}
-	else if(iter>9999){
-		k=5;
-		i[4]=iter%10;
-		iter=iter/10;
-		i[3]=iter%10;
-		iter=iter/10;
-		i[2]=iter%10;
-		iter=iter/10;
-		i[1]=iter%10;
-		iter=iter/10;
-		i[0]=iter%10;
-	}
-	else if(iter>999){
-		k=4;
-		i[3]=iter%10;
-		iter=iter/10;
-		i[2]=iter%10;
-		iter=iter/10;
-		i[1]=iter%10;
-		iter=iter/10;
-		i[0]=iter%10;
-	}
-	else if(iter>99){
-		k=3;
-		i[2]=iter%10;
-		iter=iter/10;
-		i[1]=iter%10;
-		iter=iter/10;
-		i[0]=iter%10;
-	}
-	else if(iter>9){
-		k=2;
-		i[1]=iter%10;
-		iter=iter/10;
-		i[0]=iter%10;
-	}
-	else if(iter>0){
-		k=1;
-		i[0]=iter;
-	}
-	for(int j=0;j<k;j++){
-		path+=(i[j]+48);
-	}
-	return path;
-}
-void md::connection0(){
-//some prepared imformation
-    int i,it,ion;
-	ion=0;
-        nfrozen=0;
-	for(it=0;it<ntype;it++){
-		for(i=0;i<na[it];i++){
-                        allmass[ion]=ucell.atoms[it].mass/6.02/9.109*1e5;
-	         	ionmbl[ion]=ucell.atoms[it].mbl[i];
-                        if (ionmbl[ion].x==0||ionmbl[ion].y==0||ionmbl[ion].z==0) nfrozen++;
-
-			ion++;
-		}
-	}
-	return;
-}
-	
-void md::connection1(){
-//call each atom's position
-	int i,it,ion;
-	ion=0;
-	for(it=0;it<ntype;it++){
-		for(i=0;i<na[it];i++){
-                        cartNoWrap[ion]=ucell.atoms[it].tau[i]*ucell.lat0;
-	        	taudirac[ion]=ucell.atoms[it].taud[i];
-			ion++;
-		}
-	}
-	return;
-}
-void md::connection2(){
-//to refresh the atoms' positions
-	int i,it,ion;
-	ion=0;
-	for(it=0;it<ntype;it++){
-		for(i=0;i<na[it];i++){
-                        ucell.atoms[it].tau[i]=cartNoWrap[ion]/ucell.lat0;
-	            	ucell.atoms[it].taud[i]=taudirac[ion];
-			ion++;
-		}
-	}
-	return;
-}
-void md::callforce(){
-//to call the force of each atom
-	int ion;
-	Force_Stress_LCAO FSL;
-	FSL.allocate (); 
-	matrix forceMd;
-	FSL.getForceStress(FORCE, STRESS, TEST_FORCE, TEST_STRESS, forceMd, stress_lcao);
-	for(ion=0;ion<numIon;ion++){
-			force[ion].x=forceMd(ion,0)/2.0;
-			force[ion].y=forceMd(ion,1)/2.0;
-			force[ion].z=forceMd(ion,2)/2.0;
-                       // cout<<force[ion].x/0.0195<<" "<<force[ion].y/0.0195<<" "<<force[ion].z/0.0195<<endl;
-	}
-#ifdef __MPI //2015-10-01, xiaohui
-        atom_arrange::delete_vector( SEARCH_RADIUS );
-#endif //2015-10-01, xiaohui
-	return;
-}
-void md::moveatoms(int step){
-// intent to move atoms to the center of the cell
-//intent to ensure that no atom is out of the cell
-      int ia;
-      int ii,it, id,j,ion;
-      Vector3<double> avgdc;
-      double mass,totMass;
-
-  //!! >> FUNCTION << !!
- 
-// 2015/9/18 no move of center of mass
-/* 
-        avgdc.x = 0;
-        avgdc.y=0;
-        avgdc.z=0;
-        totMass = 0;
-
-  // sum up the total momentum of all the atoms
-  // in the system. (Px,Py,Pz) =\sum m_{i}*(vx,vy,vz)_{i}
-        for(ion=0;ion<numIon;ion++){
-                         mass = allmass[ion];
-                         avgdc += mass * taudirac[ion];
-                         totMass = totMass + mass;
-        }
-	//cout<<"avgdc: "<<avgdc<<endl;
-        avgdc = avgdc / totMass;
-        avgdc.x += -0.5;
-        avgdc.y += -0.5;
-        avgdc.z += -0.5;
-       // cout<<" Average dirac coords : "<<avgdc.x<<" "<< avgdc.y<<" "<< avgdc.z<<endl;
-        for(ion=0;ion<numIon;ion++){ 
-                        taudirac[ion] = taudirac[ion] - avgdc;
-        }
-      }*/
-      for(ia=0;ia<numIon;ia++){
-//changed by a constrain box
-                 
-       if(taudirac[ia].x<0) taudirac[ia].x += 1.0;
-                        if(taudirac[ia].y<0) taudirac[ia].y += 1.0;
-                        if(taudirac[ia].z<0) taudirac[ia].z += 1.0;
-                        if(taudirac[ia].x>=1.0) taudirac[ia].x -= 1.0;
-                        if(taudirac[ia].y>=1.0) taudirac[ia].y -= 1.0;
-                        if(taudirac[ia].z>=1.0) taudirac[ia].z -= 1.0;
-/*
-                        if(taudirac[ia].x<0){
-                                    taudirac[ia].x *= -1.0;
-                                    vel[ia].x *=-1.0;
-                        }
-                        if(taudirac[ia].y<0) {
-                                    taudirac[ia].y *= -1.0;
-                                    vel[ia].y *=-1.0;
-                        }
-                        if(taudirac[ia].z<0) {
-                                    taudirac[ia].z *= -1.0;
-                                    vel[ia].z *=-1.0;
-                        }
-                        if(taudirac[ia].x>=1.0){
-                                    taudirac[ia].x = 2.0-taudirac[ia].x;
-                                    vel[ia].x *=-1.0;
-                        }
-                        if(taudirac[ia].y>=1.0){
-                                    taudirac[ia].y = 2.0-taudirac[ia].y;
-                                    vel[ia].y *=-1.0;
-                        }
-                        if(taudirac[ia].z>=1.0){
-                                    taudirac[ia].z = 2.0-taudirac[ia].z;
-                                    vel[ia].z *=-1.0;
-                        }
-*/
-                        if(taudirac[ia].x<0 || taudirac[ia].y<0
-                                || taudirac[ia].z<0 ||
-                                taudirac[ia].x>=1.0 ||
-                                taudirac[ia].y>=1.0 ||
-                                taudirac[ia].z>=1.0)
-                        {
-                                cout << " ion: "<< ia  << endl;
-                                cout << "d=" << taudirac[ia].x << " " <<
-                                taudirac[ia].y << " " << taudirac[ia].z << endl;
-                                WARNING_QUIT("md::moveatoms","the movement of atom is larger than the length of cell.");
-                        }
-
-                        cartNoWrap[ia] = taudirac[ia] * ionlatvec*ucell.lat0;
-         }
-         return;
-}
-void md::printpos(string file,int iter){
-//intend to output the positions of atoms to ordered file
-          bool pass;
-          pass = 0;
-
-
-          if (dumpmdfreq==1||iter==1||( dumpmdfreq > 1&& iter%dumpmdfreq==0 ) )
-             pass =1;
-          if (!pass) return;
-
-          string file1=file+".xyz";
-          string file2=file+".cif";
-
-          //xiaohui add 'OUT_LEVEL', 2015-09-16
-          if(OUT_LEVEL == "i"||OUT_LEVEL == "ie") ucell.print_tau();
-          if(OUT_LEVEL == "i"||OUT_LEVEL == "ie") ucell.print_cell_xyz(file1);
-          ucell.print_cell_cif(file2);
-          stringstream ss;
-
-          ss << global_out_dir << "STRU_MD";
-
-          //zhengdy modify 2015-05-06, outputfile "STRU_Restart"
-          ucell.print_stru_file(ss.str(),2);
- 
-          return;
-}
-void md::md_release(){
-          
-          if(force!=NULL)delete []force;
-          if(ionmbl!=NULL)delete []ionmbl;
-          if(cartNoWrap!=NULL)delete []cartNoWrap;
-          if(vel!=NULL)delete []vel;
-          if(taudirac!=NULL)delete []taudirac;
-          if(allmass!=NULL)delete []allmass;
-          
-          return;
-}
-void md::scalevel(){
-   int i;
-   double ke=GetAtomKE();
-   if(ke>1e-9)
-      for(i=0;i<numIon;i++){
-          vel[i]*=sqrt(3*(numIon-nfrozen)*temperature/ke/2);
-      }
-   return;
-}
-/*void md::PDF(int step){ 
-   int i,j,k;
-   double dis,c;
-   c=(ionlatvec.e11+ionlatvec.e12+ionlatvec.e13)*ucell.lat0;
-   if(step==1)for(i=0;i<100;i++){
-            rdf[i]=0;
-       }
-   for(i=0;i<numIon;i++){
-       for(j=i;j<numIon;j++){
-          for(k=0;k<100;k++){
-              dis=sqrt((cartNoWrap[i]-cartNoWrap[j])*(cartNoWrap[i]-cartNoWrap[j]));
-              if(dis<((double)(k+1)*c/(double)100)&&dis>((double)k*c/(double)100))rdf[k]=rdf[k]+1;
-          }
-       }
-   }
-   return;
-}
-void md::printRDF(int step){
-   int i;
-   ofstream file;
-//   if(step==1){
-//        file.open("RDFnow.dat");
- //       file.close();
-  // }
-   string it="RDF";
-   if(fixTemperature!=0)
-        if(step/fixTemperature!=0) intTurnTostring(step/fixTemperature,it);
-   it+=".dat";
-   if(step%dumpmdfreq==1){
-        file.open(it.c_str());
-        file<<step<<endl;
-        for(i=0;i<60;i++){
-           file<<"number "<<i<<" RDF: "<<(double)rdf[i]<<endl;
-        }
-        file.close();
-   }
-   return;
-}
-void md::RDF(){
-     int i,j,k;
-     double dis;
-     for(i=0;i<60;i++){
-            rdf[i]=0;
-     }
-     for(i=0;i<numIon;i++){
-        if(cartNoWrap[i].x-11.342>0&&cartNoWrap[i].x+11.342<0)
-            for(j=0;j<numIon;j++){
-                for(k=0;k<60;k++){
-                     dis=sqrt((cartNoWrap[i]-cartNoWrap[j])*(cartNoWrap[i]-cartNoWrap[j]));
-                     if(dis<((double)k/5.29+0.2835)&&dis>((double)k/5.29+0.0945))rdf[k]=rdf[k]+1;
-                }
-            }
-     }
-     return;
-}*/
diff --git a/ABACUS.develop/source/src_pw/md.h b/ABACUS.develop/source/src_pw/md.h
deleted file mode 100644
index 3aee8f41d5..0000000000
--- a/ABACUS.develop/source/src_pw/md.h
+++ /dev/null
@@ -1,111 +0,0 @@
-#ifndef MD_H
-#define MD_H
-
-#include<iostream>
-#include<string>
-#include<fstream>
-#include<stdlib.h>
-#include<time.h>
-#include<cmath>
-#include"tools.h"
-#include"global.h"
-#include"../src_lcao/FORCE_STRESS.h"
-#include"unitcell_pseudo.h"
-
-using namespace std;
-
-//const double fundamentalTime = 2.418884326505e-17; 
-
-//-------------------------------------------------------------------
-// mohan reconstruction note: 2021-02-07
-// 1) explain every variable and every function
-// 2) if too many, divide the class into more files 
-// 3) leave a space between different functions
-// 4) write the comments above the variables or functions
-//-------------------------------------------------------------------
-
-class md
-{
-	public:
-
-	md(int n=1);
-
-	void md_allocate();
-	void initMD();
-	bool RestartMD();
-	void mstout(int step);
-	void RemoveMovementOfCenterOfMass();
-	double GetAtomKE();
-	void MakeIntCoeff();
-	void InitVelocity();
-	void MonitorMeanSquareDisplacement(int step);
-	void OutputIonVelocityAndStats(int iter);
-	void OutputMDHeader();
-	void OutputMDGeom(int iter);
-	void ReadNewTemp(int step);
-	string intTurnTostring(int iter,string path);
-	void connection0();
-	void connection1();
-	void connection2();
-	void callforce();
-	void moveatoms(int step);
-    void printpos(string file,int iter);
-    void md_release();
-    void scalevel();
-    void RDF();
-    void printRDF(int step);
-    void PDF(int step);
-
-	// variables:
-
-	int mdtype; // 1: NVT:Nose Hoover, 2:NVT:Nose Hoover add velocity scaling 3: NPT, 0: NVE 
-	double xLogS; // thermostat's position
-	double vLogS;// thermostat's speed
-	Vector3<double> *vel;//   velocity of each atom, unit is a.u.
-	Matrix3 vBoxG;        // barostat's velocity
-	Matrix3 ionlatvec;
-	double  tauthermo;  // The thermostat fluctuation period, for Al fcc at 300K ~ 1000fs
-	double taubaro; // The barostat   fluctuation period, for Al fcc at 300K ~ 500 fs
-	double Qmass;  // Inertia of extended system variable
-	double dt ; // Time increment (hbar/E_hartree)
-	double temperature;  //Temperature (in Hartree, 1 Hartree ~ 3E5 K)
-	bool doMSD;   // compute the mean square displacement.
-	bool doMSDAtom;        // compute the mean square displacement for each atom.
-	double msdstartTime;   //beyond this time, msd is going to be computed
-	double *msd    ;
-	double diffuCoeff ;
-	int numIon;
-	int ntype;
-	int* na;
-	int nfrozen;
-	double *allmass;
-	Vector3<double> *force;
-	Vector3<int> *ionmbl;
-	Vector3<double> *taudirac;           //dirac coord of atoms
-	Vector3<double> *cartNoWrap;        //cartensian coord of atoms, *not* wrapped
-	Vector3<double> *msdcoords0;        
-	Vector3<double> *msdcoordst;
-	int nResn;   // Coarseness of Trotter factorization of Nose-Hoover thermostat
-	int nYosh;  // Order of thermostat scheme corresponding with "intCoeff"
-	double *intCoeff;
-	int outputstressperiod;// The period to output stress
-	int dumpmdfreq;       // The period to dump MD information for monitoring and restarting MD
-	int fixTemperature;       // The period to read in new temperature during MD run.
-	string mdoutputpath;         // output directory of md files: .ion .vel
-	int rstMD;// 1 : restart MD, vel.restart and ion.restart files will be read
-	// 0 : not restart from ion and vel files
-
-	double ediff;           //parameter for constrain 
-	double ediffg;          //parameter for constrain
-	int startStep;// the first step number of MD
-	bool velRescale;// whether to rescale velocities at the first step
-	ofstream out;
-	double mdenergy;
-	const static double fundamentalTime;
-	double rdf[100];
-	int step_rst; 
-
-	matrix stress_lcao;
-};
-
-#endif
diff --git a/ABACUS.develop/source/src_pw/mdNVE.cpp b/ABACUS.develop/source/src_pw/mdNVE.cpp
deleted file mode 100644
index ffb3a199f5..0000000000
--- a/ABACUS.develop/source/src_pw/mdNVE.cpp
+++ /dev/null
@@ -1,436 +0,0 @@
-#include"mdNVE.h"
-
-/*
-using namespace std;
-class mdNVE:public md{
-public:
-	mdNVE(int n=1):md(1){};
-        void	runNVE();
-	void RescaleVelocity(double KE,double newKE);
-	double Conserved(double KE, double PE);
-	double MAXVALF();
-private:
-	                      //>> INTERNAL VARIABLES <<//
-  int nfrozen;               // Number of frozen atoms
-   int ii,k,it,ion;
-   int i;                       // Counters to parse through the ion table.
-   int step;                    // step number we are on
-   //int a.nat;                  // total number of atoms
-   int nStep;                   // Total number of steps
-   int errorFlag;
-
-  //vector3<double> velocity;           // Velocity of ions (atomic units)
-  Vector3<double> fracStep;           // Step (due to velocity) in fractional coordinates
-  double mass;   // atom mass, temperay var
-   double maxStep;     // Largest movement by an ion
-   double hamiltonian,conservedE,maxForce;
-   double msd;                  //  mean square displacement
-   double diffuCoeff;           // diffusion coefficient
-   double newKE;
-   double twiceKE;     // Kinetic energy x 2
-   double oldEtot;        // old energy
-
-  //TYPE(stopwatch):: &           // Timer objects
-  //  watch, &
-  //  watch2
-  
-  //string parameterInfo;      // Line containing values for Qmass, dt, temperature for output to header
-
-  double tempNow;
-
-};
-*/
-
-void mdNVE::runNVE(int step1){
-//-------------------------------------------------------------------------------
-//
-// REFERENCES:
-//   Martyna et. al., "Explicit reversible integrators for extended system
-//      dynamics," Mol. Phys. 87 (1996), 1117.
-//   Allen and Tildesley, "Computer Simulation of Liquids" (predictor-corrector)
-//   Frenkel and Smit, "Understanding Molecular Simulation".
-//
-//------------------------------------------------------------------------------
-// REVISION LOG:
-//     04-17-2013: revised from NVT code.
-//
-//------------------------------------------------------------------------------
-
-
-
- // Optimizer // A subroutine that is called to optimize the electron density relative to the ion positions
-
-  
-   // rho              // Electron density, realspace
-   // energy                    // Energy from electronic optimization
-  //Vector3<double>  forces;           // Total forces.  
-                     // First is ion number, second direction (1,2,3 for x,y,z)
-
-  
-
-                        //>> INITIALIZATION <<//
-  TITLE("mdNVE","runNVE");
-   timer::tick("mdNVE","runNVE",'C');
-  step=step1+step_rst;
-  if(step==1)cout<<" (NVE) Start."<<endl;
-  else {
-      cout<<" (NVE) step: "<<step<<endl;
-      }
-  //watch2 = TimerStart()
-  //watch = TimerStart()
-  //numIon = sizeof(cell%ionTable);
-  
-  cout<<" (NVE) Fronzen Atom Number : "<<nfrozen<<endl;
-
-
-  // Calculate the number of total MD steps.
- // nStep = int(timeTot/dt);
-  //f.zero_out();//forces = 0;
-
-  // Allocate arrays for velocities, cartesian coordinates.
-  //velocity=new Vector3<double>[numIon];
-  //cartNoWrap=new Vector3<double>[numIon];
-  
-
-  // Set up extra output to ion optimizer / MD header
-  cout<<"Time interval       : "<< dt*fundamentalTime*1e15<< " (fs)"<<endl;
-  cout<<"Target temperature  : "<< temperature/K_BOLTZMAN_AU<< " (K)"<<endl;
-  //cout<<"Number of NVE steps : "<< nStep<<endl;
-
-  //if (rankGlobal==0)  GeometryMinimizerReportHeader(parameterInfo);
-  if(step==1){
-      for(k=0;k<numIon;k++){
-        if(ionmbl[k].x==0)vel[k].x=0;
-        if(ionmbl[k].y==0)vel[k].y=0;
-        if(ionmbl[k].z==0)vel[k].z=0;
-      }
-      if(nfrozen==0)RemoveMovementOfCenterOfMass();
-      scalevel();
-  }
-  moveatoms(step);
-  twiceKE=GetAtomKE();
-  twiceKE = twiceKE * 2;
-
-  tempNow = twiceKE/(3*double(numIon-nfrozen))/K_BOLTZMAN_AU;
-  cout<<" start temperature = "<< tempNow/K_BOLTZMAN_AU<< " (k)"<<endl;
-
-
-
-
-
-  // Set up forces 
-  callforce();
-	if(STRESS)
-	{
-		cout<<"output Pressure for check!"<<endl;
-		double press;
-		for(int i=0;i<3;i++)
-			press += stress_lcao(i,i)/3;
-		press += twiceKE/3/ucell.omega; //output virtual press = 2/3*Ek/V + sum(sigma[i][i])/3
-		double unit_transform = RYDBERG_SI / pow(BOHR_RADIUS_SI,3) * 1.0e-8 ;
-		cout<<"Virtual Pressure is "<<press*unit_transform<<" Kbar "<<endl;
-	}
-  
-  for(k=0;k<numIon;k++){
-        if(ionmbl[k].x==0)force[k].x=0;
-        if(ionmbl[k].y==0)force[k].y=0;
-        if(ionmbl[k].z==0)force[k].z=0;
-  }
-
-  //cout<<"begin maxForce"<<endl;
-  maxForce = MAXVALF();
-  cout<<"maxForce: "<<sqrt(maxForce)<<endl; 
-  mdenergy=en.etot/2;
-  conservedE = Conserved(twiceKE/2, mdenergy);
- // oldEtot=mdenergy;
-
-  // Print intial output
-  //GeometryMinimizerReportSteps(startStep, TimerStop(watch), forces, maxForce, &
-  //                                  0, .TRUE., hamiltonian, .TRUE., &
-  //                                  temperature=temperature);
-
-
-  cout<< "NVE_STEP "<<" "<<"SystemEnergy"<<" "<< "Conserved"<<" "<< "DeltaE"<<" "<< "Temperature"<<endl;
-
-  // Output the message to the screen.
-  cout<<step<<" "<< mdenergy<<" "<< conservedE<<" "<< mdenergy-oldEtot<<" "<< twiceKE/(3*double(numIon-nfrozen))/K_BOLTZMAN_AU<<endl;
- // oldEtot=mdenergy;
-                        //>> FUNCTION BODY <<//
-
-  //xiaohui move this line 2015-09-15
-  //cout<<"ready part is ready!"<<" "<<endl;
-
-  // Loop for MD step 
- // for(step = startStep+1;step< nStep+1;step++){
-
-    cout<<"------------------------------------------------------------------------------"<<endl;
-    cout<< "MD(NVE) STEP "<< step<<endl;
-    cout<<"------------------------------------------------------------------------------"<<endl;
-    if (!MY_RANK){
-         out<<"------------------------------------------------------------------------------"<<endl;
-         out<< "MD(NVE) STEP "<< step<<endl;
-         out<<"------------------------------------------------------------------------------"<<endl;
-         ofs_running<<"------------------------------------------------------------------------------"<<endl;
-         ofs_running<< "MD(NVE) STEP "<< step<<endl;
-         ofs_running<<"------------------------------------------------------------------------------"<<endl;
-    }
-   // callforce();
-  /*  cout<<"force to test:"<<endl;
-    cout<<force[0].x<<" "<<force[0].y<<" "<<force[0].z<<endl;
-    cout<<force[1].x<<" "<<force[1].y<<" "<<force[1].z<<endl;
-
-    cout<<"coord to test:"<<endl;
-    cout<<cartNoWrap[0].x<<" "<<cartNoWrap[0].y<<" "<<cartNoWrap[0].z<<endl;
-    cout<<cartNoWrap[1].x<<" "<<cartNoWrap[1].y<<" "<<cartNoWrap[1].z<<endl;
-    // Calculate the Mean-Square-Displacement.
-    cout<<"vel to test"<<endl;
-    cout<<vel[0].x<<" "<<vel[0].y<<" "<<vel[0].z<<endl;
-    cout<<vel[1].x<<" "<<vel[1].y<<" "<<vel[1].z<<endl;*/
-    //MonitorMeanSquareDisplacement(step, msd,diffuCoeff);
-  
-    // (1) 1st step of Verlet-Velocity
-    // New velocity are obtained
-	for(k=0;k<numIon;k++){       
-           mass = allmass[k];
-     //2014/6/28
-           vel[k] = vel[k] + force[k]/mass*dt/2.0;
-       	} 
-    // (2) 2nd step of Verlet-Velocity 
-    // Update the Non-Wrapped cartesion coordinate
-       twiceKE = GetAtomKE();
-       twiceKE = 2 * twiceKE;
-       if(step!=1||rstMD==1)for(k=0;k<numIon;k++){       
-           mass = allmass[k];
-     //2014/6/28
-           vel[k] = vel[k] + force[k]/mass*dt/2.0;
-       	} 
-        if(nfrozen==0) RemoveMovementOfCenterOfMass();  
-       for(i=0;i<numIon;i++){
-           cartNoWrap[i] = vel[i]*dt + cartNoWrap[i];
-		
-	}
-    // Calculate the maximal velocities.
-    // The step in fractional coordinates
-       maxStep = 0;
-       for( i = 0;i< numIon;i++){
-		if((pow(vel[i].x,2.0)+pow(vel[i].y,2.0)+pow(vel[i].z,2.0))>maxStep)
-      maxStep = pow(vel[i].x,2.0)+pow(vel[i].y,2.0)+pow(vel[i].z,2.0);
-      Mathzone::Cartesian_to_Direct(vel[i].x*dt/ucell.lat0,vel[i].y*dt/ucell.lat0,vel[i].z*dt/ucell.lat0,
-                                    ionlatvec.e11,ionlatvec.e12,ionlatvec.e13,
-                                    ionlatvec.e21,ionlatvec.e22,ionlatvec.e23,
-                                    ionlatvec.e31,ionlatvec.e32,ionlatvec.e33,
-                                    fracStep.x,fracStep.y,fracStep.z);
-//      fracStep = ionlatvec.Inverse()* vel[i]*dt/ucell.lat0;
-      taudirac[i] = taudirac[i] + fracStep;
-    }
-   // cout<<"dx: "<<fracStep.x<<" "<<fracStep.y<<" "<<fracStep.z<<endl;
-    moveatoms(step);
-    string t("md_pos_");
-    t=intTurnTostring(step,t);
-	connection2();
-    printpos(t,step);
-    maxStep = sqrt(maxStep)*dt;
-/*    if(step1){
-         PDF(step1);
-         printRDF(step1);
-    }*/
-
-    // (3.1) OFDFT optimizer: in order to calculate the force
-    // Refresh the ionpositions, including the Ewald term
-    // and the local pseudopotentials term.
-   // RefreshIonPositions(grids); // This must be called every time the 
-                                          // ion positions are changed
-
-    // (3.2)
-    // Get the new charge density in the new ion positions.
-
-    // (3.3) Calculate the forces. 
-    //CalculateForces(rho, forces);
-	//xiaohui add 2014-06-13
-	//self_consistent(step);
-/*	Force_LCAO FL;
-	FL.allocate (); 
-	FL.start_force();
-
-	int ion;
-	for(ion=0;ion<numIon;ion++){
-			force[ion].x=FL.fcs(ion,1)/2;
-			force[ion].y=FL.fcs(ion,2)/2;
-			force[ion].z=FL.fcs(ion,3)/2;
-	}
-
-    //callforce();
-	for(i=0;i<numIon;i++){
-        if(ionmbl[k].x==0)force[k].x=0;
-        if(ionmbl[k].y==0)force[k].y=0;
-        if(ionmbl[k].z==0)force[k].z=0;
-	}
-    maxForce = MAXVALF();*/
-	// (3.4) 3rd step to Update velocity
-  //2014/6/28
-  /*  for( ii=0;ii<numIon;ii++){        
-      mass =  allmass[ii];
-      vel[ii] = vel[ii] + force[ii]/mass*dt/2.0;	  
-    }*/
-
-    // calculate the kinetic energy of ions
-
-    //------------------------------------------------------
-    // this part should be not necessary in NVE
-    // mohan 2013-05-09
-    // calculate the conserved quantity during MD 
-    //hamiltonian = Conserved(twiceKE/2._DP, energy(1))
-    //newKE = 0.5*twiceKE-(hamiltonian-conservedE) 
-    //newKE = 2.0*newKE;
-
-    // updates velocity in order to adjust temperature
-    //CALL RescaleVelocity(twiceKE, newKE, numIon, nfrozen)
-    //------------------------------------------------------
-
-    // calculate the conserved quantity during MD 
-    hamiltonian = Conserved(twiceKE/2, mdenergy);
-
-// kinetic, external, coulombic, exchange-correlation,
-// ion-ion, Thomas-Fermi, von Weiszacker and
-// third term Wang-Teter, WGC,
-
-
-
-    // For monitoring MD or restarting MD
-  //   OutputMDGeom(step);
-     cout<< setprecision (9)<<hamiltonian<<" "<< setprecision (9)<<twiceKE/2<<endl;
-     if (doMSD ==1){ 
-          OutputIonVelocityAndStats(step);
-     }
-     MonitorMeanSquareDisplacement(step);
-    // output the density
-    //OutputDen(step, grids(0)%rho);
-
-
-
-
-    // Print output
-    //GeometryMinimizerReportSteps(step, TimerStop(watch), &
-    //                                  forces, maxForce, maxStep, &
-    //                                  .TRUE., hamiltonian, .TRUE., &
-     //                                 temperature=twiceKE/(3*double(numIon-nfrozen)))
-
-   // watch = TimerStart();
-    
-    // Output the message to the screen.
-    if (!MY_RANK){ 
-         out<<step<<" "<< mdenergy<<" "<< hamiltonian<<" "<< mdenergy-oldEtot<<" "<< twiceKE/(3*double(numIon-nfrozen))/K_BOLTZMAN_AU<<endl;
-         out<<step<<" "<< mdenergy<<" "<< conservedE<<" "<< mdenergy-oldEtot<<" "<< twiceKE/(3*double(numIon-nfrozen))/K_BOLTZMAN_AU<<endl;
-         ofs_running<<step<<" "<< mdenergy<<" "<< hamiltonian<<" "<< mdenergy-oldEtot<<" "<< twiceKE/(3*double(numIon-nfrozen))/K_BOLTZMAN_AU<<endl;
-         ofs_running<<step<<" "<< mdenergy<<" "<< conservedE<<" "<< mdenergy-oldEtot<<" "<< twiceKE/(3*double(numIon-nfrozen))/K_BOLTZMAN_AU<<endl;
-    }
-    oldEtot=mdenergy;
-    mstout(step);
-
-//	}
-  //----------------------------------------
-  // End of NVE MD loop
-  //----------------------------------------
-
-  //GeometryMinimizerReportFooter(TimerStop(watch2));
-
-  //cout<<"(NVE):this step finished."<<endl;
-/*
-  if(force!=NULL)delete(force);
-  if(ionmbl!=NULL)delete(ionmbl);
-  if(cartNoWrap!=NULL)delete(cartNoWrap);
-  if(vel!=NULL)delete(vel);
- // if(msdcoords0!=NULL) delete(msdcoords0);
- // if(msdcoordst!=NULL) delete(msdcoordst);
-  if(taudirac!=NULL)delete(taudirac);
-  if(allmass!=NULL)delete(allmass);*/
-     cout<<"(NVE): this step finished."<<endl;
-     if (!MY_RANK){
-         out<<"(NVE): this step finished."<<endl;
-         ofs_running<<"(NVE): this step finished."<<endl;
-     }
-//2015-09-25, xiaohui
-#ifdef __MPI
-     if(step%100==0){
-        MPI_Bcast(vel,numIon*3,MPI_DOUBLE,0,MPI_COMM_WORLD);
-        MPI_Bcast(cartNoWrap,numIon*3,MPI_DOUBLE,0,MPI_COMM_WORLD);
-        MPI_Bcast(taudirac,numIon*3,MPI_DOUBLE,0,MPI_COMM_WORLD);
-        MPI_Bcast(&xLogS,1,MPI_DOUBLE,0,MPI_COMM_WORLD);
-        MPI_Bcast(&vLogS,1,MPI_DOUBLE,0,MPI_COMM_WORLD);
-     }
-#endif
-
-     timer::tick("mdNVE","runNVE",'C');
-     return;
-}
-
-
-void mdNVE::RescaleVelocity(double KE,double newKE){
-//---------------------------------------------------------------------------
-// DESCRIPTION:
-//   This function propagates the Nose-Hoover thermostat extended-system
-//   variables
-//----------------------------------------------------------------------------
-                        //>> INTERNAL VARIABLES <<//
-
-  double scaling,tempNow, tempTarget;
-
-
-                        //>> INITIALIZATION <<//
-
-  tempNow = KE/(3*double(numIon-nfrozen))/K_BOLTZMAN_AU;
-  tempTarget = newKE/(3*double(numIon-nfrozen))/K_BOLTZMAN_AU;
-
-  out<< "Temperature now     : "<< tempNow<< " (K)"<<endl;
-  out<< "New temperature  : "<< temperature/K_BOLTZMAN_AU<< " (K)"<<endl;
-
-  scaling = sqrt(tempTarget/tempNow);
-                        //>> FUNCTION BODY <<//
-
-  out<<" scaling factor of velocity : "<< scaling<<endl;
- for(int i=0;i<numIon;i++) vel[i] = vel[i]*scaling;
-  twiceKE = twiceKE*scaling*scaling;
-  return;
-}
-
-double mdNVE::Conserved(double KE, double PE){
-//---------------------------------------------------------------------------
-// DESCRIPTION:
-//   This function calculates the conserved quantity for the NVE system. 
-//----------------------------------------------------------------------------
-
-   // KE   // Kinetic energy of particles (due to their velocity)
-   // PE   // Potential energy of particles (energy calculated using OFDFT)
-
-  double Conserved; // The conserved quantity
-
-                        //>> INITIALIZATION <<//
-                        //>> FUNCTION BODY <<//
-
-   Conserved = KE + PE ;
-   
-   if (!MY_RANK){                 
-       out<< "NVE Conservation     : "<< Conserved<<" (Hartree)"<<endl;
-       out<< "NVE Temperature      : "<< twiceKE/(3*double(numIon-nfrozen))/K_BOLTZMAN_AU<<" (K)"<<endl;
-       out<< "NVE Kinetic energy   : "<< KE<<" (Hartree)"<<endl;
-       out<< "NVE Potential energy : "<< PE<<" (Hartree)"<<endl;
-       ofs_running<< "NVE Conservation     : "<< Conserved<<" (Hartree)"<<endl;
-       ofs_running<< "NVE Temperature      : "<< twiceKE/(3*double(numIon-nfrozen))/K_BOLTZMAN_AU<<" (K)"<<endl;
-       ofs_running<< "NVE Kinetic energy   : "<< KE<<" (Hartree)"<<endl;
-       ofs_running<< "NVE Potential energy : "<< PE<<" (Hartree)"<<endl;
-   }
-   return Conserved;
-}
-double mdNVE::MAXVALF(){
-	//cout<<"enter in MAXVALF"<<endl;
-	double max,force0;
-	int i,j;
-	max=0;
-	for(i=0;i<numIon;i++){
-		force0=0;
-		force0+=pow(force[i].x,2)+pow(force[i].y,2)+pow(force[i].z,2);
-		if(max<force0)max=force0;
-	}
-	return max;
-}
-
diff --git a/ABACUS.develop/source/src_pw/mdNVE.h b/ABACUS.develop/source/src_pw/mdNVE.h
deleted file mode 100644
index 317d226886..0000000000
--- a/ABACUS.develop/source/src_pw/mdNVE.h
+++ /dev/null
@@ -1,44 +0,0 @@
-#ifndef MDNVE_H
-#define MDNVE_H
-#include"md.h"
-
-using namespace std;
-class mdNVE:public md{
-public:
-	mdNVE(int n=1):md(1){};
-	void runNVE(int step1);
-	void RescaleVelocity(double KE,double newKE);
-	double Conserved(double KE, double PE);
-	double MAXVALF();
-
-private:
-	                      //>> INTERNAL VARIABLES <<//
-//  int nfrozen;               // Number of frozen atoms
-   int ii,k,it,ion;
-   int i;                       // Counters to parse through the ion table.
-   int step;                    // step number we are on
-   //int a.nat;                  // total number of atoms
-   int nStep;                   // Total number of steps
-   int errorFlag;
-
-  //vector3<double> velocity;           // Velocity of ions (atomic units)
-  Vector3<double> fracStep;           // Step (due to velocity) in fractional coordinates
-  double mass;   // atom mass, temperay var
-   double maxStep;     // Largest movement by an ion
-   double hamiltonian,conservedE,maxForce;
-   double msd;                  //  mean square displacement
-   double diffuCoeff;           // diffusion coefficient
-   double newKE;
-   double twiceKE;     // Kinetic energy x 2
-   double oldEtot;        // old energy
-
-  //TYPE(stopwatch):: &           // Timer objects
-  //  watch, &
-  //  watch2
-  
-  //string parameterInfo;      // Line containing values for Qmass, dt, temperature for output to header
-
-  double tempNow;
-
-};
-#endif
diff --git a/ABACUS.develop/source/src_pw/mdNVT.cpp b/ABACUS.develop/source/src_pw/mdNVT.cpp
deleted file mode 100644
index 42878cbc9d..0000000000
--- a/ABACUS.develop/source/src_pw/mdNVT.cpp
+++ /dev/null
@@ -1,456 +0,0 @@
-#include "mdNVT.h"
-
-const double EXP=2.71828;
-/*
-using namespace std;
-class mdnvt:public md{
-public:	mdnvt(int n=1):md(1){};
-    void runnvt();
-	void NHIntegrator();
-	double NHhamiltonian(double KE,double PE);
-	double MAXVALF();
-private:
-	int  nfrozen;//                 ! Number of frozen atoms
-  int  ii,k,it,ion;//
-  int  i;//                       ! Counters to parse through the ion table.
-  int  step;//                    ! step number we are on
-  //int  numIon;//                  ! total number of atoms
-  int  nStep;//                   ! Total number of steps
-  int  errorFlag;  
-  double  mass;            // atom mass, temperay var
-  double  maxStep;       // Largest movement by an ion
-  //double  xLogS;         // position of thermostat
-  //double  vLogS;         // vel of thermostat
-  double  hamiltonian;   // conserved energy
-  double  maxForce;      
-  //double  msd;           // mean square displacement
-  //double  diffuCoeff;      // diffusion coefficient
-  double  twiceKE;         // Kinetic energy x 2
-  double  oldEtot;           // old energy
-};
-*/
-void mdnvt::runnvt(int step1){
-//------------------------------------------------------------------------------
-// DESCRIPTION:
-//  This subroutine integrates motion for molecular dynamics using the
-//  Nose-Hoover thermostat coupled with a predictor-corrector integrator.
-// REFERENCES:
-//  Martyna et. al., "Explicit reversible integrators for extended system
-//     dynamics," Mol. Phys. 87 (1996), 1117.
-//  Allen and Tildesley, "Computer Simulation of Liquids" (predictor-corrector)
-//  Frenkel and Smit, "Understanding Molecular Simulation".
-//------------------------------------------------------------------------------
-//REVISION LOG:
-//    10/20/2008 : Adapted from subroutines in IonOptimizers and symplectic
-//                 algorithm from Martyna et. al. with help from Frenkel/Smit.
-//
-//------------------------------------------------------------------------------
-//  Vector4 rho;   // Electron density, realspace
-//  Vector energy;      // Energy from electronic optimization
-  //Vector3 forces;// Total forces.  Final index is 0 for total force,
-  // First is ion number, second direction (1,2,3 for x,y,z)
-
-  
-  //REAL(KIND=DP), ALLOCATABLE, DIMENSION(:,:) :: vel // vel of ions (atomic units)
-  //double intCoeff[];
-  //REAL(KIND=DP), DIMENSION(3) :: fracStep // Step (due to vel) in fractional coordinates
-
-  //string message;
-  //string parameterInfo;      // Line containing values for Qmass, dt, temperature
-                                          // for output to header
-
-
-                      //  !>> INITIALIZATION <<!
-  TITLE("mdnvt","runnvt");
-   timer::tick("mdnvt","runnvt",'C');
-  step=step1+step_rst;
-  if (!MY_RANK)out<<"step: " << step <<endl;
-  if (step!=1)ReadNewTemp( step );
-  intCoeff=new double[nYosh];
-
-  //xiaohui modify 2015-09-15
-  //if(step==1)cout<<" (NVT) Start."<<endl;
-  //else {
-  //       cout<<" (NVT) step: "<<step<<endl;
-  //}
-  //// Calculate # of frozen atoms 
-  //cout<<" (NVT) Fronzen Atom Number : "<<nfrozen<<endl;
-
-
-  MakeIntCoeff();
-
-  //xiaohui modify 2015-09-15
-  //// Set up extra output to ion optimizer / MD header
-  //cout<<" Qmass for NVT       : "<< Qmass<<" (a.u.)"<<endl;
-  //cout<< " Time interval       : "<< dt*fundamentalTime*1e15<< " (fs)"<<endl;
-
-  //xiaohui add 2015-09-15
-  cout << " ---------------------------------------------------------" << endl;
-
-  cout<< " Target temperature  : "<< temperature/K_BOLTZMAN_AU<< " (K)"<<endl;
- // cout<< " Number of NVT steps : "<< nStep<<endl;
-
-  //if (rankGlobal==0)  GeometryMinimizerReportHeader(parameterInfo);
-
-  if(step==1||step%fixTemperature==1){
-      for(k=0;k<numIon;k++){
-        if(ionmbl[k].x==0)vel[k].x=0;
-        if(ionmbl[k].y==0)vel[k].y=0;
-        if(ionmbl[k].z==0)vel[k].z=0;
-      }
-      if(nfrozen==0)RemoveMovementOfCenterOfMass(); //added 20170626
-      scalevel();
-  }
-
-  // get the kinetic energy
-  twiceKE = GetAtomKE();
-  twiceKE = twiceKE * 2;
- /* vLogS=0;
-  xLogS=0;*/
-   //Set up forces 
-  //Optimizer;                        //Optimize the density
-  callforce();
-	if(STRESS)
-	{
-		cout<<"output Pressure for check!"<<endl;
-		double press;
-		for(int i=0;i<3;i++)
-			press += stress_lcao(i,i)/3;
-		press += twiceKE/3/ucell.omega; //output virtual press = 2/3 *Ek/V + sum(sigma[i][i])/3
-		double unit_transform = RYDBERG_SI / pow(BOHR_RADIUS_SI,3) * 1.0e-8 ;
-		cout<<"Virtual Pressure is "<<press*unit_transform<<" Kbar "<<endl;
-	}
-  for(k=0;k<numIon;k++){
-     if(ionmbl[k].x==0)force[k].x=0;
-     if(ionmbl[k].y==0)force[k].y=0;
-     if(ionmbl[k].z==0)force[k].z=0;
-  }
-  maxForce = MAXVALF();
-  mdenergy=en.etot/2;
-  hamiltonian = NHhamiltonian(twiceKE/2, mdenergy);
-  //xiaohui move this line 2015-09-15
-  //cout<<"ready part is ready!"<<endl;
-
-  // Print intial output
-  //GeometryMinimizerReportSteps(startStep, TimerStop(watch), forces, maxForce,0,1, hamiltonian,1,temperature)
-  //----------------------------------------------
-  // big loop
-  //-----------------------------------------------
-  cout<<" "<<std::left<<setw(12)<<"MD_STEP"<<std::left<<setw(12)<< "SystemE"<<std::left<<setw(12)<< "Conserved"<<std::left<<setw(12)<< "DeltaE"<<std::left<<setw(12)<< "Temperature"<<endl;
-
-  cout<<" "<<std::left<<setw(12)<<step<<std::left<<setw(12)<< mdenergy<<std::left<<setw(12)<< hamiltonian<<std::left<<setw(12)<< mdenergy-oldEtot<<std::left<<setw(12)<<twiceKE/(3*double(numIon-nfrozen))/K_BOLTZMAN_AU<<endl;
- 
-  oldEtot=mdenergy;
-
- 
-  // Loop for MD step 
-
-    /*xiaohui move 2015-09-15
-    //cout<<" --------------------------------------------------"<<endl;
-    cout << " ---------------------------------------------------------" << endl;
-    cout<<" Molecular Dynamics (NVT) STEP "<< step<<endl;
-    */
-
-    //cout<< " --------------------------------------------------"<<endl;
-    //cout << " ---------------------------------------------------------" << endl;
-    if (!MY_RANK){
-          out<<" --------------------------------------------------"<<endl;
-          out<<" Molecular Dynamics (NVT) STEP "<< step<<endl;
-          out<< " --------------------------------------------------"<<endl;
-          ofs_running<<" --------------------------------------------------"<<endl;
-          ofs_running<<" Molecular Dynamics (NVT) STEP "<< step<<endl;
-          ofs_running<< " --------------------------------------------------"<<endl;
-    }
-    // (1) Calculate the Mean-Square-Displacement.
-    if(step==1&&rstMD==0){ 
-    // (2) First themorstat step: updates vel, twiceKE, xLogS, vLogS
-   // MonitorMeanSquareDisplacement(step);
-         NHIntegrator();
-       
-    // (3) New vel obtained (Verlet-like step)
-        for( ii=0;ii<numIon;ii++){ 
-            mass = allmass[ii];
-            vel[ii] = vel[ii] + force[ii]/mass*dt/2.0;
-	}
-    }
-    else {
-        for( ii=0;ii<numIon;ii++){ 
-            mass = allmass[ii];
-            vel[ii] = vel[ii] + force[ii]/mass*dt/2.0;
-        }
-        twiceKE=GetAtomKE();
-        twiceKE = 2 * twiceKE;
-        NHIntegrator();
-        twiceKE=GetAtomKE();
-        twiceKE = 2 * twiceKE;
-        hamiltonian = NHhamiltonian(twiceKE/2, mdenergy);
-        NHIntegrator();
-        for( ii=0;ii<numIon;ii++){ 
-            mass = allmass[ii];
-            vel[ii] = vel[ii] + force[ii]/mass*dt/2.0;
-        }
-    }
-    // (4) Update the Non-Wrapped cartesion coordinates
-        if(nfrozen==0)RemoveMovementOfCenterOfMass();  
-        if(mdtype==2)scalevel();//choose Nose Hoover plus velocity scaling method
-	for( ii=0;ii<numIon;ii++){ 
-            cartNoWrap[ii] = vel[ii]*dt +cartNoWrap[ii];
-	}
-    
-        MonitorMeanSquareDisplacement(step);
-    // Calculate the maximal velocities.
-    // The step in fractional coordinates
-    maxStep = 0;
-    Vector3<double> fracStep;
-   	for( ii=0;ii<numIon;ii++){ 
-	    	if((pow(vel[ii].x,2)+pow(vel[ii].y,2)+pow(vel[ii].z,2))>maxStep)
-                maxStep = pow(vel[ii].x,2)+pow(vel[ii].y,2)+pow(vel[ii].z,2);
-                Mathzone::Cartesian_to_Direct(vel[ii].x*dt/ucell.lat0,vel[ii].y*dt/ucell.lat0,vel[ii].z*dt/ucell.lat0,
-                                    ionlatvec.e11,ionlatvec.e12,ionlatvec.e13,
-                                    ionlatvec.e21,ionlatvec.e22,ionlatvec.e23,
-                                    ionlatvec.e31,ionlatvec.e32,ionlatvec.e33,
-                                    fracStep.x,fracStep.y,fracStep.z);
-
-//            fracStep = ionlatvec.Inverse()*vel[ii]*dt/ucell.lat0; 
-            taudirac[ii] = taudirac[ii] + fracStep;
-	}
-
-    moveatoms(step);
-    string t("md_pos_");
-    t=intTurnTostring(step,t);
-    connection2();
-    printpos(t,step);
-    maxStep = sqrt(maxStep)*dt;
-/*   if(step1){
-         PDF(step1);
-         printRDF(step1);
-    }*/
-
-    // (5) OFDFT optimizer
-    // Refresh the ionpositions, including the Ewald term
-    // and the local pseudopotentials term.
-    //RefreshIonPositions(grids); // This must be called every time the 
-                                    // ion positions are changed
-
-    // Get the new charge density in the new ion positions.
-     //Optimizer;                        //Optimize the density
-    //CalculateStress( energy, stress);
-    //PrintStress(stress);
-
-
-    // (8)
-
-
-    
-    // (9) Second thermostat step: updates vel, twiceKE, xLogS, vLogS
-// cout<<"testNVT: x "<<xLogS<<" v "<<vLogS<<endl;
-
-    // For monitoring MD or restarting MD   
-  //   OutputMDGeom(step);
-
-     //xiaohui move this line, 2015-09-15
-     //cout<<" ekin " << setprecision (9)<<twiceKE/2<<" etot "<< setprecision (9)<<hamiltonian<<endl;
-
-     if (doMSD ==1){
-     OutputIonVelocityAndStats(step);
-}
-
-    //OutputDen(step, grids(0)%rho);
-
-    // (10) Conserved quantity during MD 
-
-// kinetic, external, coulombic, exchange-correlation,
-// ion-ion, Thomas-Fermi, von Weiszacker and
-// third term Wang-Teter, WGC,
-    if (!MY_RANK){
-         out<<setw(15)<<"maxForce=     "<<setw(15)<<"maxstep=      "<<setw(15)<<"step=     "<<endl;
-         out<<setw(15)<<maxForce<<setw(15)<<maxStep<<setw(15)<<step<<endl;
-
-    //watch = TimerStart()
-    
-    // Change temperature if needed
-
-    // Output the message to the screen.
-         out<<step<<" "<< mdenergy<<" "<< hamiltonian<<" "<< mdenergy-oldEtot<<" "<< twiceKE/(3*double(numIon-nfrozen))/K_BOLTZMAN_AU<<endl;
-    }
-    oldEtot=mdenergy;
-    mstout(step);
-    double thisRDF=0;
-    if(step%fixTemperature==0){
-        for(i=42;i<48;i++){
-           thisRDF+=rdf[i];
-        }
-  //      cout<<"RDFthirdpeak "<<step/fixTemperature<<" "<<thisRDF<<endl;
-    }
-
- //----------------------------------------
-  // End of NVT MD loop
-  //----------------------------------------
-
-  // GeometryMinimizerReportFooter(TimerStop(watch2));
-
-  //xiaohui move these two lines, 2015-09-15
-  //cout<<"(NVT): finished."<<endl;
-  //if (!MY_RANK)  ofs_running<<"(NVT): finished."<<endl;
-
-  //delete(cartNoWrap);
-  //delete(vel);
-  delete []intCoeff;
-//2015-09-25, xiaohui
-#ifdef __MPI
-  if(step%100==0){
-    MPI_Bcast(vel,numIon*3,MPI_DOUBLE,0,MPI_COMM_WORLD);
-    MPI_Bcast(cartNoWrap,numIon*3,MPI_DOUBLE,0,MPI_COMM_WORLD);
-    MPI_Bcast(taudirac,numIon*3,MPI_DOUBLE,0,MPI_COMM_WORLD);
-    MPI_Bcast(&xLogS,1,MPI_DOUBLE,0,MPI_COMM_WORLD);
-    MPI_Bcast(&vLogS,1,MPI_DOUBLE,0,MPI_COMM_WORLD);
-  }
-#endif
-  //IF (ALLOCATED(msd_coords0)) DEALLOCATE(msd_coords0)
-  //IF (ALLOCATED(msd_coordst)) DEALLOCATE(msd_coordst)
-   timer::tick("mdnvt","runnvt",'C');
-  return;
-}
-
-
-void mdnvt::NHIntegrator(){
-//---------------------------------------------------------------------------
-// DESCRIPTION:
-//   This function propagates the Nose-Hoover thermostat extended-system
-//   variables using multiple step technique.
-//   Note that only thermostat variables are updated, including
-//   thermostat positions (xLogS),
-//   thermostat velocities (vLogS),
-//   thermostat acceleration (gLogS),
-//   but not including any variables for atoms. 
-//--------------------------------------------------------------------
-
-  double freedom,scaling,gLogS;
-
-  int iResn, iYosh;
-  double wdt2, wdt4;
-
-                    //    !>> INITIALIZATION <<!
-
-  
-
-  freedom = 3*double(numIon-nfrozen);
-  scaling = 1.0;
-
-  // calculate the force of thermostat
-  gLogS = (twiceKE - freedom*temperature) / Qmass;
-//  cout<<"testNVT: x "<<xLogS<<" v "<<vLogS<<" g "<<gLogS<<" ekin "<<twiceKE<<" temperature " <<temperature/K_BOLTZMAN_AU<<endl;
-                        //!>> FUNCTION BODY <<!
-
-//  WRITE(outputUnit,*) " freedom=",freedom
-//  WRITE(outputUnit,*) " gLogS=",gLogS
-//  WRITE(outputUnit,*) " twiceKE=",twiceKE
-//  WRITE(outputUnit,*) " temperature=",temperature
-  
-  for( iResn = 0;iResn< nResn;iResn++){
-	  for( iYosh = 0;iYosh< nYosh;iYosh++){
-      wdt2 = intCoeff[iYosh]*dt/(2*double(nResn));
-      wdt4 = intCoeff[iYosh]*dt/(4*double(nResn));
-
-      // calculate the first half of integration.
-      vLogS = vLogS + gLogS*wdt4;
-      scaling = scaling * pow(EXP,(-wdt2*vLogS));
-
-      // In some cases, the exponential part will be NaN,
-      // Try to give warning at this case, added by mohan
-     // if( check_isnan(scaling) ){
-      //  out<< " iResn="<< iResn<< " iYosh="<< iYosh<<endl;
-      //  out<< " gLogS="<< gLogS<<endl;
-      //  out<< " vLogS="<< vLogS<<endl;
-      //  out<<"Scaling in NHIntegrator is NaN !"<<endl;
-      //  STOP
-      //}
-
-      gLogS = (scaling*scaling*twiceKE - freedom*temperature) / Qmass;
-      // Calculate xLogS only as a check of code
-      xLogS = xLogS + vLogS*wdt2;
-
-      // calculate the second half of integration
-      vLogS = vLogS + gLogS*wdt4;
-	  }
-  }
-  
-  if (!MY_RANK)out<< " scaling factor of vel : "<< scaling<<endl;
-  for(i=0;i<numIon;i++){
-    vel[i] = vel[i]*scaling;
-  }
-  twiceKE = twiceKE*scaling*scaling;
-  return;
-}
-
-
-double mdnvt::NHhamiltonian(double KE,double PE){
-//---------------------------------------------------------------------------
-// DESCRIPTION:
-//   This function calculates the conserved quantity for the Nose-Hoover
-//   system.
-//----------------------------------------------------------------------------
- 
-
-//  REAL(kind=DP), INTENT(IN) :: &
-//    KE, &    ! Kinetic energy of particles (due to their vel)
-//    PE, &    ! Potential energy of particles (energy calculated using OFDFT)
-//    xLogS, & ! "position" of constraint variable
-//    vLogS    ! "vel" of constraint variable
-
-  double NHhamiltonian0; // The conserved quantity
-
-                       // !>> INTERIOR VARIABLES <<!
-
-  NHhamiltonian0 = KE + PE + 0.5*vLogS*vLogS*Qmass +(3*double(numIon-nfrozen))*temperature*xLogS;
-   
-   if (!MY_RANK){
-   out<< " "<<endl;
-   out<< " --------------------------------------------------"<<endl;
-   out<< " SUMMARY OF NVT CALCULATION"<<endl;
-   out<<" --------------------------------------------------"<<endl;
-   out<<" NVT Conservation     : "<<setw(10)<< NHhamiltonian0*2<<" (Rydberg)"<<endl;
-   out<<" NVT Temperature      : "<<setw(10)<< twiceKE/(3*double(numIon-nfrozen))/K_BOLTZMAN_AU<<" (K)"<<endl;
-   out<<" NVT Kinetic energy   : "<<setw(10)<< KE*2<<" (Rydberg)"<<endl;
-   out<<" NVT Potential energy : "<<setw(10)<< PE*2<<" (Rydberg)"<<endl;
-   out<< " Thermostat Position : "<<setw(10)<< xLogS<<endl;
-   out<< " Thermostat vel : "<< vLogS<<endl;
-   ofs_running<< " "<<endl;
-   ofs_running<< " --------------------------------------------------"<<endl;
-   ofs_running<< " SUMMARY OF NVT CALCULATION"<<endl;
-   ofs_running<<" --------------------------------------------------"<<endl;
-   ofs_running<<" NVT Conservation     : "<<setw(10)<< NHhamiltonian0*2<<" (Rydberg)"<<endl;
-   ofs_running<<" NVT Temperature      : "<<setw(10)<< twiceKE/(3*double(numIon-nfrozen))/K_BOLTZMAN_AU<<" (K)"<<endl;
-   ofs_running<<" NVT Kinetic energy   : "<<setw(10)<< KE*2<<" (Rydberg)"<<endl;
-   ofs_running<<" NVT Potential energy : "<<setw(10)<< PE*2<<" (Rydberg)"<<endl;
-   ofs_running<< " Thermostat Position : "<<setw(10)<< xLogS<<endl;
-   ofs_running<< " Thermostat vel : "<< vLogS<<endl;
-
-   }
-   return NHhamiltonian0;
-}
-
-// check 
-//LOGICAL FUNCTION check_isnan(a) 
-//REAL(KIND=DP) a 
-//if (a.ne.a) then 
-//check_isnan = .true. 
-//else 
-//check_isnan = .false. 
-//end if 
-//return 
-//END FUNCTION check_isnan
-double mdnvt::MAXVALF(){
-	double max,force0;
-	int i,j;
-	max=0;
-	for(i=0;i<numIon;i++){
-		force0=0;
-		for(j=0;j<3;j++){
-			force0+=pow(force[i].x,2)+pow(force[i].y,2)+pow(force[i].z,2);
-		}
-		if(max<force0)max=force0;
-	}
-	return max;
-}
-
diff --git a/ABACUS.develop/source/src_pw/mdNVT.h b/ABACUS.develop/source/src_pw/mdNVT.h
deleted file mode 100644
index db5d874f94..0000000000
--- a/ABACUS.develop/source/src_pw/mdNVT.h
+++ /dev/null
@@ -1,32 +0,0 @@
-#ifndef MDNVT_H
-#define MDNVT_H
-
-#include "md.h"
-
-using namespace std;
-class mdnvt:public md{
-public:	mdnvt(int n=1):md(1){};
-    void runnvt(int step1);
-	void NHIntegrator();
-	double NHhamiltonian(double KE,double PE);
-	double MAXVALF();
-private:
-//	int  nfrozen;//                 ! Number of frozen atoms
-  int  ii,k,it,ion;//
-  int  i;//                       ! Counters to parse through the ion table.
-  int  step;//                    ! step number we are on
-  //int  numIon;//                  ! total number of atoms
-  int  nStep;//                   ! Total number of steps
-  int  errorFlag;  
-  double  mass;            // atom mass, temperay var
-  double  maxStep;       // Largest movement by an ion
-  //double  xLogS;         // position of thermostat
-  //double  vLogS;         // vel of thermostat
-  double  hamiltonian;   // conserved energy
-  double  maxForce;      
-  //double  msd;           // mean square displacement
-  //double  diffuCoeff;      // diffusion coefficient
-  double  twiceKE;         // Kinetic energy x 2
-  double  oldEtot;           // old energy
-};
-#endif
diff --git a/ABACUS.develop/source/src_pw/unitcell.cpp b/ABACUS.develop/source/src_pw/unitcell.cpp
index 518cc9dbab..30e3d33695 100644
--- a/ABACUS.develop/source/src_pw/unitcell.cpp
+++ b/ABACUS.develop/source/src_pw/unitcell.cpp
@@ -293,6 +293,22 @@ void UnitCell::update_pos_tau(const double* pos)
     return;
 }
 
+void UnitCell::update_pos_taud(const Vector3<double>* posd_in)
+{
+    int iat = 0;
+    for(int it = 0;it < this->ntype;it++)
+    {
+        Atom* atom = &this->atoms[it];
+        for(int ia =0;ia< atom->na;ia++)
+        {
+            this->atoms[it].taud[ia] += posd_in[ia];
+            iat++;
+        }
+    }
+    assert(iat == this->nat);
+    this->periodic_boundary_adjustment();
+}
+
 void UnitCell::periodic_boundary_adjustment()
 {
     //----------------------------------------------
diff --git a/ABACUS.develop/source/src_pw/unitcell.h b/ABACUS.develop/source/src_pw/unitcell.h
index 9d60a01af7..d634e20fce 100644
--- a/ABACUS.develop/source/src_pw/unitcell.h
+++ b/ABACUS.develop/source/src_pw/unitcell.h
@@ -67,6 +67,7 @@ class UnitCell
     const double& getNelec(void)const {return electrons_number;}
 
     void update_pos_tau(const double* pos);
+    void update_pos_taud(const Vector3<double>* posd_in);
     void periodic_boundary_adjustment();
     void bcast_atoms_tau();
     void save_cartesian_position(double* pos)const;