feature : incorporate efield into dipole.cpp

Author: YuLiu98

source/Makefile.Objects

@@ -295,7 +295,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_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

source/module_esolver/esolver_ks_lcao.cpp

@@ -14,6 +14,7 @@
 #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"
@@ -315,6 +316,7 @@ namespace ModuleESolver
         if (iter == 1)
         {
             GlobalC::CHR.set_new_e_iteration(true);
+            Dipole::first = true;
         }
         else
         {

source/module_esolver/esolver_ks_pw.cpp

@@ -21,6 +21,7 @@
 #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------------------
@@ -195,8 +196,15 @@ 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);
-        else GlobalC::CHR.set_new_e_iteration(false);
+        if (iter == 1) 
+        {
+            GlobalC::CHR.set_new_e_iteration(true);
+            Dipole::first = true;
+        }
+        else
+        {
+            GlobalC::CHR.set_new_e_iteration(false);
+        }
 
         if (GlobalV::FINAL_SCF && iter == 1)
         {

source/module_surchem/dipole.cpp

@@ -4,10 +4,15 @@
 #include "../module_base/global_variable.h"
 #include "../src_parallel/parallel_reduce.h"
 
-int Dipole::dir = 2;
-double Dipole::dipole_energy = 0.0;
-double Dipole::max_pos = 0.5;
-double Dipole::de_reg = 0.1;
+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;
 
 Dipole::Dipole(){}
 
@@ -16,41 +21,90 @@ Dipole::~Dipole(){}
 //=======================================================
 // calculate dipole potential in surface calculations
 //=======================================================
-ModuleBase::matrix Dipole::v_dipole(const UnitCell &cell, 
-                                    PW_Basis &pwb, 
-                                    const int &nspin, 
-                                    const double *const *const rho)
+ModuleBase::matrix Dipole::add_efield(const UnitCell &cell, 
+                                        PW_Basis &pwb, 
+                                        const int &nspin, 
+                                        const double *const *const rho)
 {
-    ModuleBase::TITLE("Dipole", "v_dipole");
-    ModuleBase::timer::tick("Dipole", "v_dipole");
+    ModuleBase::TITLE("Dipole", "add_efield");
+    ModuleBase::timer::tick("Dipole", "add_efield");
 
-    double h_inv;    // inverse of height
-    if(dir == 0)
+    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;
+
+    double latvec;    // latvec along the efield direction
+    if(edir == 0)
     {
-        h_inv = sqrt(cell.G.e11 * cell.G.e11 + cell.G.e12 * cell.G.e12 + cell.G.e13 * cell.G.e13);
+        bvec[0] = cell.G.e11;
+        bvec[1] = cell.G.e12; 
+        bvec[2] = cell.G.e13; 
+        latvec = cell.a1.norm();
     }
-    else if(dir == 1)
+    else if(edir == 1)
     {
-        h_inv = sqrt(cell.G.e21 * cell.G.e21 + cell.G.e22 * cell.G.e22 + cell.G.e23 * cell.G.e23);
+        bvec[0] = cell.G.e21;
+        bvec[1] = cell.G.e22; 
+        bvec[2] = cell.G.e23; 
+        latvec = cell.a2.norm();
     }
-    else if(dir = 2)
+    else if(edir = 2)
+    {
+        bvec[0] = cell.G.e31;
+        bvec[1] = cell.G.e32; 
+        bvec[2] = cell.G.e33; 
+        latvec = cell.a3.norm();
+    }
+    else
+    {
+        ModuleBase::WARNING_QUIT("Dipole::ion_dipole", "direction is wrong!");
+    }
+    bmod = sqrt(pow(bvec[0],2) + pow(bvec[1],2) + pow(bvec[2],2));
+
+    double ion_dipole = 0;
+    double elec_dipole = 0;
+
+    if(GlobalV::DIPOLE)
     {
-        h_inv = sqrt(cell.G.e31 * cell.G.e31 + cell.G.e32 * cell.G.e32 + cell.G.e33 * cell.G.e33);
+        ion_dipole = cal_ion_dipole(cell, bmod);
+        elec_dipole = cal_elec_dipole(cell, pwb, nspin, rho, bmod);
+        tot_dipole = ion_dipole - elec_dipole;
+
+        // energy correction
+        etotefield = - ModuleBase::e2 * (eamp - 0.5 * tot_dipole) * tot_dipole * cell.omega / ModuleBase::FOUR_PI;
     }
     else
     {
-        ModuleBase::WARNING_QUIT("Dipole::ion_dipole", "dipole direction is wrong!");
+        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;
     }
 
-    double ion_dipole = cal_ion_dipole(cell, h_inv);
-    double elec_dipole = cal_elec_dipole(cell, pwb, nspin, rho, h_inv);
-    double tot_dipole = ion_dipole - elec_dipole;
+    const double length = (1.0 - eopreg) * latvec * cell.lat0;
+    const double vamp = ModuleBase::e2 * (eamp - tot_dipole) * length;
 
-    // dipole energy correction
-    dipole_energy = 0.5 * ModuleBase::e2 * tot_dipole * tot_dipole * cell.omega / ModuleBase::FOUR_PI;
+    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( abs(eamp) > 0.0) 
+    {
+        ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, "Amplitute of Efield (Hartree)", eamp);
+    }
+    ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, "Potential amplitute (Ry)", vamp);
+    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)
@@ -63,80 +117,37 @@ ModuleBase::matrix Dipole::v_dipole(const UnitCell &cell,
         double z = (double)k / pwb.ncz;
         ModuleBase::Vector3<double> pos(x, y, z);
 
-        double saw = saw_function(max_pos, de_reg, pos[dir]);
+        double saw = saw_function(emaxpos, eopreg, pos[edir]);
 
         for (int is = 0; is < nspin0; is++)
         {
             v(is, ir) = saw;
         }
     }
 
