Skip to content
Merged
Show file tree
Hide file tree
Changes from 7 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
3 changes: 2 additions & 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,8 @@ 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)
UnitCell* cell_tmp = const_cast<UnitCell*>(&cell);
sf->setup_structure_factor(cell_tmp, 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