diff --git a/CMakeLists.txt b/CMakeLists.txt index 8d732a0195..cd506c2230 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -325,7 +325,7 @@ target_link_libraries(${ABACUS_BIN_NAME} neighbor orb io - ions + relax lcao gint parallel diff --git a/docs/install.md b/docs/install.md index 77a81f155c..4b76a9f87a 100644 --- a/docs/install.md +++ b/docs/install.md @@ -285,7 +285,7 @@ OPTS = ${INCLUDES} -Ofast -traceback -std=c++14 -simd -march=native -xHost -m64 - src_external - src_global - src_io -- src_ions +- module_relaxation - src_lcao - src_parallel - src_pdiag diff --git a/source/CMakeLists.txt b/source/CMakeLists.txt index 8cd202fbdf..c573b98311 100644 --- a/source/CMakeLists.txt +++ b/source/CMakeLists.txt @@ -15,7 +15,7 @@ add_subdirectory(module_xc) add_subdirectory(module_esolver) add_subdirectory(module_gint) add_subdirectory(src_io) -add_subdirectory(src_ions) +add_subdirectory(module_relaxation) add_subdirectory(src_lcao) add_subdirectory(src_parallel) add_subdirectory(src_pdiag) diff --git a/source/Makefile b/source/Makefile index 2872780da9..5043526890 100644 --- a/source/Makefile +++ b/source/Makefile @@ -27,7 +27,7 @@ VPATH=./src_global\ :./module_gint\ :./src_pw\ :./src_lcao\ -:./src_ions\ +:./module_relaxation\ :./src_io\ :./src_parallel\ :./src_pdiag\ diff --git a/source/Makefile.Objects b/source/Makefile.Objects index 202b46c648..36799e6333 100644 --- a/source/Makefile.Objects +++ b/source/Makefile.Objects @@ -16,8 +16,7 @@ OBJS_MAIN=driver.o\ input.o \ write_input.o\ input_conv.o\ -run_pw.o\ -run_lcao.o\ +driver_run.o OBJS_PW=xc_3.o \ vdwd2.o\ @@ -339,7 +338,8 @@ write_wfc_realspace.o\ magnetism.o\ optical.o\ run_md_pw.o\ -ions.o \ +ions.o\ +relaxation.o\ ions_move_methods.o\ ions_move_bfgs.o\ ions_move_cg.o\ diff --git a/source/input_conv.cpp b/source/input_conv.cpp index 56c2abab05..1ff0f0a621 100644 --- a/source/input_conv.cpp +++ b/source/input_conv.cpp @@ -12,7 +12,7 @@ #include "src_io/epsilon0_pwscf.h" #include "src_io/epsilon0_vasp.h" #include "src_io/optical.h" -#include "src_ions/ions_move_basic.h" +#include "module_relaxation/ions_move_basic.h" #include "src_pw/global.h" #include "src_pw/occupy.h" #ifdef __EXX diff --git a/source/input_update.cpp b/source/input_update.cpp index ebfae1fd67..077569b00f 100644 --- a/source/input_update.cpp +++ b/source/input_update.cpp @@ -7,7 +7,7 @@ #include "module_base/global_function.h" #include "module_base/global_variable.h" #include "input.h" -#include "src_ions/ions_move_basic.h" +#include "module_relaxation/ions_move_basic.h" #include "src_io/optical.h" #ifdef __LCAO #include "src_lcao/FORCE_STRESS.h" diff --git a/source/module_cell/unitcell.h b/source/module_cell/unitcell.h index 63ab3de228..9913d9165f 100644 --- a/source/module_cell/unitcell.h +++ b/source/module_cell/unitcell.h @@ -76,6 +76,7 @@ class UnitCell //I'm doing a bad thing here! Will change later bool ionic_position_updated = false; //whether the ionic position has been updated + bool cell_parameter_updated = false; //whether the cell parameters are updated private: ModuleBase::Matrix3 stress; //calculate stress on the cell diff --git a/source/module_deepks/test/CMakeLists.txt b/source/module_deepks/test/CMakeLists.txt index d9b8947f29..0702aa30bc 100644 --- a/source/module_deepks/test/CMakeLists.txt +++ b/source/module_deepks/test/CMakeLists.txt @@ -7,7 +7,7 @@ get_target_property(ABACUS_LINK_LIBRARIES ${ABACUS_BIN_NAME} LINK_LIBRARIES) target_link_libraries( test_deepks base cell symmetry md surchem xc_ - neighbor orb io ions gint lcao parallel mrrr pdiag pw ri driver esolver hsolver psi elecstate hamilt planewave + neighbor orb io relax gint lcao parallel mrrr pdiag pw ri driver esolver hsolver psi elecstate hamilt planewave pthread deepks ${ABACUS_LINK_LIBRARIES} diff --git a/source/module_esolver/esolver_ks.cpp b/source/module_esolver/esolver_ks.cpp index c9afb48df4..8df0dc4b46 100644 --- a/source/module_esolver/esolver_ks.cpp +++ b/source/module_esolver/esolver_ks.cpp @@ -82,6 +82,24 @@ namespace ModuleESolver // mohan add 2021-01-30 Print_Info::setup_parameters(ucell, GlobalC::kv); + if(GlobalV::BASIS_TYPE=="pw" || GlobalV::CALCULATION=="ienvelope") + { + //Envelope function is calculated as lcao_in_pw + //new plane wave basis + #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->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(); + this->print_wfcfft(inp, GlobalV::ofs_running); + } // 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, diff --git a/source/module_esolver/esolver_ks_pw.cpp b/source/module_esolver/esolver_ks_pw.cpp index 52c72ca808..d13448ca07 100644 --- a/source/module_esolver/esolver_ks_pw.cpp +++ b/source/module_esolver/esolver_ks_pw.cpp @@ -69,21 +69,6 @@ namespace ModuleESolver void ESolver_KS_PW::Init_GlobalC(Input& inp, UnitCell_pseudo& cell) { - //new plane wave basis -#ifdef __MPI - this->pw_wfc->initmpi(GlobalV::NPROC_IN_POOL, GlobalV::RANK_IN_POOL, POOL_WORLD); -#endif - this->pw_wfc->initgrids(cell.lat0, cell.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(); - ESolver_KS::print_wfcfft(inp, GlobalV::ofs_running); - this->psi = GlobalC::wf.allocate(GlobalC::kv.nks); // cout<nrxx<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, this->psi); + 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 diff --git a/source/src_ions/CMakeLists.txt b/source/module_relaxation/CMakeLists.txt similarity index 90% rename from source/src_ions/CMakeLists.txt rename to source/module_relaxation/CMakeLists.txt index 16f0d07c15..835d3e68fa 100644 --- a/source/src_ions/CMakeLists.txt +++ b/source/module_relaxation/CMakeLists.txt @@ -1,5 +1,5 @@ add_library( - ions + relax OBJECT bfgs_basic.cpp ions.cpp @@ -12,4 +12,5 @@ add_library( lattice_change_cg.cpp lattice_change_methods.cpp variable_cell.cpp + relaxation.cpp ) diff --git a/source/src_ions/bfgs_basic.cpp b/source/module_relaxation/bfgs_basic.cpp similarity index 100% rename from source/src_ions/bfgs_basic.cpp rename to source/module_relaxation/bfgs_basic.cpp diff --git a/source/src_ions/bfgs_basic.h b/source/module_relaxation/bfgs_basic.h similarity index 100% rename from source/src_ions/bfgs_basic.h rename to source/module_relaxation/bfgs_basic.h diff --git a/source/module_relaxation/ions.cpp b/source/module_relaxation/ions.cpp new file mode 100644 index 0000000000..faba21f59d --- /dev/null +++ b/source/module_relaxation/ions.cpp @@ -0,0 +1,114 @@ +#include "ions.h" +#include "../src_pw/global.h" // use chr. +#include "../src_io/print_info.h" +#include "variable_cell.h" // mohan add 2021-02-01 +#include "src_io/write_wfc_realspace.h" + +void Ions::opt_ions(ModuleESolver::ESolver *p_esolver) +{ + ModuleBase::TITLE("Ions","opt_ions"); + ModuleBase::timer::tick("Ions","opt_ions"); + + if(GlobalV::OUT_LEVEL=="i") + { + std::cout << std::setprecision(12); + std::cout<< " " << std::setw(7)<< "ISTEP" + <istep = 1; + int force_step = 1; // pengfei Li 2018-05-14 + int stress_step = 1; + bool stop= false; + + while(istep <= GlobalV::RELAX_NMAX && !stop) + { + time_t estart = time(NULL); + + if(GlobalV::OUT_LEVEL=="ie") + { + Print_Info::print_screen(stress_step, force_step, istep); + } + + // mohan added eiter to count for the electron iteration number, 2021-01-28 + p_esolver->Run(istep-1,GlobalC::ucell); + + time_t eend = time(NULL); + time_t fstart = time(NULL); + + if (GlobalV::CALCULATION=="scf" || GlobalV::CALCULATION=="relax" || GlobalV::CALCULATION=="cell-relax" || GlobalV::CALCULATION.substr(0,3)=="sto") + { + //I'm considering putting force and stress + //as part of ucell and use ucell to pass information + //back and forth between esolver and relaxation + //but I'll use force and stress explicitly here for now + + //calculate and gather all parts of total ionic forces + ModuleBase::matrix force; + if(GlobalV::CAL_FORCE) + { + p_esolver->cal_Force(force); + } + //calculate and gather all parts of stress + ModuleBase::matrix stress; + if(GlobalV::CAL_STRESS) + { + p_esolver->cal_Stress(stress); + } + stop = this->relaxation(force, stress, istep, force_step, stress_step); // pengfei Li 2018-05-14 + } + time_t fend = time(NULL); + + if(GlobalV::OUT_LEVEL=="i") + { + double etime_min = difftime(eend, estart)/60.0; + double ftime_min = difftime(fend, fstart)/60.0; + std::stringstream ss; + ss << GlobalV::RELAX_METHOD << istep; + + std::cout << " " << std::setw(7) << ss.str() + << std::setw(5) << p_esolver->getniter() + << std::setw(15) << std::setprecision(6) << GlobalC::en.etot * ModuleBase::Ry_to_eV + << std::setw(15) << IMM.get_ediff() * ModuleBase::Ry_to_eV + << std::setprecision(3) + << std::setw(15) << IMM.get_largest_grad() * ModuleBase::Ry_to_eV / 0.529177 + << std::setw(15) << IMM.get_trust_radius() + << std::setw(8) << IMM.get_update_iter() + << std::setprecision(2) << std::setw(11) << etime_min + << std::setw(11) << ftime_min << std::endl; + } + + ++istep; + + } + + if(GlobalV::OUT_LEVEL=="i") + { + std::cout << " ION DYNAMICS FINISHED :)" << std::endl; + } + + ModuleBase::timer::tick("Ions","opt_ions_pw"); + return; +} \ No newline at end of file diff --git a/source/src_ions/ions.h b/source/module_relaxation/ions.h similarity index 67% rename from source/src_ions/ions.h rename to source/module_relaxation/ions.h index 04d8589e47..6bc1b0591a 100644 --- a/source/src_ions/ions.h +++ b/source/module_relaxation/ions.h @@ -10,6 +10,12 @@ #include "lattice_change_methods.h" #include "module_esolver/esolver.h" + +//The workflow opt_ions should be moved outside module_relaxation +//since the latter is intended to perform the sole task of +//creating the next step structure based on force and stress +//according to some relaxation algorithm +//However, it will remain this way until the ucell class and MD module are sorted out class Ions { @@ -36,13 +42,13 @@ class Ions Lattice_Change_Methods LCM; //seperate force_stress function first - bool after_scf(ModuleESolver::ESolver *p_esolver,const int &istep, int &force_step, int &stress_step); + bool relaxation(ModuleBase::matrix force,ModuleBase::matrix stress,const int &istep, int &force_step, int &stress_step); bool if_do_relax(); bool if_do_cellrelax(); bool do_relax(const int& istep, int& jstep, const ModuleBase::matrix& ionic_force, const double& total_energy); bool do_cellrelax(const int& istep, const ModuleBase::matrix& stress, const double& total_energy); void reset_after_relax(const int& istep); - void reset_after_cellrelax(int& force_step, int& stress_step, ModuleESolver::ESolver *p_esolver); + void reset_after_cellrelax(int& force_step, int& stress_step); void update_pot(void); diff --git a/source/src_ions/ions_move_basic.cpp b/source/module_relaxation/ions_move_basic.cpp similarity index 100% rename from source/src_ions/ions_move_basic.cpp rename to source/module_relaxation/ions_move_basic.cpp diff --git a/source/src_ions/ions_move_basic.h b/source/module_relaxation/ions_move_basic.h similarity index 100% rename from source/src_ions/ions_move_basic.h rename to source/module_relaxation/ions_move_basic.h diff --git a/source/src_ions/ions_move_bfgs.cpp b/source/module_relaxation/ions_move_bfgs.cpp similarity index 100% rename from source/src_ions/ions_move_bfgs.cpp rename to source/module_relaxation/ions_move_bfgs.cpp diff --git a/source/src_ions/ions_move_bfgs.h b/source/module_relaxation/ions_move_bfgs.h similarity index 100% rename from source/src_ions/ions_move_bfgs.h rename to source/module_relaxation/ions_move_bfgs.h diff --git a/source/src_ions/ions_move_cg.cpp b/source/module_relaxation/ions_move_cg.cpp similarity index 100% rename from source/src_ions/ions_move_cg.cpp rename to source/module_relaxation/ions_move_cg.cpp diff --git a/source/src_ions/ions_move_cg.h b/source/module_relaxation/ions_move_cg.h similarity index 100% rename from source/src_ions/ions_move_cg.h rename to source/module_relaxation/ions_move_cg.h diff --git a/source/src_ions/ions_move_methods.cpp b/source/module_relaxation/ions_move_methods.cpp similarity index 100% rename from source/src_ions/ions_move_methods.cpp rename to source/module_relaxation/ions_move_methods.cpp diff --git a/source/src_ions/ions_move_methods.h b/source/module_relaxation/ions_move_methods.h similarity index 100% rename from source/src_ions/ions_move_methods.h rename to source/module_relaxation/ions_move_methods.h diff --git a/source/src_ions/ions_move_sd.cpp b/source/module_relaxation/ions_move_sd.cpp similarity index 100% rename from source/src_ions/ions_move_sd.cpp rename to source/module_relaxation/ions_move_sd.cpp diff --git a/source/src_ions/ions_move_sd.h b/source/module_relaxation/ions_move_sd.h similarity index 100% rename from source/src_ions/ions_move_sd.h rename to source/module_relaxation/ions_move_sd.h diff --git a/source/src_ions/lattice_change_basic.cpp b/source/module_relaxation/lattice_change_basic.cpp similarity index 100% rename from source/src_ions/lattice_change_basic.cpp rename to source/module_relaxation/lattice_change_basic.cpp diff --git a/source/src_ions/lattice_change_basic.h b/source/module_relaxation/lattice_change_basic.h similarity index 100% rename from source/src_ions/lattice_change_basic.h rename to source/module_relaxation/lattice_change_basic.h diff --git a/source/src_ions/lattice_change_cg.cpp b/source/module_relaxation/lattice_change_cg.cpp similarity index 100% rename from source/src_ions/lattice_change_cg.cpp rename to source/module_relaxation/lattice_change_cg.cpp diff --git a/source/src_ions/lattice_change_cg.h b/source/module_relaxation/lattice_change_cg.h similarity index 100% rename from source/src_ions/lattice_change_cg.h rename to source/module_relaxation/lattice_change_cg.h diff --git a/source/src_ions/lattice_change_methods.cpp b/source/module_relaxation/lattice_change_methods.cpp similarity index 100% rename from source/src_ions/lattice_change_methods.cpp rename to source/module_relaxation/lattice_change_methods.cpp diff --git a/source/src_ions/lattice_change_methods.h b/source/module_relaxation/lattice_change_methods.h similarity index 100% rename from source/src_ions/lattice_change_methods.h rename to source/module_relaxation/lattice_change_methods.h diff --git a/source/module_relaxation/relaxation.cpp b/source/module_relaxation/relaxation.cpp new file mode 100644 index 0000000000..1a26c9c4a6 --- /dev/null +++ b/source/module_relaxation/relaxation.cpp @@ -0,0 +1,124 @@ +#include "ions.h" +#include "../src_pw/global.h" // use chr. +#include "../src_io/print_info.h" +#include "variable_cell.h" // mohan add 2021-02-01 + +// The interface for relaxation +bool Ions::relaxation(ModuleBase::matrix force, ModuleBase::matrix stress, const int &istep, int &force_step, int &stress_step) +{ + ModuleBase::TITLE("Ions","after_scf"); + + // should not do it this way, will change after the refactor of ucell class + GlobalC::ucell.ionic_position_updated = false; + GlobalC::ucell.cell_parameter_updated = false; + + //stop in last step + if(istep==GlobalV::RELAX_NMAX) + { + return 1; + } + //choose what to do next + if(GlobalV::CALCULATION!="cell-relax") force_step = istep; + if(this->if_do_relax()) + { + //do relax calculation and generate next structure + bool converged = 0; + converged = this->do_relax(istep, force_step, force, GlobalC::en.etot); + if(!converged) + { + this->reset_after_relax(istep); + GlobalC::ucell.ionic_position_updated = true; + return converged; + } + else if(GlobalV::CALCULATION!="cell-relax") + { + return converged; + } + } + if(this->if_do_cellrelax()) + { + //do cell relax calculation and generate next structure + bool converged = 0; + converged = this->do_cellrelax(stress_step, stress, GlobalC::en.etot); + if(!converged) + { + GlobalC::ucell.cell_parameter_updated = true; + this->reset_after_cellrelax(force_step, stress_step); + } + return converged; + } + + return 1; +} + +bool Ions::if_do_relax() +{ + ModuleBase::TITLE("Ions","if_do_relax"); + if(GlobalV::CALCULATION=="relax"||GlobalV::CALCULATION=="cell-relax") + { + if(!GlobalC::ucell.if_atoms_can_move()) + { + ModuleBase::WARNING("Ions","No atom is allowed to move!"); + return 0; + } +// if(!IMM.get_converged()) return 1; + else + { + assert(GlobalV::CAL_FORCE==1); + return 1; + } + } + else return 0; +} +bool Ions::if_do_cellrelax() +{ + ModuleBase::TITLE("Ions","if_do_cellrelax"); + if(GlobalV::CALCULATION=="cell-relax") + { + if(!GlobalC::ucell.if_cell_can_change()) + { + ModuleBase::WARNING("Ions", "Lattice vectors are not allowed to change!"); + return 0; + } + else if(GlobalC::ucell.if_atoms_can_move()&&!IMM.get_converged()) + { + GlobalV::ofs_running<<"Note: Need to wait for atomic relaxation first!"; + return 0; + } + else + { + assert(GlobalV::CAL_STRESS==1); + return 1; + } + } + else return 0; +} +bool Ions::do_relax(const int& istep, int& jstep, const ModuleBase::matrix& ionic_force, const double& total_energy) +{ + ModuleBase::TITLE("Ions","do_relax"); + IMM.cal_movement(istep, jstep, ionic_force, total_energy); + ++jstep; + return IMM.get_converged(); +} +bool Ions::do_cellrelax(const int& istep, const ModuleBase::matrix& stress, const double& total_energy) +{ + ModuleBase::TITLE("Ions","do_cellrelax"); + LCM.cal_lattice_change(istep, stress, total_energy); + return LCM.get_converged(); +} +void Ions::reset_after_relax(const int& istep) +{ + ModuleBase::TITLE("Ions","reset_after_relax"); + GlobalV::ofs_running << " Setup the structure factor in plane wave basis." << std::endl; + GlobalC::sf.setup_structure_factor(&GlobalC::ucell,GlobalC::rhopw); +} + +void Ions::reset_after_cellrelax(int& f_step, int& s_step) +{ + ModuleBase::TITLE("Ions","reset_after_cellrelax"); + Variable_Cell::init_after_vc(); + GlobalC::pot.init_pot(s_step, GlobalC::sf.strucFac); //LiuXh add 20180619 + + f_step = 1; + ++s_step; +} diff --git a/source/src_ions/variable_cell.cpp b/source/module_relaxation/variable_cell.cpp similarity index 78% rename from source/src_ions/variable_cell.cpp rename to source/module_relaxation/variable_cell.cpp index 6860bd63dd..190d48d94b 100644 --- a/source/src_ions/variable_cell.cpp +++ b/source/module_relaxation/variable_cell.cpp @@ -5,7 +5,7 @@ Variable_Cell::Variable_Cell(){} Variable_Cell::~Variable_Cell(){} -void Variable_Cell::init_after_vc(ModuleESolver::ESolver *p_esolver) +void Variable_Cell::init_after_vc() { ModuleBase::TITLE("Variable_Cell","init_after_vc"); @@ -29,15 +29,6 @@ void Variable_Cell::init_after_vc(ModuleESolver::ESolver *p_esolver) GlobalC::rhopw->collect_uniqgg(); GlobalC::sf.setup_structure_factor(&GlobalC::ucell,GlobalC::rhopw); - if(GlobalV::BASIS_TYPE=="pw") - { - 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, p_esolver->psi); - GlobalC::wf.init_at_1(); - } - GlobalV::ofs_running << " Setup the Vl+Vh+Vxc according to new structure factor and new charge." << std::endl; //================================= // initalize local pseudopotential diff --git a/source/src_ions/variable_cell.h b/source/module_relaxation/variable_cell.h similarity index 83% rename from source/src_ions/variable_cell.h rename to source/module_relaxation/variable_cell.h index 6fc5185e90..909263b794 100644 --- a/source/src_ions/variable_cell.h +++ b/source/module_relaxation/variable_cell.h @@ -18,7 +18,7 @@ class Variable_Cell Variable_Cell(); ~Variable_Cell(); - static void init_after_vc(ModuleESolver::ESolver *p_esolver); //LiuXh add 20180515 + static void init_after_vc(); //LiuXh add 20180515 }; diff --git a/source/src_ions/ions.cpp b/source/src_ions/ions.cpp deleted file mode 100644 index 7987758ae1..0000000000 --- a/source/src_ions/ions.cpp +++ /dev/null @@ -1,225 +0,0 @@ -#include "ions.h" -#include "../src_pw/global.h" // use chr. -#include "../src_io/print_info.h" -#include "variable_cell.h" // mohan add 2021-02-01 -#include "src_io/write_wfc_realspace.h" - -void Ions::opt_ions(ModuleESolver::ESolver *p_esolver) -{ - ModuleBase::TITLE("Ions","opt_ions"); - ModuleBase::timer::tick("Ions","opt_ions"); - - if(GlobalV::OUT_LEVEL=="i") - { - std::cout << std::setprecision(12); - std::cout<< " " << std::setw(7)<< "ISTEP" - <istep = 1; - int force_step = 1; // pengfei Li 2018-05-14 - int stress_step = 1; - bool stop= false; - - while(istep <= GlobalV::RELAX_NMAX && !stop) - { - time_t estart = time(NULL); - - if(GlobalV::OUT_LEVEL=="ie") - { - Print_Info::print_screen(stress_step, force_step, istep); - } - - // mohan added eiter to count for the electron iteration number, 2021-01-28 - p_esolver->Run(istep-1,GlobalC::ucell); - - time_t eend = time(NULL); - time_t fstart = time(NULL); - - if (GlobalV::CALCULATION=="scf" || GlobalV::CALCULATION=="relax" || GlobalV::CALCULATION=="cell-relax" || GlobalV::CALCULATION.substr(0,3)=="sto") - { - stop = this->after_scf(p_esolver, istep, force_step, stress_step); // pengfei Li 2018-05-14 - } - time_t fend = time(NULL); - - if(GlobalV::OUT_LEVEL=="i") - { - double etime_min = difftime(eend, estart)/60.0; - double ftime_min = difftime(fend, fstart)/60.0; - std::stringstream ss; - ss << GlobalV::RELAX_METHOD << istep; - - std::cout << " " << std::setw(7) << ss.str() - << std::setw(5) << p_esolver->getniter() - << std::setw(15) << std::setprecision(6) << GlobalC::en.etot * ModuleBase::Ry_to_eV - << std::setw(15) << IMM.get_ediff() * ModuleBase::Ry_to_eV - << std::setprecision(3) - << std::setw(15) << IMM.get_largest_grad() * ModuleBase::Ry_to_eV / 0.529177 - << std::setw(15) << IMM.get_trust_radius() - << std::setw(8) << IMM.get_update_iter() - << std::setprecision(2) << std::setw(11) << etime_min - << std::setw(11) << ftime_min << std::endl; - } - - ++istep; - - } - - if(GlobalV::OUT_LEVEL=="i") - { - std::cout << " ION DYNAMICS FINISHED :)" << std::endl; - } - - ModuleBase::timer::tick("Ions","opt_ions_pw"); - return; -} - -bool Ions::after_scf(ModuleESolver::ESolver *p_esolver, const int &istep, int &force_step, int &stress_step) -{ - ModuleBase::TITLE("Ions","after_scf"); - - // should not do it this way, will change later - GlobalC::ucell.ionic_position_updated = false; - - //calculate and gather all parts of total ionic forces - ModuleBase::matrix force; - if(GlobalV::CAL_FORCE) - { - p_esolver->cal_Force(force); - } - //calculate and gather all parts of stress - ModuleBase::matrix stress; - if(GlobalV::CAL_STRESS) - { - p_esolver->cal_Stress(stress); - } - //stop in last step - if(istep==GlobalV::RELAX_NMAX) - { - return 1; - } - //choose what to do next - if(GlobalV::CALCULATION!="cell-relax") force_step = istep; - if(this->if_do_relax()) - { - //do relax calculation and generate next structure - bool converged = 0; - converged = this->do_relax(istep, force_step, force, GlobalC::en.etot); - if(!converged) - { - this->reset_after_relax(istep); - GlobalC::ucell.ionic_position_updated = true; - return converged; - } - else if(GlobalV::CALCULATION!="cell-relax") - { - return converged; - } - } - if(this->if_do_cellrelax()) - { - //do cell relax calculation and generate next structure - bool converged = 0; - converged = this->do_cellrelax(stress_step, stress, GlobalC::en.etot); - if(!converged) this->reset_after_cellrelax(force_step, stress_step, p_esolver); - return converged; - } - - return 1; -} - -bool Ions::if_do_relax() -{ - ModuleBase::TITLE("Ions","if_do_relax"); - if(GlobalV::CALCULATION=="relax"||GlobalV::CALCULATION=="cell-relax") - { - if(!GlobalC::ucell.if_atoms_can_move()) - { - ModuleBase::WARNING("Ions","No atom is allowed to move!"); - return 0; - } -// if(!IMM.get_converged()) return 1; - else - { - assert(GlobalV::CAL_FORCE==1); - return 1; - } - } - else return 0; -} -bool Ions::if_do_cellrelax() -{ - ModuleBase::TITLE("Ions","if_do_cellrelax"); - if(GlobalV::CALCULATION=="cell-relax") - { - if(!GlobalC::ucell.if_cell_can_change()) - { - ModuleBase::WARNING("Ions", "Lattice vectors are not allowed to change!"); - return 0; - } - else if(GlobalC::ucell.if_atoms_can_move()&&!IMM.get_converged()) - { - GlobalV::ofs_running<<"Note: Need to wait for atomic relaxation first!"; - return 0; - } - else - { - assert(GlobalV::CAL_STRESS==1); - return 1; - } - } - else return 0; -} -bool Ions::do_relax(const int& istep, int& jstep, const ModuleBase::matrix& ionic_force, const double& total_energy) -{ - ModuleBase::TITLE("Ions","do_relax"); - IMM.cal_movement(istep, jstep, ionic_force, total_energy); - ++jstep; - return IMM.get_converged(); -} -bool Ions::do_cellrelax(const int& istep, const ModuleBase::matrix& stress, const double& total_energy) -{ - ModuleBase::TITLE("Ions","do_cellrelax"); - LCM.cal_lattice_change(istep, stress, total_energy); - return LCM.get_converged(); -} -void Ions::reset_after_relax(const int& istep) -{ - ModuleBase::TITLE("Ions","reset_after_relax"); - GlobalV::ofs_running << " Setup the structure factor in plane wave basis." << std::endl; - GlobalC::sf.setup_structure_factor(&GlobalC::ucell,GlobalC::rhopw); -} - -void Ions::reset_after_cellrelax(int& f_step, int& s_step, ModuleESolver::ESolver *p_esolver) -{ - ModuleBase::TITLE("Ions","reset_after_cellrelax"); - Variable_Cell::init_after_vc(p_esolver); - GlobalC::pot.init_pot(s_step, GlobalC::sf.strucFac); //LiuXh add 20180619 - - //GlobalV::ofs_running << " Setup the new wave functions?" << std::endl; //LiuXh add 20180619 - //GlobalC::wf.wfcinit(p_esolver->psi); //LiuXh add 20180619 - f_step = 1; - ++s_step; -} diff --git a/source/src_lcao/local_orbital_charge.h b/source/src_lcao/local_orbital_charge.h index bdb631f924..0b42f0f06d 100644 --- a/source/src_lcao/local_orbital_charge.h +++ b/source/src_lcao/local_orbital_charge.h @@ -109,7 +109,7 @@ class Local_Orbital_Charge int lgd_last;// sub-FFT-mesh orbitals number in previous step. int lgd_now;// sub-FFT-mesh orbitals number in this step. - int nnrg_last;// sub-FFT-mesh orbtials number in previous step, with k. + int nnrg_last = 0;// sub-FFT-mesh orbtials number in previous step, with k. int nnrg_now; // sub-FFT-mesh orbitals number in this step, with k. // add by yshen on 9/22/2014 diff --git a/source/src_lcao/run_md_lcao.cpp b/source/src_lcao/run_md_lcao.cpp index 64257b5cf3..9f51933fd0 100644 --- a/source/src_lcao/run_md_lcao.cpp +++ b/source/src_lcao/run_md_lcao.cpp @@ -9,7 +9,7 @@ #include "../src_io/write_HS.h" #include "../src_io/cal_r_overlap_R.h" #include "../src_io/print_info.h" -#include "../src_ions/variable_cell.h" // mohan add 2021-02-01 +#include "../module_relaxation/variable_cell.h" // mohan add 2021-02-01 #include "../src_ri/exx_abfs.h" #include "../src_ri/exx_opt_orb.h" #include "../module_neighbor/sltk_atom_arrange.h" @@ -102,7 +102,7 @@ void Run_MD_LCAO::opt_ions(ModuleESolver::ESolver* p_esolver) if (cellchange) { - Variable_Cell::init_after_vc(p_esolver); + Variable_Cell::init_after_vc(); } // reset local potential diff --git a/source/src_pw/Makefile b/source/src_pw/Makefile index 119d21dd09..5354849dfd 100644 --- a/source/src_pw/Makefile +++ b/source/src_pw/Makefile @@ -72,7 +72,7 @@ VPATH=../src_global\ :../module_symmetry\ :../src_parallel\ :../src_io\ -:../src_ions\ +:../module_relaxation\ :../module_md\ :../module_symmetry\ :../\ diff --git a/source/src_pw/Makefile.Objects b/source/src_pw/Makefile.Objects index 3f0ba70178..a9035ba7ca 100644 --- a/source/src_pw/Makefile.Objects +++ b/source/src_pw/Makefile.Objects @@ -50,7 +50,8 @@ OBJS_IONS=MD_basic.o\ MD_thermo.o\ MD_fire.o\ MD_func.o\ -ions.o \ +ions.o\ +relaxation.o\ ions_move_methods.o\ ions_move_bfgs.o\ ions_move_cg.o\ diff --git a/source/src_pw/Makefile.parallel b/source/src_pw/Makefile.parallel index 551acf9176..f7c6ae8113 100644 --- a/source/src_pw/Makefile.parallel +++ b/source/src_pw/Makefile.parallel @@ -58,7 +58,7 @@ VPATH=../src_global\ :../module_symmetry\ :../src_parallel\ :../src_io\ -:../src_ions\ +:../module_relaxation\ :../module_md\ :../module_symmetry\ :../\ diff --git a/source/src_pw/Makefile.serial b/source/src_pw/Makefile.serial index 682fed75ca..1bd6369b4b 100644 --- a/source/src_pw/Makefile.serial +++ b/source/src_pw/Makefile.serial @@ -69,7 +69,7 @@ VPATH=../src_global\ :../module_cell\ :../src_parallel\ :../src_io\ -:../src_ions\ +:../module_relaxation\ :../module_md\ :../module_symmetry\ :../\ diff --git a/source/src_pw/global.h b/source/src_pw/global.h index 59e4af0c4c..85b864ce4e 100644 --- a/source/src_pw/global.h +++ b/source/src_pw/global.h @@ -4,7 +4,7 @@ #include "../module_base/global_function.h" #include "../module_base/global_variable.h" #include "../src_io/restart.h" -#include "../src_ions/ions.h" +#include "../module_relaxation/ions.h" #include "../src_lcao/exx_lip.h" #include "VNL_in_pw.h" #include "charge_broyden.h" diff --git a/source/src_pw/potential.cpp b/source/src_pw/potential.cpp index ea69b27ffd..cde1ad785f 100644 --- a/source/src_pw/potential.cpp +++ b/source/src_pw/potential.cpp @@ -229,10 +229,6 @@ void Potential::init_pot(const int &istep, // number of ionic steps GlobalC::restart.info_load.load_charge_finish = true; } } - else - { - // the extrapolation part moves to ions.cpp. - } // renormalize the charge density GlobalC::CHR.renormalize_rho(); diff --git a/source/src_pw/run_md_pw.cpp b/source/src_pw/run_md_pw.cpp index 4cfe754c1f..bf91c4b9aa 100644 --- a/source/src_pw/run_md_pw.cpp +++ b/source/src_pw/run_md_pw.cpp @@ -1,6 +1,6 @@ #include "run_md_pw.h" #include "global.h" // use chr. -#include "../src_ions/variable_cell.h" // mohan add 2021-02-01 +#include "../module_relaxation/variable_cell.h" // mohan add 2021-02-01 #include "../module_md/MD_func.h" #include "../module_md/FIRE.h" #include "../module_md/NVE.h" @@ -97,7 +97,8 @@ void Run_MD_PW::md_ions_pw(ModuleESolver::ESolver *p_esolver) if(cellchange) { - Variable_Cell::init_after_vc(p_esolver); + GlobalC::ucell.cell_parameter_updated = true; + Variable_Cell::init_after_vc(); } // reset local potential and initial wave function diff --git a/tests/integrate/170_PW_MD_2O/result.ref b/tests/integrate/170_PW_MD_2O/result.ref index 803dfee9c6..f91c08a25e 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.7989767495471 -etotperatomref -105.8994883748 -totalforceref 0.673488 -totalstressref 387.295738 +etotref -211.7989767457971 +etotperatomref -105.8994883729 +totalforceref 0.673564 +totalstressref 387.288753 totaltimeref +0.40399 diff --git a/tests/integrate/270_NO_MD_2O/result.ref b/tests/integrate/270_NO_MD_2O/result.ref index 362736b507..4b8bc71414 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.5879950287044 -etotperatomref -105.7939975144 -totalforceref 1.954326 -totalstressref 595.789229 +etotref -211.5879950293724 +etotperatomref -105.7939975147 +totalforceref 1.954342 +totalstressref 595.789919 totaltimeref +5.996