Skip to content

Commit 0339c1c

Browse files
author
wenfei-li
committed
gint : implement new method for gamma forcde & stress;
also fix a bug for old method
1 parent 27fcd05 commit 0339c1c

File tree

9 files changed

+339
-148
lines changed

9 files changed

+339
-148
lines changed

‎source/src_lcao/FORCE_gamma.cpp‎

Lines changed: 4 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -50,10 +50,10 @@ void Force_LCAO_gamma::ftable_gamma (
5050

5151
this->cal_ftvnl_dphi(loc.dm_gamma, isforce, isstress, ftvnl_dphi, stvnl_dphi);
5252
this->calFvnlDbeta(loc.dm_gamma, isforce, isstress, fvnl_dbeta, svnl_dbeta, GlobalV::vnl_method);
53-
this->cal_fvl_dphi(loc.dm_gamma, isforce, isstress, fvl_dphi, svl_dphi);
54-
55-
//this->cal_fvl_dphi_new(loc.DM, isforce, isstress, fvl_dphi, svl_dphi);
56-
53+
54+
//this->cal_fvl_dphi(loc.dm_gamma, isforce, isstress, fvl_dphi, svl_dphi);
55+
this->cal_fvl_dphi_new(loc.DM, isforce, isstress, fvl_dphi, svl_dphi);
56+
5757
//caoyu add for DeePKS
5858
#ifdef __DEEPKS
5959
if (GlobalV::deepks_scf)

‎source/src_lcao/FORCE_gamma_vl.cpp‎

Lines changed: 9 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -174,7 +174,15 @@ void Force_LCAO_gamma::cal_fvl_dphi(
174174
ModuleBase::GlobalFunc::ZEROS (this->UHM->LM->DHloc_fixed_x, this->ParaV->nloc);
175175
ModuleBase::GlobalFunc::ZEROS (this->UHM->LM->DHloc_fixed_y, this->ParaV->nloc);
176176
ModuleBase::GlobalFunc::ZEROS (this->UHM->LM->DHloc_fixed_z, this->ParaV->nloc);
177-
177+
if(GlobalV::CAL_STRESS)
178+
{
179+
ModuleBase::GlobalFunc::ZEROS (this->UHM->LM->DHloc_fixed_11, this->ParaV->nloc);
180+
ModuleBase::GlobalFunc::ZEROS (this->UHM->LM->DHloc_fixed_12, this->ParaV->nloc);
181+
ModuleBase::GlobalFunc::ZEROS (this->UHM->LM->DHloc_fixed_13, this->ParaV->nloc);
182+
ModuleBase::GlobalFunc::ZEROS (this->UHM->LM->DHloc_fixed_22, this->ParaV->nloc);
183+
ModuleBase::GlobalFunc::ZEROS (this->UHM->LM->DHloc_fixed_23, this->ParaV->nloc);
184+
ModuleBase::GlobalFunc::ZEROS (this->UHM->LM->DHloc_fixed_33, this->ParaV->nloc);
185+
}
178186
for(int ir=0; ir<GlobalC::pw.nrxx; ++ir)
179187
{
180188
GlobalC::pot.vr_eff1[ir] = GlobalC::pot.vr_eff(GlobalV::CURRENT_SPIN, ir);

‎source/src_lcao/FORCE_gamma_vl_new.cpp‎

Lines changed: 14 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -14,6 +14,7 @@ void Force_LCAO_gamma::cal_fvl_dphi_new(
1414
int istep = 1;
1515
GlobalC::pot.init_pot(istep, GlobalC::pw.strucFac);
1616
fvl_dphi.zero_out();
17+
svl_dphi.zero_out();
1718
for(int is=0; is<GlobalV::NSPIN; ++is)
1819
{
1920
GlobalV::CURRENT_SPIN = is;
@@ -22,7 +23,19 @@ void Force_LCAO_gamma::cal_fvl_dphi_new(
2223
GlobalC::pot.vr_eff1[ir] = GlobalC::pot.vr_eff(GlobalV::CURRENT_SPIN, ir);
2324
}
2425

25-
this->UHM->GG.cal_force_new(DM_in[is], GlobalC::pot.vr_eff1, fvl_dphi);
26+
this->UHM->GG.cal_force_new(DM_in[is], GlobalC::pot.vr_eff1, fvl_dphi, svl_dphi, isforce, isstress);
27+
}
28+
29+
if(isstress)
30+
{
31+
for(int i=0;i<3;i++)
32+
{
33+
for(int j=0;j<3;j++)
34+
{
35+
if(i<j) svl_dphi(j,i) = svl_dphi(i,j);
36+
svl_dphi(i,j) /= -GlobalC::ucell.omega;
37+
}
38+
}
2639
}
2740
ModuleBase::timer::tick("Force_LCAO_gamma","cal_fvl_dphi_new");
2841
}

‎source/src_lcao/gint_gamma.h‎

Lines changed: 23 additions & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -35,7 +35,9 @@ class Gint_Gamma : public Grid_Base_Beta
3535

3636
// (3) calcualte the forces related to grid
3737
void cal_force( const double*const vlocal);
38-
void cal_force_new(double** DM_in, const double*const vlocal, ModuleBase::matrix& force);
38+
void cal_force_new(double** DM_in, const double*const vlocal,
39+
ModuleBase::matrix& force, ModuleBase::matrix& stress,
40+
const bool is_force, const bool is_stress);
3941

4042
// (4) calcualte the envelope function
4143
void cal_env(const double* wfc, double* rho);
@@ -107,7 +109,9 @@ class Gint_Gamma : public Grid_Base_Beta
107109
// for calculatin of < dphi_i | Vlocal | phi_j > for foce calculation
108110
// on regular FFT real space grid.
109111
void gamma_force(const double*const vlocal) const;
110-
void gamma_force_new(const double*const*const DM, const double*const vlocal, ModuleBase::matrix& force);
112+
void gamma_force_new(const double*const*const DM, const double*const vlocal,
113+
ModuleBase::matrix& force, ModuleBase::matrix& stress,
114+
const bool is_force, const bool is_stress);
111115

112116
void cal_meshball_vlocal(
113117
const int na_grid, // how many atoms on this (i,j,k) grid
@@ -137,17 +141,28 @@ class Gint_Gamma : public Grid_Base_Beta
137141
void cal_meshball_force(
138142
const int grid_index,
139143
const int na_grid, // how many atoms on this (i,j,k) grid
140-
const int LD_pool,
141-
const int*const block_iw, // block_iw[na_grid], index of wave functions for each block
142144
const int*const block_size, // block_size[na_grid], number of columns of a band
143145
const int*const block_index, // block_index[na_grid+1], count total number of atomis orbitals
144-
const bool*const*const cal_flag, // cal_flag[GlobalC::pw.bxyz][na_grid], whether the atom-grid distance is larger than cutoff
145-
const double*const*const psir_vlbr3, // psir_vlbr3[GlobalC::pw.bxyz][LD_pool]
146+
const double*const*const psir_vlbr3_DM, // psir_vlbr3[GlobalC::pw.bxyz][LD_pool]
146147
const double*const*const dpsir_x, // psir_vlbr3[GlobalC::pw.bxyz][LD_pool]
147148
const double*const*const dpsir_y, // psir_vlbr3[GlobalC::pw.bxyz][LD_pool]
148149
const double*const*const dpsir_z, // psir_vlbr3[GlobalC::pw.bxyz][LD_pool]
149-
const double*const*const DM,
150-
ModuleBase::matrix &force);
150+
ModuleBase::matrix &force
151+
);
152+
153+
void cal_meshball_stress(
154+
const int na_grid, // how many atoms on this (i,j,k) grid
155+
const int*const block_index, // block_index[na_grid+1], count total number of atomis orbitals
156+
const double*const*const psir_vlbr3_DM,
157+
const double*const*const dpsir_xx,
158+
const double*const*const dpsir_xy,
159+
const double*const*const dpsir_xz,
160+
const double*const*const dpsir_yy,
161+
const double*const*const dpsir_yz,
162+
const double*const*const dpsir_zz,
163+
ModuleBase::matrix &stress
164+
);
165+
151166

152167
// extract the local potentials.
153168
// vldr3[GlobalC::pw.bxyz]

‎source/src_lcao/gint_gamma_fvl.cpp‎

Lines changed: 12 additions & 50 deletions
Original file line numberDiff line numberDiff line change
@@ -186,56 +186,18 @@ inline void cal_psir_ylm_dphi(
186186
const int ll = atom->iw2l[iw];
187187
const int idx_lm = atom->iw2_ylm[iw];
188188

189-
//special case for distance -> 0
190-
//Problems Remained
191-
//You have to add this two lines
192-
/*if (distance < 1e-9)
193-
{
194-
// if l==0 for localized orbital,
195-
// here is how we choose the 3D wave function psir_ylm,
196-
// and the derivative at r=0 point.
197-
if (ll == 0)
198-
{
199-
// psir_ylm is the three dimensional localized wave functions
200-
// (n,l,m), which is the multiply between 1D wave function 'tmp' and
201-
// spherical harmonic function rly.
202-
p_psir_ylm[iw] = tmp * rly[idx_lm];
203-
// the derivative of psir_ylm with respect to atom position R,
204-
// it's a std::vector, so it has x,y,z components.
205-
p_dphix[iw] = 0.0;
206-
p_dphiy[iw] = 0.0;
207-
p_dphiz[iw] = 0.0;
208-
}
209-
// if l>0, how do we choose 3D localized wave function
210-
// and the derivative.
211-
else
212-
{
213-
const Numerical_Orbital_Lm &philn = GlobalC::ORB.Phi[it].PhiLN(atom->iw2l[iw], atom->iw2n[iw]);
214-
215-
double Zty = philn.zty;
216-
p_psir_ylm[iw] = Zty * rly[idx_lm];
217-
p_dphix[iw] = Zty * grly[idx_lm][0];
218-
p_dphiy[iw] = Zty * grly[idx_lm][1];
219-
p_dphiz[iw] = Zty * grly[idx_lm][2];
220-
} // if (ll == 0)
221-
}*/
222-
// if r >0, here is how we choose the 3D wave function psir_ylm,
223-
// and the derivative at r=0 point.
224-
//else
225-
//{
226-
const double rl = pow(distance, ll);
227-
228-
// 3D wave functions
229-
p_psir_ylm[iw] = tmp * rly[idx_lm] / rl;
230-
231-
// derivative of wave functions with respect to atom positions.
232-
const double tmpdphi_rly = (dtmp - tmp * ll / distance) / rl * rly[idx_lm] / distance;
233-
const double tmprl = tmp/rl;
234-
235-
p_dphix[iw] = tmpdphi_rly * dr[0] + tmprl * grly[idx_lm][0];
236-
p_dphiy[iw] = tmpdphi_rly * dr[1] + tmprl * grly[idx_lm][1];
237-
p_dphiz[iw] = tmpdphi_rly * dr[2] + tmprl * grly[idx_lm][2];
238-
//}// if (distance < 1e-9)
189+
const double rl = pow(distance, ll);
190+
191+
// 3D wave functions
192+
p_psir_ylm[iw] = tmp * rly[idx_lm] / rl;
193+
194+
// derivative of wave functions with respect to atom positions.
195+
const double tmpdphi_rly = (dtmp - tmp * ll / distance) / rl * rly[idx_lm] / distance;
196+
const double tmprl = tmp/rl;
197+
198+
p_dphix[iw] = tmpdphi_rly * dr[0] + tmprl * grly[idx_lm][0];
199+
p_dphiy[iw] = tmpdphi_rly * dr[1] + tmprl * grly[idx_lm][1];
200+
p_dphiz[iw] = tmpdphi_rly * dr[2] + tmprl * grly[idx_lm][2];
239201
} // iw
240202
}// ib
241203
}//!id //finish loop of calc pre-info for each adjacent atom

0 commit comments

Comments
 (0)