Add function of adding compensating charge

ddhhss · ddhhss · commit ed1beb8c7992 · 2022-05-30T15:08:32.000+08:00
diff --git a/source/Makefile.Objects b/source/Makefile.Objects
@@ -240,6 +240,7 @@ OBJS_ESOLVER=esolver.o\
 
 OBJS_SURCHEM=H_correction_pw.o\
     cal_epsilon.o\
+    surchem.o\
     cal_pseudo.o\
     cal_totn.o\
     cal_vcav.o\
@@ -415,7 +416,6 @@ pzt2s.o\
 pdtrsm.o\
 pzhtrsm.o\
 
-
 PDIAG_MR_0=dcopy.o\
 dlae2.o\
 dlaebz.o\
diff --git a/source/input.cpp b/source/input.cpp
@@ -420,6 +420,11 @@ void Input::Default(void)
     sigma_k = 0.6;
     nc_k = 0.00037;
 
+    comp_q = 0.0;
+    comp_l = 1.0;
+    comp_center = 0.0;
+    comp_dim = 2;
+
     return;
 }
 
@@ -1495,6 +1500,22 @@ bool Input::Read(const std::string &fn)
         {
             read_value(ifs, nc_k);
         }
+        else if (strcmp("comp_q", word) == 0)
+        {
+            read_value(ifs, comp_q);
+        }
+        else if (strcmp("comp_l", word) == 0)
+        {
+            read_value(ifs, comp_l);
+        }
+        else if (strcmp("comp_center", word) == 0)
+        {
+            read_value(ifs, comp_center);
+        }
+        else if (strcmp("comp_dim", word) == 0)
+        {
+            read_value(ifs, comp_dim);
+        }
         //----------------------------------------------------------------------------------
         else
         {
diff --git a/source/input.h b/source/input.h
@@ -399,6 +399,11 @@ class Input
     double tau;
     double sigma_k;
     double nc_k;
+    // gauss charge
+    double comp_q;
+    double comp_l;
+    double comp_center;
+    int comp_dim;
 
     //==========================================================
     // variables for test only
diff --git a/source/input_conv.cpp b/source/input_conv.cpp
@@ -14,6 +14,7 @@
 #include "src_ions/ions_move_basic.h"
 #include "src_pw/global.h"
 #include "src_pw/occupy.h"
+#include "module_surchem/surchem.h"
 #ifdef __EXX
 #include "src_ri/exx_abfs-jle.h"
 #endif
@@ -501,6 +502,11 @@ void Input_Conv::Convert(void)
     GlobalV::tau = INPUT.tau;
     GlobalV::sigma_k = INPUT.sigma_k;
     GlobalV::nc_k = INPUT.nc_k;
+
+    GlobalC::solvent_model.comp_q = INPUT.comp_q;
+    GlobalC::solvent_model.comp_l = INPUT.comp_l;
+    GlobalC::solvent_model.comp_center = INPUT.comp_center;
+    GlobalC::solvent_model.comp_dim = INPUT.comp_dim;
     ModuleBase::timer::tick("Input_Conv", "Convert");
     return;
 }
diff --git a/source/module_surchem/cal_totn.cpp b/source/module_surchem/cal_totn.cpp
@@ -1,8 +1,94 @@
 #include "surchem.h"
 
