Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion source/module_elecstate/test/charge_extra_test.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -87,7 +87,7 @@ Structure_Factor::Structure_Factor()
Structure_Factor::~Structure_Factor()
{
}
void Structure_Factor::setup_structure_factor(UnitCell* Ucell, const ModulePW::PW_Basis* rho_basis)
void Structure_Factor::setup_structure_factor(const UnitCell* Ucell, const ModulePW::PW_Basis* rho_basis)
{
}

Expand Down
2 changes: 1 addition & 1 deletion source/module_hamilt_general/module_surchem/cal_pseudo.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -8,7 +8,7 @@ void surchem::gauss_charge(const UnitCell& cell,
complex<double>* N,
Structure_Factor* sf)
{
sf->setup_structure_factor(&GlobalC::ucell, rho_basis); // call strucFac(ntype,ngmc)
sf->setup_structure_factor(&cell, rho_basis); // call strucFac(ntype,ngmc)
for (int it = 0; it < cell.ntype; it++)
{
double RCS = GetAtom.atom_RCS[cell.atoms[it].ncpp.psd];
Expand Down
17 changes: 10 additions & 7 deletions source/module_hamilt_general/module_surchem/cal_vcav.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -2,7 +2,10 @@
#include "module_hamilt_general/module_xc/xc_functional.h"
#include "surchem.h"

void lapl_rho(const std::complex<double>* rhog, double* lapn, const ModulePW::PW_Basis* rho_basis)
void lapl_rho(const double& tpiba2,
const std::complex<double>* rhog,
double* lapn,
const ModulePW::PW_Basis* rho_basis)
{
std::complex<double> *gdrtmpg = new std::complex<double>[rho_basis->npw];
ModuleBase::GlobalFunc::ZEROS(lapn, rho_basis->nrxx);
Expand All @@ -22,9 +25,9 @@ void lapl_rho(const std::complex<double>* rhog, double* lapn, const ModulePW::PW
// bring the gdr from G --> R
rho_basis->recip2real(aux, aux);
// remember to multily 2pi/a0, which belongs to G vectors.
for (int ir = 0; ir < rho_basis->nrxx; ir++) {
lapn[ir] -= aux[ir].real() * GlobalC::ucell.tpiba2;
}
for (int ir = 0; ir < rho_basis->nrxx; ir++)
lapn[ir] -= aux[ir].real() * tpiba2;

}

delete[] gdrtmpg;
Expand Down Expand Up @@ -70,7 +73,7 @@ void surchem::createcavity(const UnitCell& ucell,
ModuleBase::GlobalFunc::ZEROS(lapn, rho_basis->nrxx);

// nabla n
XC_Functional::grad_rho(PS_TOTN, nablan, rho_basis, GlobalC::ucell.tpiba);
XC_Functional::grad_rho(PS_TOTN, nablan, rho_basis, ucell.tpiba);

// |\nabla n |^2 = nablan_2
for (int ir = 0; ir < rho_basis->nrxx; ir++)
Expand All @@ -79,7 +82,7 @@ void surchem::createcavity(const UnitCell& ucell,
}

// Laplacian of n
lapl_rho(PS_TOTN, lapn, rho_basis);
lapl_rho(ucell.tpiba2,PS_TOTN, lapn, rho_basis);

//-------------------------------------------------------------
// add -Lap(n)/|\nabla n| to vwork and copy \sqrt(|\nabla n|^2)
Expand Down Expand Up @@ -131,7 +134,7 @@ void surchem::createcavity(const UnitCell& ucell,

// \nabla(1 / |\nabla n|), ggn in real space
ModuleBase::Vector3<double> *ggn = new ModuleBase::Vector3<double>[rho_basis->nrxx];
XC_Functional::grad_rho(inv_gn, ggn, rho_basis, GlobalC::ucell.tpiba);
XC_Functional::grad_rho(inv_gn, ggn, rho_basis, ucell.tpiba);

//-------------------------------------------------------------
// add -(\nabla n . \nabla(1/ |\nabla n|)) to Vcav in real space
Expand Down
5 changes: 3 additions & 2 deletions source/module_hamilt_general/module_surchem/cal_vel.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -17,6 +17,7 @@ void shape_gradn(const double* PS_TOTN_real, const ModulePW::PW_Basis* rho_basis
}

void eps_pot(const double* PS_TOTN_real,
const double& tpiba,
const complex<double>* phi,
const ModulePW::PW_Basis* rho_basis,
double* d_eps,
Expand All @@ -36,7 +37,7 @@ void eps_pot(const double* PS_TOTN_real,
double *phisq = new double[rho_basis->nrxx];

// nabla phi
XC_Functional::grad_rho(phi, nabla_phi, rho_basis, GlobalC::ucell.tpiba);
XC_Functional::grad_rho(phi, nabla_phi, rho_basis, tpiba);

for (int ir = 0; ir < rho_basis->nrxx; ir++)
{
Expand Down Expand Up @@ -118,7 +119,7 @@ ModuleBase::matrix surchem::cal_vel(const UnitCell& cell,
this->Ael *= cell.omega / rho_basis->nxyz;

// the 2nd item of tmp_Vel
eps_pot(PS_TOTN_real, Sol_phi, rho_basis, epsilon, epspot);
eps_pot(PS_TOTN_real, cell.tpiba, Sol_phi, rho_basis, epsilon, epspot);

for (int i = 0; i < rho_basis->nrxx; i++)
{
Expand Down
20 changes: 10 additions & 10 deletions source/module_hamilt_general/module_surchem/sol_force.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -32,7 +32,7 @@ void force_cor_one(const UnitCell& cell,
vloc_at[ig] = vloc(it, rho_basis->ig2igg[ig]) * phase;
if(rho_basis->ig_gge0 == ig)
{
N[ig] = GlobalC::ucell.atoms[it].ncpp.zv / GlobalC::ucell.omega;
N[ig] = cell.atoms[it].ncpp.zv / cell.omega;
}
else
{
Expand All @@ -48,9 +48,9 @@ void force_cor_one(const UnitCell& cell,
forcesol(iat, 2) += rho_basis->gcar[ig][2] * imag(conj(delta_phi_g[ig]) * N[ig]);
}

forcesol(iat, 0) *= (GlobalC::ucell.tpiba * GlobalC::ucell.omega);
forcesol(iat, 1) *= (GlobalC::ucell.tpiba * GlobalC::ucell.omega);
forcesol(iat, 2) *= (GlobalC::ucell.tpiba * GlobalC::ucell.omega);
forcesol(iat, 0) *= (cell.tpiba * cell.omega);
forcesol(iat, 1) *= (cell.tpiba * cell.omega);
forcesol(iat, 2) *= (cell.tpiba * cell.omega);
//unit Ry/Bohr
forcesol(iat, 0) *= 2 ;
forcesol(iat, 1) *= 2 ;
Expand Down Expand Up @@ -128,9 +128,9 @@ void force_cor_two(const UnitCell& cell, const ModulePW::PW_Basis* rho_basis, Mo
forcesol(iat, 2) -= rho_basis->gcar[ig][2] * imag(conj(Vcav_g[ig]+Vel_g[ig]) * n_pseudo[ig]);
}

forcesol(iat, 0) *= (GlobalC::ucell.tpiba * GlobalC::ucell.omega);
forcesol(iat, 1) *= (GlobalC::ucell.tpiba * GlobalC::ucell.omega);
forcesol(iat, 2) *= (GlobalC::ucell.tpiba * GlobalC::ucell.omega);
forcesol(iat, 0) *= (cell.tpiba * cell.omega);
forcesol(iat, 1) *= (cell.tpiba * cell.omega);
forcesol(iat, 2) *= (cell.tpiba * cell.omega);
//eV/Ang
forcesol(iat, 0) *= 2 ;
forcesol(iat, 1) *= 2 ;
Expand All @@ -157,17 +157,17 @@ void surchem::cal_force_sol(const UnitCell& cell,
ModuleBase::TITLE("surchem", "cal_force_sol");
ModuleBase::timer::tick("surchem", "cal_force_sol");

int nat = GlobalC::ucell.nat;
int nat = cell.nat;
ModuleBase::matrix force1(nat, 3);
ModuleBase::matrix force2(nat, 3);

force_cor_one(cell, rho_basis, vloc, force1);
force_cor_two(cell, rho_basis,force2);

int iat = 0;
for (int it = 0;it < GlobalC::ucell.ntype;it++)
for (int it = 0;it < cell.ntype;it++)
{
for (int ia = 0;ia < GlobalC::ucell.atoms[it].na;ia++)
for (int ia = 0;ia < cell.atoms[it].na;ia++)
{
for(int ipol = 0; ipol < 3; ipol++)
{
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -26,6 +26,7 @@ class cal_pseudo_test : public testing::Test
{
protected:
surchem solvent_model;
UnitCell ucell;
};

TEST_F(cal_pseudo_test, gauss_charge)
Expand All @@ -41,7 +42,7 @@ TEST_F(cal_pseudo_test, gauss_charge)
double wfcecut;
bool gamma_only;

Setcell::setupcell(GlobalC::ucell);
Setcell::setupcell(ucell);

wfcecut = 80;
gamma_only = false;
Expand All @@ -60,7 +61,7 @@ TEST_F(cal_pseudo_test, gauss_charge)
#endif
// pwtest.initgrids(lat0,latvec,wfcecut);

GlobalC::rhopw->initgrids(GlobalC::ucell.lat0, GlobalC::ucell.latvec, wfcecut);
GlobalC::rhopw->initgrids(ucell.lat0, ucell.latvec, wfcecut);

GlobalC::rhopw->initparameters(gamma_only, wfcecut, distribution_type, xprime);
GlobalC::rhopw->setuptransform();
Expand All @@ -75,7 +76,7 @@ TEST_F(cal_pseudo_test, gauss_charge)
Structure_Factor sf;
sf.nbspline = -1;

solvent_model.gauss_charge(GlobalC::ucell, GlobalC::rhopw, N, &sf);
solvent_model.gauss_charge(ucell, GlobalC::rhopw, N, &sf);

EXPECT_NEAR(N[14].real(), 0.002, 1e-9);
EXPECT_NEAR(N[16].real(), -0.001573534, 1e-9);
Expand All @@ -88,7 +89,7 @@ TEST_F(cal_pseudo_test, cal_pseudo)
std::string precision_flag, device_flag;
precision_flag = "double";
device_flag = "cpu";

Setcell::setupcell(ucell);
ModulePW::PW_Basis pwtest(device_flag, precision_flag);
GlobalC::rhopw = &pwtest;
ModuleBase::Matrix3 latvec;
Expand All @@ -115,7 +116,7 @@ TEST_F(cal_pseudo_test, cal_pseudo)
#endif
// pwtest.initgrids(lat0,latvec,wfcecut);

GlobalC::rhopw->initgrids(GlobalC::ucell.lat0, GlobalC::ucell.latvec, wfcecut);
GlobalC::rhopw->initgrids(ucell.lat0, ucell.latvec, wfcecut);
GlobalC::rhopw->initparameters(gamma_only, wfcecut, distribution_type, xprime);
GlobalC::rhopw->setuptransform();
GlobalC::rhopw->collect_local_pw();
Expand All @@ -135,7 +136,7 @@ TEST_F(cal_pseudo_test, cal_pseudo)
}

complex<double>* PS_TOTN = new complex<double>[npw];
solvent_model.cal_pseudo(GlobalC::ucell, GlobalC::rhopw, Porter_g, PS_TOTN, &sf);
solvent_model.cal_pseudo(ucell, GlobalC::rhopw, Porter_g, PS_TOTN, &sf);

EXPECT_NEAR(PS_TOTN[16].real(), 0.098426466, 1e-9);
EXPECT_NEAR(PS_TOTN[14].real(), 0.102, 1e-9);
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -25,6 +25,7 @@ class cal_totn_test : public testing::Test
{
protected:
surchem solvent_model;
UnitCell ucell;
};

TEST_F(cal_totn_test, cal_totn)
Expand All @@ -40,7 +41,7 @@ TEST_F(cal_totn_test, cal_totn)
double wfcecut;
bool gamma_only;

Setcell::setupcell(GlobalC::ucell);
Setcell::setupcell(ucell);

wfcecut = 80;
gamma_only = false;
Expand All @@ -59,7 +60,7 @@ TEST_F(cal_totn_test, cal_totn)
#endif
// pwtest.initgrids(lat0,latvec,wfcecut);

GlobalC::rhopw->initgrids(GlobalC::ucell.lat0, GlobalC::ucell.latvec, wfcecut);
GlobalC::rhopw->initgrids(ucell.lat0, ucell.latvec, wfcecut);

GlobalC::rhopw->initparameters(gamma_only, wfcecut, distribution_type, xprime);
GlobalC::rhopw->setuptransform();
Expand Down Expand Up @@ -87,7 +88,7 @@ TEST_F(cal_totn_test, cal_totn)
vloc[i] = 0.1;
}

solvent_model.cal_totn(GlobalC::ucell, GlobalC::rhopw, Porter_g, N, TOTN, vloc);
solvent_model.cal_totn(ucell, GlobalC::rhopw, Porter_g, N, TOTN, vloc);

EXPECT_NEAR(TOTN[0].real(), -0.0999496256, 1e-10);
EXPECT_NEAR(TOTN[0].imag(), -1.299621928166352e-7, 1e-10);
Expand All @@ -103,7 +104,8 @@ TEST_F(cal_totn_test, induced_charge)
std::string precision_flag, device_flag;
precision_flag = "double";
device_flag = "cpu";

Setcell::setupcell(ucell);

ModulePW::PW_Basis pwtest(device_flag, precision_flag);
GlobalC::rhopw = &pwtest;
ModuleBase::Matrix3 latvec;
Expand All @@ -128,15 +130,15 @@ TEST_F(cal_totn_test, induced_charge)
#endif
// pwtest.initgrids(lat0,latvec,wfcecut);

GlobalC::rhopw->initgrids(GlobalC::ucell.lat0, GlobalC::ucell.latvec, wfcecut);
GlobalC::rhopw->initgrids(ucell.lat0, ucell.latvec, wfcecut);

GlobalC::rhopw->initparameters(gamma_only, wfcecut, distribution_type, xprime);
GlobalC::rhopw->setuptransform();
GlobalC::rhopw->collect_local_pw();
GlobalC::rhopw->collect_uniqgg();

double fac;
fac = ModuleBase::e2 * ModuleBase::FOUR_PI / (GlobalC::ucell.tpiba2 * GlobalC::rhopw->gg[0]);
fac = ModuleBase::e2 * ModuleBase::FOUR_PI / (ucell.tpiba2 * GlobalC::rhopw->gg[0]);
complex<double> delta_phi{-2.0347933860e-05, 4.5900395826e-07};
complex<double> induced_charge;
induced_charge = -delta_phi / fac;
Expand Down
19 changes: 10 additions & 9 deletions source/module_hamilt_general/module_surchem/test/cal_vcav_test.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -30,10 +30,11 @@ class cal_vcav_test : public testing::Test
{
protected:
surchem solvent_model;
UnitCell ucell;
};
TEST_F(cal_vcav_test, lapl_rho)
{
Setcell::setupcell(GlobalC::ucell);
Setcell::setupcell(ucell);

std::string precision_flag, device_flag;
precision_flag = "double";
Expand Down Expand Up @@ -61,7 +62,7 @@ TEST_F(cal_vcav_test, lapl_rho)
#ifdef __MPI
GlobalC::rhopw->initmpi(1, 0, POOL_WORLD);
#endif
GlobalC::rhopw->initgrids(GlobalC::ucell.lat0, GlobalC::ucell.latvec, wfcecut);
GlobalC::rhopw->initgrids(ucell.lat0, ucell.latvec, wfcecut);

GlobalC::rhopw->initparameters(gamma_only, wfcecut, distribution_type, xprime);
GlobalC::rhopw->setuptransform();
Expand Down Expand Up @@ -93,7 +94,7 @@ TEST_F(cal_vcav_test, lapl_rho)
GlobalC::rhopw->recip2real(aux, aux);
for (int ir = 0; ir < nrxx; ir++)
{
lapn[ir] -= aux[ir].real() * GlobalC::ucell.tpiba2;
lapn[ir] -= aux[ir].real() * ucell.tpiba2;
}
}

Expand Down Expand Up @@ -140,7 +141,7 @@ TEST_F(cal_vcav_test, shape_gradn)

TEST_F(cal_vcav_test, createcavity)
{
Setcell::setupcell(GlobalC::ucell);
Setcell::setupcell(ucell);

std::string precision_flag, device_flag;
precision_flag = "double";
Expand Down Expand Up @@ -168,7 +169,7 @@ TEST_F(cal_vcav_test, createcavity)
#ifdef __MPI
GlobalC::rhopw->initmpi(1, 0, POOL_WORLD);
#endif
GlobalC::rhopw->initgrids(GlobalC::ucell.lat0, GlobalC::ucell.latvec, wfcecut);
GlobalC::rhopw->initgrids(ucell.lat0, ucell.latvec, wfcecut);

GlobalC::rhopw->initparameters(gamma_only, wfcecut, distribution_type, xprime);
GlobalC::rhopw->setuptransform();
Expand All @@ -190,7 +191,7 @@ TEST_F(cal_vcav_test, createcavity)
PS_TOTN[ig] = 1e-7;
}

solvent_model.createcavity(GlobalC::ucell, GlobalC::rhopw, PS_TOTN, vwork);
solvent_model.createcavity(ucell, GlobalC::rhopw, PS_TOTN, vwork);

EXPECT_NEAR(vwork[0], 4.8556305312, 1e-10);
EXPECT_NEAR(vwork[1], -2.1006480538, 1e-10);
Expand All @@ -201,7 +202,7 @@ TEST_F(cal_vcav_test, createcavity)

TEST_F(cal_vcav_test, cal_vcav)
{
Setcell::setupcell(GlobalC::ucell);
Setcell::setupcell(ucell);

std::string precision_flag, device_flag;
precision_flag = "double";
Expand Down Expand Up @@ -229,7 +230,7 @@ TEST_F(cal_vcav_test, cal_vcav)
#ifdef __MPI
GlobalC::rhopw->initmpi(1, 0, POOL_WORLD);
#endif
GlobalC::rhopw->initgrids(GlobalC::ucell.lat0, GlobalC::ucell.latvec, wfcecut);
GlobalC::rhopw->initgrids(ucell.lat0, ucell.latvec, wfcecut);

GlobalC::rhopw->initparameters(gamma_only, wfcecut, distribution_type, xprime);
GlobalC::rhopw->setuptransform();
Expand All @@ -252,7 +253,7 @@ TEST_F(cal_vcav_test, cal_vcav)
int nspin = 2;
solvent_model.Vcav.create(nspin, nrxx);

solvent_model.cal_vcav(GlobalC::ucell, GlobalC::rhopw, PS_TOTN, nspin);
solvent_model.cal_vcav(ucell, GlobalC::rhopw, PS_TOTN, nspin);

EXPECT_NEAR(solvent_model.Vcav(0, 0), 4.8556305312, 1e-10);
EXPECT_NEAR(solvent_model.Vcav(0, 1), -2.1006480538, 1e-10);
Expand Down
Loading
Loading