Merge pull request #927 from YuLiu98/develop

feature : Dipole correction

Author: ouqi0711

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,6 +243,7 @@ OBJS_SURCHEM=H_correction_pw.o\
     cal_vel.o\
     corrected_energy.o\
     minimize_cg.o\
+    efield.o\
 
 OBJS_XC=xc_funct_corr_gga.o \
 xc_funct_corr_lda.o \
@@ -293,7 +294,6 @@ write_HS.o\
 write_HS_R.o\
 write_dm.o\
 write_wfc_realspace.o\
-efield.o \
 magnetism.o\
 optical.o\
 Cell_PW.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

@@ -12,7 +12,6 @@
 #include "src_io/epsilon0_vasp.h"
 #include "src_io/optical.h"
 #include "src_ions/ions_move_basic.h"
-#include "src_pw/efield.h"
 #include "src_pw/global.h"
 #include "src_pw/occupy.h"
 #ifdef __EXX
@@ -27,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)
 {
@@ -210,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_surchem/CMakeLists.txt

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