Merge md line of all esolver_types (#1480)

* refactor: merge md line of pw, lcao and others * refactor: update head files
deepmodeling · Nov 9, 2022 · bb888e9 · bb888e9
1 parent 39cc2fc
commit bb888e9
Show file tree

Hide file tree

Showing 29 changed files with 120 additions and 768 deletions.
diff --git a/source/Makefile.Objects b/source/Makefile.Objects
@@ -239,7 +239,7 @@ OBJS_HSOLVER_LCAO=hsolver_lcao.o\
       utils.o\
 
 OBJS_MD=MD_func.o\
-    run_md_classic.o\
+    run_md.o\
 	mdrun.o\
 	verlet.o\
     MSST.o\
@@ -407,7 +407,6 @@ OBJS_LCAO=DM_gamma.o\
       dm_2d.o\
       wavefunc_in_pw.o\
   	  H_TDDFT_pw.o\
-      run_md_lcao.o\
 
 OBJS_RI=exx_lip.o\
 
@@ -558,7 +557,6 @@ OBJS_SRCPW=H_Ewald_pw.o\
     wf_atomic.o\
     wf_igk.o\
     xc_3.o\
-    run_md_pw.o\
     hamilt.o\
     hamilt_pw_old.o\
     electrons.o

diff --git a/source/driver_run.cpp b/source/driver_run.cpp
@@ -4,9 +4,7 @@
 #include "src_io/winput.h"
 #include "module_neighbor/sltk_atom_arrange.h"
 #include "src_io/print_info.h"
-#include "src_lcao/run_md_lcao.h"
-#include "src_pw/run_md_pw.h"
-#include "module_md/run_md_classic.h"
+#include "module_md/run_md.h"
 
 // This is the driver function which defines the workflow of ABACUS calculations
 // It relies on the class Esolver, which is a class that organizes workflows of single point calculations.
@@ -52,24 +50,8 @@ void Driver::driver_run()
     //---------------------------MD/Relax-------------------------
     if(GlobalV::CALCULATION == "md")
     {
-        // In the future, I will universalize pw, lcao and classic line
-        if(GlobalV::ESOLVER_TYPE == "lj" || GlobalV::ESOLVER_TYPE == "dp")
-        {
-            Run_MD_CLASSIC run_md_classic;
-            run_md_classic.classic_md_line(GlobalC::ucell, p_esolver);
-        }
-        else if(GlobalV::BASIS_TYPE == "lcao")
-        {
-#ifdef __LCAO
-            Run_MD_LCAO run_md_lcao;
-            run_md_lcao.opt_ions(p_esolver);
-#endif
-        }
-        else
-        {
-            Run_MD_PW run_md_pw;
-            run_md_pw.md_ions_pw(p_esolver);
-        }
+        Run_MD run_md;
+        run_md.md_line(GlobalC::ucell, p_esolver);
     }
     else // scf; cell relaxation; nscf; etc
     {

diff --git a/source/module_esolver/esolver_dp.cpp b/source/module_esolver/esolver_dp.cpp
@@ -88,4 +88,11 @@ namespace ModuleESolver
         stress = dp_virial;
     }
 
+    void ESolver_DP::postprocess()
+    {
+        GlobalV::ofs_running << "\n\n --------------------------------------------" << std::endl;
+        GlobalV::ofs_running << std::setprecision(16);
+        GlobalV::ofs_running << " !FINAL_ETOT_IS " << dp_potential * ModuleBase::Ry_to_eV << " eV" << std::endl;
+        GlobalV::ofs_running << " --------------------------------------------\n\n" << std::endl;
+    }
 }
diff --git a/source/module_esolver/esolver_dp.h b/source/module_esolver/esolver_dp.h
@@ -29,6 +29,7 @@ namespace ModuleESolver
         void cal_Energy(double& etot) override;
         void cal_Force(ModuleBase::matrix& force) override;
         void cal_Stress(ModuleBase::matrix& stress) override;
+        void postprocess() override;
 
         //--------------temporary----------------------------
 #ifdef __DPMD

