feature : fix efield and rename Dipole class to be Efield class

Author: YuLiu98

docs/input-main.md

@@ -70,6 +70,10 @@
 
     [deepks_out_labels](#out-descriptor) | [deepks_descriptor_lmax](#lmax-descriptor) | [deepks_scf](#deepks-scf) | [deepks_model](#model-file)
 
+- [Electric field and dipole correction](#Electric-field-and-dipole-correction)
+
+    [efield](#efield) | [dipole](#dipole) | [edir](#edir) | [emaxpos](#emaxpos) | [eopreg](#eopreg) | [eamp](#eamp)
+
 [back to main page](../README.md)
 
 ## Structure of the file
@@ -910,6 +914,48 @@ Warning: this function is not robust enough for version 2.2.0. Please try these
 - **Description**: the path of the trained, traced NN model file (generated by deepks-kit). used when deepks_scf is set to 1.
 - **Default**: None
 
+### Electric field and dipole correction
+
+This part of variables are relevant to electric field and dipole correction
+
+#### efield
+
+- **Type**: Boolean
+- **Description**: If set to true, a saw-like potential simulating an electric field
+is added to the bare ionic potential. 
+- **Default**: false
+
+#### dipole
+
+- **Type**: Boolean
+- **Description**: If dipole == true and efield == true,  a dipole correction is also
+added to the bare ionic potential. If you want no electric field, parameter eamp should be zero. Must be used ONLY in a slab geometry for surface calculations, with the discontinuity FALLING IN THE EMPTY SPACE.
+- **Default**: false
+
+#### edir
+
+- **Type**: Integer
+- **Description**: The direction of the electric field or dipole correction is parallel to the reciprocal lattice vector, so the potential is constant in planes defined by FFT grid points, edir = 0, 1 or 2. Used only if efield == true.
+- **Default**: 2
+
+#### emaxpos
+
+- **Type**: Real
+- **Description**: Position of the maximum of the saw-like potential along crystal axis edir, within the  unit cell, 0 < emaxpos < 1. Used only if efield == true.
+- **Default**: 0.5
+
+#### eopreg
+
+- **Type**: Real
+- **Description**: Zone in the unit cell where the saw-like potential decreases, 0 < eopreg < 1. Used only if efield == true.
+- **Default**: 0.1
+
+#### eamp
+
+- **Type**: Real
+- **Description**: Amplitude of the electric field, in ***Hartree*** a.u.; 1 a.u. = 51.4220632*10^10 V/m. Used only if efield == true. The saw-like potential increases with slope eamp in the region from (emaxpos+eopreg-1) to (emaxpos), then decreases until (emaxpos+eopreg), in units of the crystal vector edir. Important: the change of slope of this potential must be located in the empty region, or else unphysical forces will result.
+- **Default**: 0.0
+
 ### Exact Exchange
 
 This part of variables are relevant when using hybrid functionals

source/Makefile.Objects

@@ -243,7 +243,7 @@ OBJS_SURCHEM=H_correction_pw.o\
     cal_vel.o\
     corrected_energy.o\
     minimize_cg.o\
-    dipole.o\
+    efield.o\
 
 OBJS_XC=xc_funct_corr_gga.o \
 xc_funct_corr_lda.o \

source/input.cpp

@@ -288,6 +288,15 @@ void Input::Default(void)
     lcao_dr = 0.01;
     lcao_rmax = 30; // (a.u.)
     //----------------------------------------------------------
+    // efield and dipole correction     Yu Liu add 2022-05-18
+    //----------------------------------------------------------
+    efield = false;
+    dipole = false;
+    edir = 2;
+    emaxpos = 0.5;
+    eopreg = 0.1;
+    eamp = 0.0;
+    //----------------------------------------------------------
     // vdw									//jiyy add 2019-08-04
     //----------------------------------------------------------
     vdw_method = "none";
@@ -1155,6 +1164,34 @@ bool Input::Read(const std::string &fn)
             read_value(ifs, mdp.md_damp);
         }
         //----------------------------------------------------------
+        // efield and dipole correction
+        // Yu Liu add 2022-05-18
+        //----------------------------------------------------------
+        else if (strcmp("efield", word) == 0)
+        {
+            read_value(ifs, efield);
+        }
+        else if (strcmp("dipole", word) == 0)
+        {
+            read_value(ifs, dipole);
+        }
+        else if (strcmp("edir", word) == 0)
+        {
+            read_value(ifs, edir);
+        }
+        else if (strcmp("emaxpos", word) == 0)
+        {
+            read_value(ifs, emaxpos);
+        }
+        else if (strcmp("eopreg", word) == 0)
+        {
+            read_value(ifs, eopreg);
+        }
+        else if (strcmp("eamp", word) == 0)
+        {
+            read_value(ifs, eamp);
+        }
+        //----------------------------------------------------------
         // tddft
         // Fuxiang He add 2016-10-26
         //----------------------------------------------------------
@@ -1997,6 +2034,13 @@ void Input::Bcast()
     Parallel_Common::bcast_double(mdp.msst_tscale);
     Parallel_Common::bcast_double(mdp.md_tfreq);
     Parallel_Common::bcast_double(mdp.md_damp);
+    // Yu Liu add 2022-05-18
+    Parallel_Common::bcast_bool(efield);
+    Parallel_Common::bcast_bool(dipole);
+    Parallel_Common::bcast_int(edir);
+    Parallel_Common::bcast_double(emaxpos);
+    Parallel_Common::bcast_double(eopreg);
+    Parallel_Common::bcast_double(eamp);
     /* 	// Peize Lin add 2014-04-07
         Parallel_Common::bcast_bool( vdwD2 );
         Parallel_Common::bcast_double( vdwD2_scaling );

source/input.h

@@ -255,6 +255,17 @@ class Input
         int md_msdstartTime;            //choose which step that msd be calculated */
     MD_parameters mdp;
 
+    //==========================================================
+    // efield and dipole correction
+    // Yu Liu add 2022-05-18
+    //==========================================================
+    bool efield;        // add electric field
+    bool dipole;        // dipole correction
+    int edir;           // the direction of the electric field or dipole correction
+    double emaxpos;     // position of the maximum of the saw-like potential along crystal axis edir
+    double eopreg;      // zone in the unit cell where the saw-like potential decreases
+    double eamp;        // amplitude of the electric field
+
     //==========================================================
     // vdw
     // Peize Lin add 2014-03-31, jiyy update 2019-08-01

source/input_conv.cpp

@@ -26,6 +26,7 @@
 #include "src_lcao/local_orbital_charge.h"
 #endif
 #include "module_base/timer.h"
+#include "module_surchem/efield.h"
 
 void Input_Conv::Convert(void)
 {
@@ -209,6 +210,16 @@ void Input_Conv::Convert(void)
         GlobalV::NPOL = 1;
     }
 
+//----------------------------------------------------------
+// Yu Liu add 2022-05-18
+//----------------------------------------------------------
+    GlobalV::EFIELD = INPUT.efield;
+    GlobalV::DIPOLE = INPUT.dipole;
+    Efield::edir = INPUT.edir;
+    Efield::emaxpos = INPUT.emaxpos;
+    Efield::eopreg = INPUT.eopreg;
+    Efield::eamp = INPUT.eamp;
+
 //----------------------------------------------------------
 // Fuxiang He add 2016-10-26
 //----------------------------------------------------------

source/module_esolver/esolver_ks_lcao.cpp

@@ -14,7 +14,6 @@
 #include "src_pw/symmetry_rho.h"
 #include "src_io/chi0_hilbert.h"
 #include "src_pw/threshold_elec.h"
-#include "module_surchem/dipole.h"
 
 #ifdef __DEEPKS
 #include "../module_deepks/LCAO_deepks.h"
@@ -316,7 +315,6 @@ namespace ModuleESolver
         if (iter == 1)
         {
             GlobalC::CHR.set_new_e_iteration(true);
-            Dipole::first = true;
         }
         else
         {

source/module_esolver/esolver_ks_pw.cpp

@@ -21,7 +21,6 @@
 #include "../src_io/chi0_hilbert.h"
 #include "../src_io/epsilon0_pwscf.h"
 #include "../src_io/epsilon0_vasp.h"
-#include "../module_surchem/dipole.h"
 //-----force-------------------
 #include "../src_pw/forces.h"
 //-----stress------------------
@@ -197,15 +196,8 @@ namespace ModuleESolver
     void ESolver_KS_PW::eachiterinit(const int istep, const int iter)
     {
         // mohan add 2010-07-16
-        if (iter == 1) 
-        {
-            GlobalC::CHR.set_new_e_iteration(true);
-            Dipole::first = true;
-        }
-        else
-        {
-            GlobalC::CHR.set_new_e_iteration(false);
-        }
+        if (iter == 1) GlobalC::CHR.set_new_e_iteration(true);
+        else GlobalC::CHR.set_new_e_iteration(false);
 
         if (GlobalV::FINAL_SCF && iter == 1)
         {

source/module_surchem/CMakeLists.txt

@@ -9,5 +9,5 @@ add_library(
     cal_vel.cpp
     corrected_energy.cpp
     minimize_cg.cpp
-    dipole.cpp
+    efield.cpp
 )

source/module_surchem/dipole.cpp renamed to source/module_surchem/efield.cpp

@@ -1,44 +1,32 @@
-#include "dipole.h"
+#include "efield.h"
 #include "../module_base/constants.h"
 #include "../module_base/timer.h"
 #include "../module_base/global_variable.h"
 #include "../src_parallel/parallel_reduce.h"
 
-double Dipole::etotefield = 0.0;
-double Dipole::tot_dipole = 0.0;
-int Dipole::edir = 2;
-double Dipole::emaxpos = 0.5;
-double Dipole::eopreg = 0.1;
-double Dipole::eamp = 0.0;
-bool Dipole::first = true;
-double Dipole::bvec[3];
-double Dipole::bmod;
+double Efield::etotefield = 0.0;
+double Efield::tot_dipole = 0.0;
+int Efield::edir;
+double Efield::emaxpos;
+double Efield::eopreg;
+double Efield::eamp;
+double Efield::bvec[3];
+double Efield::bmod;
 
-Dipole::Dipole(){}
+Efield::Efield(){}
 
-Dipole::~Dipole(){}
+Efield::~Efield(){}
 
 //=======================================================
 // calculate dipole potential in surface calculations
 //=======================================================
-ModuleBase::matrix Dipole::add_efield(const UnitCell &cell, 
+ModuleBase::matrix Efield::add_efield(const UnitCell &cell, 
                                         PW_Basis &pwb, 
                                         const int &nspin, 
                                         const double *const *const rho)
 {
-    ModuleBase::TITLE("Dipole", "add_efield");
-    ModuleBase::timer::tick("Dipole", "add_efield");
-
-    ModuleBase::matrix v(nspin, pwb.nrxx);
-
-    // efield only needs to be added on the first iteration, if dipfield
-    // is not used. note that for relax calculations it has to be added
-    // again on subsequent relax steps.
-    if(!GlobalV::DIPOLE && !first)
-    {
-        return v;
-    }
-    first = false;
+    ModuleBase::TITLE("Efield", "add_efield");
+    ModuleBase::timer::tick("Efield", "add_efield");
 
     double latvec;    // latvec along the efield direction
     if(edir == 0)
@@ -64,7 +52,7 @@ ModuleBase::matrix Dipole::add_efield(const UnitCell &cell,
     }
     else
     {
-        ModuleBase::WARNING_QUIT("Dipole::ion_dipole", "direction is wrong!");
+        ModuleBase::WARNING_QUIT("Efield::add_efield", "direction is wrong!");
     }
     bmod = sqrt(pow(bvec[0],2) + pow(bvec[1],2) + pow(bvec[2],2));
 
@@ -83,20 +71,22 @@ ModuleBase::matrix Dipole::add_efield(const UnitCell &cell,
     else
     {
         ion_dipole = cal_ion_dipole(cell, bmod);
-        tot_dipole = ion_dipole - elec_dipole;
 
         // energy correction
-        etotefield = - ModuleBase::e2 * eamp * tot_dipole * cell.omega / ModuleBase::FOUR_PI;
+        etotefield = - ModuleBase::e2 * eamp * ion_dipole * cell.omega / ModuleBase::FOUR_PI;
     }
 
     const double length = (1.0 - eopreg) * latvec * cell.lat0;
     const double vamp = ModuleBase::e2 * (eamp - tot_dipole) * length;
 
     GlobalV::ofs_running << "\n\n Adding external electric field: " << std::endl;
-    ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, "Computed dipole along edir", edir);
-    ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, "Elec. dipole (Ry a.u.)", elec_dipole);
-    ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, "Ion dipole (Ry a.u.)", ion_dipole);
-    ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, "Total dipole (Ry a.u.)", tot_dipole);
+    if(GlobalV::DIPOLE)
+    {
+        ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, "Computed dipole along edir", edir);
+        ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, "Elec. dipole (Ry a.u.)", elec_dipole);
+        ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, "Ion dipole (Ry a.u.)", ion_dipole);
+        ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, "Total dipole (Ry a.u.)", tot_dipole);
+    }
     if( abs(eamp) > 0.0) 
     {
         ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, "Amplitute of Efield (Hartree)", eamp);
