Merge pull request #1004 from YuLiu98/develop

refactor: update docs and check for MD parameters

Author: YuLiu98

docs/input-main.md

@@ -1114,7 +1114,7 @@ This part of variables are used to control the molecular dynamics calculations.
 
 - **Type**: Double
 - **Description**: This is the time step(fs) used in md simulation .
-- **Default**: 1
+- **Default**: 1.0
 
 #### md_tfirst & md_tlast
 
@@ -1178,31 +1178,31 @@ temperature will fluctuate violently; if it is too small, the temperature will t
 
 - **Type**: Real
 - **Description**: the velocity of shock wave ($\AA$/fs) for MSST.
-- **Default**: 0
+- **Default**: 0.0
 
 #### msst_vis
 
 - **Type**: Real
 - **Description**: artificial viscosity (mass/length/time) for MSST.
-- **Default**: 0
+- **Default**: 0.0
 
 #### msst_tscale
 
 - **Type**: Real
 - **Description**: reduction in initial temperature (0~1) used to compress volume in MSST.
-- **Default**: 0
+- **Default**: 0.01
 
 #### msst_qmass
 
-- **Type**: Double
+- **Type**: Real
 - **Description**: Inertia of extended system variable. Used only when md_type is 4, you should set a number which is larger than 0. Note that Qmass of NHC is set by md_tfreq.
 - **Default**: No default
 
 #### md_damp
 
 - **Type**: Real
 - **Description**: damping parameter (fs) used to add force in Langevin method.
-- **Default**: 1
+- **Default**: 1.0
 
 ### DFT+U correction
 

source/input.cpp

@@ -2358,6 +2358,26 @@ void Input::Check(void)
             ModuleBase::WARNING_QUIT("Input::Check", "temperature of MD calculation should be set!");
         if (mdp.md_tlast < 0.0)
             mdp.md_tlast = mdp.md_tfirst;
+
+        if(mdp.md_tfreq == 0)
+        {
+            mdp.md_tfreq = 1.0/40.0/mdp.md_dt;
+        }
+        if(mdp.md_restart) 
+        {
+            init_vel = 1;
+        }
+        if(mdp.md_ensolver == "LJ" || mdp.md_ensolver == "DP" || mdp.md_type == 4)
+        {
+            cal_stress = 1;
+        }
+        if(mdp.md_type == 4)
+        {
+            if(mdp.msst_qmass <= 0)
+            {
+                ModuleBase::WARNING_QUIT("Input::Check", "msst_qmass must be greater than 0!");
+            }
+        }
         // if(mdp.md_tfirst!=mdp.md_tlast)
         // {
         //     std::ifstream file1;

source/module_md/MD_func.cpp

@@ -242,7 +242,6 @@ void MD_func::force_virial(
 
 	if(mdp.md_ensolver == "LJ")
 	{
-        GlobalV::CAL_STRESS = 1;
 		bool which_method = unit_in.judge_big_cell();
 		if(which_method)
 		{
@@ -276,7 +275,6 @@ void MD_func::force_virial(
 	}
 	else if(mdp.md_ensolver == "DP")
 	{
-        GlobalV::CAL_STRESS = 1;
 		DP_potential::DP_pot(unit_in, potential, force, stress);
 	}
 #ifndef __CMD

source/module_md/MD_parameters.h

@@ -11,11 +11,11 @@ class MD_parameters
 		md_restart = 0;
 		md_type = 1;
         md_nstep = 10;
-		md_dt = 1;
-		md_tfirst = 0;
-		md_tlast = -1;
+		md_dt = 1.0;
+		md_tfirst = -1.0;
+		md_tlast = -1.0;
 		md_dumpfreq = 1;
-		md_restartfreq = 1;
+		md_restartfreq = 5;
 
 		// Classic MD 
 		md_ensolver = "FP";
@@ -25,17 +25,17 @@ class MD_parameters
 
 		// MSST
 		msst_direction = 2;
-		msst_qmass = 1;
-		msst_vel = 0;
-		msst_vis = 0;
+		msst_qmass = -1.0;
+		msst_vel = 0.0;
+		msst_vis = 0.0;
 		msst_tscale = 0.01;
 
 		// NHC
-		md_tfreq = 0;
+		md_tfreq = 0.0;
 		md_mnhc = 4;
 
 		// Langevin
-		md_damp = 1;
+		md_damp = 1.0;
 	};
     ~MD_parameters(){};
 

source/module_md/MSST.cpp

@@ -10,7 +10,6 @@ MSST::MSST(MD_parameters& MD_para_in, UnitCell_pseudo &unit_in) : Verlet(MD_para
 {
     std::cout << "MSST" << std::endl;
 
-    GlobalV::CAL_STRESS = 1;
     mdp.msst_qmass = mdp.msst_qmass / pow(ModuleBase::ANGSTROM_AU, 4) / pow(ModuleBase::AU_to_MASS, 2);
     mdp.msst_vel = mdp.msst_vel * ModuleBase::ANGSTROM_AU * ModuleBase::AU_to_FS;
     mdp.msst_vis = mdp.msst_vis / ModuleBase::AU_to_MASS / ModuleBase::ANGSTROM_AU * ModuleBase::AU_to_FS;

source/module_md/verlet.cpp

@@ -25,15 +25,7 @@ Verlet::Verlet(MD_parameters& MD_para_in, UnitCell_pseudo &unit_in):
     mdp.md_dt /= ModuleBase::AU_to_FS;
     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;
-    }
+    mdp.md_tfreq *= ModuleBase::AU_to_FS;
 
     // LJ parameters
     mdp.lj_rcut *= ModuleBase::ANGSTROM_AU;
@@ -42,7 +34,6 @@ Verlet::Verlet(MD_parameters& MD_para_in, UnitCell_pseudo &unit_in):
 
     step_ = 0;
     step_rst_ = 0;
-    if(mdp.md_restart) unit_in.init_vel = 1;
 
     MD_func::InitPos(ucell, pos);
     MD_func::InitVel(ucell, mdp.md_tfirst, allmass, frozen_freedom_, ionmbl, vel);

tests/integrate/120_PW_KP_MD_NVE/INPUT

@@ -24,5 +24,6 @@ mixing_beta     0.7
 
 md_type         0
 md_restart         0
+md_tfirst       0
 
 read_file_dir   ./

tests/integrate/220_NO_KP_MD_NVE/INPUT

@@ -24,5 +24,6 @@ mixing_beta     0.7
 
 md_type          0
 md_restart         0
+md_tfirst        0
 
 read_file_dir   ./

tests/integrate/320_NO_GO_MD_NVE/INPUT

@@ -25,6 +25,7 @@ mixing_beta     0.7
 
 md_type          0
 md_restart         0
+md_tfirst        0
 
 read_file_dir   ./
 

tests/integrate/601_NO_TDDFT_N2_occ/INPUT

@@ -29,6 +29,7 @@ scf_thr			1.0e-6
 tddft			1
 
 md_type		    0
+md_tfirst       0
 md_dt			0.01
 ocp                     1
 ocp_set                 5*2 2*0 1*2 4*0