diff --git a/source/module_esolver/esolver_ks_lcao.cpp b/source/module_esolver/esolver_ks_lcao.cpp
@@ -612,6 +612,9 @@ 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
@@ -22,6 +22,7 @@
 #endif
 #include "../src_pw/H_Ewald_pw.h"
 #include "module_vdw/vdw.h"
+#include "../module_relax/relax_old/variable_cell.h"    // liuyu 2022-11-07
 
 namespace ModuleESolver
 {
@@ -257,13 +258,28 @@ namespace ModuleESolver
         ModuleBase::TITLE("ESolver_KS_LCAO", "beforescf");
         ModuleBase::timer::tick("ESolver_KS_LCAO", "beforescf");
 
+        // 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();
+
+            if(GlobalC::ucell.cell_parameter_updated)
+            {
+                Variable_Cell::init_after_vc();
+            }
+
+            GlobalC::pot.init_pot(istep, GlobalC::sf.strucFac);
+        }
+
         if(GlobalV::CALCULATION=="relax" || GlobalV::CALCULATION=="cell-relax")
         {
             if(GlobalC::ucell.ionic_position_updated)
             {
                 GlobalV::ofs_running << " Setup the extrapolated charge." << std::endl;
                 // charge extrapolation if istep>0.
-                CE.update_istep(istep);
+                CE.update_istep();
                 CE.update_all_pos(GlobalC::ucell);
                 CE.extrapolate_charge();
                 CE.save_pos_next(GlobalC::ucell);

diff --git a/source/module_esolver/esolver_ks_pw.cpp b/source/module_esolver/esolver_ks_pw.cpp
@@ -17,6 +17,7 @@
 #include "../src_io/chi0_hilbert.h"
 #include "../src_io/epsilon0_pwscf.h"
 #include "../src_io/epsilon0_vasp.h"
+#include "../module_relax/relax_old/variable_cell.h"    // liuyu 2022-11-07
 //-----force-------------------
 #include "../src_pw/forces.h"
 //-----stress------------------
@@ -152,13 +153,28 @@ namespace ModuleESolver
     {
         ModuleBase::TITLE("ESolver_KS_PW", "beforescf");
 
+        // 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();
+
+            if(GlobalC::ucell.cell_parameter_updated)
+            {
+                Variable_Cell::init_after_vc();
+            }
+
+            GlobalC::pot.init_pot(istep, GlobalC::sf.strucFac);
+        }
+
         if(GlobalV::CALCULATION=="relax" || GlobalV::CALCULATION=="cell-relax")
         {
             if(GlobalC::ucell.ionic_position_updated)
             {
                 GlobalV::ofs_running << " Setup the extrapolated charge." << std::endl;
                 // charge extrapolation if istep>0.
-                CE.update_istep(istep);
+                CE.update_istep();
                 CE.update_all_pos(GlobalC::ucell);
                 CE.extrapolate_charge();
                 CE.save_pos_next(GlobalC::ucell);
@@ -414,6 +430,9 @@ namespace ModuleESolver
 
     void ESolver_KS_PW::afterscf(const int istep)
     {
+        // Temporary liuyu add 2022-11-07
+        CE.update_all_pos(GlobalC::ucell);
+
 #ifdef __LCAO
         if (GlobalC::chi0_hilbert.epsilon)                 // pengfei 2016-11-23
         {

diff --git a/source/module_esolver/esolver_lj.cpp b/source/module_esolver/esolver_lj.cpp
@@ -104,6 +104,14 @@ namespace ModuleESolver
         stress = lj_virial;
     }
 
+    void ESolver_LJ::postprocess()
+    {
+        GlobalV::ofs_running << "\n\n --------------------------------------------" << std::endl;
+        GlobalV::ofs_running << std::setprecision(16);
+        GlobalV::ofs_running << " !FINAL_ETOT_IS " << lj_potential * ModuleBase::Ry_to_eV << " eV" << std::endl;
+        GlobalV::ofs_running << " --------------------------------------------\n\n" << std::endl;
+    }
+
     double ESolver_LJ::LJ_energy(const double d)
     {
         double temp = pow(lj_sigma / d, 6);

diff --git a/source/module_esolver/esolver_lj.h b/source/module_esolver/esolver_lj.h
@@ -19,6 +19,7 @@ namespace ModuleESolver
         void cal_Energy(double& etot) override;
         void cal_Force(ModuleBase::matrix& force) override;
         void cal_Stress(ModuleBase::matrix& stress) override;
+        void postprocess() override;
 
 
         double LJ_energy(const double d);

diff --git a/source/module_esolver/esolver_of.cpp b/source/module_esolver/esolver_of.cpp
@@ -21,6 +21,7 @@
 //---------------------------------------------------
 #include "module_elecstate/elecstate_pw.h"
 #include "module_hamilt/hamilt_pw.h"
+#include "module_relax/relax_old/variable_cell.h"    // liuyu 2022-11-07
 
 namespace ModuleESolver
 {
@@ -202,13 +203,16 @@ void ESolver_OF::Init(Input &inp, UnitCell_pseudo &ucell)
     this->vw.set_para(this->nrxx, this->dV, GlobalV::of_vw_weight);
     this->wt.set_para(this->nrxx, this->dV, GlobalV::of_wt_alpha, GlobalV::of_wt_beta, this->nelec[0], GlobalV::of_tf_weight, GlobalV::of_vw_weight, GlobalV::of_read_kernel, GlobalV::of_kernel_file, this->pw_rho);
     ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, "INIT KEDF");
+
+    // Initialize charge extrapolation
+    CE.Init_CE();
 }
 
 void ESolver_OF::Run(int istep, UnitCell_pseudo& ucell)
 {
     ModuleBase::timer::tick("ESolver_OF", "Run");
     // get Ewald energy, initial rho and phi if necessary
-    this->beforeOpt();
+    this->beforeOpt(istep);
     this->iter = 0;
 
     while(true)
@@ -242,8 +246,23 @@ void ESolver_OF::Run(int istep, UnitCell_pseudo& ucell)
 // 
 // Calculate ewald energy, initialize the rho, phi, theta
 // 
-void ESolver_OF::beforeOpt()
+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.extrapolate_charge();
+
+        if(GlobalC::ucell.cell_parameter_updated)
+        {
+            Variable_Cell::init_after_vc();
+        }
+
+        GlobalC::pot.init_pot(istep, GlobalC::sf.strucFac);
+    }
+
     //calculate ewald energy
     H_Ewald_pw::compute_ewald(GlobalC::ucell, this->pw_rho);
 
@@ -829,6 +848,9 @@ 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_esolver/esolver_of.h b/source/module_esolver/esolver_of.h
@@ -2,6 +2,7 @@
 #include "../module_base/opt_TN.hpp"
 #include "../module_base/opt_DCsrch.h"
 #include "../module_psi/psi.h"
+#include "../src_pw/charge_extra.h"    // liuyu add 2022-11-07
 #include "./kedf_tf.h"
 #include "./kedf_vw.h"
 #include "./kedf_wt.h"
@@ -88,6 +89,9 @@ class ESolver_OF: public ESolver_FP
     KEDF_vW vw;
     KEDF_WT wt;
 
+    // charge extrapolation liuyu add 2022-11-07
+    Charge_Extra CE;
+
     // optimization methods
     ModuleBase::Opt_CG opt_cg;
     ModuleBase::Opt_TN opt_tn;
@@ -136,7 +140,7 @@ class ESolver_OF: public ESolver_FP
     double normdLdphi = 100.;
 
     // main process of OFDFT
-    void beforeOpt();
+    void beforeOpt(const int istep);
     void updateV();
     void solveV();
     void getNextDirect();

diff --git a/source/module_md/CMakeLists.txt b/source/module_md/CMakeLists.txt
@@ -2,7 +2,7 @@ add_library(
     md
     OBJECT
     MD_func.cpp
-    run_md_classic.cpp
+    run_md.cpp
 	mdrun.cpp
 	verlet.cpp
     MSST.cpp