Merge pull request #857 from YuLiu98/develop

dyzheng · web-flow · commit 65d35ba930cf · 2022-04-19T11:55:00.000+08:00
update the default value of md_tfreq
diff --git a/docs/input-main.md b/docs/input-main.md
@@ -1028,9 +1028,11 @@ This part of variables are used to control the molecular dynamics calculations.
 
 - **Type**: Real
 - **Description**:
-  - Oscillation frequency, used to determine Qmass of NHC;
-  - 1/(md_tfreq*md_dt) is collision probability in Anderson method.
-- **Default**: 1.0
+  - When md_type = 1, md_tfreq controls the frequency of the temperature oscillations during the simulation. If it is too large, the
+temperature will fluctuate violently; if it is too small, the temperature will take a very long time to equilibrate with the atomic system. 
+  - When md_type = 3, md_tfreq*md_dt is the collision probability in Anderson method.
+  - If md_tfreq is not set in INPUT, md_tfreq will be autoset to be 1/40/md_dt.
+- **Default**: 1/40/md_dt
 
 #### md_mnhc
 
diff --git a/source/module_md/FIRE.cpp b/source/module_md/FIRE.cpp
@@ -80,11 +80,11 @@ void FIRE::second_half()
     ModuleBase::timer::tick("FIRE", "second_half");
 }
 
-void FIRE::outputMD()
+void FIRE::outputMD(std::ofstream &ofs)
 {
-    Verlet::outputMD();
+    Verlet::outputMD(ofs);
 
-    GlobalV::ofs_running << " LARGEST GRAD (eV/A)  : " 
+    ofs << " LARGEST GRAD (eV/A)  : " 
         << max * ModuleBase::Hartree_to_eV * ModuleBase::ANGSTROM_AU << std::endl;
     std::cout << " LARGEST GRAD (eV/A)  : " 
         << max * ModuleBase::Hartree_to_eV * ModuleBase::ANGSTROM_AU << std::endl;
diff --git a/source/module_md/FIRE.h b/source/module_md/FIRE.h
@@ -12,7 +12,7 @@ class FIRE : public Verlet
     void setup(ModuleESolver::ESolver *p_ensolve);
     void first_half();
     void second_half();
-    void outputMD();
+    void outputMD(std::ofstream &ofs);
     void write_restart();
     void restart();
     void check_force();
diff --git a/source/module_md/Langevin.cpp b/source/module_md/Langevin.cpp
@@ -45,9 +45,9 @@ void Langevin::second_half()
     ModuleBase::timer::tick("Langevin", "second_half");
 }
 
-void Langevin::outputMD()
+void Langevin::outputMD(std::ofstream &ofs)
 {
-    Verlet::outputMD();
+    Verlet::outputMD(ofs);
 }
 
 void Langevin::write_restart()
diff --git a/source/module_md/Langevin.h b/source/module_md/Langevin.h
@@ -12,7 +12,7 @@ class Langevin : public Verlet
     void setup(ModuleESolver::ESolver *p_ensolve);
     void first_half();
     void second_half();
-    void outputMD();
+    void outputMD(std::ofstream &ofs);
     void write_restart();
     void restart();
     void post_force();
diff --git a/source/module_md/MD_parameters.h b/source/module_md/MD_parameters.h
@@ -31,7 +31,7 @@ class MD_parameters
 		msst_tscale = 0.01;
 
 		// NHC
-		md_tfreq = 1;
+		md_tfreq = 0;
 		md_mnhc = 4;
 
 		// Langevin
diff --git a/source/module_md/MSST.cpp b/source/module_md/MSST.cpp
@@ -160,9 +160,9 @@ void MSST::second_half()
     ModuleBase::timer::tick("MSST", "second_half");
 }
 
-void MSST::outputMD()
+void MSST::outputMD(std::ofstream &ofs)
 {
-    Verlet::outputMD();
+    Verlet::outputMD(ofs);
 }
 
 void MSST::write_restart()
diff --git a/source/module_md/MSST.h b/source/module_md/MSST.h
@@ -12,7 +12,7 @@ class MSST : public Verlet
     void setup(ModuleESolver::ESolver *p_ensolve);
     void first_half();
     void second_half();
-    void outputMD();
+    void outputMD(std::ofstream &ofs);
     void write_restart();
     void restart();
     double extra_term();
