Feature: enable three different precision levels for variable-cell MD

docs/advanced/input_files/input-main.md

@@ -32,7 +32,7 @@
   [exx_hybrid_alpha](#exx_hybrid_alpha) | [exx_hse_omega](#exx_hse_omega) | [exx_separate_loop](#exx_separate_loop) | [exx_hybrid_step](#exx_hybrid_step) | [exx_lambda](#exx_lambda) | [exx_pca_threshold](#exx_pca_threshold) | [exx_c_threshold](#exx_c_threshold) | [exx_v_threshold](#exx_v_threshold) | [exx_c_grad_threshold](#exx_c_grad_threshold) | [exx_v_grad_threshold](#exx_v_grad_threshold) | [exx_dm_threshold](#exx_dm_threshold) | [exx_schwarz_threshold](#exx_schwarz_threshold) | [exx_cauchy_threshold](#exx_cauchy_threshold) | [exx_cauchy_grad_threshold](#exx_cauchy_grad_threshold) | [exx_ccp_threshold](#exx_ccp_threshold) | [exx_ccp_rmesh_times](#exx_ccp_rmesh_times) | [exx_distribute_type](#exx_distribute_type) | [exx_opt_orb_lmax](#exx_opt_orb_lmax) | [exx_opt_orb_ecut](#exx_opt_orb_ecut) | [exx_opt_orb_tolerence](#exx_opt_orb_tolerence) | [exx_real_number](#exx_real_number)
 - [Molecular dynamics](#molecular-dynamics)
 
-  [md_type](#md_type) | [md_thermostat](#md_thermostat) | [md_nstep](#md_nstep) | [md_restart](#md_restart) | [md_dt](#md_dt) | [md_tfirst, md_tlast](#md_tfirst-md_tlast) | [md_dumpfreq](#md_dumpfreq) | [md_restartfreq](#md_restartfreq) | [md_seed](#md_seed) | [md_tfreq](#md_tfreq) | [md_tchain](#md_tchain) | [md_pmode](#md_pmode) | [md_pcouple](#md_pcouple) | [md_pfirst, md_plast](#md_pfirst-md_plast) | [md_pfreq](#md_pfreq) | [md_pchain](#md_pchain) | [dump_force](#dump_force) | [dump_vel](#dump_vel) | [dump_virial](#dump_virial) | [lj_rcut](#lj_rcut) | [lj_epsilon](#lj_epsilon) | [lj_sigma](#lj_sigma) | [pot_file](#pot_file) | [msst_direction](#msst_direction) | [msst_vel](#msst_vel) | [msst_vis](#msst_vis) | [msst_tscale](#msst_tscale) | [msst_qmass](#msst_qmass) | [md_damp](#md_damp) | [md_tolerance](#md_tolerance) | [md_nraise](#md_nraise)
+  [md_type](#md_type) | [md_thermostat](#md_thermostat) | [md_nstep](#md_nstep) | [md_restart](#md_restart) | [md_dt](#md_dt) | [md_tfirst, md_tlast](#md_tfirst-md_tlast) | [md_dumpfreq](#md_dumpfreq) | [md_restartfreq](#md_restartfreq) | [md_seed](#md_seed) | [md_prec_level](#md_prec_level) | [ref_cell_factor](#ref_cell_factor) | [md_tfreq](#md_tfreq) | [md_tchain](#md_tchain) | [md_pmode](#md_pmode) | [md_pcouple](#md_pcouple) | [md_pfirst, md_plast](#md_pfirst-md_plast) | [md_pfreq](#md_pfreq) | [md_pchain](#md_pchain) | [dump_force](#dump_force) | [dump_vel](#dump_vel) | [dump_virial](#dump_virial) | [lj_rcut](#lj_rcut) | [lj_epsilon](#lj_epsilon) | [lj_sigma](#lj_sigma) | [pot_file](#pot_file) | [msst_direction](#msst_direction) | [msst_vel](#msst_vel) | [msst_vis](#msst_vis) | [msst_tscale](#msst_tscale) | [msst_qmass](#msst_qmass) | [md_damp](#md_damp) | [md_tolerance](#md_tolerance) | [md_nraise](#md_nraise)
 - [vdW correction](#vdw-correction)
 
   [vdw_method](#vdw_method) | [vdw_s6](#vdw_s6) | [vdw_s8](#vdw_s8) | [vdw_a1](#vdw_a1) | [vdw_a2](#vdw_a2) | [vdw_d](#vdw_d) | [vdw_abc](#vdw_abc) | [vdw_C6_file](#vdw_c6_file) | [vdw_C6_unit](#vdw_c6_unit) | [vdw_R0_file](#vdw_r0_file) | [vdw_R0_unit](#vdw_r0_unit) | [vdw_cutoff_type](#vdw_cutoff_type) | [vdw_cutoff_radius](#vdw_cutoff_radius) | [vdw_radius_unit](#vdw_radius_unit) | [vdw_cutoff_period](#vdw_cutoff_period) | [vdw_cn_thr](#vdw_cn_thr) | [vdw_cn_thr_unit](#vdw_cn_thr_unit)
@@ -558,7 +558,7 @@ calculations.
   If you set gamma_only = 1, ABACUS uses gamma only, the algorithm is faster and you don't need to specify the k-points file. If you set gamma_only = 0, more than one k-point is used and the ABACUS is slower compared to the gamma only algorithm.
 
   > Note: If gamma_only is set to 1, the KPT file will be overwritten. So make sure to turn off gamma_only for multi-k calculations.
-  >
+
 - **Default**: 0
 
 ### printe
@@ -1564,7 +1564,7 @@ These variables are used to control the molecular dynamics calculations.
   - 2: NVT ensemble with Langevin method;
   - 4: MSST method;
 
-  ***Note: when md_type is set to 1, md_tfreq is required to stablize temperature. It is an empirical parameter whose value is system-dependent, ranging from 1/(40\*md_dt) to 1/(100\*md_dt). An improper choice of its value might lead to failure of job.***
+  > Note: when md_type is set to 1, md_tfreq is required to stablize temperature. It is an empirical parameter whose value is system-dependent, ranging from 1/(40\*md_dt) to 1/(100\*md_dt). An improper choice of its value might lead to failure of job.
 - **Default**: 1
 
 ### md_thermostat
@@ -1625,6 +1625,23 @@ These variables are used to control the molecular dynamics calculations.
   - md_seed >= 0: srand(md_seed) in MD initialization.
 - **Default**: -1
 
+### md_prec_level
+
+- **Type**: Integer
+- **Description**: Determine the precision level of vc-md. 
+  - 0: FFT grids do not change, only G vectors and K vectors are changed due to the change of lattice vector. This level is suitable for cases where the variation of the volume and shape is not large, and the efficiency is relatively higher.
+  - 1: A reference cell is constructed at the beginning, controlled by [ref_cell_factor](#ref_cell_factor). Then the reference cell is used to initialize FFT real-space grids and reciprocal space mesh instead of the initial cell. The cost will increase with the size of the reference cell.
+  - 2: FFT grids change per MD step. This level is suitable for cases where the variation of the volume and shape is large, such as the MSST method. However, accuracy comes at the cost of efficiency.
+  > Note: this parameter is only used in variable-cell MD!
+- **Default**: 0
+
+### ref_cell_factor
+
+- **Type**: Real
+- **Description**: 
+  Construct a reference cell bigger than the initial cell. Only used in variable-cell MD, if [md_prec_level](#md_prec_level) is set to 1. The reference cell has to be large enough so that the lattice vectors of the fluctuating cell do not exceed the reference lattice vectors during MD. Typically, 1.02 ~ 1.10 is sufficient. However, the cell fluctuations depend on the specific system and thermodynamic conditions. So users must test for a proper choice. This parameters should be used in conjunction with q2sigma, qcutz, and ecfixed. 
+- **Default**: 1.0
+
 ### md_tfreq
 
 - **Type**: Real

source/module_base/global_function.h

@@ -78,7 +78,7 @@ void OUTP(std::ofstream &ofs, const std::string &name, const T &a, const std::st
 template <class T>
 void OUT(const std::string &name,const T &a)
 {
-    std::cout << " " << std::setw(40) << name << " = " << a;
+    std::cout << " " << std::setw(40) << name << " = " << a << std::endl;
 //	std::cout << " " << name << a << std::endl;
     return;
 }

source/module_base/global_variable.cpp

@@ -58,6 +58,7 @@ int OUT_FREQ_ELEC = 0;
 int OUT_FREQ_ION = 0;
 int RELAX_NMAX = 20;
 int MD_NSTEP = 20;
+int md_prec_level = 0;
 int SCF_NMAX = 100;
 
 std::string BASIS_TYPE = "pw"; // xiaohui add 2013-09-01

source/module_base/global_variable.h

@@ -63,6 +63,7 @@ extern bool fixed_atoms;
 extern int RELAX_NMAX; // 8.3
 extern int SCF_NMAX; // 8.4
 extern int MD_NSTEP;
+extern int md_prec_level;    // liuyu 2023-03-13
 
 extern std::string BASIS_TYPE; // xiaohui add 2013-09-01
 extern std::string KS_SOLVER; // xiaohui add 2013-09-01

source/module_cell/atom_spec.cpp

@@ -13,7 +13,7 @@ Atom::Atom()
     stapos_wf = 0;
     mass = 0.0;
     tau = new ModuleBase::Vector3<double>[1];
-    tau_original = new ModuleBase::Vector3<double>[1];
+    dis = new ModuleBase::Vector3<double>[1];
     taud = new ModuleBase::Vector3<double>[1];
     vel = new ModuleBase::Vector3<double>[1];
     mag = new double[1];
@@ -32,7 +32,7 @@ Atom::Atom()
 Atom::~Atom()
 {
     delete[] tau;
-    delete[] tau_original;
+    delete[] dis;
     delete[] taud;
     delete[] vel;
     delete[] mag;
@@ -142,7 +142,7 @@ void Atom::bcast_atom(void)
     {
         assert(na!=0);
         delete[] tau;
-        delete[] tau_original;
+        delete[] dis;
 		delete[] taud;
 	    delete[] vel;
         delete[] mag;
@@ -151,7 +151,7 @@ void Atom::bcast_atom(void)
         delete[] m_loc_;
         delete[] mbl;
         tau = new ModuleBase::Vector3<double>[na];
-        tau_original = new ModuleBase::Vector3<double>[na];
+        dis = new ModuleBase::Vector3<double>[na];
 		taud = new ModuleBase::Vector3<double>[na];
 	    vel = new ModuleBase::Vector3<double>[na];
         mag = new double[na];
@@ -169,6 +169,9 @@ void Atom::bcast_atom(void)
         Parallel_Common::bcast_double( taud[i].x );
         Parallel_Common::bcast_double( taud[i].y );
         Parallel_Common::bcast_double( taud[i].z );
+        Parallel_Common::bcast_double( dis[i].x );
+        Parallel_Common::bcast_double( dis[i].y );
+        Parallel_Common::bcast_double( dis[i].z );
 	    Parallel_Common::bcast_double( vel[i].x );
 	    Parallel_Common::bcast_double( vel[i].y );
 	    Parallel_Common::bcast_double( vel[i].z );

source/module_cell/atom_spec.h

@@ -36,7 +36,7 @@ class Atom
 
     std::string label; // atomic symbol
     ModuleBase::Vector3<double> *tau;// Cartesian coordinates of each atom in this type.
-    ModuleBase::Vector3<double> *tau_original;// Cartesian coordinates of each atom in this type, but without periodic adjustment.
+    ModuleBase::Vector3<double> *dis;// direct displacements of each atom in this type in current step  liuyu modift 2023-03-22
 	ModuleBase::Vector3<double> *taud;// Direct coordinates of each atom in this type.
     ModuleBase::Vector3<double> *vel;// velocities of each atom in this type.
     ModuleBase::Vector3<double> *force; // force acting on each atom in this type.

source/module_cell/read_atoms.cpp

@@ -535,7 +535,7 @@ bool UnitCell::read_atom_positions(std::ifstream &ifpos, std::ofstream &ofs_runn
 			if (na > 0)
 			{
        			delete[] atoms[it].tau;
-				delete[] atoms[it].tau_original;
+				delete[] atoms[it].dis;
 				delete[] atoms[it].taud;
 				delete[] atoms[it].vel;
        			delete[] atoms[it].mbl;
@@ -544,7 +544,7 @@ bool UnitCell::read_atom_positions(std::ifstream &ifpos, std::ofstream &ofs_runn
                 delete[] atoms[it].angle2;
                 delete[] atoms[it].m_loc_;
        			atoms[it].tau = new ModuleBase::Vector3<double>[na];
-				atoms[it].tau_original = new ModuleBase::Vector3<double>[na];
+				atoms[it].dis = new ModuleBase::Vector3<double>[na];
        			atoms[it].taud = new ModuleBase::Vector3<double>[na];
 				atoms[it].vel = new ModuleBase::Vector3<double>[na];
        			atoms[it].mbl = new ModuleBase::Vector3<int>[na];
@@ -803,7 +803,7 @@ bool UnitCell::read_atom_positions(std::ifstream &ifpos, std::ofstream &ofs_runn
 						atoms[it].mbl[ia] = 0.0;
 						atoms[it].mbl[ia].print();
 					}
-					atoms[it].tau_original[ia] = atoms[it].tau[ia];
+                    atoms[it].dis[ia].set(0, 0, 0);
 				}//endj
 			}// end na
 			//reset some useless parameters

source/module_cell/test/prepare_unitcell.h

@@ -147,7 +147,7 @@ class UcellTestPrepare
 			ucell->atoms[it].na = this->natom[it];
 			//coordinates and related physical quantities
 			delete[] ucell->atoms[it].tau;
-			delete[] ucell->atoms[it].tau_original;
+			delete[] ucell->atoms[it].dis;
 			delete[] ucell->atoms[it].taud;
 			delete[] ucell->atoms[it].vel;
 			delete[] ucell->atoms[it].mag;
@@ -156,7 +156,7 @@ class UcellTestPrepare
 			delete[] ucell->atoms[it].m_loc_;
 			delete[] ucell->atoms[it].mbl;
 			ucell->atoms[it].tau = new ModuleBase::Vector3<double>[ucell->atoms[it].na];
-			ucell->atoms[it].tau_original = new ModuleBase::Vector3<double>[ucell->atoms[it].na];
+			ucell->atoms[it].dis = new ModuleBase::Vector3<double>[ucell->atoms[it].na];
 			ucell->atoms[it].taud = new ModuleBase::Vector3<double>[ucell->atoms[it].na];
 			ucell->atoms[it].vel = new ModuleBase::Vector3<double>[ucell->atoms[it].na];
 			ucell->atoms[it].mag = new double[ucell->atoms[it].na];
@@ -186,7 +186,7 @@ class UcellTestPrepare
 						ucell->latvec.e31, ucell->latvec.e32, ucell->latvec.e33,
 						ucell->atoms[it].taud[ia].x, ucell->atoms[it].taud[ia].y, ucell->atoms[it].taud[ia].z);
 				}
-				ucell->atoms[it].tau_original[ia] = ucell->atoms[it].tau[ia];
+                ucell->atoms[it].dis[ia].set(0, 0, 0);
 				if(this->init_vel)
 				{
 					ucell->atoms[it].vel[ia].x = this->velocity[this->atomic_index*3+0];

source/module_cell/test/unitcell_test.cpp

@@ -64,11 +64,6 @@ Magnetism::~Magnetism()
  *     - check_tau(): check if any "two atoms are too close"
  *   - SelectiveDynamics
  *     - if_atoms_can_move():it is true if any coordinates of any atom can move, i.e. mbl = 1
- *   - SaveCartesianPosition
- *     - save_cartesian_position(): make a copy of atomic Cartesian coordinates: tau
- *   - SaveCartesianPositionOriginal
- *     - save_cartesian_position_original(): make a copy of original atomic Cartesian coordinates: tau_original
- *       tau_original means without periodic adjustment
  *   - PeriodicBoundaryAdjustment
  *     - periodic_boundary_adjustment(): move atoms inside the unitcell after relaxation
  *   - PrintCell
@@ -84,10 +79,6 @@ Magnetism::~Magnetism()
  *     - Actually an integrated function to call UnitCell::print_cell and Atom::print_Atom
  *   - UpdateVel
  *     - update_vel(const ModuleBase::Vector3<double>* vel_in)
- *   - UpdatePosTau1
- *     - update_pos_tau(const double* pos)
- *   - UpdatePosTau2
- *     - update_pos_tau(const ModuleBase::Vector3<double>* posd_in)
  *   - CalUx
  *     - cal_ux(): calculate magnetic moments of cell
  *   - ReadOrbFile
@@ -602,68 +593,6 @@ TEST_F(UcellTest,SelectiveDynamics)
 	EXPECT_TRUE(ucell->if_atoms_can_move());
 }
 
-TEST_F(UcellTest,SaveCartesianPosition)
-{
-	UcellTestPrepare utp = UcellTestLib["C1H2-Index"];
-	GlobalV::relax_new = utp.relax_new;
-	ucell = utp.SetUcellInfo();
-	//the first realization
-	double* pos = new double[3*ucell->nat];
-	ucell->save_cartesian_position(pos);
-	//another realization
-	ModuleBase::Vector3<double>* pos1 = new ModuleBase::Vector3<double>[ucell->nat];
-	ucell->save_cartesian_position(pos1);
-	int iat = 0;
-	for(int it=0; it<utp.natom.size(); ++it)
-	{
-		for(int ia=0; ia<utp.natom[it];++ia)
-		{
-			//first
-			EXPECT_DOUBLE_EQ(pos[3*iat],ucell->atoms[it].tau[ia].x*ucell->lat0);
-			EXPECT_DOUBLE_EQ(pos[3*iat+1],ucell->atoms[it].tau[ia].y*ucell->lat0);
-			EXPECT_DOUBLE_EQ(pos[3*iat+2],ucell->atoms[it].tau[ia].z*ucell->lat0);
-			//second
-			EXPECT_DOUBLE_EQ(pos1[iat].x,ucell->atoms[it].tau[ia].x*ucell->lat0);
-			EXPECT_DOUBLE_EQ(pos1[iat].y,ucell->atoms[it].tau[ia].y*ucell->lat0);
-			EXPECT_DOUBLE_EQ(pos1[iat].z,ucell->atoms[it].tau[ia].z*ucell->lat0);
-			++iat;
-		}
-	}
-	delete[] pos;
-	delete[] pos1;
-}
-
-TEST_F(UcellTest,SaveCartesianPositionOriginal)
-{
-	UcellTestPrepare utp = UcellTestLib["C1H2-Index"];
-	GlobalV::relax_new = utp.relax_new;
-	ucell = utp.SetUcellInfo();
-	//the first realization
-	double* pos = new double[3*ucell->nat];
-	ucell->save_cartesian_position_original(pos);
-	//another realization
-	ModuleBase::Vector3<double>* pos1 = new ModuleBase::Vector3<double>[ucell->nat];
-	ucell->save_cartesian_position_original(pos1);
-	int iat = 0;
-	for(int it=0; it<utp.natom.size(); ++it)
-	{
-		for(int ia=0; ia<utp.natom[it];++ia)
-		{
-			//first
-			EXPECT_DOUBLE_EQ(pos[3*iat],ucell->atoms[it].tau_original[ia].x*ucell->lat0);
-			EXPECT_DOUBLE_EQ(pos[3*iat+1],ucell->atoms[it].tau_original[ia].y*ucell->lat0);
-			EXPECT_DOUBLE_EQ(pos[3*iat+2],ucell->atoms[it].tau[ia].z*ucell->lat0);
-			//second
-			EXPECT_DOUBLE_EQ(pos1[iat].x,ucell->atoms[it].tau_original[ia].x*ucell->lat0);
-			EXPECT_DOUBLE_EQ(pos1[iat].y,ucell->atoms[it].tau_original[ia].y*ucell->lat0);
-			EXPECT_DOUBLE_EQ(pos1[iat].z,ucell->atoms[it].tau_original[ia].z*ucell->lat0);
-			++iat;
-		}
-	}
-	delete[] pos;
-	delete[] pos1;
-}
-
 TEST_F(UcellDeathTest,PeriodicBoundaryAdjustment1)
 {
 	UcellTestPrepare utp = UcellTestLib["C1H2-PBA"];
@@ -842,54 +771,6 @@ TEST_F(UcellTest,UpdateVel)
 	delete[] vel_in;
 }
 
-TEST_F(UcellTest,UpdatePosTau1)
-{
-	UcellTestPrepare utp = UcellTestLib["C1H2-Index"];
-	GlobalV::relax_new = utp.relax_new;
-	ucell = utp.SetUcellInfo();
-	double* pos_in = new double[ucell->nat*3];
-	ucell->set_iat2itia();
-	for(int iat=0; iat<ucell->nat; ++iat)
-	{
-		pos_in[iat*3] = 0.01*ucell->lat0;
-		pos_in[iat*3+1] = 0.01*ucell->lat0;
-		pos_in[iat*3+2] = 0.01*ucell->lat0;
-	}
-	ucell->update_pos_tau(pos_in);
-	for(int iat=0; iat<ucell->nat; ++iat)
-	{
-		int it, ia;
-		ucell->iat2iait(iat,&ia,&it);
-		if(ucell->atoms[it].mbl[ia].x != 0) EXPECT_DOUBLE_EQ(ucell->atoms[it].tau[ia].x,0.01);
-		if(ucell->atoms[it].mbl[ia].y != 0) EXPECT_DOUBLE_EQ(ucell->atoms[it].tau[ia].y,0.01);
-		if(ucell->atoms[it].mbl[ia].z != 0) EXPECT_DOUBLE_EQ(ucell->atoms[it].tau[ia].z,0.01);
-	}
-	delete[] pos_in;
-}
-
-TEST_F(UcellTest,UpdatePosTau2)
-{
-	UcellTestPrepare utp = UcellTestLib["C1H2-Index"];
-	GlobalV::relax_new = utp.relax_new;
-	ucell = utp.SetUcellInfo();
-	ModuleBase::Vector3<double>* pos_in = new ModuleBase::Vector3<double>[ucell->nat];
-	ucell->set_iat2itia();
-	for(int iat=0; iat<ucell->nat; ++iat)
-	{
-		pos_in[iat].set(0.01*ucell->lat0,0.01*ucell->lat0,0.01*ucell->lat0);
-	}
-	ucell->update_pos_tau(pos_in);
-	for(int iat=0; iat<ucell->nat; ++iat)
-	{
-		int it, ia;
-		ucell->iat2iait(iat,&ia,&it);
-		if(ucell->atoms[it].mbl[ia].x != 0) EXPECT_DOUBLE_EQ(ucell->atoms[it].tau[ia].x,0.01);
-		if(ucell->atoms[it].mbl[ia].y != 0) EXPECT_DOUBLE_EQ(ucell->atoms[it].tau[ia].y,0.01);
-		if(ucell->atoms[it].mbl[ia].z != 0) EXPECT_DOUBLE_EQ(ucell->atoms[it].tau[ia].z,0.01);
-	}
-	delete[] pos_in;
-}
-
 TEST_F(UcellTest,CalUx1)
 {
 	UcellTestPrepare utp = UcellTestLib["C1H2-Read"];