From 5f186a2cba6adb6b16be7cfca3a0db28c1451be8 Mon Sep 17 00:00:00 2001 From: YuLiu98 Date: Fri, 10 Mar 2023 11:04:51 +0800 Subject: [PATCH 01/17] Refactor: enable call Parallel_Grid::init() again --- .../hamilt_pwdft/parallel_grid.cpp | 46 ++++++++----------- .../hamilt_pwdft/parallel_grid.h | 10 ++-- 2 files changed, 22 insertions(+), 34 deletions(-) diff --git a/source/module_hamilt_pw/hamilt_pwdft/parallel_grid.cpp b/source/module_hamilt_pw/hamilt_pwdft/parallel_grid.cpp index c1e30daf18..df0735dd7c 100644 --- a/source/module_hamilt_pw/hamilt_pwdft/parallel_grid.cpp +++ b/source/module_hamilt_pw/hamilt_pwdft/parallel_grid.cpp @@ -3,13 +3,12 @@ Parallel_Grid::Parallel_Grid() { - this->allocate = false; + this->allocate = false; this->allocate_final_scf = false; //LiuXh add 20180619 } Parallel_Grid::~Parallel_Grid() { - //if(this->allocate) //LiuXh modify 20180619 if(this->allocate || this->allocate_final_scf) //LiuXh add 20180619 { for(int ip=0; ipncxy = ncx * ncy; this->ncxyz = ncxy * ncz; + // enable to call this function again liuyu 2023-03-10 + if(this->allocate) + { + for(int ip=0; ipallocate = false; + } + // (2) assert(allocate==false); assert(GlobalV::KPAR > 0); @@ -81,8 +97,6 @@ void Parallel_Grid::init( this->numz = new int*[GlobalV::KPAR]; this->startz = new int*[GlobalV::KPAR]; this->whichpro = new int*[GlobalV::KPAR]; - this->numdata = new int*[GlobalV::KPAR]; - this->startdata = new int*[GlobalV::KPAR]; for(int ip=0; ipnumz[ip] = new int[nproc]; this->startz[ip] = new int[nproc]; this->whichpro[ip] = new int[this->ncz]; - this->numdata[ip] = new int[nproc]; - this->startdata[ip] = new int[nproc]; ModuleBase::GlobalFunc::ZEROS(this->numz[ip], nproc); ModuleBase::GlobalFunc::ZEROS(this->startz[ip], nproc); ModuleBase::GlobalFunc::ZEROS(this->whichpro[ip], this->ncz); - ModuleBase::GlobalFunc::ZEROS(this->numdata[ip], nproc); - ModuleBase::GlobalFunc::ZEROS(this->startdata[ip], nproc); } this->allocate = true; @@ -164,20 +174,6 @@ void Parallel_Grid::z_distribution(void) // { // GlobalV::ofs_running << "\n iz=" << iz << " whichpro=" << whichpro[ip][iz]; // } - - //(4) - for(int proc=0; procnumz = new int*[GlobalV::KPAR]; this->startz = new int*[GlobalV::KPAR]; this->whichpro = new int*[GlobalV::KPAR]; - this->numdata = new int*[GlobalV::KPAR]; - this->startdata = new int*[GlobalV::KPAR]; for(int ip=0; ipnumz[ip] = new int[nproc]; this->startz[ip] = new int[nproc]; this->whichpro[ip] = new int[this->ncz]; - this->numdata[ip] = new int[nproc]; - this->startdata[ip] = new int[nproc]; ModuleBase::GlobalFunc::ZEROS(this->numz[ip], nproc); ModuleBase::GlobalFunc::ZEROS(this->startz[ip], nproc); ModuleBase::GlobalFunc::ZEROS(this->whichpro[ip], this->ncz); - ModuleBase::GlobalFunc::ZEROS(this->numdata[ip], nproc); - ModuleBase::GlobalFunc::ZEROS(this->startdata[ip], nproc); } this->allocate_final_scf = true; diff --git a/source/module_hamilt_pw/hamilt_pwdft/parallel_grid.h b/source/module_hamilt_pw/hamilt_pwdft/parallel_grid.h index 4cb7f79c79..222c304baa 100644 --- a/source/module_hamilt_pw/hamilt_pwdft/parallel_grid.h +++ b/source/module_hamilt_pw/hamilt_pwdft/parallel_grid.h @@ -31,12 +31,10 @@ class Parallel_Grid void z_distribution(void); - int *nproc_in_pool; - int **numz; - int **startz; - int **whichpro; - int **numdata; - int **startdata; + int *nproc_in_pool = nullptr; + int **numz = nullptr; + int **startz = nullptr; + int **whichpro = nullptr; int ncx; int ncy; From 665f0f3bf722a146dd46261a1661e76978088f3a Mon Sep 17 00:00:00 2001 From: YuLiu98 Date: Sun, 12 Mar 2023 10:12:09 +0800 Subject: [PATCH 02/17] Feature: add init_after_vc in esolver --- source/module_esolver/esolver_fp.cpp | 14 ++++++++ source/module_esolver/esolver_fp.h | 1 + source/module_esolver/esolver_ks.cpp | 25 +++++++++++++ source/module_esolver/esolver_ks.h | 2 ++ source/module_esolver/esolver_ks_pw.cpp | 35 +++++++++++++++++++ source/module_esolver/esolver_ks_pw.h | 1 + .../hamilt_pwdft/wavefunc.cpp | 11 +++--- 7 files changed, 85 insertions(+), 4 deletions(-) diff --git a/source/module_esolver/esolver_fp.cpp b/source/module_esolver/esolver_fp.cpp index ed8dd44cbe..05ac583c06 100644 --- a/source/module_esolver/esolver_fp.cpp +++ b/source/module_esolver/esolver_fp.cpp @@ -39,6 +39,20 @@ namespace ModuleESolver } + void ESolver_FP::init_after_vc(Input& inp, UnitCell& cell) + { + if (inp.nx * inp.ny * inp.nz == 0) + this->pw_rho->initgrids(cell.lat0, cell.latvec, inp.ecutrho); + else + this->pw_rho->initgrids(cell.lat0, cell.latvec, inp.nx, inp.ny, inp.nz); + + this->pw_rho->initparameters(false, inp.ecutrho); + this->pw_rho->setuptransform(); + this->pw_rho->collect_local_pw(); + this->pw_rho->collect_uniqgg(); + this->print_rhofft(inp, GlobalV::ofs_running); + } + void ESolver_FP::print_rhofft(Input&inp, ofstream &ofs) { std::cout << " UNIFORM GRID DIM : " << GlobalC::rhopw->nx << " * " << GlobalC::rhopw->ny << " * " << GlobalC::rhopw->nz << std::endl; diff --git a/source/module_esolver/esolver_fp.h b/source/module_esolver/esolver_fp.h index 208bde8a8b..8c90ca0708 100644 --- a/source/module_esolver/esolver_fp.h +++ b/source/module_esolver/esolver_fp.h @@ -14,6 +14,7 @@ namespace ModuleESolver ESolver_FP(); virtual ~ESolver_FP(); virtual void Init(Input& inp, UnitCell& cell) override; + virtual void init_after_vc(Input& inp, UnitCell& cell); // liuyu add 2023-03-09 // Hamilt* phamilt; elecstate::ElecState* pelec = nullptr; diff --git a/source/module_esolver/esolver_ks.cpp b/source/module_esolver/esolver_ks.cpp index b1626c0ff3..65c0c901cd 100644 --- a/source/module_esolver/esolver_ks.cpp +++ b/source/module_esolver/esolver_ks.cpp @@ -102,6 +102,31 @@ namespace ModuleESolver CE.Init_CE(); } + template + void ESolver_KS::init_after_vc(Input& inp, UnitCell& ucell) + { + ESolver_FP::init_after_vc(inp,ucell); + + // symm_flag == 0 in md calculation, thus this part is annotated + // if (ModuleSymmetry::Symmetry::symm_flag == 1) + // { + // GlobalC::symm.analy_sys(ucell, GlobalV::ofs_running); + // ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, "SYMMETRY"); + // } + + // reset cartesian coordinates due to the change of lattice + GlobalC::kv.set_after_vc(GlobalC::symm, GlobalV::global_kpoint_card, GlobalV::NSPIN, GlobalC::ucell.G, GlobalC::ucell.latvec); + ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, "INIT K-POINTS"); + + // initialize the real-space uniform grid for FFT and parallel + // distribution of plane waves + GlobalC::Pgrid.init(GlobalC::rhopw->nx, GlobalC::rhopw->ny, GlobalC::rhopw->nz, GlobalC::rhopw->nplane, + GlobalC::rhopw->nrxx, GlobalC::bigpw->nbz, GlobalC::bigpw->bz); // mohan add 2010-07-22, update 2011-05-04 + + // Calculate Structure factor + GlobalC::sf.setup_structure_factor(&GlobalC::ucell, GlobalC::rhopw); + } + template void ESolver_KS::hamilt2density(const int istep, const int iter, const FPTYPE ethr) { diff --git a/source/module_esolver/esolver_ks.h b/source/module_esolver/esolver_ks.h index 9dc378a122..3345f60d60 100644 --- a/source/module_esolver/esolver_ks.h +++ b/source/module_esolver/esolver_ks.h @@ -28,6 +28,8 @@ namespace ModuleESolver int out_freq_elec;// frequency for output virtual void Init(Input& inp, UnitCell& cell) override; + virtual void init_after_vc(Input& inp, UnitCell& cell) override; // liuyu add 2023-03-09 + virtual void Run(const int istep, UnitCell& cell) override; // calculate electron density from a specific Hamiltonian diff --git a/source/module_esolver/esolver_ks_pw.cpp b/source/module_esolver/esolver_ks_pw.cpp index 6f28773a8c..76224c9ba5 100644 --- a/source/module_esolver/esolver_ks_pw.cpp +++ b/source/module_esolver/esolver_ks_pw.cpp @@ -87,6 +87,7 @@ namespace ModuleESolver template void ESolver_KS_PW::Init_GlobalC(Input& inp, UnitCell& cell) { + if(this->psi != nullptr) delete this->psi; this->psi = GlobalC::wf.allocate(GlobalC::kv.nks); // cout<nrxx< + void ESolver_KS_PW::init_after_vc(Input& inp, UnitCell& ucell) + { + ESolver_KS::init_after_vc(inp,ucell); + + this->pw_wfc->initgrids(ucell.lat0, ucell.latvec, GlobalC::rhopw->nx, GlobalC::rhopw->ny, GlobalC::rhopw->nz); + this->pw_wfc->initparameters(false, inp.ecutwfc, GlobalC::kv.nks, GlobalC::kv.kvec_d.data()); +#ifdef __MPI + if(INPUT.pw_seed > 0) MPI_Allreduce(MPI_IN_PLACE, &this->pw_wfc->ggecut, 1, MPI_DOUBLE, MPI_MAX , MPI_COMM_WORLD); + //qianrui add 2021-8-13 to make different kpar parameters can get the same results +#endif + this->pw_wfc->setuptransform(); + for(int ik = 0 ; ik < GlobalC::kv.nks; ++ik) GlobalC::kv.ngk[ik] = this->pw_wfc->npwk[ik]; + this->pw_wfc->collect_local_pw(); + + delete this->phsol; this->phsol = new hsolver::HSolverPW(GlobalC::wfcpw); + + delete this->pelec; this->pelec = new elecstate::ElecStatePW( GlobalC::wfcpw, &(this->chr), (K_Vectors*)(&(GlobalC::kv))); + + this->pelec->charge->allocate(GlobalV::NSPIN, GlobalC::rhopw->nrxx, GlobalC::rhopw->npw); + + delete this->pelec->pot; + this->pelec->pot = new elecstate::Potential( + GlobalC::rhopw, + &GlobalC::ucell, + &(GlobalC::ppcell.vloc), + &(GlobalC::sf.strucFac), + &(GlobalC::en.etxc), + &(GlobalC::en.vtxc)); + + //temporary + this->Init_GlobalC(inp,ucell); + } + template void ESolver_KS_PW::beforescf(int istep) { diff --git a/source/module_esolver/esolver_ks_pw.h b/source/module_esolver/esolver_ks_pw.h index f7212d5038..babda11fe5 100644 --- a/source/module_esolver/esolver_ks_pw.h +++ b/source/module_esolver/esolver_ks_pw.h @@ -18,6 +18,7 @@ namespace ModuleESolver ESolver_KS_PW(); ~ESolver_KS_PW(); void Init(Input& inp, UnitCell& cell) override; + void init_after_vc(Input& inp, UnitCell& cell) override; void cal_Energy(double& etot) override; void cal_Force(ModuleBase::matrix& force) override; void cal_Stress(ModuleBase::matrix& stress) override; diff --git a/source/module_hamilt_pw/hamilt_pwdft/wavefunc.cpp b/source/module_hamilt_pw/hamilt_pwdft/wavefunc.cpp index ed30f1c1a3..63788e6a67 100644 --- a/source/module_hamilt_pw/hamilt_pwdft/wavefunc.cpp +++ b/source/module_hamilt_pw/hamilt_pwdft/wavefunc.cpp @@ -19,7 +19,7 @@ wavefunc::~wavefunc() { std::cout << " ~wavefunc()" << std::endl; } - if(this->irindex == nullptr) + if(this->irindex != nullptr) { delete[] this->irindex; this->irindex=nullptr; @@ -62,7 +62,8 @@ psi::Psi>* wavefunc::allocate(const int nks) } else if(GlobalV::BASIS_TYPE!="pw") { - this->evc = new ModuleBase::ComplexMatrix [nks2]; + if(this->evc != nullptr) delete[] this->evc; + this->evc = new ModuleBase::ComplexMatrix [nks2]; for (int ik = 0; ik < nks2; ik++) { this->evc[ik].create(GlobalV::NBANDS, npwx * GlobalV::NPOL); @@ -73,7 +74,8 @@ psi::Psi>* wavefunc::allocate(const int nks) if((GlobalV::BASIS_TYPE=="lcao" || GlobalV::BASIS_TYPE=="lcao_in_pw") || winput::out_spillage==2) {//for lcao_in_pw - this->wanf2 = new ModuleBase::ComplexMatrix [nks2]; + if(this->wanf2 != nullptr) delete[] this->wanf2; + this->wanf2 = new ModuleBase::ComplexMatrix [nks2]; for (int ik = 0; ik < nks2; ik++) { this->wanf2[ik].create(GlobalV::NLOCAL, npwx * GlobalV::NPOL); @@ -520,7 +522,8 @@ void wavefunc::wfcinit_k(psi::Psi>* psi_in) if(GlobalV::BASIS_TYPE=="pw") { - this->irindex = new int [GlobalC::wfcpw->fftnxy]; + if(this->irindex != nullptr) delete[] this->irindex; + this->irindex = new int [GlobalC::wfcpw->fftnxy]; GlobalC::wfcpw->getfftixy2is(this->irindex); #if defined(__CUDA) || defined(__ROCM) if (GlobalV::device_flag == "gpu") { From 845f5c1c173fddb0e6b14752d6ed0c4e116e6259 Mon Sep 17 00:00:00 2001 From: YuLiu98 Date: Mon, 13 Mar 2023 10:08:53 +0800 Subject: [PATCH 03/17] Feature: add a para to control md prec level --- docs/advanced/input_files/input-main.md | 10 +++++++++- source/module_io/input.cpp | 5 +++++ source/module_io/test/INPUT | 1 + source/module_io/test/input_test.cpp | 2 ++ source/module_io/write_input.cpp | 1 + source/module_md/MD_parameters.h | 2 ++ 6 files changed, 20 insertions(+), 1 deletion(-) diff --git a/docs/advanced/input_files/input-main.md b/docs/advanced/input_files/input-main.md index 03ff396e5a..bbd29e9c01 100644 --- a/docs/advanced/input_files/input-main.md +++ b/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) | [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) @@ -1611,6 +1611,14 @@ 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. Note that this parameter is only used in variable-cell 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: 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. +- **Default**: 0 + ### md_tfreq - **Type**: Real diff --git a/source/module_io/input.cpp b/source/module_io/input.cpp index 8e3a3f3c5c..c9b85c698f 100644 --- a/source/module_io/input.cpp +++ b/source/module_io/input.cpp @@ -1353,6 +1353,10 @@ bool Input::Read(const std::string &fn) { read_value(ifs, mdp.md_seed); } + else if (strcmp("md_prec_level", word) == 0) + { + read_value(ifs, mdp.md_prec_level); + } else if (strcmp("md_restart", word) == 0) { read_value(ifs, mdp.md_restart); @@ -2852,6 +2856,7 @@ void Input::Bcast() Parallel_Common::bcast_int(mdp.md_dumpfreq); Parallel_Common::bcast_int(mdp.md_restartfreq); Parallel_Common::bcast_int(mdp.md_seed); + Parallel_Common::bcast_int(mdp.md_prec_level); Parallel_Common::bcast_bool(mdp.md_restart); Parallel_Common::bcast_double(mdp.lj_rcut); Parallel_Common::bcast_double(mdp.lj_epsilon); diff --git a/source/module_io/test/INPUT b/source/module_io/test/INPUT index 7024fba0ed..b4e8dc8922 100644 --- a/source/module_io/test/INPUT +++ b/source/module_io/test/INPUT @@ -177,6 +177,7 @@ md_tlast -1 #temperature last md_dumpfreq 1 #The period to dump MD information md_restartfreq 5 #The period to output MD restart information md_seed -1 #random seed for MD +md_prec_level 1 #precision level for vc-md md_restart 0 #whether restart lj_rcut 8.5 #cutoff radius of LJ potential lj_epsilon 0.01032 #the value of epsilon for LJ potential diff --git a/source/module_io/test/input_test.cpp b/source/module_io/test/input_test.cpp index 782c9b4025..2d5baccd05 100644 --- a/source/module_io/test/input_test.cpp +++ b/source/module_io/test/input_test.cpp @@ -333,6 +333,7 @@ TEST_F(InputTest, Default) EXPECT_EQ(INPUT.mdp.md_restart,0); EXPECT_EQ(INPUT.mdp.md_restartfreq,5); EXPECT_EQ(INPUT.mdp.md_seed,-1); + EXPECT_EQ(INPUT.mdp.md_prec_level,0); EXPECT_EQ(INPUT.mdp.md_tchain,1); EXPECT_DOUBLE_EQ(INPUT.mdp.md_tfirst,-1); EXPECT_DOUBLE_EQ(INPUT.mdp.md_tfreq,0); @@ -665,6 +666,7 @@ TEST_F(InputTest, Read) EXPECT_EQ(INPUT.mdp.md_restart,0); EXPECT_EQ(INPUT.mdp.md_restartfreq,5); EXPECT_EQ(INPUT.mdp.md_seed,-1); + EXPECT_EQ(INPUT.mdp.md_prec_level,1); EXPECT_EQ(INPUT.mdp.md_tchain,1); EXPECT_DOUBLE_EQ(INPUT.mdp.md_tfirst,-1); EXPECT_DOUBLE_EQ(INPUT.mdp.md_tfreq,0); diff --git a/source/module_io/write_input.cpp b/source/module_io/write_input.cpp index 35e45f6d75..2b5f881901 100644 --- a/source/module_io/write_input.cpp +++ b/source/module_io/write_input.cpp @@ -247,6 +247,7 @@ ModuleBase::GlobalFunc::OUTP(ofs, "out_bandgap", out_bandgap, "if true, print ou ModuleBase::GlobalFunc::OUTP(ofs,"md_dumpfreq",mdp.md_dumpfreq,"The period to dump MD information"); ModuleBase::GlobalFunc::OUTP(ofs,"md_restartfreq",mdp.md_restartfreq,"The period to output MD restart information"); ModuleBase::GlobalFunc::OUTP(ofs,"md_seed",mdp.md_seed,"random seed for MD"); + ModuleBase::GlobalFunc::OUTP(ofs,"md_prec_level",mdp.md_prec_level,"precision level for vc-md"); ModuleBase::GlobalFunc::OUTP(ofs,"md_restart",mdp.md_restart,"whether restart"); ModuleBase::GlobalFunc::OUTP(ofs,"lj_rcut",mdp.lj_rcut,"cutoff radius of LJ potential"); ModuleBase::GlobalFunc::OUTP(ofs,"lj_epsilon",mdp.lj_epsilon,"the value of epsilon for LJ potential"); diff --git a/source/module_md/MD_parameters.h b/source/module_md/MD_parameters.h index 7a5a3b16c9..b108d7589b 100644 --- a/source/module_md/MD_parameters.h +++ b/source/module_md/MD_parameters.h @@ -18,6 +18,7 @@ class MD_parameters md_dumpfreq = 1; md_restartfreq = 5; md_seed = -1; + md_prec_level = 0; // Classic MD lj_rcut = 8.5; @@ -61,6 +62,7 @@ class MD_parameters int md_dumpfreq; // The period to dump MD information int md_restartfreq; // The period to output MD restart information int md_seed; // random seed for MD + int md_prec_level; // precision level for vc-md; 0: do not reinit FFT grid, 1: reinit FFT grid per step // Classic MD // liuyu 2021-07-30 double lj_rcut; // cutoff radius of LJ potential (\AA) From fb06439f62efefce1e821ccb38b9c7b3006e52e8 Mon Sep 17 00:00:00 2001 From: YuLiu98 Date: Mon, 13 Mar 2023 11:44:48 +0800 Subject: [PATCH 04/17] Feature: enable pw basis vc-md --- source/module_base/global_variable.cpp | 1 + source/module_base/global_variable.h | 1 + source/module_elecstate/elecstate.cpp | 2 +- .../module_elecstate/module_charge/charge.cpp | 63 +++++++++++-------- .../module_elecstate/module_charge/charge.h | 22 +++---- source/module_esolver/esolver_ks_pw.cpp | 34 ++++++---- source/module_md/MSST.cpp | 3 +- source/module_md/Nose_Hoover.cpp | 6 ++ source/module_md/run_md.cpp | 1 - 9 files changed, 81 insertions(+), 52 deletions(-) diff --git a/source/module_base/global_variable.cpp b/source/module_base/global_variable.cpp index b41ef7120b..e8984c0916 100644 --- a/source/module_base/global_variable.cpp +++ b/source/module_base/global_variable.cpp @@ -57,6 +57,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 diff --git a/source/module_base/global_variable.h b/source/module_base/global_variable.h index 54ec3a720d..bb6cddcd30 100644 --- a/source/module_base/global_variable.h +++ b/source/module_base/global_variable.h @@ -62,6 +62,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 diff --git a/source/module_elecstate/elecstate.cpp b/source/module_elecstate/elecstate.cpp index f7b16f20d0..38ffbcd49e 100644 --- a/source/module_elecstate/elecstate.cpp +++ b/source/module_elecstate/elecstate.cpp @@ -378,7 +378,7 @@ void ElecState::init_scf(const int istep, const ModuleBase::ComplexMatrix& struc // (2) other effective potentials need charge density, // choose charge density from ionic step 0. //-------------------------------------------------------------------- - if (istep == 0) + if (istep == 0 || GlobalV::MD_PREC_LEVEL == 1) { this->charge->init_rho(); } diff --git a/source/module_elecstate/module_charge/charge.cpp b/source/module_elecstate/module_charge/charge.cpp index 5130157683..4827ff7268 100644 --- a/source/module_elecstate/module_charge/charge.cpp +++ b/source/module_elecstate/module_charge/charge.cpp @@ -41,40 +41,49 @@ Charge::Charge() Charge::~Charge() { - //if(allocate_rho) //LiuXh modify 20180619 - if(allocate_rho || allocate_rho_final_scf) //LiuXh add 20180619 - { - for(int i=0; idestroy(); } +void Charge::destroy() +{ + if(allocate_rho || allocate_rho_final_scf) //LiuXh add 20180619 + { + for(int i=0; idestroy(); + allocate_rho = false; + } + assert(allocate_rho == false); // mohan add 2021-02-20 diff --git a/source/module_elecstate/module_charge/charge.h b/source/module_elecstate/module_charge/charge.h index 7a51436a86..c738ea796a 100644 --- a/source/module_elecstate/module_charge/charge.h +++ b/source/module_elecstate/module_charge/charge.h @@ -29,21 +29,19 @@ class Charge // NAME : rhog_core [ngm], the core charge in reciprocal space //========================================================== - double **rho; - double **rho_save; + double **rho = nullptr; + double **rho_save = nullptr; - std::complex **rhog; - std::complex **rhog_save; + std::complex **rhog = nullptr; + std::complex **rhog_save = nullptr; - double **kin_r; // kinetic energy density in real space, for meta-GGA - double **kin_r_save; // kinetic energy density in real space, for meta-GGA - // wenfei 2021-07-28 + double **kin_r = nullptr; // kinetic energy density in real space, for meta-GGA + double **kin_r_save = nullptr; // kinetic energy density in real space, for meta-GGA + // wenfei 2021-07-28 - double *rho_core; - std::complex *rhog_core; + double *rho_core = nullptr; + std::complex *rhog_core = nullptr; - double *start_mag_type; - double *start_mag_atom; int prenspin = 1; void init_rho(); @@ -87,6 +85,8 @@ class Charge private: double sum_rho(void) const; + void destroy(); // free arrays liuyu 2023-03-12 + bool allocate_rho; bool allocate_rho_final_scf; // LiuXh add 20180606 diff --git a/source/module_esolver/esolver_ks_pw.cpp b/source/module_esolver/esolver_ks_pw.cpp index 76224c9ba5..6ed49a80a8 100644 --- a/source/module_esolver/esolver_ks_pw.cpp +++ b/source/module_esolver/esolver_ks_pw.cpp @@ -238,12 +238,20 @@ namespace ModuleESolver this->CE.save_pos_next(GlobalC::ucell); this->CE.extrapolate_charge(this->pelec->charge); - if(GlobalC::ucell.cell_parameter_updated) + // different precision level for vc-md + if(GlobalC::ucell.cell_parameter_updated && GlobalV::MD_PREC_LEVEL) + { + this->init_after_vc(INPUT, GlobalC::ucell); + } + else { Variable_Cell::init_after_vc(); + GlobalC::wfcpw->initgrids(GlobalC::ucell.lat0, GlobalC::ucell.latvec, GlobalC::wfcpw->nx, GlobalC::wfcpw->ny, GlobalC::wfcpw->nz); + GlobalC::wfcpw->initparameters(false, INPUT.ecutwfc, GlobalC::kv.nks, GlobalC::kv.kvec_d.data()); + GlobalC::wfcpw->collect_local_pw(); + GlobalC::wf.init_after_vc(GlobalC::kv.nks); + GlobalC::wf.init_at_1(); } - - //this->pelec->init_scf(istep, GlobalC::sf.strucFac); } if(GlobalV::CALCULATION=="relax" || GlobalV::CALCULATION=="cell-relax") @@ -262,15 +270,17 @@ namespace ModuleESolver // the new charge density. //this->pelec->init_scf( istep, GlobalC::sf.strucFac ); } + + if(GlobalC::ucell.cell_parameter_updated) + { + GlobalC::wfcpw->initgrids(GlobalC::ucell.lat0, GlobalC::ucell.latvec, GlobalC::wfcpw->nx, GlobalC::wfcpw->ny, GlobalC::wfcpw->nz); + GlobalC::wfcpw->initparameters(false, INPUT.ecutwfc, GlobalC::kv.nks, GlobalC::kv.kvec_d.data()); + GlobalC::wfcpw->collect_local_pw(); + GlobalC::wf.init_after_vc(GlobalC::kv.nks); + GlobalC::wf.init_at_1(); + } } - if(GlobalC::ucell.cell_parameter_updated) - { - GlobalC::wfcpw->initgrids(GlobalC::ucell.lat0, GlobalC::ucell.latvec, GlobalC::wfcpw->nx, GlobalC::wfcpw->ny, GlobalC::wfcpw->nz); - GlobalC::wfcpw->initparameters(false, INPUT.ecutwfc, GlobalC::kv.nks, GlobalC::kv.kvec_d.data()); - GlobalC::wfcpw->collect_local_pw(); - GlobalC::wf.init_after_vc(GlobalC::kv.nks); - GlobalC::wf.init_at_1(); - } + //init Hamilt, this should be allocated before each scf loop //Operators in HamiltPW should be reallocated once cell changed //delete Hamilt if not first scf @@ -578,6 +588,7 @@ namespace ModuleESolver void ESolver_KS_PW::cal_Force(ModuleBase::matrix& force) { Forces ff; + if (this->__kspw_psi != nullptr) this->__kspw_psi = nullptr; if (this->__kspw_psi == nullptr) { this->__kspw_psi = GlobalV::precision_flag == "single" ? new psi::Psi, Device>(this->kspw_psi[0]) : @@ -590,6 +601,7 @@ namespace ModuleESolver void ESolver_KS_PW::cal_Stress(ModuleBase::matrix& stress) { Stress_PW ss(this->pelec); + if (this->__kspw_psi != nullptr) this->__kspw_psi = nullptr; if (this->__kspw_psi == nullptr) { this->__kspw_psi = GlobalV::precision_flag == "single" ? new psi::Psi, Device>(this->kspw_psi[0]) : diff --git a/source/module_md/MSST.cpp b/source/module_md/MSST.cpp index 039c60e500..a6c76b84f3 100644 --- a/source/module_md/MSST.cpp +++ b/source/module_md/MSST.cpp @@ -8,7 +8,8 @@ MSST::MSST(MD_parameters& MD_para_in, UnitCell &unit_in) : MDrun(MD_para_in, unit_in) { - std::cout << "MSST" << std::endl; + GlobalV::MD_PREC_LEVEL = mdp.md_prec_level; + ucell.cell_parameter_updated = true; 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; diff --git a/source/module_md/Nose_Hoover.cpp b/source/module_md/Nose_Hoover.cpp index 7ba9ee0cf4..44388f0e4c 100644 --- a/source/module_md/Nose_Hoover.cpp +++ b/source/module_md/Nose_Hoover.cpp @@ -17,6 +17,12 @@ Nose_Hoover::Nose_Hoover(MD_parameters& MD_para_in, UnitCell &unit_in) : MDrun(M ModuleBase::WARNING_QUIT("Nose_Hoover", " md_tfirst must be larger than 0 in NHC !!! "); } + if(mdp.md_pmode != "none") + { + GlobalV::MD_PREC_LEVEL = mdp.md_prec_level; + ucell.cell_parameter_updated = true; + } + // init NPT related variables for(int i=0; i<6; ++i) { diff --git a/source/module_md/run_md.cpp b/source/module_md/run_md.cpp index 3209df4310..18a340d501 100644 --- a/source/module_md/run_md.cpp +++ b/source/module_md/run_md.cpp @@ -39,7 +39,6 @@ void Run_MD::md_line(UnitCell &unit_in, ModuleESolver::ESolver *p_esolver) else if(INPUT.mdp.md_type == 4) { mdrun = new MSST(INPUT.mdp, unit_in); - unit_in.cell_parameter_updated = true; } // md cycle From e66c3bbaee04468c5bae1470e4b650d5f38d6e6c Mon Sep 17 00:00:00 2001 From: YuLiu98 Date: Mon, 13 Mar 2023 16:30:02 +0800 Subject: [PATCH 05/17] Fix: fix original vc-md --- source/module_esolver/esolver_ks_pw.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/source/module_esolver/esolver_ks_pw.cpp b/source/module_esolver/esolver_ks_pw.cpp index 6ed49a80a8..4d5db99379 100644 --- a/source/module_esolver/esolver_ks_pw.cpp +++ b/source/module_esolver/esolver_ks_pw.cpp @@ -243,7 +243,7 @@ namespace ModuleESolver { this->init_after_vc(INPUT, GlobalC::ucell); } - else + else if(GlobalC::ucell.cell_parameter_updated) { Variable_Cell::init_after_vc(); GlobalC::wfcpw->initgrids(GlobalC::ucell.lat0, GlobalC::ucell.latvec, GlobalC::wfcpw->nx, GlobalC::wfcpw->ny, GlobalC::wfcpw->nz); From 3bfb5558d9a20321bc95b6144dc41037a61039c0 Mon Sep 17 00:00:00 2001 From: YuLiu98 Date: Mon, 13 Mar 2023 20:08:28 +0800 Subject: [PATCH 06/17] Fix: reset ig_gge0 when reinit FFT --- source/module_base/global_function.h | 2 +- source/module_esolver/esolver_fp.cpp | 1 - source/module_pw/pw_basis.cpp | 2 ++ 3 files changed, 3 insertions(+), 2 deletions(-) diff --git a/source/module_base/global_function.h b/source/module_base/global_function.h index fc908c4df5..fa77c76b26 100644 --- a/source/module_base/global_function.h +++ b/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 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; } diff --git a/source/module_esolver/esolver_fp.cpp b/source/module_esolver/esolver_fp.cpp index 05ac583c06..f603d91721 100644 --- a/source/module_esolver/esolver_fp.cpp +++ b/source/module_esolver/esolver_fp.cpp @@ -50,7 +50,6 @@ namespace ModuleESolver this->pw_rho->setuptransform(); this->pw_rho->collect_local_pw(); this->pw_rho->collect_uniqgg(); - this->print_rhofft(inp, GlobalV::ofs_running); } void ESolver_FP::print_rhofft(Input&inp, ofstream &ofs) diff --git a/source/module_pw/pw_basis.cpp b/source/module_pw/pw_basis.cpp index 6f7b71e0b7..1bbaab28fc 100644 --- a/source/module_pw/pw_basis.cpp +++ b/source/module_pw/pw_basis.cpp @@ -108,6 +108,7 @@ void PW_Basis::getstartgr() void PW_Basis::collect_local_pw() { if(this->npw <= 0) return; + this->ig_gge0 = -1; delete[] this->gg; this->gg = new double[this->npw]; delete[] this->gdirect; this->gdirect = new ModuleBase::Vector3[this->npw]; delete[] this->gcar; this->gcar = new ModuleBase::Vector3[this->npw]; @@ -143,6 +144,7 @@ void PW_Basis::collect_local_pw() void PW_Basis::collect_uniqgg() { if(this->npw <= 0) return; + this->ig_gge0 = -1; delete[] this->ig2igg; this->ig2igg = new int [this->npw]; int *sortindex = new int [this->npw]; double *tmpgg = new double [this->npw]; From 2ca2a60aa50f91b46c2d6a08c19bb28f8d1d2595 Mon Sep 17 00:00:00 2001 From: YuLiu98 Date: Fri, 17 Mar 2023 14:55:12 +0800 Subject: [PATCH 07/17] Feature: update esolver_lcao --- source/module_esolver/esolver_dp.cpp | 5 ++++ source/module_esolver/esolver_ks_lcao.cpp | 25 +++++++++++++++++++ source/module_esolver/esolver_ks_lcao.h | 1 + .../module_esolver/esolver_ks_lcao_elec.cpp | 19 +++++++++----- source/module_esolver/esolver_ks_pw.cpp | 24 ++++++++++-------- 5 files changed, 57 insertions(+), 17 deletions(-) diff --git a/source/module_esolver/esolver_dp.cpp b/source/module_esolver/esolver_dp.cpp index c2d6e1f9bc..31dba2332e 100644 --- a/source/module_esolver/esolver_dp.cpp +++ b/source/module_esolver/esolver_dp.cpp @@ -1,4 +1,5 @@ #include "esolver_dp.h" +#include "module_base/timer.h" namespace ModuleESolver @@ -17,6 +18,9 @@ namespace ModuleESolver void ESolver_DP::Run(const int istep, UnitCell& ucell) { + ModuleBase::TITLE("ESolver_DP", "Run"); + ModuleBase::timer::tick("ESolver_DP", "Run"); + cell[0] = ucell.latvec.e11 * ucell.lat0_angstrom; cell[1] = ucell.latvec.e12 * ucell.lat0_angstrom; cell[2] = ucell.latvec.e13 * ucell.lat0_angstrom; @@ -71,6 +75,7 @@ namespace ModuleESolver #else ModuleBase::WARNING_QUIT("DP_pot", "Please recompile with -D__DPMD !"); #endif + ModuleBase::timer::tick("ESolver_DP", "Run"); } void ESolver_DP::cal_Energy(double& etot) diff --git a/source/module_esolver/esolver_ks_lcao.cpp b/source/module_esolver/esolver_ks_lcao.cpp index 578cdd7a6c..17af48e6a9 100644 --- a/source/module_esolver/esolver_ks_lcao.cpp +++ b/source/module_esolver/esolver_ks_lcao.cpp @@ -195,6 +195,31 @@ void ESolver_KS_LCAO::Init(Input& inp, UnitCell& ucell) } } +void ESolver_KS_LCAO::init_after_vc(Input& inp, UnitCell& ucell) +{ + ESolver_KS::init_after_vc(inp, ucell); + + delete this->pelec; + this->pelec = new elecstate::ElecStateLCAO(&(chr), &(GlobalC::kv), GlobalC::kv.nks, &(this->LOC), &(this->UHM), &(this->LOWF)); + + GlobalC::ppcell.init_vloc(GlobalC::ppcell.vloc, GlobalC::rhopw); + + this->pelec->charge->allocate(GlobalV::NSPIN, GlobalC::rhopw->nrxx, GlobalC::rhopw->npw); + + if(this->pelec->pot != nullptr) + { + delete this->pelec->pot; + this->pelec->pot = new elecstate::Potential( + GlobalC::rhopw, + &GlobalC::ucell, + &(GlobalC::ppcell.vloc), + &(GlobalC::sf.strucFac), + &(GlobalC::en.etxc), + &(GlobalC::en.vtxc) + ); + } +} + void ESolver_KS_LCAO::cal_Energy(double& etot) { etot = GlobalC::en.etot; diff --git a/source/module_esolver/esolver_ks_lcao.h b/source/module_esolver/esolver_ks_lcao.h index 86eaf10233..23c6e15275 100644 --- a/source/module_esolver/esolver_ks_lcao.h +++ b/source/module_esolver/esolver_ks_lcao.h @@ -18,6 +18,7 @@ namespace ModuleESolver ~ESolver_KS_LCAO(); void Init(Input& inp, UnitCell& cell) override; + void init_after_vc(Input& inp, UnitCell& cell) override; void cal_Energy(double& etot) override; void cal_Force(ModuleBase::matrix& force) override; diff --git a/source/module_esolver/esolver_ks_lcao_elec.cpp b/source/module_esolver/esolver_ks_lcao_elec.cpp index 793a815aa4..052bb92daf 100644 --- a/source/module_esolver/esolver_ks_lcao_elec.cpp +++ b/source/module_esolver/esolver_ks_lcao_elec.cpp @@ -261,13 +261,20 @@ namespace ModuleESolver // Temporary, md and relax will merge later liuyu add 2022-11-07 if(GlobalV::CALCULATION == "md" && istep) { - CE.update_istep(); - CE.save_pos_next(GlobalC::ucell); - CE.extrapolate_charge(pelec->charge); - - if(GlobalC::ucell.cell_parameter_updated) + // different precision level for vc-md + if(GlobalC::ucell.cell_parameter_updated && GlobalV::MD_PREC_LEVEL) + { + this->init_after_vc(INPUT, GlobalC::ucell); + } + else { - Variable_Cell::init_after_vc(); + this->CE.update_istep(); + this->CE.save_pos_next(GlobalC::ucell); + this->CE.extrapolate_charge(this->pelec->charge); + if(GlobalC::ucell.cell_parameter_updated) + { + Variable_Cell::init_after_vc(); + } } } diff --git a/source/module_esolver/esolver_ks_pw.cpp b/source/module_esolver/esolver_ks_pw.cpp index 4d5db99379..e46d66dfd8 100644 --- a/source/module_esolver/esolver_ks_pw.cpp +++ b/source/module_esolver/esolver_ks_pw.cpp @@ -234,23 +234,25 @@ namespace ModuleESolver // Temporary, md and relax will merge later liuyu add 2022-11-07 if(GlobalV::CALCULATION == "md" && istep) { - this->CE.update_istep(); - this->CE.save_pos_next(GlobalC::ucell); - this->CE.extrapolate_charge(this->pelec->charge); - // different precision level for vc-md if(GlobalC::ucell.cell_parameter_updated && GlobalV::MD_PREC_LEVEL) { this->init_after_vc(INPUT, GlobalC::ucell); } - else if(GlobalC::ucell.cell_parameter_updated) + else { - Variable_Cell::init_after_vc(); - GlobalC::wfcpw->initgrids(GlobalC::ucell.lat0, GlobalC::ucell.latvec, GlobalC::wfcpw->nx, GlobalC::wfcpw->ny, GlobalC::wfcpw->nz); - GlobalC::wfcpw->initparameters(false, INPUT.ecutwfc, GlobalC::kv.nks, GlobalC::kv.kvec_d.data()); - GlobalC::wfcpw->collect_local_pw(); - GlobalC::wf.init_after_vc(GlobalC::kv.nks); - GlobalC::wf.init_at_1(); + this->CE.update_istep(); + this->CE.save_pos_next(GlobalC::ucell); + this->CE.extrapolate_charge(this->pelec->charge); + if(GlobalC::ucell.cell_parameter_updated) + { + Variable_Cell::init_after_vc(); + GlobalC::wfcpw->initgrids(GlobalC::ucell.lat0, GlobalC::ucell.latvec, GlobalC::wfcpw->nx, GlobalC::wfcpw->ny, GlobalC::wfcpw->nz); + GlobalC::wfcpw->initparameters(false, INPUT.ecutwfc, GlobalC::kv.nks, GlobalC::kv.kvec_d.data()); + GlobalC::wfcpw->collect_local_pw(); + GlobalC::wf.init_after_vc(GlobalC::kv.nks); + GlobalC::wf.init_at_1(); + } } } From 72dc9078a95cb64c0951b8a4f217a21368330f5e Mon Sep 17 00:00:00 2001 From: YuLiu98 Date: Wed, 22 Mar 2023 10:22:02 +0800 Subject: [PATCH 08/17] Feature: enable reference cell for vc-md --- docs/advanced/input_files/input-main.md | 19 ++++++++++++----- source/module_base/global_variable.cpp | 2 +- source/module_base/global_variable.h | 2 +- source/module_elecstate/elecstate.cpp | 2 +- source/module_esolver/esolver_fp.cpp | 4 ++-- source/module_esolver/esolver_ks.cpp | 2 +- source/module_esolver/esolver_ks_lcao.cpp | 3 --- .../module_esolver/esolver_ks_lcao_elec.cpp | 5 +++-- source/module_esolver/esolver_ks_pw.cpp | 8 +++---- source/module_io/input.cpp | 21 +++++++++++++++++++ source/module_io/input.h | 1 + source/module_io/test/input_test.cpp | 2 ++ source/module_io/test/support/INPUT | 1 + source/module_io/write_input.cpp | 3 ++- source/module_md/MD_parameters.h | 2 +- source/module_md/MSST.cpp | 1 - source/module_md/Nose_Hoover.cpp | 1 - 17 files changed, 54 insertions(+), 25 deletions(-) diff --git a/docs/advanced/input_files/input-main.md b/docs/advanced/input_files/input-main.md index fb654addb1..c4a0828146 100644 --- a/docs/advanced/input_files/input-main.md +++ b/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_prec_level](#md_prec_level) | [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 @@ -1554,7 +1554,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 @@ -1618,11 +1618,20 @@ These variables are used to control the molecular dynamics calculations. ### md_prec_level - **Type**: Integer -- **Description**: Determine the precision level of vc-md. Note that this parameter is only used in variable-cell MD! +- **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: 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. + - 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 diff --git a/source/module_base/global_variable.cpp b/source/module_base/global_variable.cpp index e8984c0916..6afa72f75b 100644 --- a/source/module_base/global_variable.cpp +++ b/source/module_base/global_variable.cpp @@ -57,7 +57,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 md_prec_level = 0; int SCF_NMAX = 100; std::string BASIS_TYPE = "pw"; // xiaohui add 2013-09-01 diff --git a/source/module_base/global_variable.h b/source/module_base/global_variable.h index bb6cddcd30..c79125436c 100644 --- a/source/module_base/global_variable.h +++ b/source/module_base/global_variable.h @@ -62,7 +62,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 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 diff --git a/source/module_elecstate/elecstate.cpp b/source/module_elecstate/elecstate.cpp index 38ffbcd49e..4c8aecb585 100644 --- a/source/module_elecstate/elecstate.cpp +++ b/source/module_elecstate/elecstate.cpp @@ -378,7 +378,7 @@ void ElecState::init_scf(const int istep, const ModuleBase::ComplexMatrix& struc // (2) other effective potentials need charge density, // choose charge density from ionic step 0. //-------------------------------------------------------------------- - if (istep == 0 || GlobalV::MD_PREC_LEVEL == 1) + if (istep == 0 || GlobalV::md_prec_level == 2) { this->charge->init_rho(); } diff --git a/source/module_esolver/esolver_fp.cpp b/source/module_esolver/esolver_fp.cpp index f603d91721..1c281bb65a 100644 --- a/source/module_esolver/esolver_fp.cpp +++ b/source/module_esolver/esolver_fp.cpp @@ -27,9 +27,9 @@ namespace ModuleESolver if (this->classname == "ESolver_OF") this->pw_rho->setfullpw(inp.of_full_pw, inp.of_full_pw_dim); // Initalize the plane wave basis set if (inp.nx * inp.ny * inp.nz == 0) - this->pw_rho->initgrids(cell.lat0, cell.latvec, inp.ecutrho); + this->pw_rho->initgrids(inp.ref_cell_factor * cell.lat0, cell.latvec, inp.ecutrho); else - this->pw_rho->initgrids(cell.lat0, cell.latvec, inp.nx, inp.ny, inp.nz); + this->pw_rho->initgrids(inp.ref_cell_factor * cell.lat0, cell.latvec, inp.nx, inp.ny, inp.nz); this->pw_rho->initparameters(false, inp.ecutrho); this->pw_rho->setuptransform(); diff --git a/source/module_esolver/esolver_ks.cpp b/source/module_esolver/esolver_ks.cpp index ebddc3407d..852147f21d 100644 --- a/source/module_esolver/esolver_ks.cpp +++ b/source/module_esolver/esolver_ks.cpp @@ -79,7 +79,7 @@ namespace ModuleESolver #ifdef __MPI this->pw_wfc->initmpi(GlobalV::NPROC_IN_POOL, GlobalV::RANK_IN_POOL, POOL_WORLD); #endif - this->pw_wfc->initgrids(ucell.lat0, ucell.latvec, GlobalC::rhopw->nx, GlobalC::rhopw->ny, GlobalC::rhopw->nz); + this->pw_wfc->initgrids(inp.ref_cell_factor * ucell.lat0, ucell.latvec, GlobalC::rhopw->nx, GlobalC::rhopw->ny, GlobalC::rhopw->nz); this->pw_wfc->initparameters(false, inp.ecutwfc, GlobalC::kv.nks, GlobalC::kv.kvec_d.data()); #ifdef __MPI if(INPUT.pw_seed > 0) MPI_Allreduce(MPI_IN_PLACE, &this->pw_wfc->ggecut, 1, MPI_DOUBLE, MPI_MAX , MPI_COMM_WORLD); diff --git a/source/module_esolver/esolver_ks_lcao.cpp b/source/module_esolver/esolver_ks_lcao.cpp index 3da148257a..35a418dede 100644 --- a/source/module_esolver/esolver_ks_lcao.cpp +++ b/source/module_esolver/esolver_ks_lcao.cpp @@ -788,9 +788,6 @@ void ESolver_KS_LCAO::eachiterfinish(int iter) void ESolver_KS_LCAO::afterscf(const int istep) { - // Temporary liuyu add 2022-11-07 - CE.update_all_pos(GlobalC::ucell); - // if (this->conv_elec || iter == GlobalV::SCF_NMAX) // { //-------------------------------------- diff --git a/source/module_esolver/esolver_ks_lcao_elec.cpp b/source/module_esolver/esolver_ks_lcao_elec.cpp index 052bb92daf..829b554d5c 100644 --- a/source/module_esolver/esolver_ks_lcao_elec.cpp +++ b/source/module_esolver/esolver_ks_lcao_elec.cpp @@ -262,16 +262,17 @@ namespace ModuleESolver if(GlobalV::CALCULATION == "md" && istep) { // different precision level for vc-md - if(GlobalC::ucell.cell_parameter_updated && GlobalV::MD_PREC_LEVEL) + if(GlobalC::ucell.cell_parameter_updated && GlobalV::md_prec_level == 2) { this->init_after_vc(INPUT, GlobalC::ucell); } else { this->CE.update_istep(); + this->CE.update_all_pos(GlobalC::ucell); this->CE.save_pos_next(GlobalC::ucell); this->CE.extrapolate_charge(this->pelec->charge); - if(GlobalC::ucell.cell_parameter_updated) + if(GlobalC::ucell.cell_parameter_updated && GlobalV::md_prec_level == 0) { Variable_Cell::init_after_vc(); } diff --git a/source/module_esolver/esolver_ks_pw.cpp b/source/module_esolver/esolver_ks_pw.cpp index a1924a7d11..fcc4b577cf 100644 --- a/source/module_esolver/esolver_ks_pw.cpp +++ b/source/module_esolver/esolver_ks_pw.cpp @@ -235,16 +235,17 @@ namespace ModuleESolver if(GlobalV::CALCULATION == "md" && istep) { // different precision level for vc-md - if(GlobalC::ucell.cell_parameter_updated && GlobalV::MD_PREC_LEVEL) + if(GlobalC::ucell.cell_parameter_updated && GlobalV::md_prec_level == 2) { this->init_after_vc(INPUT, GlobalC::ucell); } else { this->CE.update_istep(); + this->CE.update_all_pos(GlobalC::ucell); this->CE.save_pos_next(GlobalC::ucell); this->CE.extrapolate_charge(this->pelec->charge); - if(GlobalC::ucell.cell_parameter_updated) + if(GlobalC::ucell.cell_parameter_updated && GlobalV::md_prec_level == 0) { Variable_Cell::init_after_vc(); GlobalC::wfcpw->initgrids(GlobalC::ucell.lat0, GlobalC::ucell.latvec, GlobalC::wfcpw->nx, GlobalC::wfcpw->ny, GlobalC::wfcpw->nz); @@ -520,9 +521,6 @@ namespace ModuleESolver template void ESolver_KS_PW::afterscf(const int istep) { - // Temporary liuyu add 2022-11-07 - this->CE.update_all_pos(GlobalC::ucell); - for (int is = 0; is < GlobalV::NSPIN; is++) { std::stringstream ssc; diff --git a/source/module_io/input.cpp b/source/module_io/input.cpp index 3dcc3e22c1..02e0ad6361 100644 --- a/source/module_io/input.cpp +++ b/source/module_io/input.cpp @@ -185,6 +185,7 @@ void Input::Default(void) search_pbc = true; symmetry = 0; init_vel = false; + ref_cell_factor = 1.0; symmetry_prec = 1.0e-5; // LiuXh add 2021-08-12, accuracy for symmetry cal_force = 0; dump_force = true; @@ -841,6 +842,10 @@ bool Input::Read(const std::string &fn) { read_bool(ifs, init_vel); } + else if (strcmp("ref_cell_factor", word) == 0) + { + read_value(ifs, ref_cell_factor); + } else if (strcmp("symmetry_prec", word) == 0) // LiuXh add 2021-08-12, accuracy for symmetry { read_value(ifs, symmetry_prec); @@ -2565,6 +2570,11 @@ void Input::Default_2(void) // jiyy add 2019-08-04 { cal_stress = 1; } + + if(mdp.md_type == 4 || (mdp.md_type == 1 && mdp.md_pmode != "none")) + { + GlobalV::md_prec_level = mdp.md_prec_level; + } } else if (calculation == "cell-relax") // mohan add 2011-11-04 { @@ -2651,6 +2661,11 @@ void Input::Default_2(void) // jiyy add 2019-08-04 { bessel_descriptor_ecut = std::to_string(ecutwfc); } + + if (GlobalV::md_prec_level != 1) + { + ref_cell_factor = 1.0; + } } #ifdef __MPI void Input::Bcast() @@ -2717,6 +2732,7 @@ void Input::Bcast() Parallel_Common::bcast_double(search_radius); Parallel_Common::bcast_int(symmetry); Parallel_Common::bcast_bool(init_vel); // liuyu 2021-07-14 + Parallel_Common::bcast_int(ref_cell_factor); Parallel_Common::bcast_double(symmetry_prec); // LiuXh add 2021-08-12, accuracy for symmetry Parallel_Common::bcast_bool(cal_force); Parallel_Common::bcast_bool(dump_force); @@ -3133,6 +3149,11 @@ void Input::Check(void) ModuleBase::WARNING_QUIT("Input", "gate field cannot be used with efield if dip_cor_flag=false !"); } + if(ref_cell_factor < 1.0) + { + ModuleBase::WARNING_QUIT("Input", "ref_cell_factor must not be less than 1.0"); + } + //---------------------------------------------------------- // main parameters / electrons / spin ( 1/16 ) //---------------------------------------------------------- diff --git a/source/module_io/input.h b/source/module_io/input.h index 2a8087cd29..f292ef5086 100644 --- a/source/module_io/input.h +++ b/source/module_io/input.h @@ -47,6 +47,7 @@ class Input bool init_vel; // read velocity from STRU or not liuyu 2021-07-14 bool dump_vel; // output atomic velocities into the file MD_dump or not. liuyu 2023-03-01 + int ref_cell_factor; // construct a reference cell bigger than the initial cell liuyu 2023-03-21 /* symmetry level: -1, no symmetry at all; diff --git a/source/module_io/test/input_test.cpp b/source/module_io/test/input_test.cpp index 4a2774a444..fc4ad0b2da 100644 --- a/source/module_io/test/input_test.cpp +++ b/source/module_io/test/input_test.cpp @@ -79,6 +79,7 @@ TEST_F(InputTest, Default) EXPECT_TRUE(INPUT.search_pbc); EXPECT_EQ(INPUT.symmetry,0); EXPECT_FALSE(INPUT.init_vel); + EXPECT_DOUBLE_EQ(INPUT.ref_cell_factor,1.0); EXPECT_DOUBLE_EQ(INPUT.symmetry_prec,1.0e-5); EXPECT_EQ(INPUT.cal_force,0); EXPECT_TRUE(INPUT.dump_force); @@ -669,6 +670,7 @@ TEST_F(InputTest, Read) EXPECT_EQ(INPUT.mdp.md_restartfreq,5); EXPECT_EQ(INPUT.mdp.md_seed,-1); EXPECT_EQ(INPUT.mdp.md_prec_level,1); + EXPECT_DOUBLE_EQ(INPUT.ref_cell_factor,1.2); EXPECT_EQ(INPUT.mdp.md_tchain,1); EXPECT_DOUBLE_EQ(INPUT.mdp.md_tfirst,-1); EXPECT_DOUBLE_EQ(INPUT.mdp.md_tfreq,0); diff --git a/source/module_io/test/support/INPUT b/source/module_io/test/support/INPUT index 0a13529c04..074b8b06f5 100644 --- a/source/module_io/test/support/INPUT +++ b/source/module_io/test/support/INPUT @@ -179,6 +179,7 @@ md_dumpfreq 1 #The period to dump MD information md_restartfreq 5 #The period to output MD restart information md_seed -1 #random seed for MD md_prec_level 1 #precision level for vc-md +ref_cell_factor 1.2 #construct a reference cell bigger than the initial cell md_restart 0 #whether restart lj_rcut 8.5 #cutoff radius of LJ potential lj_epsilon 0.01032 #the value of epsilon for LJ potential diff --git a/source/module_io/write_input.cpp b/source/module_io/write_input.cpp index 6b7432e4a9..54a62b2db7 100644 --- a/source/module_io/write_input.cpp +++ b/source/module_io/write_input.cpp @@ -248,7 +248,8 @@ ModuleBase::GlobalFunc::OUTP(ofs, "out_bandgap", out_bandgap, "if true, print ou ModuleBase::GlobalFunc::OUTP(ofs,"md_dumpfreq",mdp.md_dumpfreq,"The period to dump MD information"); ModuleBase::GlobalFunc::OUTP(ofs,"md_restartfreq",mdp.md_restartfreq,"The period to output MD restart information"); ModuleBase::GlobalFunc::OUTP(ofs,"md_seed",mdp.md_seed,"random seed for MD"); - ModuleBase::GlobalFunc::OUTP(ofs,"md_prec_level",mdp.md_prec_level,"precision level for vc-md"); + ModuleBase::GlobalFunc::OUTP(ofs,"md_prec_level",GlobalV::md_prec_level,"precision level for vc-md"); + ModuleBase::GlobalFunc::OUTP(ofs,"ref_cell_factor",ref_cell_factor,"construct a reference cell bigger than the initial cell"); ModuleBase::GlobalFunc::OUTP(ofs,"md_restart",mdp.md_restart,"whether restart"); ModuleBase::GlobalFunc::OUTP(ofs,"lj_rcut",mdp.lj_rcut,"cutoff radius of LJ potential"); ModuleBase::GlobalFunc::OUTP(ofs,"lj_epsilon",mdp.lj_epsilon,"the value of epsilon for LJ potential"); diff --git a/source/module_md/MD_parameters.h b/source/module_md/MD_parameters.h index b108d7589b..549c0baf2a 100644 --- a/source/module_md/MD_parameters.h +++ b/source/module_md/MD_parameters.h @@ -62,7 +62,7 @@ class MD_parameters int md_dumpfreq; // The period to dump MD information int md_restartfreq; // The period to output MD restart information int md_seed; // random seed for MD - int md_prec_level; // precision level for vc-md; 0: do not reinit FFT grid, 1: reinit FFT grid per step + int md_prec_level; // precision level for vc-md; 0: do not reinit FFT grid, 1: reference cell, 2: reinit FFT grid per step // Classic MD // liuyu 2021-07-30 double lj_rcut; // cutoff radius of LJ potential (\AA) diff --git a/source/module_md/MSST.cpp b/source/module_md/MSST.cpp index a6c76b84f3..2a156d0a59 100644 --- a/source/module_md/MSST.cpp +++ b/source/module_md/MSST.cpp @@ -8,7 +8,6 @@ MSST::MSST(MD_parameters& MD_para_in, UnitCell &unit_in) : MDrun(MD_para_in, unit_in) { - GlobalV::MD_PREC_LEVEL = mdp.md_prec_level; ucell.cell_parameter_updated = true; mdp.msst_qmass = mdp.msst_qmass / pow(ModuleBase::ANGSTROM_AU, 4) / pow(ModuleBase::AU_to_MASS, 2); diff --git a/source/module_md/Nose_Hoover.cpp b/source/module_md/Nose_Hoover.cpp index 44388f0e4c..d213a3c1c7 100644 --- a/source/module_md/Nose_Hoover.cpp +++ b/source/module_md/Nose_Hoover.cpp @@ -19,7 +19,6 @@ Nose_Hoover::Nose_Hoover(MD_parameters& MD_para_in, UnitCell &unit_in) : MDrun(M if(mdp.md_pmode != "none") { - GlobalV::MD_PREC_LEVEL = mdp.md_prec_level; ucell.cell_parameter_updated = true; } From e115cc4a3d8bf9afe7db08703d83bca915eb5471 Mon Sep 17 00:00:00 2001 From: YuLiu98 Date: Wed, 22 Mar 2023 13:13:45 +0800 Subject: [PATCH 09/17] Refactor: delete useless functions and variables --- source/module_cell/unitcell.cpp | 87 +++++-------------- source/module_cell/unitcell.h | 2 - .../module_charge/charge_extra.cpp | 10 ++- .../module_charge/charge_extra.h | 1 + source/module_esolver/esolver_ks_pw.cpp | 2 +- source/module_md/FIRE.cpp | 11 ++- source/module_md/Langevin.cpp | 11 ++- source/module_md/MD_func.cpp | 15 ---- source/module_md/MD_func.h | 2 - source/module_md/MSST.cpp | 17 +++- source/module_md/Nose_Hoover.cpp | 12 +-- source/module_md/mdrun.cpp | 12 +-- source/module_md/mdrun.h | 2 +- source/module_md/test/MD_func_test.cpp | 22 ----- source/module_relax/relax_old/bfgs_basic.cpp | 18 +--- source/module_relax/relax_old/bfgs_basic.h | 8 +- .../relax_old/ions_move_basic.cpp | 36 ++++---- .../module_relax/relax_old/ions_move_basic.h | 4 +- .../module_relax/relax_old/ions_move_bfgs.cpp | 8 +- .../module_relax/relax_old/ions_move_cg.cpp | 17 +--- source/module_relax/relax_old/ions_move_cg.h | 1 - .../module_relax/relax_old/ions_move_sd.cpp | 13 +-- source/module_relax/relax_old/ions_move_sd.h | 1 - 23 files changed, 111 insertions(+), 201 deletions(-) diff --git a/source/module_cell/unitcell.cpp b/source/module_cell/unitcell.cpp index faeb70d6a1..62e6fe3e9c 100644 --- a/source/module_cell/unitcell.cpp +++ b/source/module_cell/unitcell.cpp @@ -293,72 +293,6 @@ void UnitCell::set_iat2itia(void) return; } -void UnitCell::update_pos_tau(const double* pos) -{ - int iat = 0; - for (int it = 0; it < this->ntype; it++) - { - Atom* atom = &this->atoms[it]; - for (int ia = 0; ia < atom->na; ia++) - { - if (atom->mbl[ia].x != 0) - { - atom->tau_original[ia].x += (pos[3 * iat] / this->lat0 - atom->tau[ia].x); - atom->tau[ia].x = pos[3 * iat] / this->lat0; - } - if (atom->mbl[ia].y != 0) - { - atom->tau_original[ia].y += (pos[3 * iat + 1] / this->lat0 - atom->tau[ia].y); - atom->tau[ia].y = pos[3 * iat + 1] / this->lat0; - } - if (atom->mbl[ia].z != 0) - { - atom->tau_original[ia].z += (pos[3 * iat + 2] / this->lat0 - atom->tau[ia].z); - atom->tau[ia].z = pos[3 * iat + 2] / this->lat0; - } - - // the direct coordinates also need to be updated. - atom->taud[ia] = atom->tau[ia] * this->GT; - iat++; - } - } - assert(iat == this->nat); - return; -} - -void UnitCell::update_pos_tau(const ModuleBase::Vector3* 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) - { - if (atom->mbl[ia].x != 0) - { - atom->tau_original[ia].x += (posd_in[iat].x / this->lat0 - atom->tau[ia].x); - atom->tau[ia].x = posd_in[iat].x / this->lat0; - } - if (atom->mbl[ia].y != 0) - { - atom->tau_original[ia].y += (posd_in[iat].y / this->lat0 - atom->tau[ia].y); - atom->tau[ia].y = posd_in[iat].y / this->lat0; - } - if (atom->mbl[ia].z != 0) - { - atom->tau_original[ia].z += (posd_in[iat].z / this->lat0 - atom->tau[ia].z); - atom->tau[ia].z = posd_in[iat].z / this->lat0; - } - - // the direct coordinates also need to be updated. - atom->taud[ia] = atom->tau[ia] * this->GT; - iat++; - } - } - assert(iat == this->nat); - return; -} - // Note : note here we are not keeping track of 'tau_original', namely // the Cartesian coordinate before periodic adjustment // The reason is that this is only used in relaxation @@ -384,6 +318,27 @@ void UnitCell::update_pos_taud(double* posd_in) this->bcast_atoms_tau(); } +// posd_in is atomic displacements here liuyu 2023-03-22 +void UnitCell::update_pos_taud(const ModuleBase::Vector3* posd_in) +{ + for (int it = 0; it < this->ntype; it++) + { + Atom* atom = &this->atoms[it]; + for (int ia = 0; ia < atom->na; ia++) + { + for ( int ik = 0; ik < 3; ++ik) + { + if (atom->mbl[ia][ik]) + { + atom->taud[ia][ik] += posd_in[ia][ik]; + } + } + } + } + this->periodic_boundary_adjustment(); + this->bcast_atoms_tau(); +} + void UnitCell::update_vel(const ModuleBase::Vector3* vel_in) { int iat = 0; diff --git a/source/module_cell/unitcell.h b/source/module_cell/unitcell.h index a832e3d407..f04c8b28fb 100644 --- a/source/module_cell/unitcell.h +++ b/source/module_cell/unitcell.h @@ -180,8 +180,6 @@ class UnitCell void print_cell_xyz(const std::string &fn)const; void print_cell_cif(const std::string &fn)const; - void update_pos_tau(const double* pos); - void update_pos_tau(const ModuleBase::Vector3* posd_in); void update_pos_taud(const ModuleBase::Vector3* posd_in); void update_pos_taud(double* posd_in); void update_vel(const ModuleBase::Vector3* vel_in); diff --git a/source/module_elecstate/module_charge/charge_extra.cpp b/source/module_elecstate/module_charge/charge_extra.cpp index f54d83b8b6..c8877cb74e 100644 --- a/source/module_elecstate/module_charge/charge_extra.cpp +++ b/source/module_elecstate/module_charge/charge_extra.cpp @@ -64,6 +64,7 @@ void Charge_Extra::Init_CE() } natom = GlobalC::ucell.nat; + omega_old = GlobalC::ucell.omega; pos_old1 = new ModuleBase::Vector3[natom]; pos_old2 = new ModuleBase::Vector3[natom]; @@ -100,7 +101,6 @@ void Charge_Extra::extrapolate_charge(Charge* chr) rho_extr = min(istep, pot_order); if(rho_extr == 0) { - // if(cellchange) scale(); GlobalC::sf.setup_structure_factor(&GlobalC::ucell, GlobalC::rhopw); GlobalV::ofs_running << " charge density from previous step !" << std::endl; return; @@ -125,6 +125,10 @@ void Charge_Extra::extrapolate_charge(Charge* chr) for(int ir=0; irnrxx; ir++) { chr->rho[is][ir] -= rho_atom[is][ir]; + if(GlobalC::ucell.cell_parameter_updated) + { + chr->rho[is][ir] *= omega_old; + } } } @@ -215,6 +219,10 @@ void Charge_Extra::extrapolate_charge(Charge* chr) { for(int ir=0; irnrxx; ir++) { + if(GlobalC::ucell.cell_parameter_updated) + { + chr->rho[is][ir] /= GlobalC::ucell.omega; + } chr->rho[is][ir] += rho_atom[is][ir]; } } diff --git a/source/module_elecstate/module_charge/charge_extra.h b/source/module_elecstate/module_charge/charge_extra.h index b913da7a18..1ec16b5b6e 100644 --- a/source/module_elecstate/module_charge/charge_extra.h +++ b/source/module_elecstate/module_charge/charge_extra.h @@ -33,6 +33,7 @@ class Charge_Extra int natom; int pot_order; int rho_extr; + double omega_old; // the old volume of the last step // for the second-order extrapolation ModuleBase::Vector3 *pos_old1 = nullptr; diff --git a/source/module_esolver/esolver_ks_pw.cpp b/source/module_esolver/esolver_ks_pw.cpp index fcc4b577cf..19bd501e1f 100644 --- a/source/module_esolver/esolver_ks_pw.cpp +++ b/source/module_esolver/esolver_ks_pw.cpp @@ -241,8 +241,8 @@ namespace ModuleESolver } else { - this->CE.update_istep(); this->CE.update_all_pos(GlobalC::ucell); + this->CE.update_istep(); this->CE.save_pos_next(GlobalC::ucell); this->CE.extrapolate_charge(this->pelec->charge); if(GlobalC::ucell.cell_parameter_updated && GlobalV::md_prec_level == 0) diff --git a/source/module_md/FIRE.cpp b/source/module_md/FIRE.cpp index c0054edf62..6876671dc7 100644 --- a/source/module_md/FIRE.cpp +++ b/source/module_md/FIRE.cpp @@ -57,9 +57,14 @@ void FIRE::first_half() { if(ionmbl[i][k]) { - pos[i][k] += vel[i][k]*mdp.md_dt; + pos[i][k] = vel[i][k] * mdp.md_dt / ucell.lat0; + } + else + { + pos[i][k] = 0; } } + pos[i] = pos[i] * ucell.GT; } } #ifdef __MPI @@ -67,9 +72,7 @@ void FIRE::first_half() MPI_Bcast(vel , ucell.nat*3,MPI_DOUBLE,0,MPI_COMM_WORLD); #endif - ucell.update_pos_tau(pos); - ucell.periodic_boundary_adjustment(); - MD_func::InitPos(ucell, pos); + ucell.update_pos_taud(pos); ModuleBase::timer::tick("FIRE", "first_half"); } diff --git a/source/module_md/Langevin.cpp b/source/module_md/Langevin.cpp index 5dd36e8de4..5d473c93b5 100644 --- a/source/module_md/Langevin.cpp +++ b/source/module_md/Langevin.cpp @@ -42,9 +42,14 @@ void Langevin::first_half() if(ionmbl[i][k]) { vel[i][k] += 0.5 * (force[i][k] + fictitious_force[i][k]) * mdp.md_dt / allmass[i]; - pos[i][k] += vel[i][k] * mdp.md_dt; + pos[i][k] = vel[i][k] * mdp.md_dt / ucell.lat0; + } + else + { + pos[i][k] = 0; } } + pos[i] = pos[i] * ucell.GT; } } @@ -53,9 +58,7 @@ void Langevin::first_half() MPI_Bcast(vel , ucell.nat*3,MPI_DOUBLE,0,MPI_COMM_WORLD); #endif - ucell.update_pos_tau(pos); - ucell.periodic_boundary_adjustment(); - MD_func::InitPos(ucell, pos); + ucell.update_pos_taud(pos); ModuleBase::timer::tick("Langevin", "first_half"); } diff --git a/source/module_md/MD_func.cpp b/source/module_md/MD_func.cpp index e113b7b634..909f95ac70 100644 --- a/source/module_md/MD_func.cpp +++ b/source/module_md/MD_func.cpp @@ -191,21 +191,6 @@ void MD_func::InitVel( std::cout << "--------------------------------- INITVEL DONE ------------------------------------" << std::endl; } -void MD_func::InitPos( - const UnitCell &unit_in, - ModuleBase::Vector3* pos) -{ - int ion=0; - for(int it=0;it* pos); static void InitVel( const UnitCell &unit_in, diff --git a/source/module_md/MSST.cpp b/source/module_md/MSST.cpp index 2a156d0a59..c6af464f94 100644 --- a/source/module_md/MSST.cpp +++ b/source/module_md/MSST.cpp @@ -114,15 +114,25 @@ void MSST::first_half() // propagate atom positions 1 time step for(int i=0; i *pos; // atom position + ModuleBase::Vector3 *pos; // atom displacements liuyu modify 2023-03-22 ModuleBase::Vector3 *vel; // atom velocity ModuleBase::Vector3 *ionmbl; // atom is frozen or not ModuleBase::Vector3 *force; // force of each atom diff --git a/source/module_md/test/MD_func_test.cpp b/source/module_md/test/MD_func_test.cpp index a01817cd63..a988fdf8f6 100644 --- a/source/module_md/test/MD_func_test.cpp +++ b/source/module_md/test/MD_func_test.cpp @@ -15,9 +15,6 @@ * - MD_func::gaussrand * - genarate Gaussian random number * - * - MD_func::InitPos - * - initialize atomic positions - * * - MD_func::RandomVel * - initialize atomic velocity randomly * @@ -84,24 +81,6 @@ TEST_F(MD_func_test, gaussrand) EXPECT_DOUBLE_EQ(MD_func::gaussrand(), 0.60805637857480721); } -TEST_F(MD_func_test, initpos) -{ - MD_func::InitPos(ucell, pos); - - EXPECT_DOUBLE_EQ(pos[0].x, 0.0); - EXPECT_DOUBLE_EQ(pos[0].y, 0.0); - EXPECT_DOUBLE_EQ(pos[0].z, 0.0); - EXPECT_DOUBLE_EQ(pos[1].x, 5.2); - EXPECT_DOUBLE_EQ(pos[1].y, 5.2); - EXPECT_DOUBLE_EQ(pos[1].z, 0.0); - EXPECT_DOUBLE_EQ(pos[2].x, 5.1); - EXPECT_DOUBLE_EQ(pos[2].y, 0.0); - EXPECT_DOUBLE_EQ(pos[2].z, 5.0); - EXPECT_DOUBLE_EQ(pos[3].x, 0.0); - EXPECT_DOUBLE_EQ(pos[3].y, 5.3); - EXPECT_DOUBLE_EQ(pos[3].z, 5.0); -} - TEST_F(MD_func_test, randomvel) { ucell.init_vel = 0; @@ -174,7 +153,6 @@ TEST_F(MD_func_test, compute_stress) TEST_F(MD_func_test, MDdump) { - MD_func::InitPos(ucell, pos); MD_func::MDdump(0, ucell, INPUT, virial, force, vel); std::ifstream ifs("MD_dump"); std::string output_str; diff --git a/source/module_relax/relax_old/bfgs_basic.cpp b/source/module_relax/relax_old/bfgs_basic.cpp index 086c393df8..0e58d3ab0e 100644 --- a/source/module_relax/relax_old/bfgs_basic.cpp +++ b/source/module_relax/relax_old/bfgs_basic.cpp @@ -8,8 +8,6 @@ double BFGS_Basic::relax_bfgs_w2 = -1.0; // defalut is 0.05 BFGS_Basic::BFGS_Basic() { - pos = nullptr; - pos_p = nullptr; grad = nullptr; grad_p = nullptr; move = nullptr; @@ -20,8 +18,6 @@ BFGS_Basic::BFGS_Basic() BFGS_Basic::~BFGS_Basic() { - delete[] pos; - delete[] pos_p; delete[] grad; delete[] grad_p; delete[] move; @@ -32,23 +28,17 @@ void BFGS_Basic::allocate_basic(void) { assert(dim>0); - delete[] pos; - delete[] pos_p; delete[] grad; delete[] grad_p; delete[] move; delete[] move_p; - pos = new double[dim]; - pos_p = new double [dim]; grad = new double[dim]; grad_p = new double [dim]; move = new double [dim]; move_p = new double [dim]; - ModuleBase::GlobalFunc::ZEROS(pos, dim); ModuleBase::GlobalFunc::ZEROS(grad, dim); - ModuleBase::GlobalFunc::ZEROS(pos_p, dim); ModuleBase::GlobalFunc::ZEROS(grad_p, dim); ModuleBase::GlobalFunc::ZEROS(move, dim); ModuleBase::GlobalFunc::ZEROS(move_p, dim); @@ -71,9 +61,8 @@ void BFGS_Basic::update_inverse_hessian(void) for(int i=0;ipos[i] - this->pos_p[i]; // mohan update 2010-07-27 - s[i] = this->check_move( pos[i], pos_p[i] ); + s[i] = this->check_move( move_p[i] ); s[i] *= GlobalC::ucell.lat0; y[i] = this->grad[i] - this->grad_p[i]; @@ -201,7 +190,6 @@ void BFGS_Basic::save_bfgs(void) this->save_flag = true; for(int i=0;ipos_p[i] = this->pos[i]; this->grad_p[i] = this->grad[i]; this->move_p[i] = this->move[i]; } @@ -390,12 +378,12 @@ void BFGS_Basic::compute_trust_radius(void) return; } -double BFGS_Basic::check_move(const double &pos, const double &pos_p) +double BFGS_Basic::check_move(const double &move_p) { // this must be careful. // unit is GlobalC::ucell.lat0. assert(GlobalC::ucell.lat0>0.0); - const double direct_move = (pos - pos_p)/GlobalC::ucell.lat0; + const double direct_move = move_p / GlobalC::ucell.lat0; double shortest_move = direct_move; for(int cell=-1; cell<=1; ++cell) { diff --git a/source/module_relax/relax_old/bfgs_basic.h b/source/module_relax/relax_old/bfgs_basic.h index 6a9207245a..478fd822fc 100644 --- a/source/module_relax/relax_old/bfgs_basic.h +++ b/source/module_relax/relax_old/bfgs_basic.h @@ -29,11 +29,9 @@ class BFGS_Basic void reset_hessian(void); void save_bfgs(void); - double* pos; // std::vector containing 3N coordinates of the system ( x ) - double* grad; //std::vector containing 3N components of ( grad( V(x) ) ) - double* move; // pos = pos_p + move. + double* grad; // std::vector containing 3N components of ( grad( V(x) ) ) + double* move; // movements - double* pos_p; // p: previous double* grad_p; // p: previous double* move_p; @@ -49,7 +47,7 @@ class BFGS_Basic // to the minimum value at the previous step // mohan add 2010-07-27 - double check_move(const double &pos, const double &pos_p); + double check_move(const double &move_p); private: bool wolfe_flag; diff --git a/source/module_relax/relax_old/ions_move_basic.cpp b/source/module_relax/relax_old/ions_move_basic.cpp index a9f7edb0ea..0970d8e4a5 100644 --- a/source/module_relax/relax_old/ions_move_basic.cpp +++ b/source/module_relax/relax_old/ions_move_basic.cpp @@ -20,23 +20,19 @@ double Ions_Move_Basic::best_xxx = 1.0; int Ions_Move_Basic::out_stru=0; -void Ions_Move_Basic::setup_gradient(double* pos, double *grad, const ModuleBase::matrix &force) +void Ions_Move_Basic::setup_gradient(double *grad, const ModuleBase::matrix &force) { ModuleBase::TITLE("Ions_Move_Basic","setup_gradient"); assert(GlobalC::ucell.ntype>0); - assert(pos!=NULL); assert(grad!=NULL); assert(dim == 3*GlobalC::ucell.nat); - ModuleBase::GlobalFunc::ZEROS(pos, dim); ModuleBase::GlobalFunc::ZEROS(grad, dim); // (1) init gradient - // the unit of pos: Bohr. // the unit of force: Ry/Bohr. // the unit of gradient: - GlobalC::ucell.save_cartesian_position(pos); int iat=0; for(int it = 0;it < GlobalC::ucell.ntype;it++) { @@ -64,12 +60,11 @@ void Ions_Move_Basic::setup_gradient(double* pos, double *grad, const ModuleBase return; } -void Ions_Move_Basic::move_atoms(double *move, double *pos) +void Ions_Move_Basic::move_atoms(double *move) { ModuleBase::TITLE("Ions_Move_Basic","move_atoms"); assert(move!=NULL); - assert(pos!=NULL); //------------------------ // for test only @@ -97,15 +92,23 @@ void Ions_Move_Basic::move_atoms(double *move, double *pos) } const double move_threshold = 1.0e-10; - const int total_freedom = GlobalC::ucell.nat * 3; - for(int i =0;i move_threshold ) - { - pos[i] += move[i]; - } - } - GlobalC::ucell.update_pos_tau(pos); + ModuleBase::Vector3 *move_ion = new ModuleBase::Vector3 [GlobalC::ucell.nat]; + for (int i = 0; i < GlobalC::ucell.nat; ++i) + { + for (int k = 0; k < 3; ++k) + { + if( abs(move[3*i + k]) > move_threshold ) + { + move_ion[i][k] = move[3*i + k] / GlobalC::ucell.lat0; + } + else + { + move_ion[i][k] = 0; + } + } + move_ion[i] = move_ion[i] * GlobalC::ucell.GT; + } + GlobalC::ucell.update_pos_taud(move_ion); GlobalC::ucell.periodic_boundary_adjustment(); @@ -122,6 +125,7 @@ void Ions_Move_Basic::move_atoms(double *move, double *pos) { GlobalC::ucell.print_cell_cif("STRU_NOW.cif"); } + delete[] move_ion; return; } diff --git a/source/module_relax/relax_old/ions_move_basic.h b/source/module_relax/relax_old/ions_move_basic.h index 32cbade730..2cca011cd3 100644 --- a/source/module_relax/relax_old/ions_move_basic.h +++ b/source/module_relax/relax_old/ions_move_basic.h @@ -26,12 +26,12 @@ namespace Ions_Move_Basic //---------------------------------------------------------------------------- // setup the gradient, all the same for any geometry optimization methods. //---------------------------------------------------------------------------- - void setup_gradient(double *pos, double* grad, const ModuleBase::matrix &force); + void setup_gradient(double* grad, const ModuleBase::matrix &force); //---------------------------------------------------------------------------- // move the atom positions, considering the periodic boundary condition. //---------------------------------------------------------------------------- - void move_atoms(double* move, double *pos); + void move_atoms(double* move); //---------------------------------------------------------------------------- // check the converged conditions ( if largest gradient is smaller than diff --git a/source/module_relax/relax_old/ions_move_bfgs.cpp b/source/module_relax/relax_old/ions_move_bfgs.cpp index 43031e1e68..cd11ea760f 100644 --- a/source/module_relax/relax_old/ions_move_bfgs.cpp +++ b/source/module_relax/relax_old/ions_move_bfgs.cpp @@ -40,7 +40,7 @@ void Ions_Move_BFGS::start(const ModuleBase::matrix& force, const double &energy // istep must be set eariler. // use force to setup gradient. - Ions_Move_Basic::setup_gradient(this->pos, this->grad, force); + Ions_Move_Basic::setup_gradient(this->grad, force); // use energy_in and istep to setup etot and etot_old. Ions_Move_Basic::setup_etot(energy_in, 0); // use gradient and etot and etot_old to check @@ -72,7 +72,7 @@ void Ions_Move_BFGS::start(const ModuleBase::matrix& force, const double &energy // even if the energy is higher, we save the information. this->save_bfgs(); - Ions_Move_Basic::move_atoms(move, pos); + Ions_Move_Basic::move_atoms(move); } return; } @@ -109,14 +109,13 @@ void Ions_Move_BFGS::restart_bfgs(void) { // mohan add 2010-07-26. // there must be one of the two has the correct sign and value. - this->move_p[i] = this->check_move(pos[i], pos_p[i])/trust_radius_old; + this->move_p[i] = this->check_move(move_p[i])/trust_radius_old; //std::cout << " " << std::setw(20) << move_p[i] << std::setw(20) << dpmin << std::endl; } } else { // bfgs initialization - ModuleBase::GlobalFunc::ZEROS(pos_p, dim); ModuleBase::GlobalFunc::ZEROS(grad_p, dim); ModuleBase::GlobalFunc::ZEROS(move_p, dim); @@ -234,7 +233,6 @@ void Ions_Move_BFGS::bfgs_routine(void) etot = etot_p; for(int i=0;ipos[i] = pos_p[i]; this->grad[i] = grad_p[i]; } diff --git a/source/module_relax/relax_old/ions_move_cg.cpp b/source/module_relax/relax_old/ions_move_cg.cpp index 7cf0bd7271..efcb5483a7 100644 --- a/source/module_relax/relax_old/ions_move_cg.cpp +++ b/source/module_relax/relax_old/ions_move_cg.cpp @@ -23,7 +23,6 @@ double Ions_Move_CG::RELAX_CG_THR =-1.0; //default is 0.5 Ions_Move_CG::Ions_Move_CG() { - this->pos0 = nullptr; this->grad0 = nullptr; this->cg_grad0 = nullptr; this->move0 = nullptr; @@ -31,7 +30,6 @@ Ions_Move_CG::Ions_Move_CG() Ions_Move_CG::~Ions_Move_CG() { - delete[] pos0; delete[] grad0; delete[] cg_grad0; delete[] move0; @@ -41,15 +39,12 @@ void Ions_Move_CG::allocate(void) { ModuleBase::TITLE("Ions_Move_CG","allocate"); assert( dim > 0); - delete[] pos0; delete[] grad0; delete[] cg_grad0; delete[] move0; - this->pos0 = new double[dim]; this->grad0 = new double[dim]; this->cg_grad0 = new double[dim]; this->move0 = new double[dim]; - ModuleBase::GlobalFunc::ZEROS(pos0, dim); ModuleBase::GlobalFunc::ZEROS(grad0, dim); ModuleBase::GlobalFunc::ZEROS(cg_grad0, dim); ModuleBase::GlobalFunc::ZEROS(move0, dim); @@ -60,7 +55,6 @@ void Ions_Move_CG::start(const ModuleBase::matrix& force, const double& etot_in) { ModuleBase::TITLE("Ions_Move_CG","start"); assert(dim>0); - assert(pos0!=0); assert(grad0!=0); assert(cg_grad0!=0); assert(move0!=0); @@ -72,7 +66,6 @@ void Ions_Move_CG::start(const ModuleBase::matrix& force, const double& etot_in) static double xa,xb,xc,xpt,steplength,fmax; // the steepest descent method static int nbrent; - double* pos = new double[dim]; double* grad = new double[dim]; double* cg_gradn = new double[dim]; double* move =new double[dim]; @@ -82,7 +75,6 @@ void Ions_Move_CG::start(const ModuleBase::matrix& force, const double& etot_in) int flag = 0; - ModuleBase::GlobalFunc::ZEROS(pos, dim); ModuleBase::GlobalFunc::ZEROS(grad, dim); ModuleBase::GlobalFunc::ZEROS(move, dim); ModuleBase::GlobalFunc::ZEROS(cg_grad, dim); @@ -101,7 +93,7 @@ void Ions_Move_CG::start(const ModuleBase::matrix& force, const double& etot_in) nbrent = 0; } - Ions_Move_Basic::setup_gradient(pos, grad, force); + Ions_Move_Basic::setup_gradient(grad, force); // use energy_in and istep to setup etot and etot_old. Ions_Move_Basic::setup_etot(etot_in, 0); // use gradient and etot and etot_old to check @@ -142,7 +134,7 @@ void Ions_Move_CG::start(const ModuleBase::matrix& force, const double& etot_in) } setup_move(move0, cg_gradn, steplength); // move the atom position - Ions_Move_Basic::move_atoms(move0, pos); + Ions_Move_Basic::move_atoms(move0); for (int i=0;ienergy_saved = 1.0e10; this->grad_saved = nullptr; - this->pos_saved = nullptr; } Ions_Move_SD::~Ions_Move_SD() { delete[] grad_saved; - delete[] pos_saved; } void Ions_Move_SD::allocate(void) @@ -21,9 +19,7 @@ void Ions_Move_SD::allocate(void) ModuleBase::TITLE("Ions_Move_SD","allocate"); assert( dim > 0); delete[] grad_saved; - delete[] pos_saved; this->grad_saved = new double[dim]; - this->pos_saved = new double[dim]; } void Ions_Move_SD::start(const ModuleBase::matrix& force, const double& etot_in) @@ -32,12 +28,9 @@ void Ions_Move_SD::start(const ModuleBase::matrix& force, const double& etot_in) assert(dim>0); assert(grad_saved!=0); - assert(pos_saved!=0); - double* pos = new double[dim]; double* grad = new double[dim]; double* move =new double[dim]; - ModuleBase::GlobalFunc::ZEROS(pos, dim); ModuleBase::GlobalFunc::ZEROS(grad, dim); ModuleBase::GlobalFunc::ZEROS(move, dim); @@ -45,12 +38,11 @@ void Ions_Move_SD::start(const ModuleBase::matrix& force, const double& etot_in) // 0: ediff < 0 bool judgement = 0; setup_etot(etot_in, judgement); - setup_gradient(pos, grad, force); + setup_gradient(grad, force); if(istep==1 || etot_in <= energy_saved) { energy_saved = etot_in; - for(int i=0; i Date: Wed, 22 Mar 2023 13:48:11 +0800 Subject: [PATCH 10/17] Refactor: update first_half and second_half --- source/module_cell/unitcell.cpp | 88 -------------------------------- source/module_cell/unitcell.h | 4 -- source/module_md/FIRE.cpp | 39 ++------------ source/module_md/Langevin.cpp | 40 ++------------- source/module_md/MSST.cpp | 22 +------- source/module_md/Nose_Hoover.cpp | 56 ++------------------ source/module_md/Nose_Hoover.h | 6 --- source/module_md/mdrun.cpp | 55 +++++++++++++------- source/module_md/mdrun.h | 10 ++++ source/module_md/verlet.cpp | 5 +- 10 files changed, 61 insertions(+), 264 deletions(-) diff --git a/source/module_cell/unitcell.cpp b/source/module_cell/unitcell.cpp index 62e6fe3e9c..27ee32bbda 100644 --- a/source/module_cell/unitcell.cpp +++ b/source/module_cell/unitcell.cpp @@ -408,94 +408,6 @@ void UnitCell::bcast_atoms_tau() #endif } -void UnitCell::save_cartesian_position(double* pos) const -{ - int iat = 0; - for (int it = 0; it < this->ntype; it++) - { - Atom* atom = &this->atoms[it]; - for (int ia = 0; ia < atom->na; ia++) - { - pos[3 * iat] = atom->tau[ia].x * this->lat0; - pos[3 * iat + 1] = atom->tau[ia].y * this->lat0; - pos[3 * iat + 2] = atom->tau[ia].z * this->lat0; - iat++; - } - } - assert(iat == this->nat); - return; -} - -void UnitCell::save_cartesian_position(ModuleBase::Vector3* pos) const -{ - int iat = 0; - for (int it = 0; it < this->ntype; ++it) - { - Atom* atom = &this->atoms[it]; - for (int ia = 0; ia < atom->na; ++ia) - { - pos[iat] = atom->tau_original[ia] * this->lat0; - iat++; - } - } - assert(iat == this->nat); - return; -} - -void UnitCell::save_cartesian_position_original(double* pos) const -{ - int iat = 0; - for (int it = 0; it < this->ntype; it++) - { - Atom* atom = &this->atoms[it]; - for (int ia = 0; ia < atom->na; ia++) - { - pos[3 * iat] = atom->tau_original[ia].x * this->lat0; - pos[3 * iat + 1] = atom->tau_original[ia].y * this->lat0; - pos[3 * iat + 2] = atom->tau_original[ia].z * this->lat0; - iat++; - } - } - assert(iat == this->nat); - return; -} - -void UnitCell::save_cartesian_position_original(ModuleBase::Vector3* pos) const -{ - int iat = 0; - for (int it = 0; it < this->ntype; ++it) - { - Atom* atom = &this->atoms[it]; - for (int ia = 0; ia < atom->na; ++ia) - { - pos[iat] = atom->tau_original[ia] * this->lat0; - iat++; - } - } - assert(iat == this->nat); - return; -} - -/* -bool UnitCell::judge_big_cell(void) const -{ - double diameter = 2 * GlobalV::SEARCH_RADIUS; - - double dis_x = this->omega / cross(a1 * lat0, a2 * lat0).norm(); - double dis_y = this->omega / cross(a2 * lat0, a3 * lat0).norm(); - double dis_z = this->omega / cross(a3 * lat0, a1 * lat0).norm(); - - if (dis_x > diameter && dis_y > diameter && dis_z > diameter) - { - return 1; - } - else - { - return 0; - } -} -*/ - #ifndef __CMD void UnitCell::cal_ux() { diff --git a/source/module_cell/unitcell.h b/source/module_cell/unitcell.h index f04c8b28fb..453bab1a99 100644 --- a/source/module_cell/unitcell.h +++ b/source/module_cell/unitcell.h @@ -185,10 +185,6 @@ class UnitCell void update_vel(const ModuleBase::Vector3* vel_in); void periodic_boundary_adjustment(); void bcast_atoms_tau(); - void save_cartesian_position(double* pos)const; - void save_cartesian_position(ModuleBase::Vector3* pos)const; - void save_cartesian_position_original(double* pos)const; - void save_cartesian_position_original(ModuleBase::Vector3* pos)const; bool judge_big_cell(void)const; void update_stress(ModuleBase::matrix &scs); //updates stress diff --git a/source/module_md/FIRE.cpp b/source/module_md/FIRE.cpp index 6876671dc7..906cf6ffde 100644 --- a/source/module_md/FIRE.cpp +++ b/source/module_md/FIRE.cpp @@ -36,43 +36,12 @@ void FIRE::first_half() { ModuleBase::TITLE("FIRE", "first_half"); ModuleBase::timer::tick("FIRE", "first_half"); - if(GlobalV::MY_RANK == 0) - { - for(int i=0; i* force) { - for(int i=0; i* force); + /** * @brief output MD information such as energy, temperature, and pressure * @param ofs determine the output files diff --git a/source/module_md/verlet.cpp b/source/module_md/verlet.cpp index 9598712ed7..4073dc6ce4 100644 --- a/source/module_md/verlet.cpp +++ b/source/module_md/verlet.cpp @@ -21,7 +21,8 @@ void Verlet::first_half() ModuleBase::TITLE("Verlet", "first_half"); ModuleBase::timer::tick("Verlet", "first_half"); - MDrun::first_half(); + MDrun::update_vel(force); + MDrun::update_pos(); ModuleBase::timer::tick("Verlet", "first_half"); } @@ -31,7 +32,7 @@ void Verlet::second_half() ModuleBase::TITLE("Verlet", "second_half"); ModuleBase::timer::tick("Verlet", "second_half"); - MDrun::second_half(); + MDrun::update_vel(force); apply_thermostat(); ModuleBase::timer::tick("Verlet", "second_half"); From 12f224d394514fd90c1bb87369b7c1dbdedda0df Mon Sep 17 00:00:00 2001 From: YuLiu98 Date: Wed, 22 Mar 2023 14:10:07 +0800 Subject: [PATCH 11/17] Fix: fix bug in update_vel() --- source/module_md/Langevin.cpp | 19 ++++++++++--------- source/module_md/Langevin.h | 2 +- source/module_md/mdrun.cpp | 2 +- 3 files changed, 12 insertions(+), 11 deletions(-) diff --git a/source/module_md/Langevin.cpp b/source/module_md/Langevin.cpp index c11833118a..490ce05b55 100644 --- a/source/module_md/Langevin.cpp +++ b/source/module_md/Langevin.cpp @@ -8,12 +8,12 @@ Langevin::Langevin(MD_parameters& MD_para_in, UnitCell &unit_in) : MDrun(MD_para // convert to a.u. unit mdp.md_damp /= ModuleBase::AU_to_FS; - fictitious_force = new ModuleBase::Vector3 [ucell.nat]; + total_force = new ModuleBase::Vector3 [ucell.nat]; } Langevin::~Langevin() { - delete []fictitious_force; + delete []total_force; } void Langevin::setup(ModuleESolver::ESolver *p_ensolve) @@ -33,7 +33,7 @@ void Langevin::first_half() ModuleBase::TITLE("Langevin", "first_half"); ModuleBase::timer::tick("Langevin", "first_half"); - MDrun::update_vel(force + fictitious_force); + MDrun::update_vel(total_force); MDrun::update_pos(); ModuleBase::timer::tick("Langevin", "first_half"); @@ -45,7 +45,7 @@ void Langevin::second_half() ModuleBase::timer::tick("Langevin", "second_half"); post_force(); - MDrun::update_vel(force + fictitious_force); + MDrun::update_vel(total_force); ModuleBase::timer::tick("Langevin", "second_half"); } @@ -67,21 +67,22 @@ void Langevin::restart() void Langevin::post_force() { - double t_target = MD_func::target_temp(step_ + step_rst_, mdp.md_tfirst, mdp.md_tlast); - if(GlobalV::MY_RANK==0) { + double t_target = MD_func::target_temp(step_ + step_rst_, mdp.md_tfirst, mdp.md_tlast); + ModuleBase::Vector3 fictitious_force; for(int i=0; i *fictitious_force; // Langevin fictitious_force + ModuleBase::Vector3 *total_force; // total force = true force + Langevin fictitious_force }; diff --git a/source/module_md/mdrun.cpp b/source/module_md/mdrun.cpp index e8ac4e2030..a4a14a4929 100644 --- a/source/module_md/mdrun.cpp +++ b/source/module_md/mdrun.cpp @@ -100,7 +100,7 @@ void MDrun::update_pos() ucell.update_pos_taud(pos); } -void Nose_Hoover::update_vel(const ModuleBase::Vector3* force) +void MDrun::update_vel(const ModuleBase::Vector3* force) { if(GlobalV::MY_RANK == 0) { From a2d25ee560bd2e8ab7af2efba95a36574e123171 Mon Sep 17 00:00:00 2001 From: YuLiu98 Date: Wed, 22 Mar 2023 14:50:51 +0800 Subject: [PATCH 12/17] Test: update unit tests --- source/module_cell/test/unitcell_test.cpp | 119 ------------------ source/module_md/test/FIRE_test.cpp | 48 ++++---- source/module_md/test/Langevin_test.cpp | 48 ++++---- source/module_md/test/MSST_test.cpp | 48 ++++---- source/module_md/test/NHC_test.cpp | 48 ++++---- source/module_md/test/verlet_test.cpp | 144 +++++++++++----------- 6 files changed, 168 insertions(+), 287 deletions(-) diff --git a/source/module_cell/test/unitcell_test.cpp b/source/module_cell/test/unitcell_test.cpp index 4a94a4b718..6365d9a324 100644 --- a/source/module_cell/test/unitcell_test.cpp +++ b/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* vel_in) - * - UpdatePosTau1 - * - update_pos_tau(const double* pos) - * - UpdatePosTau2 - * - update_pos_tau(const ModuleBase::Vector3* 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* pos1 = new ModuleBase::Vector3[ucell->nat]; - ucell->save_cartesian_position(pos1); - int iat = 0; - for(int it=0; itatoms[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* pos1 = new ModuleBase::Vector3[ucell->nat]; - ucell->save_cartesian_position_original(pos1); - int iat = 0; - for(int it=0; itatoms[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; iatnat; ++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; iatnat; ++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* pos_in = new ModuleBase::Vector3[ucell->nat]; - ucell->set_iat2itia(); - for(int iat=0; iatnat; ++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; iatnat; ++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"]; diff --git a/source/module_md/test/FIRE_test.cpp b/source/module_md/test/FIRE_test.cpp index 02e876d2a1..16f7493950 100644 --- a/source/module_md/test/FIRE_test.cpp +++ b/source/module_md/test/FIRE_test.cpp @@ -73,18 +73,18 @@ TEST_F(FIRE_test, first_half) { mdrun->first_half(); - EXPECT_NEAR(mdrun->pos[0].x, 9.995455294044568, doublethreshold); - EXPECT_NEAR(mdrun->pos[0].y, 0.003264683323249327, doublethreshold); - EXPECT_NEAR(mdrun->pos[0].z, 9.9994784290476932, doublethreshold); - EXPECT_NEAR(mdrun->pos[1].x, 5.2052136746814073, doublethreshold); - EXPECT_NEAR(mdrun->pos[1].y, 5.19405131115556, doublethreshold); - EXPECT_NEAR(mdrun->pos[1].z, 0.0035886062145122004, doublethreshold); - EXPECT_NEAR(mdrun->pos[2].x, 5.0947593079696478, doublethreshold); - EXPECT_NEAR(mdrun->pos[2].y, 0.00048706739346586134, doublethreshold); - EXPECT_NEAR(mdrun->pos[2].z, 4.9979870945593055, doublethreshold); - EXPECT_NEAR(mdrun->pos[3].x, 0.0045717233044145923, doublethreshold); - EXPECT_NEAR(mdrun->pos[3].y, 5.3021969381277305, doublethreshold); - EXPECT_NEAR(mdrun->pos[3].z, 4.9989458701784848, doublethreshold); + EXPECT_NEAR(mdrun->pos[0].x, -0.00045447059554315662, doublethreshold); + EXPECT_NEAR(mdrun->pos[0].y, 0.00032646833232493271, doublethreshold); + EXPECT_NEAR(mdrun->pos[0].z, -5.215709523063016e-05, doublethreshold); + EXPECT_NEAR(mdrun->pos[1].x, 0.0005213674681407162, doublethreshold); + EXPECT_NEAR(mdrun->pos[1].y, -0.00059486888444406608, doublethreshold); + EXPECT_NEAR(mdrun->pos[1].z, 0.00035886062145122004, doublethreshold); + EXPECT_NEAR(mdrun->pos[2].x, -0.00052406920303529794, doublethreshold); + EXPECT_NEAR(mdrun->pos[2].y, 4.8706739346586155e-05, doublethreshold); + EXPECT_NEAR(mdrun->pos[2].z, -0.00020129054406946794, doublethreshold); + EXPECT_NEAR(mdrun->pos[3].x, 0.00045717233044145918, doublethreshold); + EXPECT_NEAR(mdrun->pos[3].y, 0.00021969381277291936, doublethreshold); + EXPECT_NEAR(mdrun->pos[3].z, -0.00010541298215149392, doublethreshold); EXPECT_NEAR(mdrun->vel[0].x, -0.00010993118004167345, doublethreshold); EXPECT_NEAR(mdrun->vel[0].y, 7.8968913216100539e-05, doublethreshold); @@ -105,18 +105,18 @@ TEST_F(FIRE_test, second_half) mdrun->first_half(); mdrun->second_half(); - EXPECT_NEAR(mdrun->pos[0].x, 9.995455294044568, doublethreshold); - EXPECT_NEAR(mdrun->pos[0].y, 0.003264683323249327, doublethreshold); - EXPECT_NEAR(mdrun->pos[0].z, 9.9994784290476932, doublethreshold); - EXPECT_NEAR(mdrun->pos[1].x, 5.2052136746814073, doublethreshold); - EXPECT_NEAR(mdrun->pos[1].y, 5.19405131115556, doublethreshold); - EXPECT_NEAR(mdrun->pos[1].z, 0.0035886062145122004, doublethreshold); - EXPECT_NEAR(mdrun->pos[2].x, 5.0947593079696478, doublethreshold); - EXPECT_NEAR(mdrun->pos[2].y, 0.00048706739346586134, doublethreshold); - EXPECT_NEAR(mdrun->pos[2].z, 4.9979870945593055, doublethreshold); - EXPECT_NEAR(mdrun->pos[3].x, 0.0045717233044145923, doublethreshold); - EXPECT_NEAR(mdrun->pos[3].y, 5.3021969381277305, doublethreshold); - EXPECT_NEAR(mdrun->pos[3].z, 4.9989458701784848, doublethreshold); + EXPECT_NEAR(mdrun->pos[0].x, -0.00045447059554315662, doublethreshold); + EXPECT_NEAR(mdrun->pos[0].y, 0.00032646833232493271, doublethreshold); + EXPECT_NEAR(mdrun->pos[0].z, -5.215709523063016e-05, doublethreshold); + EXPECT_NEAR(mdrun->pos[1].x, 0.0005213674681407162, doublethreshold); + EXPECT_NEAR(mdrun->pos[1].y, -0.00059486888444406608, doublethreshold); + EXPECT_NEAR(mdrun->pos[1].z, 0.00035886062145122004, doublethreshold); + EXPECT_NEAR(mdrun->pos[2].x, -0.00052406920303529794, doublethreshold); + EXPECT_NEAR(mdrun->pos[2].y, 4.8706739346586155e-05, doublethreshold); + EXPECT_NEAR(mdrun->pos[2].z, -0.00020129054406946794, doublethreshold); + EXPECT_NEAR(mdrun->pos[3].x, 0.00045717233044145918, doublethreshold); + EXPECT_NEAR(mdrun->pos[3].y, 0.00021969381277291936, doublethreshold); + EXPECT_NEAR(mdrun->pos[3].z, -0.00010541298215149392, doublethreshold); EXPECT_NEAR(mdrun->vel[0].x, -0.00010978976887416819, doublethreshold); EXPECT_NEAR(mdrun->vel[0].y, 7.9202349471957007e-05, doublethreshold); diff --git a/source/module_md/test/Langevin_test.cpp b/source/module_md/test/Langevin_test.cpp index d6db645e32..c9cd7f06c4 100644 --- a/source/module_md/test/Langevin_test.cpp +++ b/source/module_md/test/Langevin_test.cpp @@ -73,18 +73,18 @@ TEST_F(Langevin_test, first_half) { mdrun->first_half(); - EXPECT_NEAR(mdrun->pos[0].x, 9.9957116654640092, doublethreshold); - EXPECT_NEAR(mdrun->pos[0].y, 0.0016393608896004908, doublethreshold); - EXPECT_NEAR(mdrun->pos[0].z, 0.0049409894499896556, doublethreshold); - EXPECT_NEAR(mdrun->pos[1].x, 5.2079697932877451, doublethreshold); - EXPECT_NEAR(mdrun->pos[1].y, 5.1985329235797453, doublethreshold); - EXPECT_NEAR(mdrun->pos[1].z, 0.0045070523389717319, doublethreshold); - EXPECT_NEAR(mdrun->pos[2].x, 5.0930848914994087, doublethreshold); - EXPECT_NEAR(mdrun->pos[2].y, 0.0011838145470033955, doublethreshold); - EXPECT_NEAR(mdrun->pos[2].z, 4.9932869712840313, doublethreshold); - EXPECT_NEAR(mdrun->pos[3].x, 9.9993642687852322, doublethreshold); - EXPECT_NEAR(mdrun->pos[3].y, 5.2974098983662827, doublethreshold); - EXPECT_NEAR(mdrun->pos[3].z, 5.0029326569457702, doublethreshold); + EXPECT_NEAR(mdrun->pos[0].x, -0.00042883345359910814, doublethreshold); + EXPECT_NEAR(mdrun->pos[0].y, 0.00016393608896004904, doublethreshold); + EXPECT_NEAR(mdrun->pos[0].z, 0.00049409894499896569, doublethreshold); + EXPECT_NEAR(mdrun->pos[1].x, 0.00079697932877452634, doublethreshold); + EXPECT_NEAR(mdrun->pos[1].y, -0.00014670764202547791, doublethreshold); + EXPECT_NEAR(mdrun->pos[1].z, 0.0004507052338971732, doublethreshold); + EXPECT_NEAR(mdrun->pos[2].x, -0.00069151085005912606, doublethreshold); + EXPECT_NEAR(mdrun->pos[2].y, 0.00011838145470033956, doublethreshold); + EXPECT_NEAR(mdrun->pos[2].z, -0.00067130287159685429, doublethreshold); + EXPECT_NEAR(mdrun->pos[3].x, -6.3573121476911994e-05, doublethreshold); + EXPECT_NEAR(mdrun->pos[3].y, -0.00025901016337184232, doublethreshold); + EXPECT_NEAR(mdrun->pos[3].z, 0.00029326569457701463, doublethreshold); EXPECT_NEAR(mdrun->vel[0].x, -0.00010372985195918919, doublethreshold); EXPECT_NEAR(mdrun->vel[0].y, 3.9654243613399205e-05, doublethreshold); @@ -105,18 +105,18 @@ TEST_F(Langevin_test, second_half) mdrun->first_half(); mdrun->second_half(); - EXPECT_NEAR(mdrun->pos[0].x, 9.9933045979909725, doublethreshold); - EXPECT_NEAR(mdrun->pos[0].y, 0.00033862365219131471, doublethreshold); - EXPECT_NEAR(mdrun->pos[0].z, 9.9954281801131337, doublethreshold); - EXPECT_NEAR(mdrun->pos[1].x, 5.1986310958164275, doublethreshold); - EXPECT_NEAR(mdrun->pos[1].y, 5.2027340532086015, doublethreshold); - EXPECT_NEAR(mdrun->pos[1].z, 9.9987348662023798, doublethreshold); - EXPECT_NEAR(mdrun->pos[2].x, 5.0973799076212742, doublethreshold); - EXPECT_NEAR(mdrun->pos[2].y, 0.003868919168865627, doublethreshold); - EXPECT_NEAR(mdrun->pos[2].z, 5.0001845767835942, doublethreshold); - EXPECT_NEAR(mdrun->pos[3].x, 9.9999789728863995, doublethreshold); - EXPECT_NEAR(mdrun->pos[3].y, 5.3031968974372354, doublethreshold); - EXPECT_NEAR(mdrun->pos[3].z, 4.9996952920372832, doublethreshold); + EXPECT_NEAR(mdrun->pos[0].x, -0.00066954020090275205, doublethreshold); + EXPECT_NEAR(mdrun->pos[0].y, 3.3862365219131354e-05, doublethreshold); + EXPECT_NEAR(mdrun->pos[0].z, -0.00045718198868662484, doublethreshold); + EXPECT_NEAR(mdrun->pos[1].x, -0.0001368904183573199, doublethreshold); + EXPECT_NEAR(mdrun->pos[1].y, 0.00027340532086011393, doublethreshold); + EXPECT_NEAR(mdrun->pos[1].z, -0.00012651337976204397, doublethreshold); + EXPECT_NEAR(mdrun->pos[2].x, -0.00026200923787255071, doublethreshold); + EXPECT_NEAR(mdrun->pos[2].y, 0.00038689191688656276, doublethreshold); + EXPECT_NEAR(mdrun->pos[2].z, 1.8457678359430833e-05, doublethreshold); + EXPECT_NEAR(mdrun->pos[3].x, -2.1027113600492346e-06, doublethreshold); + EXPECT_NEAR(mdrun->pos[3].y, 0.000319689743723507, doublethreshold); + EXPECT_NEAR(mdrun->pos[3].z, -3.0470796271690045e-05, doublethreshold); EXPECT_NEAR(mdrun->vel[0].x, -8.2630969616448438e-05, doublethreshold); EXPECT_NEAR(mdrun->vel[0].y, 0.0001366029202159129, doublethreshold); diff --git a/source/module_md/test/MSST_test.cpp b/source/module_md/test/MSST_test.cpp index 77037f4d51..80dcb32a9a 100644 --- a/source/module_md/test/MSST_test.cpp +++ b/source/module_md/test/MSST_test.cpp @@ -98,18 +98,18 @@ TEST_F(MSST_test, first_half) EXPECT_NEAR(ucell.latvec.e33, 9.9959581179144905, doublethreshold); EXPECT_NEAR(ucell.omega, 999.59581179144902, doublethreshold); - EXPECT_NEAR(mdrun->pos[0].x, 9.9945728176928519, doublethreshold); - EXPECT_NEAR(mdrun->pos[0].y, 0.0029442816868202821, doublethreshold); - EXPECT_NEAR(mdrun->pos[0].z, 9.9953814995781549, doublethreshold); - EXPECT_NEAR(mdrun->pos[1].x, 5.2062875654254492, doublethreshold); - EXPECT_NEAR(mdrun->pos[1].y, 5.1939646253791674, doublethreshold); - EXPECT_NEAR(mdrun->pos[1].z, 0.0039673531204680191, doublethreshold); - EXPECT_NEAR(mdrun->pos[2].x, 5.0944925283175158, doublethreshold); - EXPECT_NEAR(mdrun->pos[2].y, 9.9998842371692671, doublethreshold); - EXPECT_NEAR(mdrun->pos[2].z, 4.9957537084439876, doublethreshold); - EXPECT_NEAR(mdrun->pos[3].x, 0.0046470885642230716, doublethreshold); - EXPECT_NEAR(mdrun->pos[3].y, 5.3032068557647492, doublethreshold); - EXPECT_NEAR(mdrun->pos[3].z, 4.9968136746863676, doublethreshold); + EXPECT_NEAR(mdrun->pos[0].x, -0.00054271823071484467, doublethreshold); + EXPECT_NEAR(mdrun->pos[0].y, 0.00029442816868202821, doublethreshold); + EXPECT_NEAR(mdrun->pos[0].z, -5.7685149290774873e-05, doublethreshold); + EXPECT_NEAR(mdrun->pos[1].x, 0.00062875654254500096, doublethreshold); + EXPECT_NEAR(mdrun->pos[1].y, -0.00060353746208327032, doublethreshold); + EXPECT_NEAR(mdrun->pos[1].z, 0.0003968957326219519, doublethreshold); + EXPECT_NEAR(mdrun->pos[2].x, -0.00055074716824834991, doublethreshold); + EXPECT_NEAR(mdrun->pos[2].y, -1.1576283073263842e-05, doublethreshold); + EXPECT_NEAR(mdrun->pos[2].z, -0.00022262503373940808, doublethreshold); + EXPECT_NEAR(mdrun->pos[3].x, 0.00046470885642230719, doublethreshold); + EXPECT_NEAR(mdrun->pos[3].y, 0.00032068557647491725, doublethreshold); + EXPECT_NEAR(mdrun->pos[3].z, -0.00011658554959218052, doublethreshold); EXPECT_NEAR(mdrun->vel[0].x, -0.00013127726219846624, doublethreshold); EXPECT_NEAR(mdrun->vel[0].y, 7.121876825065284e-05, doublethreshold); @@ -143,18 +143,18 @@ TEST_F(MSST_test, second_half) EXPECT_NEAR(ucell.latvec.e33, 9.9959581179144905, doublethreshold); EXPECT_NEAR(ucell.omega, 999.59581179144902, doublethreshold); - EXPECT_NEAR(mdrun->pos[0].x, 9.9945728176928519, doublethreshold); - EXPECT_NEAR(mdrun->pos[0].y, 0.0029442816868202821, doublethreshold); - EXPECT_NEAR(mdrun->pos[0].z, 9.9953814995781549, doublethreshold); - EXPECT_NEAR(mdrun->pos[1].x, 5.2062875654254492, doublethreshold); - EXPECT_NEAR(mdrun->pos[1].y, 5.1939646253791674, doublethreshold); - EXPECT_NEAR(mdrun->pos[1].z, 0.0039673531204680191, doublethreshold); - EXPECT_NEAR(mdrun->pos[2].x, 5.0944925283175158, doublethreshold); - EXPECT_NEAR(mdrun->pos[2].y, 9.9998842371692671, doublethreshold); - EXPECT_NEAR(mdrun->pos[2].z, 4.9957537084439876, doublethreshold); - EXPECT_NEAR(mdrun->pos[3].x, 0.0046470885642230716, doublethreshold); - EXPECT_NEAR(mdrun->pos[3].y, 5.3032068557647492, doublethreshold); - EXPECT_NEAR(mdrun->pos[3].z, 4.9968136746863676, doublethreshold); + EXPECT_NEAR(mdrun->pos[0].x, -0.00054271823071484467, doublethreshold); + EXPECT_NEAR(mdrun->pos[0].y, 0.00029442816868202821, doublethreshold); + EXPECT_NEAR(mdrun->pos[0].z, -5.7685149290774873e-05, doublethreshold); + EXPECT_NEAR(mdrun->pos[1].x, 0.00062875654254500096, doublethreshold); + EXPECT_NEAR(mdrun->pos[1].y, -0.00060353746208327032, doublethreshold); + EXPECT_NEAR(mdrun->pos[1].z, 0.0003968957326219519, doublethreshold); + EXPECT_NEAR(mdrun->pos[2].x, -0.00055074716824834991, doublethreshold); + EXPECT_NEAR(mdrun->pos[2].y, -1.1576283073263842e-05, doublethreshold); + EXPECT_NEAR(mdrun->pos[2].z, -0.00022262503373940808, doublethreshold); + EXPECT_NEAR(mdrun->pos[3].x, 0.00046470885642230719, doublethreshold); + EXPECT_NEAR(mdrun->pos[3].y, 0.00032068557647491725, doublethreshold); + EXPECT_NEAR(mdrun->pos[3].z, -0.00011658554959218052, doublethreshold); EXPECT_NEAR(mdrun->vel[0].x, -0.00013113585103096098, doublethreshold); EXPECT_NEAR(mdrun->vel[0].y, 7.1452204506509308e-05, doublethreshold); diff --git a/source/module_md/test/NHC_test.cpp b/source/module_md/test/NHC_test.cpp index 139aea91b5..bdd98e0f36 100644 --- a/source/module_md/test/NHC_test.cpp +++ b/source/module_md/test/NHC_test.cpp @@ -74,18 +74,18 @@ TEST_F(NHC_test, first_half) { mdrun->first_half(); - EXPECT_NEAR(mdrun->pos[0].x, 20.786480524287398, doublethreshold); - EXPECT_NEAR(mdrun->pos[0].y, -0.17664650841074978, doublethreshold); - EXPECT_NEAR(mdrun->pos[0].z, 12.763557808832209, doublethreshold); - EXPECT_NEAR(mdrun->pos[1].x, 10.417962231096082, doublethreshold); - EXPECT_NEAR(mdrun->pos[1].y, 8.0482704580843034, doublethreshold); - EXPECT_NEAR(mdrun->pos[1].z, 0.0037631015822075175, doublethreshold); - EXPECT_NEAR(mdrun->pos[2].x, 10.326609893752673, doublethreshold); - EXPECT_NEAR(mdrun->pos[2].y, 15.398062039826284, doublethreshold); - EXPECT_NEAR(mdrun->pos[2].z, 6.3799415879158659, doublethreshold); - EXPECT_NEAR(mdrun->pos[3].x, 0.11327789743365901, doublethreshold); - EXPECT_NEAR(mdrun->pos[3].y, 8.1219079731358974, doublethreshold); - EXPECT_NEAR(mdrun->pos[3].z, 6.380946983836691, doublethreshold); + EXPECT_NEAR(mdrun->pos[0].x, -0.00023793471204889866, doublethreshold); + EXPECT_NEAR(mdrun->pos[0].y, 0.00017779705725471447, doublethreshold); + EXPECT_NEAR(mdrun->pos[0].z, -4.2849245001782489e-05, doublethreshold); + EXPECT_NEAR(mdrun->pos[1].x, 0.00024735043800072474, doublethreshold); + EXPECT_NEAR(mdrun->pos[1].y, -0.00035118366823219693, doublethreshold); + EXPECT_NEAR(mdrun->pos[1].z, 0.00029481907729065568, doublethreshold); + EXPECT_NEAR(mdrun->pos[2].x, -0.0002466402027798198, doublethreshold); + EXPECT_NEAR(mdrun->pos[2].y, -1.7332719211775936e-05, doublethreshold); + EXPECT_NEAR(mdrun->pos[2].z, -0.00016536863874868282, doublethreshold); + EXPECT_NEAR(mdrun->pos[3].x, 0.00023722447682995553, doublethreshold); + EXPECT_NEAR(mdrun->pos[3].y, 0.00019071933018949112, doublethreshold); + EXPECT_NEAR(mdrun->pos[3].z, -8.6601193540496082e-05, doublethreshold); EXPECT_NEAR(mdrun->vel[0].x, -8.2616217816642025e-05, doublethreshold); EXPECT_NEAR(mdrun->vel[0].y, 5.3600943884620578e-05, doublethreshold); @@ -106,18 +106,18 @@ TEST_F(NHC_test, second_half) mdrun->first_half(); mdrun->second_half(); - EXPECT_NEAR(mdrun->pos[0].x, 20.786455241418096, doublethreshold); - EXPECT_NEAR(mdrun->pos[0].y, -0.1766462896653096, doublethreshold); - EXPECT_NEAR(mdrun->pos[0].z, 12.763542287437854, doublethreshold); - EXPECT_NEAR(mdrun->pos[1].x, 10.417949570368947, doublethreshold); - EXPECT_NEAR(mdrun->pos[1].y, 8.0482606637809138, doublethreshold); - EXPECT_NEAR(mdrun->pos[1].z, 0.0037631021542052384, doublethreshold); - EXPECT_NEAR(mdrun->pos[2].x, 10.326597330593856, doublethreshold); - EXPECT_NEAR(mdrun->pos[2].y, 15.398043315340637, doublethreshold); - EXPECT_NEAR(mdrun->pos[2].z, 6.3799338269394124, doublethreshold); - EXPECT_NEAR(mdrun->pos[3].x, 0.11327776681410975, doublethreshold); - EXPECT_NEAR(mdrun->pos[3].y, 8.1218981009532172, doublethreshold); - EXPECT_NEAR(mdrun->pos[3].z, 6.3809392230130584, doublethreshold); + EXPECT_NEAR(mdrun->pos[0].x, -0.00023793503786683287, doublethreshold); + EXPECT_NEAR(mdrun->pos[0].y, 0.0001777972998948069, doublethreshold); + EXPECT_NEAR(mdrun->pos[0].z, -4.2849303620229072e-05, doublethreshold); + EXPECT_NEAR(mdrun->pos[1].x, 0.00024735077679297626, doublethreshold); + EXPECT_NEAR(mdrun->pos[1].y, -0.00035118414817912879, doublethreshold); + EXPECT_NEAR(mdrun->pos[1].z, 0.00029481948060782815, doublethreshold); + EXPECT_NEAR(mdrun->pos[2].x, -0.00024664053996740619, doublethreshold); + EXPECT_NEAR(mdrun->pos[2].y, -1.7332743502090304e-05, doublethreshold); + EXPECT_NEAR(mdrun->pos[2].z, -0.00016536886497560272, doublethreshold); + EXPECT_NEAR(mdrun->pos[3].x, 0.00023722480104322458, doublethreshold); + EXPECT_NEAR(mdrun->pos[3].y, 0.00019071959178664489, doublethreshold); + EXPECT_NEAR(mdrun->pos[3].z, -8.6601312012301988e-05, doublethreshold); EXPECT_NEAR(mdrun->vel[0].x, -5.6948727377433429e-05, doublethreshold); EXPECT_NEAR(mdrun->vel[0].y, 4.4909079540334415e-05, doublethreshold); diff --git a/source/module_md/test/verlet_test.cpp b/source/module_md/test/verlet_test.cpp index 733e2253f2..2456f7d46c 100644 --- a/source/module_md/test/verlet_test.cpp +++ b/source/module_md/test/verlet_test.cpp @@ -74,18 +74,18 @@ TEST_F(Verlet_test, first_half) { mdrun->first_half(); - EXPECT_NEAR(mdrun->pos[0].x, 9.9945454470992772, doublethreshold); - EXPECT_NEAR(mdrun->pos[0].y, 0.0029590658162135359, doublethreshold); - EXPECT_NEAR(mdrun->pos[0].z, 9.9994204767196599, doublethreshold); - EXPECT_NEAR(mdrun->pos[1].x, 5.2063192793031234, doublethreshold); - EXPECT_NEAR(mdrun->pos[1].y, 5.1939342598421803, doublethreshold); - EXPECT_NEAR(mdrun->pos[1].z, 0.0039873402383468889, doublethreshold); - EXPECT_NEAR(mdrun->pos[2].x, 5.0944648036273872, doublethreshold); - EXPECT_NEAR(mdrun->pos[2].y, 9.9998836061438716, doublethreshold); - EXPECT_NEAR(mdrun->pos[2].z, 4.997763438399228, doublethreshold); - EXPECT_NEAR(mdrun->pos[3].x, 0.0046704699702541427, doublethreshold); - EXPECT_NEAR(mdrun->pos[3].y, 5.303223068197739, doublethreshold); - EXPECT_NEAR(mdrun->pos[3].z, 4.9988287446427613, doublethreshold); + EXPECT_NEAR(mdrun->pos[0].x, -0.00054545529007222658, doublethreshold); + EXPECT_NEAR(mdrun->pos[0].y, 0.00029590658162135359, doublethreshold); + EXPECT_NEAR(mdrun->pos[0].z, -5.7952328034033513e-05, doublethreshold); + EXPECT_NEAR(mdrun->pos[1].x, 0.00063192793031220879, doublethreshold); + EXPECT_NEAR(mdrun->pos[1].y, -0.00060657401578200095, doublethreshold); + EXPECT_NEAR(mdrun->pos[1].z, 0.00039873402383468892, doublethreshold); + EXPECT_NEAR(mdrun->pos[2].x, -0.00055351963726126224, doublethreshold); + EXPECT_NEAR(mdrun->pos[2].y, -1.1639385612741475e-05, doublethreshold); + EXPECT_NEAR(mdrun->pos[2].z, -0.00022365616007718661, doublethreshold); + EXPECT_NEAR(mdrun->pos[3].x, 0.00046704699702541431, doublethreshold); + EXPECT_NEAR(mdrun->pos[3].y, 0.00032230681977380224, doublethreshold); + EXPECT_NEAR(mdrun->pos[3].z, -0.00011712553572388214, doublethreshold); EXPECT_NEAR(mdrun->vel[0].x, -0.00013193932519649473, doublethreshold); EXPECT_NEAR(mdrun->vel[0].y, 7.1576379239356465e-05, doublethreshold); @@ -106,18 +106,18 @@ TEST_F(Verlet_test, NVE) mdrun->first_half(); mdrun->second_half(); - EXPECT_NEAR(mdrun->pos[0].x, 9.9945454470992772, doublethreshold); - EXPECT_NEAR(mdrun->pos[0].y, 0.0029590658162135359, doublethreshold); - EXPECT_NEAR(mdrun->pos[0].z, 9.9994204767196599, doublethreshold); - EXPECT_NEAR(mdrun->pos[1].x, 5.2063192793031234, doublethreshold); - EXPECT_NEAR(mdrun->pos[1].y, 5.1939342598421803, doublethreshold); - EXPECT_NEAR(mdrun->pos[1].z, 0.0039873402383468889, doublethreshold); - EXPECT_NEAR(mdrun->pos[2].x, 5.0944648036273872, doublethreshold); - EXPECT_NEAR(mdrun->pos[2].y, 9.9998836061438716, doublethreshold); - EXPECT_NEAR(mdrun->pos[2].z, 4.997763438399228, doublethreshold); - EXPECT_NEAR(mdrun->pos[3].x, 0.0046704699702541427, doublethreshold); - EXPECT_NEAR(mdrun->pos[3].y, 5.303223068197739, doublethreshold); - EXPECT_NEAR(mdrun->pos[3].z, 4.9988287446427613, doublethreshold); + EXPECT_NEAR(mdrun->pos[0].x, -0.00054545529007222658, doublethreshold); + EXPECT_NEAR(mdrun->pos[0].y, 0.00029590658162135359, doublethreshold); + EXPECT_NEAR(mdrun->pos[0].z, -5.7952328034033513e-05, doublethreshold); + EXPECT_NEAR(mdrun->pos[1].x, 0.00063192793031220879, doublethreshold); + EXPECT_NEAR(mdrun->pos[1].y, -0.00060657401578200095, doublethreshold); + EXPECT_NEAR(mdrun->pos[1].z, 0.00039873402383468892, doublethreshold); + EXPECT_NEAR(mdrun->pos[2].x, -0.00055351963726126224, doublethreshold); + EXPECT_NEAR(mdrun->pos[2].y, -1.1639385612741475e-05, doublethreshold); + EXPECT_NEAR(mdrun->pos[2].z, -0.00022365616007718661, doublethreshold); + EXPECT_NEAR(mdrun->pos[3].x, 0.00046704699702541431, doublethreshold); + EXPECT_NEAR(mdrun->pos[3].y, 0.00032230681977380224, doublethreshold); + EXPECT_NEAR(mdrun->pos[3].z, -0.00011712553572388214, doublethreshold); EXPECT_NEAR(mdrun->vel[0].x, -0.00013179791402898947, doublethreshold); EXPECT_NEAR(mdrun->vel[0].y, 7.1809815495212933e-05, doublethreshold); @@ -139,18 +139,18 @@ TEST_F(Verlet_test, Anderson) mdrun->mdp.md_thermostat = "anderson"; mdrun->second_half(); - EXPECT_NEAR(mdrun->pos[0].x, 9.9945454470992772, doublethreshold); - EXPECT_NEAR(mdrun->pos[0].y, 0.0029590658162135359, doublethreshold); - EXPECT_NEAR(mdrun->pos[0].z, 9.9994204767196599, doublethreshold); - EXPECT_NEAR(mdrun->pos[1].x, 5.2063192793031234, doublethreshold); - EXPECT_NEAR(mdrun->pos[1].y, 5.1939342598421803, doublethreshold); - EXPECT_NEAR(mdrun->pos[1].z, 0.0039873402383468889, doublethreshold); - EXPECT_NEAR(mdrun->pos[2].x, 5.0944648036273872, doublethreshold); - EXPECT_NEAR(mdrun->pos[2].y, 9.9998836061438716, doublethreshold); - EXPECT_NEAR(mdrun->pos[2].z, 4.997763438399228, doublethreshold); - EXPECT_NEAR(mdrun->pos[3].x, 0.0046704699702541427, doublethreshold); - EXPECT_NEAR(mdrun->pos[3].y, 5.303223068197739, doublethreshold); - EXPECT_NEAR(mdrun->pos[3].z, 4.9988287446427613, doublethreshold); + EXPECT_NEAR(mdrun->pos[0].x, -0.00054545529007222658, doublethreshold); + EXPECT_NEAR(mdrun->pos[0].y, 0.00029590658162135359, doublethreshold); + EXPECT_NEAR(mdrun->pos[0].z, -5.7952328034033513e-05, doublethreshold); + EXPECT_NEAR(mdrun->pos[1].x, 0.00063192793031220879, doublethreshold); + EXPECT_NEAR(mdrun->pos[1].y, -0.00060657401578200095, doublethreshold); + EXPECT_NEAR(mdrun->pos[1].z, 0.00039873402383468892, doublethreshold); + EXPECT_NEAR(mdrun->pos[2].x, -0.00055351963726126224, doublethreshold); + EXPECT_NEAR(mdrun->pos[2].y, -1.1639385612741475e-05, doublethreshold); + EXPECT_NEAR(mdrun->pos[2].z, -0.00022365616007718661, doublethreshold); + EXPECT_NEAR(mdrun->pos[3].x, 0.00046704699702541431, doublethreshold); + EXPECT_NEAR(mdrun->pos[3].y, 0.00032230681977380224, doublethreshold); + EXPECT_NEAR(mdrun->pos[3].z, -0.00011712553572388214, doublethreshold); EXPECT_NEAR(mdrun->vel[0].x, -0.00013179791402898947, doublethreshold); EXPECT_NEAR(mdrun->vel[0].y, 7.1809815495212933e-05, doublethreshold); @@ -172,18 +172,18 @@ TEST_F(Verlet_test, Berendsen) mdrun->mdp.md_thermostat = "berendsen"; mdrun->second_half(); - EXPECT_NEAR(mdrun->pos[0].x, 9.9945454470992772, doublethreshold); - EXPECT_NEAR(mdrun->pos[0].y, 0.0029590658162135359, doublethreshold); - EXPECT_NEAR(mdrun->pos[0].z, 9.9994204767196599, doublethreshold); - EXPECT_NEAR(mdrun->pos[1].x, 5.2063192793031234, doublethreshold); - EXPECT_NEAR(mdrun->pos[1].y, 5.1939342598421803, doublethreshold); - EXPECT_NEAR(mdrun->pos[1].z, 0.0039873402383468889, doublethreshold); - EXPECT_NEAR(mdrun->pos[2].x, 5.0944648036273872, doublethreshold); - EXPECT_NEAR(mdrun->pos[2].y, 9.9998836061438716, doublethreshold); - EXPECT_NEAR(mdrun->pos[2].z, 4.997763438399228, doublethreshold); - EXPECT_NEAR(mdrun->pos[3].x, 0.0046704699702541427, doublethreshold); - EXPECT_NEAR(mdrun->pos[3].y, 5.303223068197739, doublethreshold); - EXPECT_NEAR(mdrun->pos[3].z, 4.9988287446427613, doublethreshold); + EXPECT_NEAR(mdrun->pos[0].x, -0.00054545529007222658, doublethreshold); + EXPECT_NEAR(mdrun->pos[0].y, 0.00029590658162135359, doublethreshold); + EXPECT_NEAR(mdrun->pos[0].z, -5.7952328034033513e-05, doublethreshold); + EXPECT_NEAR(mdrun->pos[1].x, 0.00063192793031220879, doublethreshold); + EXPECT_NEAR(mdrun->pos[1].y, -0.00060657401578200095, doublethreshold); + EXPECT_NEAR(mdrun->pos[1].z, 0.00039873402383468892, doublethreshold); + EXPECT_NEAR(mdrun->pos[2].x, -0.00055351963726126224, doublethreshold); + EXPECT_NEAR(mdrun->pos[2].y, -1.1639385612741475e-05, doublethreshold); + EXPECT_NEAR(mdrun->pos[2].z, -0.00022365616007718661, doublethreshold); + EXPECT_NEAR(mdrun->pos[3].x, 0.00046704699702541431, doublethreshold); + EXPECT_NEAR(mdrun->pos[3].y, 0.00032230681977380224, doublethreshold); + EXPECT_NEAR(mdrun->pos[3].z, -0.00011712553572388214, doublethreshold); EXPECT_NEAR(mdrun->vel[0].x, -0.00013179175250738632, doublethreshold); EXPECT_NEAR(mdrun->vel[0].y, 7.1806458403162173e-05, doublethreshold); @@ -205,18 +205,18 @@ TEST_F(Verlet_test, rescaling) mdrun->mdp.md_thermostat = "rescaling"; mdrun->second_half(); - EXPECT_NEAR(mdrun->pos[0].x, 9.9945454470992772, doublethreshold); - EXPECT_NEAR(mdrun->pos[0].y, 0.0029590658162135359, doublethreshold); - EXPECT_NEAR(mdrun->pos[0].z, 9.9994204767196599, doublethreshold); - EXPECT_NEAR(mdrun->pos[1].x, 5.2063192793031234, doublethreshold); - EXPECT_NEAR(mdrun->pos[1].y, 5.1939342598421803, doublethreshold); - EXPECT_NEAR(mdrun->pos[1].z, 0.0039873402383468889, doublethreshold); - EXPECT_NEAR(mdrun->pos[2].x, 5.0944648036273872, doublethreshold); - EXPECT_NEAR(mdrun->pos[2].y, 9.9998836061438716, doublethreshold); - EXPECT_NEAR(mdrun->pos[2].z, 4.997763438399228, doublethreshold); - EXPECT_NEAR(mdrun->pos[3].x, 0.0046704699702541427, doublethreshold); - EXPECT_NEAR(mdrun->pos[3].y, 5.303223068197739, doublethreshold); - EXPECT_NEAR(mdrun->pos[3].z, 4.9988287446427613, doublethreshold); + EXPECT_NEAR(mdrun->pos[0].x, -0.00054545529007222658, doublethreshold); + EXPECT_NEAR(mdrun->pos[0].y, 0.00029590658162135359, doublethreshold); + EXPECT_NEAR(mdrun->pos[0].z, -5.7952328034033513e-05, doublethreshold); + EXPECT_NEAR(mdrun->pos[1].x, 0.00063192793031220879, doublethreshold); + EXPECT_NEAR(mdrun->pos[1].y, -0.00060657401578200095, doublethreshold); + EXPECT_NEAR(mdrun->pos[1].z, 0.00039873402383468892, doublethreshold); + EXPECT_NEAR(mdrun->pos[2].x, -0.00055351963726126224, doublethreshold); + EXPECT_NEAR(mdrun->pos[2].y, -1.1639385612741475e-05, doublethreshold); + EXPECT_NEAR(mdrun->pos[2].z, -0.00022365616007718661, doublethreshold); + EXPECT_NEAR(mdrun->pos[3].x, 0.00046704699702541431, doublethreshold); + EXPECT_NEAR(mdrun->pos[3].y, 0.00032230681977380224, doublethreshold); + EXPECT_NEAR(mdrun->pos[3].z, -0.00011712553572388214, doublethreshold); EXPECT_NEAR(mdrun->vel[0].x, -0.00013178559069770653, doublethreshold); EXPECT_NEAR(mdrun->vel[0].y, 7.1803101154153484e-05, doublethreshold); @@ -238,18 +238,18 @@ TEST_F(Verlet_test, rescale_v) mdrun->mdp.md_thermostat = "rescale_v"; mdrun->second_half(); - EXPECT_NEAR(mdrun->pos[0].x, 9.9945454470992772, doublethreshold); - EXPECT_NEAR(mdrun->pos[0].y, 0.0029590658162135359, doublethreshold); - EXPECT_NEAR(mdrun->pos[0].z, 9.9994204767196599, doublethreshold); - EXPECT_NEAR(mdrun->pos[1].x, 5.2063192793031234, doublethreshold); - EXPECT_NEAR(mdrun->pos[1].y, 5.1939342598421803, doublethreshold); - EXPECT_NEAR(mdrun->pos[1].z, 0.0039873402383468889, doublethreshold); - EXPECT_NEAR(mdrun->pos[2].x, 5.0944648036273872, doublethreshold); - EXPECT_NEAR(mdrun->pos[2].y, 9.9998836061438716, doublethreshold); - EXPECT_NEAR(mdrun->pos[2].z, 4.997763438399228, doublethreshold); - EXPECT_NEAR(mdrun->pos[3].x, 0.0046704699702541427, doublethreshold); - EXPECT_NEAR(mdrun->pos[3].y, 5.303223068197739, doublethreshold); - EXPECT_NEAR(mdrun->pos[3].z, 4.9988287446427613, doublethreshold); + EXPECT_NEAR(mdrun->pos[0].x, -0.00054545529007222658, doublethreshold); + EXPECT_NEAR(mdrun->pos[0].y, 0.00029590658162135359, doublethreshold); + EXPECT_NEAR(mdrun->pos[0].z, -5.7952328034033513e-05, doublethreshold); + EXPECT_NEAR(mdrun->pos[1].x, 0.00063192793031220879, doublethreshold); + EXPECT_NEAR(mdrun->pos[1].y, -0.00060657401578200095, doublethreshold); + EXPECT_NEAR(mdrun->pos[1].z, 0.00039873402383468892, doublethreshold); + EXPECT_NEAR(mdrun->pos[2].x, -0.00055351963726126224, doublethreshold); + EXPECT_NEAR(mdrun->pos[2].y, -1.1639385612741475e-05, doublethreshold); + EXPECT_NEAR(mdrun->pos[2].z, -0.00022365616007718661, doublethreshold); + EXPECT_NEAR(mdrun->pos[3].x, 0.00046704699702541431, doublethreshold); + EXPECT_NEAR(mdrun->pos[3].y, 0.00032230681977380224, doublethreshold); + EXPECT_NEAR(mdrun->pos[3].z, -0.00011712553572388214, doublethreshold); EXPECT_NEAR(mdrun->vel[0].x, -0.00013178559069770653, doublethreshold); EXPECT_NEAR(mdrun->vel[0].y, 7.1803101154153484e-05, doublethreshold); From 0795e51a13b68bbb1681ce7b41734492313c5099 Mon Sep 17 00:00:00 2001 From: YuLiu98 Date: Wed, 22 Mar 2023 16:28:36 +0800 Subject: [PATCH 13/17] Refactor: refactor second-order method --- source/module_cell/atom_spec.cpp | 8 +-- source/module_cell/atom_spec.h | 2 +- source/module_cell/read_atoms.cpp | 6 +- source/module_cell/test/prepare_unitcell.h | 6 +- source/module_cell/unitcell.cpp | 12 +--- .../module_charge/charge_extra.cpp | 68 +++++++++---------- .../module_charge/charge_extra.h | 11 ++- source/module_esolver/esolver_ks_lcao.cpp | 3 - .../module_esolver/esolver_ks_lcao_elec.cpp | 10 +-- source/module_esolver/esolver_ks_pw.cpp | 11 +-- source/module_esolver/esolver_of.cpp | 6 +- .../module_io/test_serial/prepare_unitcell.h | 6 +- source/module_md/test/setcell.h | 2 +- .../relax_old/ions_move_basic.cpp | 33 +++++---- 14 files changed, 79 insertions(+), 105 deletions(-) diff --git a/source/module_cell/atom_spec.cpp b/source/module_cell/atom_spec.cpp index f6ee0a0029..8a0601950d 100644 --- a/source/module_cell/atom_spec.cpp +++ b/source/module_cell/atom_spec.cpp @@ -13,7 +13,7 @@ Atom::Atom() stapos_wf = 0; mass = 0.0; tau = new ModuleBase::Vector3[1]; - tau_original = new ModuleBase::Vector3[1]; + dis = new ModuleBase::Vector3[1]; taud = new ModuleBase::Vector3[1]; vel = new ModuleBase::Vector3[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[na]; - tau_original = new ModuleBase::Vector3[na]; + dis = new ModuleBase::Vector3[na]; taud = new ModuleBase::Vector3[na]; vel = new ModuleBase::Vector3[na]; mag = new double[na]; diff --git a/source/module_cell/atom_spec.h b/source/module_cell/atom_spec.h index 5b3910f4d9..4a1848a14b 100644 --- a/source/module_cell/atom_spec.h +++ b/source/module_cell/atom_spec.h @@ -36,7 +36,7 @@ class Atom std::string label; // atomic symbol ModuleBase::Vector3 *tau;// Cartesian coordinates of each atom in this type. - ModuleBase::Vector3 *tau_original;// Cartesian coordinates of each atom in this type, but without periodic adjustment. + ModuleBase::Vector3 *dis;// direct displacements of each atom in this type in current step liuyu modift 2023-03-22 ModuleBase::Vector3 *taud;// Direct coordinates of each atom in this type. ModuleBase::Vector3 *vel;// velocities of each atom in this type. ModuleBase::Vector3 *force; // force acting on each atom in this type. diff --git a/source/module_cell/read_atoms.cpp b/source/module_cell/read_atoms.cpp index dc487c40dc..4b1e2f3380 100644 --- a/source/module_cell/read_atoms.cpp +++ b/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[na]; - atoms[it].tau_original = new ModuleBase::Vector3[na]; + atoms[it].dis = new ModuleBase::Vector3[na]; atoms[it].taud = new ModuleBase::Vector3[na]; atoms[it].vel = new ModuleBase::Vector3[na]; atoms[it].mbl = new ModuleBase::Vector3[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 diff --git a/source/module_cell/test/prepare_unitcell.h b/source/module_cell/test/prepare_unitcell.h index 1ddc7aa6f7..11898c820a 100644 --- a/source/module_cell/test/prepare_unitcell.h +++ b/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[ucell->atoms[it].na]; - ucell->atoms[it].tau_original = new ModuleBase::Vector3[ucell->atoms[it].na]; + ucell->atoms[it].dis = new ModuleBase::Vector3[ucell->atoms[it].na]; ucell->atoms[it].taud = new ModuleBase::Vector3[ucell->atoms[it].na]; ucell->atoms[it].vel = new ModuleBase::Vector3[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]; diff --git a/source/module_cell/unitcell.cpp b/source/module_cell/unitcell.cpp index 27ee32bbda..316e56cee0 100644 --- a/source/module_cell/unitcell.cpp +++ b/source/module_cell/unitcell.cpp @@ -293,12 +293,6 @@ void UnitCell::set_iat2itia(void) return; } -// Note : note here we are not keeping track of 'tau_original', namely -// the Cartesian coordinate before periodic adjustment -// The reason is that this is only used in relaxation -// for which we are not using 2nd order extrapolation -// and tau_original is only relevant for 2nd order extrapolation -// and is only meaningful in the context of MD void UnitCell::update_pos_taud(double* posd_in) { int iat = 0; @@ -328,10 +322,8 @@ void UnitCell::update_pos_taud(const ModuleBase::Vector3* posd_in) { for ( int ik = 0; ik < 3; ++ik) { - if (atom->mbl[ia][ik]) - { - atom->taud[ia][ik] += posd_in[ia][ik]; - } + atom->taud[ia][ik] += posd_in[ia][ik]; + atom->dis[ia][ik] = posd_in[ia][ik]; } } } diff --git a/source/module_elecstate/module_charge/charge_extra.cpp b/source/module_elecstate/module_charge/charge_extra.cpp index c8877cb74e..4c4b0f0510 100644 --- a/source/module_elecstate/module_charge/charge_extra.cpp +++ b/source/module_elecstate/module_charge/charge_extra.cpp @@ -21,10 +21,12 @@ Charge_Extra::~Charge_Extra() delete[] delta_rho2; } - delete[] pos_old1; - delete[] pos_old2; - delete[] pos_now; - delete[] pos_next; + if(pot_order == 3) + { + delete[] dis_old1; + delete[] dis_old2; + delete[] dis_now; + } } void Charge_Extra::Init_CE() @@ -66,10 +68,12 @@ void Charge_Extra::Init_CE() natom = GlobalC::ucell.nat; omega_old = GlobalC::ucell.omega; - pos_old1 = new ModuleBase::Vector3[natom]; - pos_old2 = new ModuleBase::Vector3[natom]; - pos_now = new ModuleBase::Vector3[natom]; - pos_next = new ModuleBase::Vector3[natom]; + if(pot_order == 3) + { + dis_old1 = new ModuleBase::Vector3[natom]; + dis_old2 = new ModuleBase::Vector3[natom]; + dis_now = new ModuleBase::Vector3[natom]; + } alpha = 1.0; beta = 0.0; @@ -227,6 +231,8 @@ void Charge_Extra::extrapolate_charge(Charge* chr) } } + omega_old = GlobalC::ucell.omega; + for(int is=0; ispos_next); - return; -} - -void Charge_Extra::update_istep() -{ - this->istep++; - return; -} - -void Charge_Extra::update_all_pos(const UnitCell& ucell) -{ - for(int i=0; ipos_old2[i] = this->pos_old1[i]; - this->pos_old1[i] = this->pos_now[i]; - if(GlobalV::CALCULATION=="relax"||GlobalV::CALCULATION=="cell-relax") + int iat = 0; + for (int it = 0; it < ucell.ntype; it++) { - this->pos_now[i] = this->pos_next[i]; + Atom* atom = &ucell.atoms[it]; + for (int ia = 0; ia < atom->na; ia++) + { + dis_old2[iat] = dis_old1[iat]; + dis_old1[iat] = dis_now[iat]; + dis_now[iat] = atom->dis[ia]; + iat++; + } } + assert(iat == ucell.nat); } - if(GlobalV::CALCULATION=="md") - ucell.save_cartesian_position_original(this->pos_now); return; } diff --git a/source/module_elecstate/module_charge/charge_extra.h b/source/module_elecstate/module_charge/charge_extra.h index 1ec16b5b6e..3236367ed9 100644 --- a/source/module_elecstate/module_charge/charge_extra.h +++ b/source/module_elecstate/module_charge/charge_extra.h @@ -24,9 +24,7 @@ class Charge_Extra //Init_CE will be removed and everything put back in the constructor void Init_CE(); void extrapolate_charge(Charge* chr); - void save_pos_next(const UnitCell& ucell); - void update_istep(); - void update_all_pos(const UnitCell& ucell); + void update_all_dis(const UnitCell& ucell); private: int istep = 0; @@ -36,10 +34,9 @@ class Charge_Extra double omega_old; // the old volume of the last step // for the second-order extrapolation - ModuleBase::Vector3 *pos_old1 = nullptr; - ModuleBase::Vector3 *pos_old2 = nullptr; - ModuleBase::Vector3 *pos_now = nullptr; - ModuleBase::Vector3 *pos_next = nullptr; + ModuleBase::Vector3 *dis_old1 = nullptr; // dis_old2 = pos_old1 - pos_old2 + ModuleBase::Vector3 *dis_old2 = nullptr; // dis_old1 = pos_now - pos_old1 + ModuleBase::Vector3 *dis_now = nullptr; // dis_now = pos_next - pos_now double** delta_rho1 = nullptr; double** delta_rho2 = nullptr; diff --git a/source/module_esolver/esolver_ks_lcao.cpp b/source/module_esolver/esolver_ks_lcao.cpp index 6ab1c8ac20..92f832b4dd 100644 --- a/source/module_esolver/esolver_ks_lcao.cpp +++ b/source/module_esolver/esolver_ks_lcao.cpp @@ -830,9 +830,6 @@ void ESolver_KS_LCAO::eachiterfinish(int iter) void ESolver_KS_LCAO::afterscf(const int istep) { - // Temporary liuyu add 2022-11-07 - CE.update_all_pos(GlobalC::ucell); - if (this->LOC.out_dm1 == 1) { double** dm2d; diff --git a/source/module_esolver/esolver_ks_lcao_elec.cpp b/source/module_esolver/esolver_ks_lcao_elec.cpp index 1936947a96..e78db27de0 100644 --- a/source/module_esolver/esolver_ks_lcao_elec.cpp +++ b/source/module_esolver/esolver_ks_lcao_elec.cpp @@ -268,9 +268,7 @@ namespace ModuleESolver } else { - this->CE.update_istep(); - // this->CE.update_all_pos(GlobalC::ucell); - this->CE.save_pos_next(GlobalC::ucell); + this->CE.update_all_dis(GlobalC::ucell); this->CE.extrapolate_charge(this->pelec->charge); if(GlobalC::ucell.cell_parameter_updated && GlobalV::md_prec_level == 0) { @@ -285,10 +283,8 @@ namespace ModuleESolver { GlobalV::ofs_running << " Setup the extrapolated charge." << std::endl; // charge extrapolation if istep>0. - CE.update_istep(); - CE.update_all_pos(GlobalC::ucell); - CE.extrapolate_charge(pelec->charge); - CE.save_pos_next(GlobalC::ucell); + this->CE.update_all_dis(GlobalC::ucell); + this->CE.extrapolate_charge(this->pelec->charge); GlobalV::ofs_running << " Setup the Vl+Vh+Vxc according to new structure factor and new charge." << std::endl; // calculate the new potential accordint to diff --git a/source/module_esolver/esolver_ks_pw.cpp b/source/module_esolver/esolver_ks_pw.cpp index 0fc21834c7..a5bc31f40a 100644 --- a/source/module_esolver/esolver_ks_pw.cpp +++ b/source/module_esolver/esolver_ks_pw.cpp @@ -241,9 +241,7 @@ namespace ModuleESolver } else { - // this->CE.update_all_pos(GlobalC::ucell); - this->CE.update_istep(); - this->CE.save_pos_next(GlobalC::ucell); + this->CE.update_all_dis(GlobalC::ucell); this->CE.extrapolate_charge(this->pelec->charge); if(GlobalC::ucell.cell_parameter_updated && GlobalV::md_prec_level == 0) { @@ -263,10 +261,8 @@ namespace ModuleESolver { GlobalV::ofs_running << " Setup the extrapolated charge." << std::endl; // charge extrapolation if istep>0. - this->CE.update_istep(); - this->CE.update_all_pos(GlobalC::ucell); + this->CE.update_all_dis(GlobalC::ucell); this->CE.extrapolate_charge(this->pelec->charge); - this->CE.save_pos_next(GlobalC::ucell); GlobalV::ofs_running << " Setup the Vl+Vh+Vxc according to new structure factor and new charge." << std::endl; // calculate the new potential accordint to @@ -551,9 +547,6 @@ namespace ModuleESolver template void ESolver_KS_PW::afterscf(const int istep) { - // Temporary liuyu add 2022-11-07 - this->CE.update_all_pos(GlobalC::ucell); - if (GlobalV::out_pot == 1) // output the effective potential, sunliang 2023-03-16 { for (int is = 0; is < GlobalV::NSPIN; is++) diff --git a/source/module_esolver/esolver_of.cpp b/source/module_esolver/esolver_of.cpp index db9d7e5429..61116ba6f6 100644 --- a/source/module_esolver/esolver_of.cpp +++ b/source/module_esolver/esolver_of.cpp @@ -283,8 +283,7 @@ void ESolver_OF::beforeOpt(const int istep) // Temporary, md and relax will merge later liuyu add 2022-11-07 if(GlobalV::CALCULATION == "md" && istep) { - CE.update_istep(); - CE.save_pos_next(GlobalC::ucell); + CE.update_all_dis(GlobalC::ucell); CE.extrapolate_charge(pelec->charge); if(GlobalC::ucell.cell_parameter_updated) @@ -886,9 +885,6 @@ void ESolver_OF::printInfo() void ESolver_OF::afterOpt() { - // Temporary liuyu add 2022-11-07 - CE.update_all_pos(GlobalC::ucell); - if (this->conv) { GlobalV::ofs_running << "\n charge density convergence is achieved" << std::endl; diff --git a/source/module_io/test_serial/prepare_unitcell.h b/source/module_io/test_serial/prepare_unitcell.h index 14ecfb6a1e..0136324ec3 100644 --- a/source/module_io/test_serial/prepare_unitcell.h +++ b/source/module_io/test_serial/prepare_unitcell.h @@ -148,7 +148,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; @@ -157,7 +157,7 @@ class UcellTestPrepare delete[] ucell->atoms[it].m_loc_; delete[] ucell->atoms[it].mbl; ucell->atoms[it].tau = new ModuleBase::Vector3[ucell->atoms[it].na]; - ucell->atoms[it].tau_original = new ModuleBase::Vector3[ucell->atoms[it].na]; + ucell->atoms[it].dis = new ModuleBase::Vector3[ucell->atoms[it].na]; ucell->atoms[it].taud = new ModuleBase::Vector3[ucell->atoms[it].na]; ucell->atoms[it].vel = new ModuleBase::Vector3[ucell->atoms[it].na]; ucell->atoms[it].mag = new double[ucell->atoms[it].na]; @@ -187,7 +187,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]; diff --git a/source/module_md/test/setcell.h b/source/module_md/test/setcell.h index b8c21a758b..0b375021fa 100644 --- a/source/module_md/test/setcell.h +++ b/source/module_md/test/setcell.h @@ -56,7 +56,7 @@ class Setcell delete[] ucell.atoms[0].vel; delete[] ucell.atoms[0].mbl; ucell.atoms[0].tau = new ModuleBase::Vector3[4]; - ucell.atoms[0].tau_original = new ModuleBase::Vector3[4]; + ucell.atoms[0].dis = new ModuleBase::Vector3[4]; ucell.atoms[0].taud = new ModuleBase::Vector3[4]; ucell.atoms[0].vel = new ModuleBase::Vector3[4]; ucell.atoms[0].mbl = new ModuleBase::Vector3[4]; diff --git a/source/module_relax/relax_old/ions_move_basic.cpp b/source/module_relax/relax_old/ions_move_basic.cpp index 0970d8e4a5..2e83fedd4d 100644 --- a/source/module_relax/relax_old/ions_move_basic.cpp +++ b/source/module_relax/relax_old/ions_move_basic.cpp @@ -91,28 +91,31 @@ void Ions_Move_Basic::move_atoms(double *move) assert( iat == GlobalC::ucell.nat ); } - const double move_threshold = 1.0e-10; + int iat = 0; + const double move_threshold = 1.0e-10; ModuleBase::Vector3 *move_ion = new ModuleBase::Vector3 [GlobalC::ucell.nat]; - for (int i = 0; i < GlobalC::ucell.nat; ++i) + for(int it = 0; it < GlobalC::ucell.ntype; it++) { - for (int k = 0; k < 3; ++k) + Atom* atom = &GlobalC::ucell.atoms[it]; + for(int ia = 0; ia < atom->na; ia++) { - if( abs(move[3*i + k]) > move_threshold ) + for (int k = 0; k < 3; ++k) { - move_ion[i][k] = move[3*i + k] / GlobalC::ucell.lat0; - } - else - { - move_ion[i][k] = 0; + if( abs(move[3*iat + k]) > move_threshold && atom->mbl[ia][k]) + { + move_ion[iat][k] = move[3*iat + k] / GlobalC::ucell.lat0; + } + else + { + move_ion[iat][k] = 0; + } } + move_ion[iat] = move_ion[iat] * GlobalC::ucell.GT; + iat++; } - move_ion[i] = move_ion[i] * GlobalC::ucell.GT; } - GlobalC::ucell.update_pos_taud(move_ion); - - GlobalC::ucell.periodic_boundary_adjustment(); - - GlobalC::ucell.bcast_atoms_tau(); + assert( iat == GlobalC::ucell.nat ); + GlobalC::ucell.update_pos_taud(move_ion); //-------------------------------------------- // Print out the structure file. From 0ab74a15bed9ab20d273b9a1d6603707bb1536e9 Mon Sep 17 00:00:00 2001 From: YuLiu98 Date: Wed, 22 Mar 2023 22:19:08 +0800 Subject: [PATCH 14/17] Fix: fix nan when t_first = 0 --- source/module_elecstate/module_charge/charge_extra.cpp | 4 ++-- source/module_md/MD_func.cpp | 2 +- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/source/module_elecstate/module_charge/charge_extra.cpp b/source/module_elecstate/module_charge/charge_extra.cpp index 4c4b0f0510..8f711f48dd 100644 --- a/source/module_elecstate/module_charge/charge_extra.cpp +++ b/source/module_elecstate/module_charge/charge_extra.cpp @@ -273,7 +273,7 @@ void Charge_Extra::find_alpha_and_beta(void) a21 = a12; det = a11 * a22 - a12 * a21; - if(det < -1e-16) + if(det < -1e-20) { alpha = 0.0; beta = 0.0; @@ -281,7 +281,7 @@ void Charge_Extra::find_alpha_and_beta(void) ModuleBase::GlobalFunc::OUT(GlobalV::ofs_warning,"in find_alpha_and beta() det = ", det); } - if(det > 1e-16) + if(det > 1e-20) { alpha = (b1 * a22 - b2 * a12) / det; beta = (a11 * b2 - a21 * b1) / det; diff --git a/source/module_md/MD_func.cpp b/source/module_md/MD_func.cpp index 909f95ac70..e484554f9a 100644 --- a/source/module_md/MD_func.cpp +++ b/source/module_md/MD_func.cpp @@ -144,7 +144,7 @@ void MD_func::RandomVel( } double factor; - if(3*numIon == frozen_freedom) + if(3*numIon == frozen_freedom || temperature == 0) { factor = 0; } From 58e9be1b17ca23c74016a056872abca0fa0c0217 Mon Sep 17 00:00:00 2001 From: YuLiu98 Date: Wed, 22 Mar 2023 22:25:06 +0800 Subject: [PATCH 15/17] Test: update tests due to refactor of chg_extrap --- tests/integrate/170_PW_MD_2O/result.ref | 6 +++--- tests/integrate/270_NO_MD_2O/result.ref | 8 ++++---- 2 files changed, 7 insertions(+), 7 deletions(-) diff --git a/tests/integrate/170_PW_MD_2O/result.ref b/tests/integrate/170_PW_MD_2O/result.ref index be1016e1da..e7efce7c2b 100644 --- a/tests/integrate/170_PW_MD_2O/result.ref +++ b/tests/integrate/170_PW_MD_2O/result.ref @@ -1,5 +1,5 @@ -etotref -211.7989767295324 +etotref -211.7989767295072 etotperatomref -105.8994883648 totalforceref 0.673568 -totalstressref 387.288868 -totaltimeref +1.2656 +totalstressref 387.288858 +totaltimeref +3.5444 diff --git a/tests/integrate/270_NO_MD_2O/result.ref b/tests/integrate/270_NO_MD_2O/result.ref index 5c43d329e3..7f084cdeca 100644 --- a/tests/integrate/270_NO_MD_2O/result.ref +++ b/tests/integrate/270_NO_MD_2O/result.ref @@ -1,5 +1,5 @@ -etotref -211.5884429564663 +etotref -211.5884429564845 etotperatomref -105.7942214782 -totalforceref 1.962895 -totalstressref 594.289483 -totaltimeref +25.661 +totalforceref 1.962896 +totalstressref 594.289512 +totaltimeref +62.39 From 821b2cfbc7694c0aa280ee6beb87c1bd40dc4f96 Mon Sep 17 00:00:00 2001 From: YuLiu98 Date: Wed, 22 Mar 2023 23:36:45 +0800 Subject: [PATCH 16/17] Fix: mpi bug in atom.bcast() --- source/module_cell/atom_spec.cpp | 3 +++ source/module_elecstate/module_charge/charge_extra.cpp | 3 +++ source/module_io/test/prepare_unitcell.h | 6 +++--- tests/integrate/170_PW_MD_2O/result.ref | 6 +++--- tests/integrate/270_NO_MD_2O/result.ref | 8 ++++---- 5 files changed, 16 insertions(+), 10 deletions(-) diff --git a/source/module_cell/atom_spec.cpp b/source/module_cell/atom_spec.cpp index 8a0601950d..c21ead6663 100644 --- a/source/module_cell/atom_spec.cpp +++ b/source/module_cell/atom_spec.cpp @@ -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 ); diff --git a/source/module_elecstate/module_charge/charge_extra.cpp b/source/module_elecstate/module_charge/charge_extra.cpp index 8f711f48dd..5174514db4 100644 --- a/source/module_elecstate/module_charge/charge_extra.cpp +++ b/source/module_elecstate/module_charge/charge_extra.cpp @@ -297,6 +297,9 @@ void Charge_Extra::find_alpha_and_beta(void) } } + GlobalV::ofs_running << " alpha = " << alpha << std::endl; + GlobalV::ofs_running << " beta = " << beta << std::endl; + return; } diff --git a/source/module_io/test/prepare_unitcell.h b/source/module_io/test/prepare_unitcell.h index 14ecfb6a1e..b2a9e70e46 100644 --- a/source/module_io/test/prepare_unitcell.h +++ b/source/module_io/test/prepare_unitcell.h @@ -148,7 +148,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; @@ -157,7 +157,7 @@ class UcellTestPrepare delete[] ucell->atoms[it].m_loc_; delete[] ucell->atoms[it].mbl; ucell->atoms[it].tau = new ModuleBase::Vector3[ucell->atoms[it].na]; - ucell->atoms[it].tau_original = new ModuleBase::Vector3[ucell->atoms[it].na]; + ucell->atoms[it].dis = new ModuleBase::Vector3[ucell->atoms[it].na]; ucell->atoms[it].taud = new ModuleBase::Vector3[ucell->atoms[it].na]; ucell->atoms[it].vel = new ModuleBase::Vector3[ucell->atoms[it].na]; ucell->atoms[it].mag = new double[ucell->atoms[it].na]; @@ -187,7 +187,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]; diff --git a/tests/integrate/170_PW_MD_2O/result.ref b/tests/integrate/170_PW_MD_2O/result.ref index e7efce7c2b..be1016e1da 100644 --- a/tests/integrate/170_PW_MD_2O/result.ref +++ b/tests/integrate/170_PW_MD_2O/result.ref @@ -1,5 +1,5 @@ -etotref -211.7989767295072 +etotref -211.7989767295324 etotperatomref -105.8994883648 totalforceref 0.673568 -totalstressref 387.288858 -totaltimeref +3.5444 +totalstressref 387.288868 +totaltimeref +1.2656 diff --git a/tests/integrate/270_NO_MD_2O/result.ref b/tests/integrate/270_NO_MD_2O/result.ref index 7f084cdeca..5c43d329e3 100644 --- a/tests/integrate/270_NO_MD_2O/result.ref +++ b/tests/integrate/270_NO_MD_2O/result.ref @@ -1,5 +1,5 @@ -etotref -211.5884429564845 +etotref -211.5884429564663 etotperatomref -105.7942214782 -totalforceref 1.962896 -totalstressref 594.289512 -totaltimeref +62.39 +totalforceref 1.962895 +totalstressref 594.289483 +totaltimeref +25.661 From 311bb1134b9942ef4fb7eeb66c2d6e2ce1cd15fb Mon Sep 17 00:00:00 2001 From: YuLiu98 Date: Thu, 23 Mar 2023 09:42:59 +0800 Subject: [PATCH 17/17] Fix: wrong data type for ref_cell_factor --- source/module_io/input.cpp | 2 +- source/module_io/input.h | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/source/module_io/input.cpp b/source/module_io/input.cpp index 3de93ad6af..a1958246b3 100644 --- a/source/module_io/input.cpp +++ b/source/module_io/input.cpp @@ -2737,7 +2737,7 @@ void Input::Bcast() Parallel_Common::bcast_double(search_radius); Parallel_Common::bcast_int(symmetry); Parallel_Common::bcast_bool(init_vel); // liuyu 2021-07-14 - Parallel_Common::bcast_int(ref_cell_factor); + Parallel_Common::bcast_double(ref_cell_factor); Parallel_Common::bcast_double(symmetry_prec); // LiuXh add 2021-08-12, accuracy for symmetry Parallel_Common::bcast_bool(cal_force); Parallel_Common::bcast_bool(dump_force); diff --git a/source/module_io/input.h b/source/module_io/input.h index 6bd598265f..73be441dfd 100644 --- a/source/module_io/input.h +++ b/source/module_io/input.h @@ -47,7 +47,7 @@ class Input bool init_vel; // read velocity from STRU or not liuyu 2021-07-14 bool dump_vel; // output atomic velocities into the file MD_dump or not. liuyu 2023-03-01 - int ref_cell_factor; // construct a reference cell bigger than the initial cell liuyu 2023-03-21 + double ref_cell_factor; // construct a reference cell bigger than the initial cell liuyu 2023-03-21 /* symmetry level: -1, no symmetry at all;