diff --git a/source/module_md/NVE.cpp b/source/module_md/NVE.cpp
@@ -36,9 +36,9 @@ void NVE::second_half()
     ModuleBase::timer::tick("NVE", "second_half");
 }
 
-void NVE::outputMD()
+void NVE::outputMD(std::ofstream &ofs)
 {
-    Verlet::outputMD();
+    Verlet::outputMD(ofs);
 }
 
 void NVE::write_restart()
diff --git a/source/module_md/NVE.h b/source/module_md/NVE.h
@@ -12,7 +12,7 @@ class NVE : public Verlet
     void setup(ModuleESolver::ESolver *p_ensolve);
     void first_half();
     void second_half();
-    void outputMD();
+    void outputMD(std::ofstream &ofs);
     void write_restart();
     void restart();
 
diff --git a/source/module_md/NVT_ADS.cpp b/source/module_md/NVT_ADS.cpp
@@ -7,9 +7,6 @@
 
 NVT_ADS::NVT_ADS(MD_parameters& MD_para_in, UnitCell_pseudo &unit_in) : Verlet(MD_para_in, unit_in)
 {
-    // convert to a.u. unit
-    mdp.md_tfreq *= ModuleBase::AU_to_FS;
-
     nraise = mdp.md_tfreq * mdp.md_dt;
 }
 
