fix: update the default value of md_tfreq

Author: YuLiu98

doc/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
 

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

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)

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;
@@ -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

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;

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");