Skip to content

Commit 60a408a

Browse files
ErjieWuErjieWu
authored
Fix: Sign problem in DeePKS cal_pdm(). (#6082)
* Add HR precalc functions for DeePKS and fix some bugs. * Remove and combine some files in DeePKS. * Fix wrong function position. * Fix: sign problem in DeePKS.
1 parent 35448cb commit 60a408a

File tree

7 files changed

+54
-38
lines changed

7 files changed

+54
-38
lines changed

‎source/module_hamilt_lcao/module_deepks/deepks_force.cpp‎

Lines changed: 3 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -102,9 +102,9 @@ void DeePKS_domain::cal_f_delta(const std::vector<std::vector<TK>>& dm,
102102
int dRz = 0;
103103
if constexpr (std::is_same<TK, std::complex<double>>::value) // for multi-k
104104
{
105-
dRx = dR2.x - dR1.x;
106-
dRy = dR2.y - dR1.y;
107-
dRz = dR2.z - dR1.z;
105+
dRx = dR1.x - dR2.x;
106+
dRy = dR1.y - dR2.y;
107+
dRz = dR1.z - dR2.z;
108108
}
109109
ModuleBase::Vector3<double> dR(dRx, dRy, dRz);
110110

‎source/module_hamilt_lcao/module_deepks/deepks_fpre.cpp‎

Lines changed: 5 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -70,14 +70,15 @@ void DeePKS_domain::cal_gdmx(const int lmaxd,
7070
int dRz = 0;
7171
if constexpr (std::is_same<TK, std::complex<double>>::value)
7272
{
73-
dRx = (dR2 - dR1).x;
74-
dRy = (dR2 - dR1).y;
75-
dRz = (dR2 - dR1).z;
73+
dRx = (dR1 - dR2).x;
74+
dRy = (dR1 - dR2).y;
75+
dRz = (dR1 - dR2).z;
7676
}
7777
ModuleBase::Vector3<double> dR(dRx, dRy, dRz);
7878

7979
hamilt::AtomPair<double> dm_pair(ibt1, ibt2, dRx, dRy, dRz, &pv);
80-
dm_pair.allocate(nullptr, 1);
80+
dm_pair.allocate(nullptr, true);
81+
// not support nspin = 2 now
8182
for (int ik = 0; ik < nks; ik++)
8283
{
8384
TK kphase = TK(0);

‎source/module_hamilt_lcao/module_deepks/deepks_orbpre.cpp‎

Lines changed: 4 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -197,15 +197,15 @@ void DeePKS_domain::cal_orbital_precalc(const std::vector<TH>& dm_hl,
197197
{
198198
hamilt::AtomPair<double> dm_pair(ibt1,
199199
ibt2,
200-
(dR2 - dR1).x,
201-
(dR2 - dR1).y,
202-
(dR2 - dR1).z,
200+
(dR1 - dR2).x,
201+
(dR1 - dR2).y,
202+
(dR1 - dR2).z,
203203
&pv);
204204

205205
dm_pair.allocate(&dm_array[ik * row_size * col_size], 0);
206206

207207
const double arg
208-
= -(kvec_d[ik] * ModuleBase::Vector3<double>(dR2 - dR1)) * ModuleBase::TWO_PI;
208+
= -(kvec_d[ik] * ModuleBase::Vector3<double>(dR1 - dR2)) * ModuleBase::TWO_PI;
209209

210210
double sinp, cosp;
211211

‎source/module_hamilt_lcao/module_deepks/deepks_pdm.cpp‎

Lines changed: 33 additions & 18 deletions
Original file line numberDiff line numberDiff line change
@@ -273,36 +273,51 @@ void DeePKS_domain::cal_pdm(bool& init_pdm,
273273
}
274274
}
275275
// prepare DM from DMR
276-
std::vector<double> dm_array(row_size * col_size, 0.0);
277-
const double* dm_current = nullptr;
278276
int dRx = 0, dRy = 0, dRz = 0;
279277
if constexpr (std::is_same<TK, std::complex<double>>::value)
280278
{
281-
dRx = dR2.x - dR1.x;
282-
dRy = dR2.y - dR1.y;
283-
dRz = dR2.z - dR1.z;
279+
dRx = dR1.x - dR2.x;
280+
dRy = dR1.y - dR2.y;
281+
dRz = dR1.z - dR2.z;
284282
}
285-
// dm_k
283+
ModuleBase::Vector3<double> dR(dRx, dRy, dRz);
284+
285+
hamilt::AtomPair<double> dm_pair(ibt1, ibt2, dRx, dRy, dRz, &pv);
286+
dm_pair.allocate(nullptr, true);
286287
auto dm_k = dm->get_DMK_vector();
287-
const int nrow = pv.nrow;
288-
for (int ir = 0; ir < row_size; ir++)
288+
289+
if constexpr (std::is_same<TK, double>::value) // for gamma-only
290+
{
291+
for (int is = 0; is < dm_k.size(); is++)
292+
{
293+
if (ModuleBase::GlobalFunc::IS_COLUMN_MAJOR_KS_SOLVER(PARAM.inp.ks_solver))
294+
{
295+
dm_pair.add_from_matrix(dm_k[is].data(), pv.get_row_size(), 1.0, 1);
296+
}
297+
else
298+
{
299+
dm_pair.add_from_matrix(dm_k[is].data(), pv.get_col_size(), 1.0, 0);
300+
}
301+
}
302+
}
303+
else // for multi-k
289304
{
290-
for (int ic = 0; ic < col_size; ic++)
305+
for (int ik = 0; ik < dm_k.size(); ik++)
291306
{
292-
int iglob = (pv.atom_begin_row[ibt1] + ir) + nrow * (pv.atom_begin_col[ibt2] + ic);
293-
int iloc = ir * col_size + ic;
294-
std::complex<double> tmp = 0.0;
295-
for(int ik = 0; ik < dm_k.size(); ik++) // dm_k.size() == _nk * _nspin
307+
const double arg = -(kvec_d[ik] * dR) * ModuleBase::TWO_PI;
308+
const std::complex<double> kphase = std::complex<double>(cos(arg), sin(arg));
309+
if (ModuleBase::GlobalFunc::IS_COLUMN_MAJOR_KS_SOLVER(PARAM.inp.ks_solver))
310+
{
311+
dm_pair.add_from_matrix(dm_k[ik].data(), pv.get_row_size(), kphase, 1);
312+
}
313+
else
296314
{
297-
const double arg = (kvec_d[ik] * ModuleBase::Vector3<double>(dR1 - dR2)) * ModuleBase::TWO_PI;
298-
const std::complex<double> kphase = std::complex<double>(cos(arg), sin(arg));
299-
tmp += dm_k[ik][iglob] * kphase;
315+
dm_pair.add_from_matrix(dm_k[ik].data(), pv.get_col_size(), kphase, 0);
300316
}
301-
dm_array[iloc] += tmp.real();
302317
}
303318
}
304319

305-
dm_current = dm_array.data();
320+
const double* dm_current = dm_pair.get_pointer();
306321
// use s_2t and dm_current to get g_1dmt
307322
// dgemm_: C = alpha * A * B + beta * C
308323
// C = g_1dmt, A = dm_current, B = s_2t

‎source/module_hamilt_lcao/module_deepks/deepks_spre.cpp‎

Lines changed: 3 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -78,9 +78,9 @@ void DeePKS_domain::cal_gdmepsl(const int lmaxd,
7878
int dRz = 0;
7979
if constexpr (std::is_same<TK, std::complex<double>>::value)
8080
{
81-
dRx = (dR2 - dR1).x;
82-
dRy = (dR2 - dR1).y;
83-
dRz = (dR2 - dR1).z;
81+
dRx = (dR1 - dR2).x;
82+
dRy = (dR1 - dR2).y;
83+
dRz = (dR1 - dR2).z;
8484
}
8585
ModuleBase::Vector3<double> dR(dRx, dRy, dRz);
8686

‎source/module_hamilt_lcao/module_hcontainer/atom_pair.cpp‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -701,7 +701,7 @@ void AtomPair<T>::add_from_matrix(const std::complex<T>* hk,
701701
hk_real_pointer = (T*)hk_tmp;
702702
hk_imag_pointer = hk_real_pointer+1;
703703
BlasConnector::axpy(this->col_size, kphase.real(), hk_real_pointer, 2, hr_tmp, 1);
704-
BlasConnector::axpy(this->col_size, kphase.imag(), hk_imag_pointer, 2, hr_tmp, 1);
704+
BlasConnector::axpy(this->col_size, -kphase.imag(), hk_imag_pointer, 2, hr_tmp, 1);
705705
hk_tmp += ld_hk;
706706
hr_tmp += this->col_size;
707707
}
@@ -715,7 +715,7 @@ void AtomPair<T>::add_from_matrix(const std::complex<T>* hk,
715715
hk_real_pointer = (T*)hk_tmp;
716716
hk_imag_pointer = hk_real_pointer+1;
717717
BlasConnector::axpy(this->col_size, kphase.real(), hk_real_pointer, ld_hk_2, hr_tmp, 1);
718-
BlasConnector::axpy(this->col_size, kphase.imag(), hk_imag_pointer, ld_hk_2, hr_tmp, 1);
718+
BlasConnector::axpy(this->col_size, -kphase.imag(), hk_imag_pointer, ld_hk_2, hr_tmp, 1);
719719
hk_tmp ++;
720720
hr_tmp += this->col_size;
721721
}

‎tests/deepks/603_NO_deepks_H2O_multik/result.ref‎

Lines changed: 4 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -1,4 +1,4 @@
1-
etotref -466.8143178204656
2-
etotperatomref -155.6047726068
3-
totalforceref 10.244182
4-
totaltimeref 7.52
1+
etotref -466.8999964506086
2+
etotperatomref -155.6333321502
3+
totalforceref 10.047085
4+
totaltimeref 6.22

0 commit comments

Comments
 (0)