@@ -47,7 +44,7 @@ void NVT_ADS::second_half()
         double deviation;
         for(int i=0; i<ucell.nat; ++i)
         {
-            if(rand()/double(RAND_MAX) <= 1.0/nraise)
+            if(rand()/double(RAND_MAX) <= nraise)
             {
                 deviation = sqrt(mdp.md_tlast / allmass[i]);
                 for(int k=0; k<3; ++k)
@@ -68,9 +65,9 @@ void NVT_ADS::second_half()
     ModuleBase::timer::tick("NVT_ADS", "second_half");
 }
 
-void NVT_ADS::outputMD()
+void NVT_ADS::outputMD(std::ofstream &ofs)
 {
-    Verlet::outputMD();
+    Verlet::outputMD(ofs);
 }
 
 void NVT_ADS::write_restart()
diff --git a/source/module_md/NVT_ADS.h b/source/module_md/NVT_ADS.h
@@ -12,7 +12,7 @@ class NVT_ADS : public Verlet
     void setup(ModuleESolver::ESolver *p_ensolve);
     void first_half();
     void second_half();
-    void outputMD();
+    void outputMD(std::ofstream &ofs);
     void write_restart();
     void restart();
 
diff --git a/source/module_md/NVT_NHC.cpp b/source/module_md/NVT_NHC.cpp
@@ -7,9 +7,6 @@
 
 NVT_NHC::NVT_NHC(MD_parameters& MD_para_in, UnitCell_pseudo &unit_in) : Verlet(MD_para_in, unit_in)
 {
-    // convert to a.u. unit
-    mdp.md_tfreq *= ModuleBase::AU_to_FS;
-
     if(mdp.md_tfirst == 0)
     {
         std::cout << " md_tfirst must be larger than 0 in NHC !!! " << std::endl;
@@ -92,9 +89,9 @@ void NVT_NHC::second_half()
     ModuleBase::timer::tick("NVT_NHC", "second_half");
 }
 
-void NVT_NHC::outputMD()
+void NVT_NHC::outputMD(std::ofstream &ofs)
 {
-    Verlet::outputMD();
+    Verlet::outputMD(ofs);
 }
 
 void NVT_NHC::write_restart()
@@ -197,6 +194,12 @@ void NVT_NHC::integrate()
             }
 
             scale *= exp(-veta[0]*delta/2.0);
+            if(!isfinite(scale))
+            {
+                std::cout << " Please set a proper md_tfreq !!! " << std::endl;
+                ModuleBase::WARNING_QUIT("NVT_NHC", " Please set a proper md_tfreq !!! ");
+            }
+            
             KE = kinetic * scale * scale;
 
             // update force
diff --git a/source/module_md/NVT_NHC.h b/source/module_md/NVT_NHC.h
@@ -12,7 +12,7 @@ class NVT_NHC : public Verlet
     void setup(ModuleESolver::ESolver *p_ensolve);
     void first_half();
     void second_half();
-    void outputMD();
+    void outputMD(std::ofstream &ofs);
     void write_restart();
     void restart();
     void integrate();
diff --git a/source/module_md/run_md_classic.cpp b/source/module_md/run_md_classic.cpp
@@ -59,6 +59,8 @@ void Run_MD_CLASSIC::classic_md_line(void)
     // md cycle
     while ( (verlet->step_ + verlet->step_rst_) <= GlobalV::MD_NSTEP && !verlet->stop )
     {
+        Print_Info::print_screen(0, 0, verlet->step_ + verlet->step_rst_);
+
         if(verlet->step_ == 0)
         {
             verlet->setup(p_esolver);
@@ -79,8 +81,8 @@ void Run_MD_CLASSIC::classic_md_line(void)
 
         if((verlet->step_ + verlet->step_rst_) % verlet->mdp.md_dumpfreq == 0)
         {
-            Print_Info::print_screen(0, 0, verlet->step_ + verlet->step_rst_);
-            verlet->outputMD();
+            // Print_Info::print_screen(0, 0, verlet->step_ + verlet->step_rst_);
+            verlet->outputMD(GlobalV::ofs_running);
 
             MD_func::MDdump(verlet->step_ + verlet->step_rst_, verlet->ucell, verlet->virial, verlet->force);
         }
diff --git a/source/module_md/verlet.cpp b/source/module_md/verlet.cpp
@@ -25,6 +25,15 @@ Verlet::Verlet(MD_parameters& MD_para_in, UnitCell_pseudo &unit_in):
     mdp.md_tfirst /= ModuleBase::Hartree_to_K;
     mdp.md_tlast /= ModuleBase::Hartree_to_K;
 
+    if(mdp.md_tfreq == 0)
+    {
+        mdp.md_tfreq = 1.0/40.0/mdp.md_dt;
+    }
+    else
+    {
+        mdp.md_tfreq *= ModuleBase::AU_to_FS;
+    }
+
     // LJ parameters
     mdp.lj_rcut *= ModuleBase::ANGSTROM_AU;
 	mdp.lj_epsilon /= ModuleBase::Hartree_to_eV;
@@ -94,20 +103,10 @@ void Verlet::second_half()
     }
 }
 
-void Verlet::outputMD()
+void Verlet::outputMD(std::ofstream &ofs)
 {
     if(GlobalV::MY_RANK) return;
 
-	std::ofstream ofs;
-	if(step_ + step_rst_ == 0)
-	{
-		ofs.open("run.log", ios::trunc);
-	}
-	else
-	{
-		ofs.open("run.log", ios::app);
-	}
-
     temperature_ = 2*kinetic/(double(3*ucell.nat-frozen_freedom_))*ModuleBase::Hartree_to_K;
 
     const double unit_transform = ModuleBase::HARTREE_SI / pow(ModuleBase::BOHR_RADIUS_SI,3) * 1.0e-8;
@@ -117,7 +116,7 @@ void Verlet::outputMD()
         press += stress(i,i)/3;
     }
 
-    // std::cout << std::setprecision(6) << std::setiosflags(ios::showpos) << std::setiosflags(ios::fixed);
+    std::cout << " ------------------------------------------------------------------------------------------------" << std::endl;
     std::cout << " " << std::left << std::setw(20) << "Energy" 
             << std::left << std::setw(20) << "Potential" 
             << std::left << std::setw(20) << "Kinetic" 
@@ -128,25 +127,25 @@ void Verlet::outputMD()
             << std::left << std::setw(20) << kinetic
             << std::left << std::setw(20) << temperature_
             << std::left << std::setw(20) << press*unit_transform <<std::endl;
-	std::cout << " ------------------------------------------------------------" << std::endl;
-    MD_func::outStress(virial, stress);
-
-    if(step_ == 0)
-	{
-		ofs << std::left << std::setw(20) << "Step" 
-            << std::left << std::setw(20) << "Energy" 
-            << std::left << std::setw(20) << "Potential" 
-            << std::left << std::setw(20) << "Kinetic" 
-            << std::left << std::setw(20) << "Temperature" 
-            << std::left << std::setw(20) << "Pressure (KBAR)" <<std::endl;
-	}
-    ofs << std::left << std::setw(20) << step_ + step_rst_
-        << std::left << std::setw(20) << potential+kinetic
+	std::cout << " ------------------------------------------------------------------------------------------------" << std::endl;
+
+    ofs << std::endl;
+    ofs << std::endl;
+    ofs << " ------------------------------------------------------------------------------------------------" << std::endl;
+	ofs << " " << std::left << std::setw(20) << "Energy" 
+        << std::left << std::setw(20) << "Potential" 
+        << std::left << std::setw(20) << "Kinetic" 
+        << std::left << std::setw(20) << "Temperature" 
+        << std::left << std::setw(20) << "Pressure (KBAR)" <<std::endl;
+    ofs << " " << std::left << std::setw(20) << potential+kinetic
         << std::left << std::setw(20) << potential
         << std::left << std::setw(20) << kinetic
         << std::left << std::setw(20) << temperature_
         << std::left << std::setw(20) << press*unit_transform <<std::endl;
-    ofs.close();
+    ofs << " ------------------------------------------------------------------------------------------------" << std::endl;
+    MD_func::outStress(virial, stress);
+    ofs << std::endl;
+    ofs << std::endl;
 }
 
 void Verlet::write_restart()
