Skip to content

Commit 990e415

Browse files
dyzhengdyzheng
committed
Fix: error of force and stress after refactor
1 parent 1d234a5 commit 990e415

5 files changed

Lines changed: 71 additions & 73 deletions

File tree

‎source/source_lcao/FORCE_STRESS.cpp‎

Lines changed: 17 additions & 29 deletions
Original file line numberDiff line numberDiff line change
@@ -204,50 +204,44 @@ void Force_Stress_LCAO<T>::getForceStress(UnitCell& ucell,
204204
nullptr, kv.kvec_d, nullptr, &ucell, orb.cutoffs(), &gd,
205205
two_center_bundle.overlap_orb_beta.get());
206206
tmp_nonlocal.cal_force_stress(isforce, isstress, dmR, fvnl_dbeta, svnl_dbeta);
207+
208+
// Switch back to spin channel 0
209+
if (PARAM.inp.nspin == 2)
210+
{
211+
dmat.dm->switch_dmr(0);
212+
}
207213

208214
// Calculate local potential force/stress (vl_dphi)
209215
// This uses grid integration, not operator-based method
210216
flk.ParaV = dmat.dm->get_paraV_pointer();
211217
PulayForceStress::cal_pulay_fs(fvl_dphi, svl_dphi, *dmat.dm, ucell, pelec->pot,
212218
isforce, isstress, false /*reset dm to gint*/);
213-
214-
// Switch back to spin channel 0
215-
if (PARAM.inp.nspin == 2)
216-
{
217-
dmat.dm->switch_dmr(0);
218-
}
219219
}
220220
else if (PARAM.inp.nspin == 4)
221221
{
222-
// For nspin=4 (non-collinear), need complex DMR
223-
// Create temporary complex DMR for DM
224-
hamilt::HContainer<std::complex<double>> tmp_dmr(dmat.dm->get_DMR_pointer(1)->get_paraV());
225-
std::vector<int> ijrs = dmat.dm->get_DMR_pointer(1)->get_ijr_info();
226-
tmp_dmr.insert_ijrs(&ijrs);
227-
tmp_dmr.allocate();
228-
dmat.dm->cal_DMR_full(&tmp_dmr);
229-
230-
// Create temporary complex DMR for EDM
231-
hamilt::HContainer<std::complex<double>> tmp_edmr(edm.get_DMR_pointer(1)->get_paraV());
232-
tmp_edmr.insert_ijrs(&ijrs);
233-
tmp_edmr.allocate();
234-
edm.cal_DMR_full(&tmp_edmr);
235222

236223
// Calculate kinetic force/stress (uses DM)
237224
if (PARAM.inp.t_in_h)
238225
{
239226
hamilt::EkineticNew<hamilt::OperatorLCAO<std::complex<double>, std::complex<double>>> tmp_ekinetic(
240227
nullptr, kv.kvec_d, nullptr, &ucell, orb.cutoffs(), &gd,
241228
two_center_bundle.kinetic_orb.get());
242-
tmp_ekinetic.cal_force_stress(isforce, isstress, &tmp_dmr, ftvnl_dphi, stvnl_dphi);
229+
tmp_ekinetic.cal_force_stress(isforce, isstress, dmat.dm->get_DMR_pointer(1), ftvnl_dphi, stvnl_dphi);
243230
}
244231

245232
// Calculate overlap force/stress (uses EDM)
246233
hamilt::OverlapNew<hamilt::OperatorLCAO<std::complex<double>, std::complex<double>>> tmp_overlap(
247234
nullptr, kv.kvec_d, nullptr, nullptr, &ucell, orb.cutoffs(), &gd,
248235
two_center_bundle.overlap_orb.get());
249-
tmp_overlap.cal_force_stress(isforce, isstress, &tmp_edmr, foverlap, soverlap);
236+
tmp_overlap.cal_force_stress(isforce, isstress, edm.get_DMR_pointer(1), foverlap, soverlap);
250237

238+
// For nspin=4 (non-collinear), need complex DMR
239+
// Create temporary complex DMR for DM
240+
hamilt::HContainer<std::complex<double>> tmp_dmr(dmat.dm->get_DMR_pointer(1)->get_paraV());
241+
std::vector<int> ijrs = dmat.dm->get_DMR_pointer(1)->get_ijr_info();
242+
tmp_dmr.insert_ijrs(&ijrs);
243+
tmp_dmr.allocate();
244+
dmat.dm->cal_DMR_full(&tmp_dmr);
251245
// Calculate nonlocal force/stress (uses DM)
252246
hamilt::NonlocalNew<hamilt::OperatorLCAO<std::complex<double>, std::complex<double>>> tmp_nonlocal(
253247
nullptr, kv.kvec_d, nullptr, &ucell, orb.cutoffs(), &gd,
@@ -263,18 +257,12 @@ void Force_Stress_LCAO<T>::getForceStress(UnitCell& ucell,
263257
// MPI reduction for forces
264258
if (isforce)
265259
{
266-
Parallel_Reduce::reduce_pool(foverlap.c, foverlap.nr * foverlap.nc);
267-
Parallel_Reduce::reduce_pool(ftvnl_dphi.c, ftvnl_dphi.nr * ftvnl_dphi.nc);
268-
Parallel_Reduce::reduce_pool(fvnl_dbeta.c, fvnl_dbeta.nr * fvnl_dbeta.nc);
269260
Parallel_Reduce::reduce_pool(fvl_dphi.c, fvl_dphi.nr * fvl_dphi.nc);
270261
}
271262

272263
// MPI reduction for stresses
273264
if (isstress)
274265
{
275-
Parallel_Reduce::reduce_pool(soverlap.c, soverlap.nr * soverlap.nc);
276-
Parallel_Reduce::reduce_pool(stvnl_dphi.c, stvnl_dphi.nr * stvnl_dphi.nc);
277-
Parallel_Reduce::reduce_pool(svnl_dbeta.c, svnl_dbeta.nr * svnl_dbeta.nc);
278266
Parallel_Reduce::reduce_pool(svl_dphi.c, svl_dphi.nr * svl_dphi.nc);
279267
}
280268

@@ -599,8 +587,8 @@ void Force_Stress_LCAO<T>::getForceStress(UnitCell& ucell,
599587
//-----------------------------
600588
// this->print_force("OVERLAP FORCE",foverlap,1,ry);
601589
ModuleIO::print_force(GlobalV::ofs_running, ucell, "OVERLAP FORCE", foverlap, false);
602-
// this->print_force("TVNL_DPHI force",ftvnl_dphi,PARAM.inp.test_force);
603-
// this->print_force("VNL_DBETA force",fvnl_dbeta,PARAM.inp.test_force);
590+
ModuleIO::print_force(GlobalV::ofs_running, ucell, "TVNL_DPHI force",ftvnl_dphi,false);
591+
ModuleIO::print_force(GlobalV::ofs_running, ucell, "VNL_DBETA force",fvnl_dbeta,false);
604592
// this->print_force("T_VNL FORCE",ftvnl,1,ry);
605593
ModuleIO::print_force(GlobalV::ofs_running, ucell, "T_VNL FORCE", ftvnl, false);
606594
ModuleIO::print_force(GlobalV::ofs_running, ucell, "VL_dPHI FORCE", fvl_dphi, false);

‎source/source_lcao/module_operator_lcao/ekinetic_force_stress.hpp‎

Lines changed: 26 additions & 21 deletions
Original file line numberDiff line numberDiff line change
@@ -22,7 +22,7 @@ inline double get_real<double>(const double& val)
2222
template <typename TK, typename TR>
2323
void EkineticNew<OperatorLCAO<TK, TR>>::cal_force_stress(const bool cal_force,
2424
const bool cal_stress,
25-
const HContainer<TR>* dmR,
25+
const HContainer<double>* dmR,
2626
ModuleBase::matrix& force,
2727
ModuleBase::matrix& stress)
2828
{
@@ -71,7 +71,7 @@ void EkineticNew<OperatorLCAO<TK, TR>>::cal_force_stress(const bool cal_force,
7171
}
7272

7373
// Find density matrix for this atom pair
74-
const hamilt::BaseMatrix<TR>* dm_matrix = dmR->find_matrix(iat1, iat2, R_index[0], R_index[1], R_index[2]);
74+
const hamilt::BaseMatrix<double>* dm_matrix = dmR->find_matrix(iat1, iat2, R_index[0], R_index[1], R_index[2]);
7575
if (dm_matrix == nullptr)
7676
{
7777
continue;
@@ -90,9 +90,17 @@ void EkineticNew<OperatorLCAO<TK, TR>>::cal_force_stress(const bool cal_force,
9090
continue;
9191
}
9292

93-
const TR* dm_pointer = dm_matrix->get_pointer();
93+
const double* dm_pointer = dm_matrix->get_pointer();
9494
double olm[4] = {0, 0, 0, 0}; // value, dx, dy, dz
9595

96+
// step_trace = 0 for npol=1; ={0, 1, col_size, col_size+1} for npol=2
97+
std::vector<int> step_trace(npol * npol, 0);
98+
if (npol == 2) {
99+
step_trace[1] = 1;
100+
step_trace[2] = col_indexes.size();
101+
step_trace[3] = col_indexes.size() + 1;
102+
}
103+
96104
// Loop over orbital pairs
97105
for (int iw1l = 0; iw1l < row_indexes.size(); iw1l += npol)
98106
{
@@ -115,34 +123,36 @@ void EkineticNew<OperatorLCAO<TK, TR>>::cal_force_stress(const bool cal_force,
115123
dtau * this->ucell->lat0,
116124
&olm[0], &olm[1]);
117125

126+
// only charge should be considered
127+
double dm_current = get_real(dm_pointer[0]);
128+
118129
// Calculate force contribution
119130
if (cal_force)
120131
{
121-
double dm_real = get_real(dm_pointer[iw2l]);
122-
123132
// F = -sum(dm * dT/dr)
124133
// Factor of 2 for Hermitian matrix will be applied later
125134
for (int i = 0; i < 3; i++)
126135
{
127-
force_tmp1[i] -= dm_real * olm[i + 1];
128-
force_tmp2[i] += dm_real * olm[i + 1];
136+
force_tmp1[i] += dm_current * olm[i + 1];
137+
force_tmp2[i] -= dm_current * olm[i + 1];
129138
}
130139
}
131140

132141
// Calculate stress contribution
133142
if (cal_stress)
134143
{
135-
double dm_real = get_real(dm_pointer[iw2l]);
136-
137144
// stress_ij = sum(dm * dT/dr_i * r_j)
138-
stress_local[0] += dm_real * olm[1] * dtau.x; // xx
139-
stress_local[1] += dm_real * olm[2] * dtau.y; // yy
140-
stress_local[2] += dm_real * olm[3] * dtau.z; // zz
141-
stress_local[3] += dm_real * olm[1] * dtau.y; // xy
142-
stress_local[4] += dm_real * olm[2] * dtau.z; // yz
143-
stress_local[5] += dm_real * olm[3] * dtau.x; // zx
145+
stress_local[0] -= dm_current * olm[1] * dtau.x; // xx
146+
stress_local[1] -= dm_current * olm[1] * dtau.y; // xy
147+
stress_local[2] -= dm_current * olm[1] * dtau.z; // xz
148+
stress_local[3] -= dm_current * olm[2] * dtau.y; // yy
149+
stress_local[4] -= dm_current * olm[2] * dtau.z; // yz
150+
stress_local[5] -= dm_current * olm[3] * dtau.z; // zz
144151
}
152+
153+
dm_pointer += npol;
145154
}
155+
dm_pointer += (npol - 1) * col_indexes.size();
146156
}
147157
}
148158
}
@@ -168,11 +178,6 @@ void EkineticNew<OperatorLCAO<TK, TR>>::cal_force_stress(const bool cal_force,
168178
#ifdef __MPI
169179
Parallel_Reduce::reduce_all(force.c, force.nr * force.nc);
170180
#endif
171-
// Factor of 2 for Hermitian matrix
172-
for (int i = 0; i < force.nr * force.nc; i++)
173-
{
174-
force.c[i] *= 2.0;
175-
}
176181
}
177182

178183
if (cal_stress)
@@ -183,7 +188,7 @@ void EkineticNew<OperatorLCAO<TK, TR>>::cal_force_stress(const bool cal_force,
183188
const double weight = this->ucell->lat0 / this->ucell->omega;
184189
for (int i = 0; i < 6; i++)
185190
{
186-
stress.c[i] = stress_tmp[i] * weight * 2.0; // Factor of 2 for Hermitian
191+
stress.c[i] = stress_tmp[i] * weight;
187192
}
188193
// Rearrange to 3x3 matrix format
189194
stress.c[8] = stress.c[5]; // stress(2,2)

‎source/source_lcao/module_operator_lcao/ekinetic_new.h‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -71,7 +71,7 @@ class EkineticNew<OperatorLCAO<TK, TR>> : public OperatorLCAO<TK, TR>
7171
*/
7272
void cal_force_stress(const bool cal_force,
7373
const bool cal_stress,
74-
const HContainer<TR>* dmR,
74+
const HContainer<double>* dmR,
7575
ModuleBase::matrix& force,
7676
ModuleBase::matrix& stress);
7777

‎source/source_lcao/module_operator_lcao/overlap_force_stress.hpp‎

Lines changed: 26 additions & 21 deletions
Original file line numberDiff line numberDiff line change
@@ -22,7 +22,7 @@ inline double get_real_overlap<double>(const double& val)
2222
template <typename TK, typename TR>
2323
void OverlapNew<OperatorLCAO<TK, TR>>::cal_force_stress(const bool cal_force,
2424
const bool cal_stress,
25-
const HContainer<TR>* dmR,
25+
const HContainer<double>* dmR,
2626
ModuleBase::matrix& force,
2727
ModuleBase::matrix& stress)
2828
{
@@ -71,7 +71,7 @@ void OverlapNew<OperatorLCAO<TK, TR>>::cal_force_stress(const bool cal_force,
7171
}
7272

7373
// Find density matrix for this atom pair
74-
const hamilt::BaseMatrix<TR>* dm_matrix = dmR->find_matrix(iat1, iat2, R_index[0], R_index[1], R_index[2]);
74+
const hamilt::BaseMatrix<double>* dm_matrix = dmR->find_matrix(iat1, iat2, R_index[0], R_index[1], R_index[2]);
7575
if (dm_matrix == nullptr)
7676
{
7777
continue;
@@ -90,9 +90,17 @@ void OverlapNew<OperatorLCAO<TK, TR>>::cal_force_stress(const bool cal_force,
9090
continue;
9191
}
9292

93-
const TR* dm_pointer = dm_matrix->get_pointer();
93+
const double* dm_pointer = dm_matrix->get_pointer();
9494
double olm[4] = {0, 0, 0, 0}; // value, dx, dy, dz
9595

96+
// step_trace = 0 for npol=1; ={0, 1, col_size, col_size+1} for npol=2
97+
std::vector<int> step_trace(npol * npol, 0);
98+
if (npol == 2) {
99+
step_trace[1] = 1;
100+
step_trace[2] = col_indexes.size();
101+
step_trace[3] = col_indexes.size() + 1;
102+
}
103+
96104
// Loop over orbital pairs
97105
for (int iw1l = 0; iw1l < row_indexes.size(); iw1l += npol)
98106
{
@@ -115,34 +123,36 @@ void OverlapNew<OperatorLCAO<TK, TR>>::cal_force_stress(const bool cal_force,
115123
dtau * this->ucell->lat0,
116124
&olm[0], &olm[1]);
117125

126+
// only charge should be considered
127+
double dm_current = get_real_overlap(dm_pointer[0]);
128+
118129
// Calculate force contribution
119130
if (cal_force)
120131
{
121-
double dm_real = get_real_overlap(dm_pointer[iw2l]);
122-
123132
// F = -sum(dm * dT/dr)
124133
// Factor of 2 for Hermitian matrix will be applied later
125134
for (int i = 0; i < 3; i++)
126135
{
127-
force_tmp1[i] -= dm_real * olm[i + 1];
128-
force_tmp2[i] += dm_real * olm[i + 1];
136+
force_tmp1[i] -= dm_current * olm[i + 1];
137+
force_tmp2[i] += dm_current * olm[i + 1];
129138
}
130139
}
131140

132141
// Calculate stress contribution
133142
if (cal_stress)
134143
{
135-
double dm_real = get_real_overlap(dm_pointer[iw2l]);
136-
137144
// stress_ij = sum(dm * dT/dr_i * r_j)
138-
stress_local[0] += dm_real * olm[1] * dtau.x; // xx
139-
stress_local[1] += dm_real * olm[2] * dtau.y; // yy
140-
stress_local[2] += dm_real * olm[3] * dtau.z; // zz
141-
stress_local[3] += dm_real * olm[1] * dtau.y; // xy
142-
stress_local[4] += dm_real * olm[2] * dtau.z; // yz
143-
stress_local[5] += dm_real * olm[3] * dtau.x; // zx
145+
stress_local[0] += dm_current * olm[1] * dtau.x; // xx
146+
stress_local[1] += dm_current * olm[1] * dtau.y; // xy
147+
stress_local[2] += dm_current * olm[1] * dtau.z; // xz
148+
stress_local[3] += dm_current * olm[2] * dtau.y; // yy
149+
stress_local[4] += dm_current * olm[2] * dtau.z; // yz
150+
stress_local[5] += dm_current * olm[3] * dtau.z; // zz
144151
}
152+
153+
dm_pointer += npol;
145154
}
155+
dm_pointer += (npol - 1) * col_indexes.size();
146156
}
147157
}
148158
}
@@ -168,11 +178,6 @@ void OverlapNew<OperatorLCAO<TK, TR>>::cal_force_stress(const bool cal_force,
168178
#ifdef __MPI
169179
Parallel_Reduce::reduce_all(force.c, force.nr * force.nc);
170180
#endif
171-
// Factor of 2 for Hermitian matrix
172-
for (int i = 0; i < force.nr * force.nc; i++)
173-
{
174-
force.c[i] *= 2.0;
175-
}
176181
}
177182

178183
if (cal_stress)
@@ -183,7 +188,7 @@ void OverlapNew<OperatorLCAO<TK, TR>>::cal_force_stress(const bool cal_force,
183188
const double weight = this->ucell->lat0 / this->ucell->omega;
184189
for (int i = 0; i < 6; i++)
185190
{
186-
stress.c[i] = stress_tmp[i] * weight * 2.0; // Factor of 2 for Hermitian
191+
stress.c[i] = stress_tmp[i] * weight;
187192
}
188193
// Rearrange to 3x3 matrix format
189194
stress.c[8] = stress.c[5]; // stress(2,2)

‎source/source_lcao/module_operator_lcao/overlap_new.h‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -62,7 +62,7 @@ class OverlapNew<OperatorLCAO<TK, TR>> : public OperatorLCAO<TK, TR>
6262
*/
6363
void cal_force_stress(const bool cal_force,
6464
const bool cal_stress,
65-
const HContainer<TR>* dmR,
65+
const HContainer<double>* dmR,
6666
ModuleBase::matrix& force,
6767
ModuleBase::matrix& stress);
6868

0 commit comments

Comments
 (0)