Merge pull request #1186 from wenfei-li/develop

refactor : fully decouple esolver and relaxation ; fix bug in envelop function

CMakeLists.txt

@@ -325,7 +325,7 @@ target_link_libraries(${ABACUS_BIN_NAME}
     neighbor
     orb
     io
-    ions
+    relax
     lcao
     gint
     parallel

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

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)

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\

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\

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

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"

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

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}

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,

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<<GlobalC::rhopw->nrxx<<endl;
@@ -177,10 +162,17 @@ namespace ModuleESolver
                 GlobalV::ofs_running << " Setup the Vl+Vh+Vxc according to new structure factor and new charge." << std::endl;
                 // calculate the new potential accordint to
                 // the new charge density.
-                GlobalC::pot.init_pot( istep-1, GlobalC::sf.strucFac );
+                GlobalC::pot.init_pot( 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, 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

source/src_ions/CMakeLists.txt → 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
 )

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" 
+		<<std::setw(5)<< "NE"
+		<<std::setw(15)<< "ETOT(eV)"
+		<<std::setw(15)<< "EDIFF(eV)"
+        <<std::setw(15)<< "MAX_F(eV/A)"
+        <<std::setw(15)<< "TRADIUS(Bohr)"
+		<<std::setw(8)<< "UPDATE"
+		<<std::setw(11)<< "ETIME(MIN)"
+		<<std::setw(11)<< "FTIME(MIN)"
+        <<std::endl;
+	}
+
+    // Geometry optimization algorithm setup.
+    if (GlobalV::CALCULATION=="relax")
+    {
+        //Ions_Move_Methods
+        IMM.allocate();
+    }
+    if (GlobalV::CALCULATION=="cell-relax")
+    {
+        //Ions_Move_Methods
+        IMM.allocate();
+        // allocate arrays related to changes of lattice vectors
+        LCM.allocate();
+    } 
+
+    this->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;
+}

source/src_ions/ions.h → 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);