diff --git a/source/module_md/verlet.h b/source/module_md/verlet.h
@@ -15,7 +15,7 @@ class Verlet
     virtual void setup(ModuleESolver::ESolver *p_esolver);
     virtual void first_half();
     virtual void second_half();
-    virtual void outputMD();
+    virtual void outputMD(std::ofstream &ofs);
     virtual void write_restart();
     virtual void restart();
 
diff --git a/source/src_io/print_info.cpp b/source/src_io/print_info.cpp
@@ -37,18 +37,30 @@ void Print_Info::setup_parameters(UnitCell_pseudo &ucell, K_Vectors &kv)
 
 			std::cout << " ---------------------------------------------------------" << std::endl;
 
-			if(INPUT.mdp.md_type ==1 || INPUT.mdp.md_type==3)
-			{
-				std::cout << " ENSEMBLE                 : " << "NVT" << std::endl;
-			}
-			else if(INPUT.mdp.md_type==0)
-			{
-				std::cout << " ENSEMBLE                 : " << "NVE" << std::endl;
-			}
-			
-			std::cout << " Time interval(fs)        : " << INPUT.mdp.md_dt << std::endl;
-		}
-		std::cout << " ---------------------------------------------------------" << std::endl;
+            if(INPUT.mdp.md_type == -1)
+            {
+                std::cout << " ENSEMBLE                 : " << "FIRE" << std::endl;
+            }
+            else if(INPUT.mdp.md_type == 0)
+            {
+                std::cout << " ENSEMBLE                 : " << "NVE" << std::endl;
+            }
+            else if(INPUT.mdp.md_type == 1 || INPUT.mdp.md_type == 3)
+            {
+                std::cout << " ENSEMBLE                 : " << "NVT" << std::endl;
+            }
+            else if(INPUT.mdp.md_type==2)
+            {
+                std::cout << " ENSEMBLE                 : " << "Langevin" << std::endl;
+            }
+            else if(INPUT.mdp.md_type==4)
+            {
+                std::cout << " ENSEMBLE                 : " << "MSST" << std::endl;
+            }
+
+            std::cout << " Time interval(fs)        : " << INPUT.mdp.md_dt << std::endl;
+        }
+        std::cout << " ---------------------------------------------------------" << std::endl;
 
 
 		std::cout << " " << std::setw(8) << "SPIN"
diff --git a/source/src_io/write_input.cpp b/source/src_io/write_input.cpp
@@ -175,7 +175,7 @@ void Input::Print(const std::string &fn)const
 	ModuleBase::GlobalFunc::OUTP(ofs,"msst_vis",mdp.msst_vis,"artificial viscosity");
 	ModuleBase::GlobalFunc::OUTP(ofs,"msst_tscale",mdp.msst_tscale,"reduction in initial temperature");
 	ModuleBase::GlobalFunc::OUTP(ofs,"msst_qmass",mdp.msst_qmass,"mass of thermostat");
-	ModuleBase::GlobalFunc::OUTP(ofs,"md_tfreq",mdp.md_tfreq,"oscillation frequency, used to determine msst_qmass of NHC");
+	ModuleBase::GlobalFunc::OUTP(ofs,"md_tfreq",mdp.md_tfreq,"oscillation frequency, used to determine qmass of NHC");
 	ModuleBase::GlobalFunc::OUTP(ofs,"md_damp",mdp.md_damp,"damping parameter (time units) used to add force in Langevin method");
 
 