@@ -105,6 +95,7 @@ ModuleBase::matrix Dipole::add_efield(const UnitCell &cell,
     ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, "Total length (Bohr)", length);
 
     // dipole potential
+    ModuleBase::matrix v(nspin, pwb.nrxx);
     const int nspin0 = (nspin == 2) ? 2 : 1;
 
     for (int ir = 0; ir < pwb.nrxx; ++ir)
@@ -127,15 +118,15 @@ ModuleBase::matrix Dipole::add_efield(const UnitCell &cell,
 
     double fac = ModuleBase::e2 * (eamp - tot_dipole) * cell.lat0 / bmod;
 
-    ModuleBase::timer::tick("Dipole", "add_efield");
+    ModuleBase::timer::tick("Efield", "add_efield");
     return v * fac;
 }
 
 
 //=======================================================
 // calculate dipole density in surface calculations
 //=======================================================
-double Dipole::cal_ion_dipole(const UnitCell &cell, const double &bmod)
+double Efield::cal_ion_dipole(const UnitCell &cell, const double &bmod)
 {
     double ion_dipole = 0;
     for(int it=0; it<cell.ntype; ++it)
@@ -152,7 +143,7 @@ double Dipole::cal_ion_dipole(const UnitCell &cell, const double &bmod)
     return ion_dipole;
 }
 
-double Dipole::cal_elec_dipole(const UnitCell &cell, 
+double Efield::cal_elec_dipole(const UnitCell &cell, 
                             PW_Basis &pwb, 
                             const int &nspin, 
                             const double *const *const rho,
@@ -185,7 +176,7 @@ double Dipole::cal_elec_dipole(const UnitCell &cell,
     return elec_dipole;
 }
 
-double Dipole::saw_function(const double &a, const double &b, const double &x)
+double Efield::saw_function(const double &a, const double &b, const double &x)
 {
     assert(x>=0);
     assert(x<=1);
@@ -206,7 +197,7 @@ double Dipole::saw_function(const double &a, const double &b, const double &x)
     }
 }
 
-void Dipole::compute_force(const UnitCell &cell, ModuleBase::matrix &fdip)
+void Efield::compute_force(const UnitCell &cell, ModuleBase::matrix &fdip)
 {
     if(GlobalV::DIPOLE)
     {