-    double fac = - ModuleBase::e2 * tot_dipole * cell.lat0 / h_inv;
+    double fac = ModuleBase::e2 * (eamp - tot_dipole) * cell.lat0 / bmod;
 
-    ModuleBase::timer::tick("Dipole", "v_dipole");
+    ModuleBase::timer::tick("Dipole", "add_efield");
     return v * fac;
 }
 
 
 //=======================================================
 // calculate dipole density in surface calculations
 //=======================================================
-// double Dipole::dipole(const UnitCell &cell, 
-//                         PW_Basis &pwb, 
-//                         const double *const *const rho,
-//                         complex<double> *TOTN)
-// {
-//     double *Porter = new double[pwb.nrxx];
-
-//     for (int ir = 0; ir < pwb.nrxx; ++ir)
-//     {
-//         int i = ir / (pwb.ncy * pwb.nczp);
-//         int j = ir / pwb.nczp - i * pwb.ncy;
-//         int k = ir % pwb.nczp + pwb.nczp_start;
-//         double x = (double)i / pwb.ncx;
-//         double y = (double)j / pwb.ncy;
-//         double z = (double)k / pwb.ncz;
-//         ModuleBase::Vector3<double> pos(x, y, z);
-
-//         pos = cell.latvec * pos;
-
-//         Porter[ir] = cell.lat0 * pos[dir];
-//     }
-
-//     complex<double> *Porter_g = new complex<double>[pwb.ngmc];
-//     GlobalC::UFFT.ToReciSpace(Porter, Porter_g);
-
-//     double m = 0;
-//     for (int ig = 0; ig < pwb.ngmc; ig++)
-//     {
-//         m += (conj(TOTN[ig]) * Porter_g[ig]).real();
-//     }
-
-//     Parallel_Reduce::reduce_double_pool(m);
-
-//     // height
-//     const double h = cell.latvec.to_matrix()(dir, dir) * cell.lat0;
-
-//     delete[] Porter, Porter_g;
-
-//     return -m * h * ModuleBase::e2;
-// }
-
-double Dipole::cal_ion_dipole(const UnitCell &cell, const double &h_inv)
+double Dipole::cal_ion_dipole(const UnitCell &cell, const double &bmod)
 {
     double ion_dipole = 0;
     for(int it=0; it<cell.ntype; ++it)
     {
         double sum = 0;
         for(int ia=0; ia<cell.atoms[it].na; ++ia)
         {
-            sum += saw_function(max_pos, de_reg, cell.atoms[it].taud[ia][dir]);
+            sum += saw_function(emaxpos, eopreg, cell.atoms[it].taud[ia][edir]);
         }
         ion_dipole += sum * cell.atoms[it].zv;
     }
-    ion_dipole *= cell.lat0 / h_inv * ModuleBase::FOUR_PI / cell.omega;
-
-    // std::cout << "ion_dipole = " << ion_dipole << std::endl;
+    ion_dipole *= cell.lat0 / bmod * ModuleBase::FOUR_PI / cell.omega;
 
     return ion_dipole;
 }