diff --git a/source/src_lcao/run_md_lcao.cpp b/source/src_lcao/run_md_lcao.cpp
@@ -116,6 +116,8 @@ void Run_MD_LCAO::opt_ions(ModuleESolver::ESolver *p_esolver)
     // md cycle
     while ((verlet->step_ + verlet->step_rst_) <= GlobalV::MD_NSTEP && !verlet->stop)
     {
+        Print_Info::print_screen(0, 0, verlet->step_ + verlet->step_rst_);
+
         if(verlet->step_ == 0)
         {
             MD_func::ParaV = this->LM_md.ParaV;
@@ -159,8 +161,8 @@ void Run_MD_LCAO::opt_ions(ModuleESolver::ESolver *p_esolver)
 
         if((verlet->step_ + verlet->step_rst_) % verlet->mdp.md_dumpfreq == 0)
         {
-            Print_Info::print_screen(0, 0, verlet->step_ + verlet->step_rst_);
-            verlet->outputMD();
+            // Print_Info::print_screen(0, 0, verlet->step_ + verlet->step_rst_);
+            verlet->outputMD(GlobalV::ofs_running);
 
             MD_func::MDdump(verlet->step_ + verlet->step_rst_, verlet->ucell, verlet->virial, verlet->force);
         }
diff --git a/source/src_pw/electrons.cpp b/source/src_pw/electrons.cpp
diff --git a/source/src_pw/run_md_pw.cpp b/source/src_pw/run_md_pw.cpp

Original file line number	Diff line number	Diff line change
`@@ -80,11 +80,11 @@ void FIRE::second_half()`
`80`	`80`	`ModuleBase::timer::tick("FIRE", "second_half");`
`81`	`81`	`}`
`82`	`82`
`83`		`-void FIRE::outputMD()`
	`83`	`+void FIRE::outputMD(std::ofstream &ofs)`
`84`	`84`	`{`
`85`		`- Verlet::outputMD();`
	`85`	`+ Verlet::outputMD(ofs);`
`86`	`86`
`87`		`- GlobalV::ofs_running << " LARGEST GRAD (eV/A) : "`
	`87`	`+ ofs << " LARGEST GRAD (eV/A) : "`
`88`	`88`	`<< max * ModuleBase::Hartree_to_eV * ModuleBase::ANGSTROM_AU << std::endl;`
`89`	`89`	`std::cout << " LARGEST GRAD (eV/A) : "`
`90`	`90`	`<< max * ModuleBase::Hartree_to_eV * ModuleBase::ANGSTROM_AU << std::endl;`
Original file line number	Diff line number	Diff line change
`@@ -45,9 +45,9 @@ void Langevin::second_half()`
`45`	`45`	`ModuleBase::timer::tick("Langevin", "second_half");`
`46`	`46`	`}`
`47`	`47`
`48`		`-void Langevin::outputMD()`
	`48`	`+void Langevin::outputMD(std::ofstream &ofs)`
`49`	`49`	`{`
`50`		`- Verlet::outputMD();`
	`50`	`+ Verlet::outputMD(ofs);`
`51`	`51`	`}`
`52`	`52`
`53`	`53`	`void Langevin::write_restart()`
Original file line number	Diff line number	Diff line change
`@@ -160,9 +160,9 @@ void MSST::second_half()`
`160`	`160`	`ModuleBase::timer::tick("MSST", "second_half");`
`161`	`161`	`}`
`162`	`162`
`163`		`-void MSST::outputMD()`
	`163`	`+void MSST::outputMD(std::ofstream &ofs)`
`164`	`164`	`{`
`165`		`- Verlet::outputMD();`
	`165`	`+ Verlet::outputMD(ofs);`
`166`	`166`	`}`
`167`	`167`
`168`	`168`	`void MSST::write_restart()`
Original file line number	Diff line number	Diff line change
`@@ -36,9 +36,9 @@ void NVE::second_half()`
`36`	`36`	`ModuleBase::timer::tick("NVE", "second_half");`
`37`	`37`	`}`
`38`	`38`
`39`		`-void NVE::outputMD()`
	`39`	`+void NVE::outputMD(std::ofstream &ofs)`
`40`	`40`	`{`
`41`		`- Verlet::outputMD();`
	`41`	`+ Verlet::outputMD(ofs);`
`42`	`42`	`}`
`43`	`43`
`44`	`44`	`void NVE::write_restart()`