add calculation of Ael and Acav.

Author: ddhhss

source/module_surchem/cal_vcav.cpp

@@ -122,7 +122,7 @@ void surchem::createcavity(const UnitCell &ucell, PW_Basis &pwb, const complex<d
     // quantum surface area, integral of (gamma*A / n) * |\nabla n|
     //=term1 * sqrt_nablan_2
     //-------------------------------------------------------------
-    double qs = 0;
+    qs = 0;
 
     for (int ir = 0; ir < pwb.nrxx; ir++)
     {
@@ -136,7 +136,7 @@ void surchem::createcavity(const UnitCell &ucell, PW_Basis &pwb, const complex<d
     // cavitation energy
     //-------------------------------------------------------------
 
-    double Ael = surchem::cal_Acav(ucell, pwb, qs);
+    double Ael = cal_Acav(ucell, pwb);
 
     //  packs the real array into a complex one
     //  to G space
@@ -199,6 +199,7 @@ ModuleBase::matrix surchem::cal_vcav(const UnitCell &ucell, PW_Basis &pwb, const
             }
         }
     }
+
     delete[] tmp_Vcav;
     ModuleBase::timer::tick("surchem", "cal_vcav");
     return Vcav;

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, TOTN_real, delta_phi);
+    double Ael = cal_Ael(cell, pwb);
 
     // the 2nd item of tmp_Vel
     double *Vel2 = new double[pwb.nrxx];

source/module_surchem/corrected_energy.cpp

@@ -1,5 +1,6 @@
 #include "surchem.h"
 
+/*
 double surchem::cal_Ael(const UnitCell &cell, PW_Basis &pwb, const double *TOTN_real, const double *delta_phi_R)
 {
     double Ael = 0.0;
@@ -14,6 +15,30 @@ double surchem::cal_Ael(const UnitCell &cell, PW_Basis &pwb, const double *TOTN_
 }
 
 double surchem::cal_Acav(const UnitCell &cell, PW_Basis &pwb, double qs)
+{
+    double Acav = 0.0;
+    Acav = GlobalV::tau * qs;
+    Acav = Acav * cell.omega * 0.5 / pwb.nxyz;
+    Parallel_Reduce::reduce_double_pool(Acav);
+    cout << "Acav: " << Acav << endl;
+    return Acav;
+}
+*/
+
+double surchem::cal_Ael(const UnitCell &cell, PW_Basis &pwb)
+{
+    double Ael = 0.0;
+    for (int ir = 0; ir < pwb.nrxx; ir++)
+    {
+        Ael -= TOTN_real[ir] * delta_phi[ir];
+    }
+    Parallel_Reduce::reduce_double_pool(Ael);
+    Ael = Ael * cell.omega * 0.5 / pwb.nxyz;
+    cout << "Ael: " << Ael << endl;
+    return Ael;
+}
+
+double surchem::cal_Acav(const UnitCell &cell, PW_Basis &pwb)
 {
     double Acav = 0.0;
     Acav = GlobalV::tau * qs;

source/module_surchem/surchem.h

@@ -52,9 +52,9 @@ class surchem
                                       const complex<double> *PS_TOTN,
                                       int nspin);
 
-    double cal_Ael(const UnitCell &cell, PW_Basis &pwb, const double *TOTN_real, const double *delta_phi_R);
+    double cal_Ael(const UnitCell &cell, PW_Basis &pwb);
 
-    double cal_Acav(const UnitCell &cell, PW_Basis &pwb, double qs);
+    double cal_Acav(const UnitCell &cell, PW_Basis &pwb);
 
     void minimize_cg(const UnitCell &ucell,
                             PW_Basis &pwb,

source/src_pw/potential.cpp

@@ -363,9 +363,12 @@ ModuleBase::matrix Potential::v_of_rho(const double *const *const rho_in, const
         v += H_Hartree_pw::v_hartree(GlobalC::ucell, GlobalC::pw, GlobalV::NSPIN, rho_in);
         if (GlobalV::imp_sol)
         {
-            // surchem tmp;
-            // v += tmp.v_correction(GlobalC::ucell, GlobalC::pw, GlobalV::NSPIN, rho_in);
             v += GlobalC::solvent_model.v_correction(GlobalC::ucell, GlobalC::pw, GlobalV::NSPIN, rho_in);
+
+            // test energy outside
+            cout << endl << "energy Outside: " << endl;
+            GlobalC::solvent_model.cal_Ael(GlobalC::ucell, GlobalC::pw);
+            GlobalC::solvent_model.cal_Acav(GlobalC::ucell, GlobalC::pw);
         }
     }