@@ -145,7 +156,7 @@ double Dipole::cal_elec_dipole(const UnitCell &cell,
                             PW_Basis &pwb, 
                             const int &nspin, 
                             const double *const *const rho,
-                            const double &h_inv)
+                            const double &bmod)
 {
     double elec_dipole = 0;
     const int nspin0 = (nspin == 2) ? 2 : 1;
@@ -160,7 +171,7 @@ double Dipole::cal_elec_dipole(const UnitCell &cell,
         double z = (double)k / pwb.ncz;
         ModuleBase::Vector3<double> pos(x, y, z);
 
-        double saw = saw_function(max_pos, de_reg, pos[dir]);
+        double saw = saw_function(emaxpos, eopreg, pos[edir]);
 
         for (int is = 0; is < nspin0; is++)
         {
@@ -169,9 +180,7 @@ double Dipole::cal_elec_dipole(const UnitCell &cell,
     }
 
     Parallel_Reduce::reduce_double_pool(elec_dipole);
-    elec_dipole *= cell.lat0 / h_inv * ModuleBase::FOUR_PI / pwb.ncxyz;
-
-    // std::cout << "elec_dipole = " << elec_dipole << std::endl;
+    elec_dipole *= cell.lat0 / bmod * ModuleBase::FOUR_PI / pwb.ncxyz;
 
     return elec_dipole;
 }
@@ -195,4 +204,38 @@ double Dipole::saw_function(const double &a, const double &b, const double &x)
     {
         return 0.5 * fac - fac * (x - a) / b;
     }
+}
+
+void Dipole::compute_force(const UnitCell &cell, ModuleBase::matrix &fdip)
+{
+    if(GlobalV::DIPOLE)
+    {
+        int iat = 0;
+        for(int it=0; it<cell.ntype; ++it)
+        {
+            for(int ia=0; ia<cell.atoms[it].na; ++ia)
+            {
+                for(int jj=0; jj<3; ++jj)
+                {
+                    fdip(iat, jj) = ModuleBase::e2 * (eamp - tot_dipole) * cell.atoms[it].zv * bvec[jj] / bmod;
+                }
+                ++iat;
+            }
+        }
+    }
+    else
+    {
+        int iat = 0;
+        for(int it=0; it<cell.ntype; ++it)
+        {
+            for(int ia=0; ia<cell.atoms[it].na; ++ia)
+            {
+                for(int jj=0; jj<3; ++jj)
+                {
+                    fdip(iat, jj) = ModuleBase::e2 * eamp * cell.atoms[it].zv * bvec[jj] / bmod;
+                }
+                ++iat;
+            }
+        }
+    }
 }

source/module_surchem/dipole.h

@@ -10,15 +10,10 @@ class Dipole
     Dipole();
     ~Dipole();
 
-    static double dipole(const UnitCell &cell, 
-                            PW_Basis &pwb, 
-                            const double *const *const rho, 
-                            complex<double> *TOTN);
-
-    static ModuleBase::matrix v_dipole(const UnitCell &cell, 
-                                        PW_Basis &pwb, 
-                                        const int &nspin, 
-                                        const double *const *const rho);
+    static ModuleBase::matrix add_efield(const UnitCell &cell, 
+                                            PW_Basis &pwb, 
+                                            const int &nspin, 
+                                            const double *const *const rho);
 
     static double cal_elec_dipole(const UnitCell &cell, 
                                 PW_Basis &pwb, 
@@ -30,12 +25,19 @@ class Dipole
 
     static double saw_function(const double &a, const double &b, const double &x);
 
+    static void compute_force(const UnitCell &cell, ModuleBase::matrix &fdip);
+
 
 
-    static double dipole_energy;        // dipole energy
-    static int dir;                     // 0, 1, 2 denotes x, y, z direction for dipole correction
-    static double max_pos;              // the maximum position of the saw function
-    static double de_reg;               // the decrease region length of the saw function
+    static double etotefield;           // dipole energy
+    static double tot_dipole;           // total dipole
+    static int edir;                    // 0, 1, 2 denotes x, y, z direction for dipole correction
+    static double emaxpos;              // the maximum position of the saw function
+    static double eopreg;               // the decrease region length of the saw function
+    static double eamp;                 // field amplitude (in a.u.) (1 a.u. = 51.44 10^10 V/m)
+    static bool first;                  // efield only needs to be added on the first iteration, if dipfield is not used.
+    static double bvec[3];
+    static double bmod;
 };