+void add_comp_chg(const UnitCell &cell, PW_Basis &pwb, double q, double l, double center, complex<double> *NG, int dim)
+{
+    // x dim
+    if (dim == 0)
+    {
+        double L = cell.a1[0];
+        q /= (cell.a2[1] * cell.a3[2] * l);
+        ModuleBase::GlobalFunc::ZEROS(NG, pwb.ngmc);
+        for (int ig = pwb.gstart; ig < pwb.ngmc; ig++)
+        {
+            double GX = pwb.get_G_cartesian_projection(ig, 0);
+            double GY = pwb.get_G_cartesian_projection(ig, 1);
+            double GZ = pwb.get_G_cartesian_projection(ig, 2);
+            GX = GX * 2 * ModuleBase::PI;
+            if (GY == 0 && GZ == 0 && GX != 0)
+            {
+                NG[ig] = exp(ModuleBase::NEG_IMAG_UNIT * GX * center) * complex<double>(2.0 * q * sin(GX * l / 2.0) / (L * GX), 0.0);
+            }
+        }
+        NG[0] = complex<double>(q * l / L, 0.0);
+    }
+    // y dim
+    else if (dim == 1)
+    {
+        double L = cell.a2[1];
+        q /= (cell.a1[0] * cell.a3[2] * l);
+        ModuleBase::GlobalFunc::ZEROS(NG, pwb.ngmc);
+        for (int ig = pwb.gstart; ig < pwb.ngmc; ig++)
+        {
+            double GX = pwb.get_G_cartesian_projection(ig, 0);
+            double GY = pwb.get_G_cartesian_projection(ig, 1);
+            double GZ = pwb.get_G_cartesian_projection(ig, 2);
+            GY = GY * 2 * ModuleBase::PI;
+            if (GX == 0 && GZ == 0 && GY != 0)
+            {
+                NG[ig] = exp(ModuleBase::NEG_IMAG_UNIT * GY * center) * complex<double>(2.0 * q * sin(GY * l / 2.0) / (L * GY), 0.0);
+            }
+        }
+        NG[0] = complex<double>(q * l / L, 0.0);
+    }
+    // z dim
+    else if (dim == 2)
+    {
+        double L = cell.a3[2];
+        q /= (cell.a1[0] * cell.a2[1] * l);
+        ModuleBase::GlobalFunc::ZEROS(NG, pwb.ngmc);
+        for (int ig = pwb.gstart; ig < pwb.ngmc; ig++)
+        {
+            double GX = pwb.get_G_cartesian_projection(ig, 0);
+            double GY = pwb.get_G_cartesian_projection(ig, 1);
+            double GZ = pwb.get_G_cartesian_projection(ig, 2);
+            GZ = GZ * 2 * ModuleBase::PI;
+            if (GX == 0 && GY == 0 && GZ != 0)
+            {
+                NG[ig] = exp(ModuleBase::NEG_IMAG_UNIT * GZ * center) * complex<double>(2.0 * q * sin(GZ * l / 2.0) / (L * GZ), 0.0);
+            }
+        }
+        NG[0] = complex<double>(q * l / L, 0.0);
+    }
+}
+
+void test_sum(const UnitCell &cell, PW_Basis &pwb, double* NR)
+{
+    double sum = 0.0;
+    for (int i=0; i < pwb.nrxx; i++)
+    {
+        sum += NR[i];
+    }
+    cout << sum * cell.omega / pwb.nrxx << endl;
+}
+
 void surchem::cal_totn(const UnitCell &cell, PW_Basis &pwb,
                        const complex<double> *Porter_g, complex<double> *N,
                        complex<double> *TOTN) {
+    /*
+    // test compensating charge
+    complex<double> *comp_reci = new complex<double>[pwb.ngmc];
+    double *comp_real = new double[pwb.nrxx];
+    add_comp_chg(cell, pwb, comp_q, comp_l, comp_center, comp_reci, comp_dim);
+    GlobalC::UFFT.ToRealSpace(comp_reci, comp_real);
+    for (int ir = 0; ir < pwb.nz; ir++)
+    {
+        cout << comp_real[ir] << endl;
+    }
+    cout << "sum" << endl;
+    test_sum(cell, pwb, comp_real);
+    delete[] comp_real;
+    delete[] comp_reci;
+    */
 
     // vloc to N
     complex<double> *vloc_g = new complex<double>[pwb.ngmc];
@@ -23,11 +109,13 @@ void surchem::cal_totn(const UnitCell &cell, PW_Basis &pwb,
 
     if (GlobalV::MY_RANK == 0) {
         N[0] = Porter_g[0];
+        // cout << "Ng[0]" << N[0] << endl;
     }
 
     for (int ig = 0; ig < pwb.ngmc; ig++) {
         TOTN[ig] = N[ig] - Porter_g[ig];
     }
+
     delete[] vloc_g;
     return;
 }
diff --git a/source/module_surchem/cal_vcav.cpp b/source/module_surchem/cal_vcav.cpp
@@ -136,7 +136,7 @@ void surchem::createcavity(const UnitCell &ucell, PW_Basis &pwb, const complex<d
     // cavitation energy
     //-------------------------------------------------------------
 
-    double Ael = cal_Acav(ucell, pwb);
+    // double Ael = cal_Acav(ucell, pwb);
 
     //  packs the real array into a complex one
     //  to G space
diff --git a/source/module_surchem/cal_vel.cpp b/source/module_surchem/cal_vel.cpp
@@ -106,7 +106,7 @@ ModuleBase::matrix surchem::cal_vel(const UnitCell &cell,
     }
 
     // calculate Ael
-    double Ael = cal_Ael(cell, pwb);
+    // double Ael = cal_Ael(cell, pwb);
 
     // the 2nd item of tmp_Vel
     eps_pot(PS_TOTN_real, Sol_phi, pwb, epsilon, epspot);
diff --git a/source/module_surchem/surchem.h b/source/module_surchem/surchem.h
@@ -24,6 +24,12 @@ class surchem
     ModuleBase::matrix Vel;
     double qs;
 
+    // compensating charge params
+    double comp_q;
+    double comp_l;
+    double comp_center;
+    int comp_dim;
+    // get atom info
     atom_in GetAtom;
 
     void allocate(const int &nrxx, const int &nspin);
diff --git a/source/src_io/write_input.cpp b/source/src_io/write_input.cpp
@@ -330,6 +330,11 @@ void Input::Print(const std::string &fn) const
     ModuleBase::GlobalFunc::OUTP(ofs, "tau", tau, "the effective surface tension parameter");
     ModuleBase::GlobalFunc::OUTP(ofs, "sigma_k", sigma_k, " the width of the diffuse cavity");
     ModuleBase::GlobalFunc::OUTP(ofs, "nc_k", nc_k, " the cut-off charge density");
+
+    ModuleBase::GlobalFunc::OUTP(ofs, "comp_q", comp_q, " total charge of compensating charge");
+    ModuleBase::GlobalFunc::OUTP(ofs, "comp_l", comp_l, " total length of compensating charge");
+    ModuleBase::GlobalFunc::OUTP(ofs, "comp_center", comp_center, " center of compensating charge on dim");
+    ModuleBase::GlobalFunc::OUTP(ofs, "comp_dim", comp_dim, " dimension of compensating charge(x, y or z)");
     ofs.close();
     return;
 }
diff --git a/source/src_pw/potential.cpp b/source/src_pw/potential.cpp
@@ -360,7 +360,7 @@ ModuleBase::matrix Potential::v_of_rho(const double *const *const rho_in, const
             v += GlobalC::solvent_model.v_correction(GlobalC::ucell, GlobalC::pw, GlobalV::NSPIN, rho_in);
 
             // test energy outside
-            cout << endl << "energy Outside: " << endl;
+            cout << "energy Outside: " << endl;
             GlobalC::solvent_model.cal_Ael(GlobalC::ucell, GlobalC::pw);
             GlobalC::solvent_model.cal_Acav(GlobalC::ucell, GlobalC::pw);
         }

Original file line number	Diff line number	Diff line change
`@@ -106,7 +106,7 @@ ModuleBase::matrix surchem::cal_vel(const UnitCell &cell,`
`106`	`106`	`}`
`107`	`107`
`108`	`108`	`// calculate Ael`
`109`		`- double Ael = cal_Ael(cell, pwb);`
	`109`	`+ // double Ael = cal_Ael(cell, pwb);`
`110`	`110`
`111`	`111`	`// the 2nd item of tmp_Vel`
`112`	`112`	`eps_pot(PS_TOTN_real, Sol_phi, pwb, epsilon, epspot);`
Original file line number	Diff line number	Diff line change
`@@ -360,7 +360,7 @@ ModuleBase::matrix Potential::v_of_rho(const double const const rho_in, const`
`360`	`360`	`v += GlobalC::solvent_model.v_correction(GlobalC::ucell, GlobalC::pw, GlobalV::NSPIN, rho_in);`
`361`	`361`
`362`	`362`	`// test energy outside`
`363`		`- cout << endl << "energy Outside: " << endl;`
	`363`	`+ cout << "energy Outside: " << endl;`
`364`	`364`	`GlobalC::solvent_model.cal_Ael(GlobalC::ucell, GlobalC::pw);`
`365`	`365`	`GlobalC::solvent_model.cal_Acav(GlobalC::ucell, GlobalC::pw);`
`366`	`366`	`}`