add vcav, vcel.

Author: ddhhss

source/module_surchem/cal_vcav.cpp

@@ -176,17 +176,17 @@ ModuleBase::matrix surchem::cal_vcav(const UnitCell &ucell, PW_Basis &pwb, const
     ModuleBase::TITLE("surchem", "cal_vcav");
     ModuleBase::timer::tick("surchem", "cal_vcav");
 
-    double *Vcav = new double[pwb.nrxx];
-    ModuleBase::GlobalFunc::ZEROS(Vcav, pwb.nrxx);
+    double *tmp_Vcav = new double[pwb.nrxx];
+    ModuleBase::GlobalFunc::ZEROS(tmp_Vcav, pwb.nrxx);
 
-    createcavity(ucell, pwb, PS_TOTN, Vcav);
+    createcavity(ucell, pwb, PS_TOTN, tmp_Vcav);
 
-    ModuleBase::matrix v(nspin, pwb.nrxx);
+    ModuleBase::GlobalFunc::ZEROS(Vcav.c, nspin * pwb.nrxx);
     if (nspin == 4)
     {
         for (int ir = 0; ir < pwb.nrxx; ir++)
         {
-            v(0, ir) += Vcav[ir];
+            Vcav(0, ir) += tmp_Vcav[ir];
         }
     }
     else
@@ -195,11 +195,11 @@ ModuleBase::matrix surchem::cal_vcav(const UnitCell &ucell, PW_Basis &pwb, const
         {
             for (int ir = 0; ir < pwb.nrxx; ir++)
             {
-                v(is, ir) += Vcav[ir];
+                Vcav(is, ir) += tmp_Vcav[ir];
             }
         }
     }
-    delete[] Vcav;
+    delete[] tmp_Vcav;
     ModuleBase::timer::tick("surchem", "cal_vcav");
-    return v;
+    return Vcav;
 }

source/module_surchem/cal_vel.cpp

@@ -80,8 +80,8 @@ ModuleBase::matrix surchem::cal_vel(const UnitCell &cell,
     complex<double> *Sol_phi0 = new complex<double>[pwb.ngmc];
     int ncgsol = 0;
 
-    double *Vel = new double[pwb.nrxx];
-    ModuleBase::GlobalFunc::ZEROS(Vel, pwb.nrxx);
+    double *tmp_Vel = new double[pwb.nrxx];
+    ModuleBase::GlobalFunc::ZEROS(tmp_Vel, pwb.nrxx);
 
     // Calculate Sol_phi with epsilon.
     ncgsol = 0;
@@ -98,34 +98,35 @@ ModuleBase::matrix surchem::cal_vel(const UnitCell &cell,
     GlobalC::UFFT.ToRealSpace(Sol_phi, phi_tilda_R);
     GlobalC::UFFT.ToRealSpace(Sol_phi0, phi_tilda_R0);
 
-    // the 1st item of Vel
+    // the 1st item of tmp_Vel
     for (int i = 0; i < pwb.nrxx; i++)
     {
         delta_phi[i] = phi_tilda_R[i] - phi_tilda_R0[i];
-        Vel[i] += delta_phi[i];
+        tmp_Vel[i] += delta_phi[i];
     }
 
     // calculate Ael
     double Ael = cal_Ael(cell, pwb, TOTN_real, delta_phi);
 
-    // the 2nd item of Vel
+    // the 2nd item of tmp_Vel
     double *Vel2 = new double[pwb.nrxx];
     ModuleBase::GlobalFunc::ZEROS(Vel2, pwb.nrxx);
 
     eps_pot(PS_TOTN_real, Sol_phi, pwb, epsilon, Vel2);
 
     for (int i = 0; i < pwb.nrxx; i++)
     {
-        Vel[i] += Vel2[i];
+        tmp_Vel[i] += Vel2[i];
     }
 
-    ModuleBase::matrix v(nspin, pwb.nrxx);
+    // ModuleBase::matrix v(nspin, pwb.nrxx);
+    ModuleBase::GlobalFunc::ZEROS(Vel.c, nspin * pwb.nrxx);
 
     if (nspin == 4)
     {
         for (int ir = 0; ir < pwb.nrxx; ir++)
         {
-            v(0, ir) += Vel[ir];
+            Vel(0, ir) += tmp_Vel[ir];
         }
     }
     else
@@ -134,7 +135,7 @@ ModuleBase::matrix surchem::cal_vel(const UnitCell &cell,
         {
             for (int ir = 0; ir < pwb.nrxx; ir++)
             {
-                v(is, ir) += Vel[ir];
+                Vel(is, ir) += tmp_Vel[ir];
             }
         }
     }
@@ -145,13 +146,13 @@ ModuleBase::matrix surchem::cal_vel(const UnitCell &cell,
     delete[] B;
     delete[] epsilon;
     delete[] epsilon0;
-    delete[] Vel;
+    delete[] tmp_Vel;
     delete[] Vel2;
     // delete[] TOTN_real;
     delete[] phi_tilda_R;
     delete[] phi_tilda_R0;
     // delete[] delta_phi_R;
 
     ModuleBase::timer::tick("surchem", "cal_vel");
-    return v;
+    return Vel;
 }

source/module_surchem/surchem.cpp

@@ -9,14 +9,27 @@ surchem::surchem()
 {
     TOTN_real = new double[1];
     delta_phi = new double[1];
+    Vcav = ModuleBase::matrix();
+    Vel = ModuleBase::matrix();
+    qs = 0;
 }
 
-void surchem::allocate(const int &nrxx)
+void surchem::allocate(const int &nrxx, const int &nspin)
 {
+    assert(nrxx > 0);
+    assert(nspin > 0);
+
     delete[] TOTN_real;
     delete[] delta_phi;
+
     TOTN_real = new double[nrxx];
     delta_phi = new double[nrxx];
+    Vcav.create(nspin, nrxx);
+    Vel.create(nspin, nrxx);
+
+    ModuleBase::GlobalFunc::ZEROS(delta_phi, nrxx);
+    ModuleBase::GlobalFunc::ZEROS(TOTN_real, nrxx);
+    return;
 }
 
 surchem::~surchem()

source/module_surchem/surchem.h

@@ -19,10 +19,13 @@ class surchem
 
     double *TOTN_real;
     double *delta_phi;
+    ModuleBase::matrix Vcav;
+    ModuleBase::matrix Vel;
+    double qs;
 
     atom_in GetAtom;
 
-    void allocate(const int &nrxx);
+    void allocate(const int &nrxx, const int &nspin);
 
     void cal_epsilon(PW_Basis &pwb, const double *PS_TOTN_real, double *epsilon, double *epsilon0);
 

source/src_pw/potential.cpp

@@ -63,7 +63,7 @@ void Potential::allocate(const int nrxx)
 
     if (GlobalV::imp_sol)
     {
-        GlobalC::solvent_model.allocate(nrxx);
+        GlobalC::solvent_model.allocate(nrxx, GlobalV::NSPIN);
     }
 
     return;