Skip to content

Commit 40c7882

Browse files
A-006A-006
committed
change ucell in pwdft/stress_func_ewa.cpp
1 parent 405ac9b commit 40c7882

File tree

6 files changed

+29
-25
lines changed

6 files changed

+29
-25
lines changed

‎source/module_hamilt_lcao/hamilt_lcaodft/FORCE_STRESS.cpp‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1018,7 +1018,7 @@ void Force_Stress_LCAO<T>::calStressPwPart(const UnitCell& ucell,
10181018
//--------------------------------------------------------
10191019
// ewald stress: use plane wave only.
10201020
//--------------------------------------------------------
1021-
sc_pw.stress_ewa(sigmaewa, rhopw, 0); // remain problem
1021+
sc_pw.stress_ewa(ucell,sigmaewa, rhopw, 0); // remain problem
10221022

10231023
//--------------------------------------------------------
10241024
// stress due to core correlation.

‎source/module_hamilt_pw/hamilt_ofdft/of_stress_pw.cpp‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -66,7 +66,7 @@ void OF_Stress_PW::cal_stress(ModuleBase::matrix& sigmatot,
6666
stress_har(sigmahar, this->rhopw, true, pelec->charge);
6767

6868
// ewald contribution
69-
stress_ewa(sigmaewa, this->rhopw, true);
69+
stress_ewa(ucell,sigmaewa, this->rhopw, true);
7070

7171
// xc contribution: add gradient corrections(non diagonal)
7272
for (int i = 0; i < 3; i++)

‎source/module_hamilt_pw/hamilt_pwdft/stress_func.h‎

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -77,7 +77,8 @@ class Stress_Func
7777

7878
// 3) the stress from the ewald term (ion-ion intraction under
7979
// periodic boundary conditions).
80-
void stress_ewa(ModuleBase::matrix& sigma,
80+
void stress_ewa(const UnitCell& ucell,
81+
ModuleBase::matrix& sigma,
8182
ModulePW::PW_Basis* rho_basis,
8283
const bool is_pw); // ewald part in PW or LCAO basis
8384

‎source/module_hamilt_pw/hamilt_pwdft/stress_func_ewa.cpp‎

Lines changed: 23 additions & 20 deletions
Original file line numberDiff line numberDiff line change
@@ -11,15 +11,18 @@
1111

1212
//calcualte the Ewald stress term in PW and LCAO
1313
template<typename FPTYPE, typename Device>
14-
void Stress_Func<FPTYPE, Device>::stress_ewa(ModuleBase::matrix& sigma, ModulePW::PW_Basis* rho_basis, const bool is_pw)
14+
void Stress_Func<FPTYPE, Device>::stress_ewa(const UnitCell& ucell,
15+
ModuleBase::matrix& sigma,
16+
ModulePW::PW_Basis* rho_basis,
17+
const bool is_pw)
1518
{
1619
ModuleBase::TITLE("Stress_Func","stress_ewa");
1720
ModuleBase::timer::tick("Stress_Func","stress_ewa");
1821

1922
FPTYPE charge=0;
20-
for(int it=0; it < GlobalC::ucell.ntype; it++)
23+
for(int it=0; it < ucell.ntype; it++)
2124
{
22-
charge = charge + GlobalC::ucell.atoms[it].ncpp.zv * GlobalC::ucell.atoms[it].na;
25+
charge = charge + ucell.atoms[it].ncpp.zv * ucell.atoms[it].na;
2326
}
2427
//choose alpha in order to have convergence in the sum over G
2528
//upperbound is a safe upper bound for the error ON THE ENERGY
@@ -30,7 +33,7 @@ void Stress_Func<FPTYPE, Device>::stress_ewa(ModuleBase::matrix& sigma, ModulePW
3033
alpha-=0.1;
3134
if(alpha==0.0)
3235
ModuleBase::WARNING_QUIT("stres_ew", "optimal alpha not found");
33-
upperbound =ModuleBase::e2 * pow(charge,2) * sqrt( 2 * alpha / (ModuleBase::TWO_PI)) * erfc(sqrt(GlobalC::ucell.tpiba2 * rho_basis->ggecut / 4.0 / alpha));
36+
upperbound =ModuleBase::e2 * pow(charge,2) * sqrt( 2 * alpha / (ModuleBase::TWO_PI)) * erfc(sqrt(ucell.tpiba2 * rho_basis->ggecut / 4.0 / alpha));
3437
}
3538
while(upperbound>1e-7);
3639

@@ -40,7 +43,7 @@ void Stress_Func<FPTYPE, Device>::stress_ewa(ModuleBase::matrix& sigma, ModulePW
4043
const int ig0 = rho_basis->ig_gge0;
4144
if( ig0 >= 0)
4245
{
43-
sdewald = (ModuleBase::TWO_PI) * ModuleBase::e2 / 4.0 / alpha * pow(charge/GlobalC::ucell.omega,2);
46+
sdewald = (ModuleBase::TWO_PI) * ModuleBase::e2 / 4.0 / alpha * pow(charge/ucell.omega,2);
4447
}
4548
else
4649
{
@@ -79,27 +82,27 @@ void Stress_Func<FPTYPE, Device>::stress_ewa(ModuleBase::matrix& sigma, ModulePW
7982
for(; ig < ig_end; ig++)
8083
{
8184
if(ig == ig0) continue;
82-
g2 = rho_basis->gg[ig]* GlobalC::ucell.tpiba2;
85+
g2 = rho_basis->gg[ig]* ucell.tpiba2;
8386
g2a = g2 /4.0/alpha;
8487
rhostar=std::complex<FPTYPE>(0.0,0.0);
85-
for(int it=0; it < GlobalC::ucell.ntype; it++)
88+
for(int it=0; it < ucell.ntype; it++)
8689
{
87-
for(int i=0; i<GlobalC::ucell.atoms[it].na; i++)
90+
for(int i=0; i<ucell.atoms[it].na; i++)
8891
{
89-
arg = (rho_basis->gcar[ig] * GlobalC::ucell.atoms[it].tau[i]) * (ModuleBase::TWO_PI);
92+
arg = (rho_basis->gcar[ig] * ucell.atoms[it].tau[i]) * (ModuleBase::TWO_PI);
9093
FPTYPE sinp, cosp;
9194
ModuleBase::libm::sincos(arg, &sinp, &cosp);
92-
rhostar = rhostar + std::complex<FPTYPE>(GlobalC::ucell.atoms[it].ncpp.zv * cosp,GlobalC::ucell.atoms[it].ncpp.zv * sinp);
95+
rhostar = rhostar + std::complex<FPTYPE>(ucell.atoms[it].ncpp.zv * cosp,ucell.atoms[it].ncpp.zv * sinp);
9396
}
9497
}
95-
rhostar /= GlobalC::ucell.omega;
98+
rhostar /= ucell.omega;
9699
sewald = fact* (ModuleBase::TWO_PI) * ModuleBase::e2 * ModuleBase::libm::exp(-g2a) / g2 * pow(std::abs(rhostar),2);
97100
local_sdewald -= sewald;
98101
for(int l=0;l<3;l++)
99102
{
100103
for(int m=0;m<l+1;m++)
101104
{
102-
local_sigma(l, m) += sewald * GlobalC::ucell.tpiba2 * 2.0 * rho_basis->gcar[ig][l] * rho_basis->gcar[ig][m] / g2 * (g2a + 1);
105+
local_sigma(l, m) += sewald * ucell.tpiba2 * 2.0 * rho_basis->gcar[ig][l] * rho_basis->gcar[ig][m] / g2 * (g2a + 1);
103106
}
104107
}
105108
}
@@ -123,25 +126,25 @@ void Stress_Func<FPTYPE, Device>::stress_ewa(ModuleBase::matrix& sigma, ModulePW
123126

124127
FPTYPE sqa = sqrt(alpha);
125128
FPTYPE sq8a_2pi = sqrt(8 * alpha / (ModuleBase::TWO_PI));
126-
rmax = 4.0/sqa/GlobalC::ucell.lat0;
129+
rmax = 4.0/sqa/ucell.lat0;
127130

128131
// collapse it, ia, jt, ja loop into a single loop
129132
long long ijat, ijat_end;
130133
int it, i, jt, j;
131-
ModuleBase::TASK_DIST_1D(num_threads, thread_id, (long long)GlobalC::ucell.nat * GlobalC::ucell.nat, ijat, ijat_end);
134+
ModuleBase::TASK_DIST_1D(num_threads, thread_id, (long long)ucell.nat * ucell.nat, ijat, ijat_end);
132135
ijat_end = ijat + ijat_end;
133-
GlobalC::ucell.ijat2iaitjajt(ijat, &i, &it, &j, &jt);
136+
ucell.ijat2iaitjajt(ijat, &i, &it, &j, &jt);
134137

135138
while (ijat < ijat_end)
136139
{
137140
//calculate tau[na]-tau[nb]
138-
d_tau = GlobalC::ucell.atoms[it].tau[i] - GlobalC::ucell.atoms[jt].tau[j];
141+
d_tau = ucell.atoms[it].tau[i] - ucell.atoms[jt].tau[j];
139142
//generates nearest-neighbors shells
140-
H_Ewald_pw::rgen(d_tau, rmax, irr, GlobalC::ucell.latvec, GlobalC::ucell.G, r, r2, nrm);
143+
H_Ewald_pw::rgen(d_tau, rmax, irr, ucell.latvec, ucell.G, r, r2, nrm);
141144
for(int nr=0 ; nr<nrm ; nr++)
142145
{
143-
rr=sqrt(r2[nr]) * GlobalC::ucell.lat0;
144-
fac = -ModuleBase::e2/2.0/GlobalC::ucell.omega*pow(GlobalC::ucell.lat0,2)*GlobalC::ucell.atoms[it].ncpp.zv * GlobalC::ucell.atoms[jt].ncpp.zv / pow(rr,3) * (erfc(sqa*rr)+rr * sq8a_2pi * ModuleBase::libm::exp(-alpha * pow(rr,2)));
146+
rr=sqrt(r2[nr]) * ucell.lat0;
147+
fac = -ModuleBase::e2/2.0/ucell.omega*pow(ucell.lat0,2)*ucell.atoms[it].ncpp.zv * ucell.atoms[jt].ncpp.zv / pow(rr,3) * (erfc(sqa*rr)+rr * sq8a_2pi * ModuleBase::libm::exp(-alpha * pow(rr,2)));
145148
for(int l=0; l<3; l++)
146149
{
147150
for(int m=0; m<l+1; m++)
@@ -155,7 +158,7 @@ void Stress_Func<FPTYPE, Device>::stress_ewa(ModuleBase::matrix& sigma, ModulePW
155158
}//end nr
156159

157160
++ijat;
158-
GlobalC::ucell.step_jajtiait(&j, &jt, &i, &it);
161+
ucell.step_jajtiait(&j, &jt, &i, &it);
159162
}
160163

161164
delete[] r;

‎source/module_hamilt_pw/hamilt_pwdft/stress_pw.cpp‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -69,7 +69,7 @@ void Stress_PW<FPTYPE, Device>::cal_stress(ModuleBase::matrix& sigmatot,
6969
this->stress_har(sigmahar, rho_basis, 1, pelec->charge);
7070

7171
// ewald contribution
72-
this->stress_ewa(sigmaewa, rho_basis, 1);
72+
this->stress_ewa(ucell,sigmaewa, rho_basis, 1);
7373

7474
// xc contribution: add gradient corrections(non diagonal)
7575
for (int i = 0; i < 3; i++)

‎source/module_hamilt_pw/hamilt_stodft/sto_stress_pw.cpp‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -42,7 +42,7 @@ void Sto_Stress_PW<FPTYPE, Device>::cal_stress(ModuleBase::matrix& sigmatot,
4242
this->stress_har(sigmahar, rho_basis, true, chr);
4343

4444
// ewald contribution
45-
this->stress_ewa(sigmaewa, rho_basis, true);
45+
this->stress_ewa(ucell_in,sigmaewa, rho_basis, true);
4646

4747
// xc contribution: add gradient corrections(non diagonal)
4848
for (int i = 0; i < 3; ++i)

0 commit comments

Comments
 (0)