Skip to content

Commit 28df43d

Browse files
ErjieWuErjieWu
authored
Refactor: Replace nlm_save in DeePKS by HContainer object phialpha. (#5766)
* Add set_size() for BaseMatrix. * Refactor: replace nlm_save and nlm_save_k in DeePKS by psialpha. * Add cal_stress value check for deepks_out_unittest. * Modify the integrate test result.ref to reasonable results. * Change contributeHR in deepks_lcao into template. * Update functions used in deepks/test. * Change all 'psi' into 'phi' in DeePKS. * Update deepks_lcao.h
1 parent 0cdab49 commit 28df43d

38 files changed

+1962
-2251
lines changed

‎docs/advanced/input_files/input-main.md‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -2060,7 +2060,7 @@ Warning: this function is not robust enough for the current version. Please try
20602060
- **Type**: int
20612061
- **Availability**: numerical atomic orbital basis
20622062
- **Description**: Include V_delta label for DeePKS training. When `deepks_out_labels` is true and `deepks_v_delta` > 0, ABACUS will output h_base.npy, v_delta.npy and h_tot.npy(h_tot=h_base+v_delta).
2063-
Meanwhile, when `deepks_v_delta` equals 1, ABACUS will also output v_delta_precalc.npy, which is used to calculate V_delta during DeePKS training. However, when the number of atoms grows, the size of v_delta_precalc.npy will be very large. In this case, it's recommended to set `deepks_v_delta` as 2, and ABACUS will output psialpha.npy and grad_evdm.npy but not v_delta_precalc.npy. These two files are small and can be used to calculate v_delta_precalc in the procedure of training DeePKS.
2063+
Meanwhile, when `deepks_v_delta` equals 1, ABACUS will also output v_delta_precalc.npy, which is used to calculate V_delta during DeePKS training. However, when the number of atoms grows, the size of v_delta_precalc.npy will be very large. In this case, it's recommended to set `deepks_v_delta` as 2, and ABACUS will output phialpha.npy and grad_evdm.npy but not v_delta_precalc.npy. These two files are small and can be used to calculate v_delta_precalc in the procedure of training DeePKS.
20642064
- **Default**: 0
20652065

20662066
### deepks_out_unittest

‎source/Makefile.Objects‎

Lines changed: 2 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -189,13 +189,12 @@ OBJS_CELL=atom_pseudo.o\
189189
check_atomic_stru.o\
190190

191191
OBJS_DEEPKS=LCAO_deepks.o\
192-
deepks_fgamma.o\
193-
deepks_fk.o\
192+
deepks_force.o\
194193
LCAO_deepks_odelta.o\
195194
LCAO_deepks_io.o\
196195
LCAO_deepks_mpi.o\
197196
LCAO_deepks_pdm.o\
198-
LCAO_deepks_psialpha.o\
197+
LCAO_deepks_phialpha.o\
199198
LCAO_deepks_torch.o\
200199
LCAO_deepks_vdelta.o\
201200
deepks_hmat.o\

‎source/module_esolver/lcao_before_scf.cpp‎

Lines changed: 6 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -205,17 +205,19 @@ void ESolver_KS_LCAO<TK, TR>::before_scf(UnitCell& ucell, const int istep)
205205
}
206206

207207
#ifdef __DEEPKS
208-
// for each ionic step, the overlap <psi|alpha> must be rebuilt
208+
// for each ionic step, the overlap <phi|alpha> must be rebuilt
209209
// since it depends on ionic positions
210210
if (PARAM.globalv.deepks_setorb)
211211
{
212212
const Parallel_Orbitals* pv = &this->pv;
213-
// build and save <psi(0)|alpha(R)> at beginning
214-
GlobalC::ld.build_psialpha(PARAM.inp.cal_force, ucell, orb_, this->gd, *(two_center_bundle_.overlap_orb_alpha));
213+
// allocate <phi(0)|alpha(R)>, phialpha is different every ion step, so it is allocated here
214+
GlobalC::ld.allocate_phialpha(PARAM.inp.cal_force, ucell, orb_, this->gd);
215+
// build and save <phi(0)|alpha(R)> at beginning
216+
GlobalC::ld.build_phialpha(PARAM.inp.cal_force, ucell, orb_, this->gd, *(two_center_bundle_.overlap_orb_alpha));
215217

216218
if (PARAM.inp.deepks_out_unittest)
217219
{
218-
GlobalC::ld.check_psialpha(PARAM.inp.cal_force, ucell, orb_, this->gd);
220+
GlobalC::ld.check_phialpha(PARAM.inp.cal_force, ucell, orb_, this->gd);
219221
}
220222
}
221223
#endif

‎source/module_esolver/lcao_others.cpp‎

Lines changed: 6 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -211,17 +211,19 @@ void ESolver_KS_LCAO<TK, TR>::others(UnitCell& ucell, const int istep)
211211
}
212212

213213
#ifdef __DEEPKS
214-
// for each ionic step, the overlap <psi|alpha> must be rebuilt
214+
// for each ionic step, the overlap <phi|alpha> must be rebuilt
215215
// since it depends on ionic positions
216216
if (PARAM.globalv.deepks_setorb)
217217
{
218218
const Parallel_Orbitals* pv = &this->pv;
219-
// build and save <psi(0)|alpha(R)> at beginning
220-
GlobalC::ld.build_psialpha(PARAM.inp.cal_force, ucell, orb_, this->gd, *(two_center_bundle_.overlap_orb_alpha));
219+
// allocate <phi(0)|alpha(R)>, phialpha is different every ion step, so it is allocated here
220+
GlobalC::ld.allocate_phialpha(PARAM.inp.cal_force, ucell, orb_, this->gd);
221+
// build and save <phi(0)|alpha(R)> at beginning
222+
GlobalC::ld.build_phialpha(PARAM.inp.cal_force, ucell, orb_, this->gd, *(two_center_bundle_.overlap_orb_alpha));
221223

222224
if (PARAM.inp.deepks_out_unittest)
223225
{
224-
GlobalC::ld.check_psialpha(PARAM.inp.cal_force, ucell, orb_, this->gd);
226+
GlobalC::ld.check_phialpha(PARAM.inp.cal_force, ucell, orb_, this->gd);
225227
}
226228
}
227229
#endif

‎source/module_hamilt_lcao/hamilt_lcaodft/FORCE_STRESS.cpp‎

Lines changed: 11 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -12,12 +12,11 @@
1212
#include "module_hamilt_general/module_surchem/surchem.h" //sunml add 2022-08-10
1313
#include "module_hamilt_general/module_vdw/vdw.h"
1414
#include "module_parameter/parameter.h"
15-
#ifdef __DEEPKS
1615
#include "module_elecstate/elecstate_lcao.h"
16+
#ifdef __DEEPKS
1717
#include "module_hamilt_lcao/module_deepks/LCAO_deepks.h" //caoyu add for deepks 2021-06-03
1818
#include "module_hamilt_lcao/module_deepks/LCAO_deepks_io.h" // mohan add 2024-07-22
1919
#endif
20-
#include "module_elecstate/elecstate_lcao.h"
2120
#include "module_hamilt_lcao/hamilt_lcaodft/operator_lcao/dftu_lcao.h"
2221
#include "module_hamilt_lcao/hamilt_lcaodft/operator_lcao/dspin_lcao.h"
2322
#include "module_hamilt_lcao/hamilt_lcaodft/operator_lcao/nonlocal_new.h"
@@ -511,7 +510,14 @@ void Force_Stress_LCAO<T>::getForceStress(UnitCell& ucell,
511510
{
512511
const std::vector<std::vector<double>>& dm_gamma
513512
= dynamic_cast<const elecstate::ElecStateLCAO<double>*>(pelec)->get_DM()->get_DMK_vector();
514-
GlobalC::ld.cal_gdmx(dm_gamma, ucell, orb, gd, kv.get_nks(), kv.kvec_d, isstress);
513+
GlobalC::ld.cal_gdmx(dm_gamma,
514+
ucell,
515+
orb,
516+
gd,
517+
kv.get_nks(),
518+
kv.kvec_d,
519+
GlobalC::ld.phialpha,
520+
isstress);
515521
}
516522
else
517523
{
@@ -520,7 +526,8 @@ void Force_Stress_LCAO<T>::getForceStress(UnitCell& ucell,
520526
->get_DM()
521527
->get_DMK_vector();
522528

523-
GlobalC::ld.cal_gdmx(dm_k, ucell, orb, gd, kv.get_nks(), kv.kvec_d, isstress);
529+
GlobalC::ld
530+
.cal_gdmx(dm_k, ucell, orb, gd, kv.get_nks(), kv.kvec_d, GlobalC::ld.phialpha, isstress);
524531
}
525532
if (PARAM.inp.deepks_out_unittest)
526533
{

‎source/module_hamilt_lcao/hamilt_lcaodft/FORCE_gamma.cpp‎

Lines changed: 69 additions & 49 deletions
Original file line numberDiff line numberDiff line change
@@ -1,19 +1,19 @@
11
#include "FORCE.h"
22
#include "module_base/memory.h"
3-
#include "module_parameter/parameter.h"
43
#include "module_base/parallel_reduce.h"
54
#include "module_base/timer.h"
65
#include "module_cell/module_neighbor/sltk_grid_driver.h"
76
#include "module_hamilt_pw/hamilt_pwdft/global.h"
7+
#include "module_parameter/parameter.h"
88
#ifdef __DEEPKS
99
#include "module_hamilt_lcao/module_deepks/LCAO_deepks.h" //caoyu add for deepks on 20210813
1010
#include "module_hamilt_lcao/module_deepks/LCAO_deepks_io.h"
1111
#endif
1212
#include "module_cell/module_neighbor/sltk_grid_driver.h" //GridD
1313
#include "module_elecstate/elecstate_lcao.h"
1414
#include "module_hamilt_lcao/hamilt_lcaodft/LCAO_domain.h"
15-
#include "module_io/write_HS.h"
1615
#include "module_hamilt_lcao/hamilt_lcaodft/pulay_force_stress.h"
16+
#include "module_io/write_HS.h"
1717

1818
template <>
1919
void Force_LCAO<double>::allocate(const UnitCell& ucell,
@@ -147,28 +147,28 @@ void Force_LCAO<double>::finish_ftable(ForceStressArrays& fsr)
147147
return;
148148
}
149149

150-
//template <>
151-
//void Force_LCAO<double>::test(Parallel_Orbitals& pv, double* mm, const std::string& name)
150+
// template <>
151+
// void Force_LCAO<double>::test(Parallel_Orbitals& pv, double* mm, const std::string& name)
152152
//{
153-
// std::cout << "\n PRINT " << name << std::endl;
154-
// std::cout << std::setprecision(6) << std::endl;
155-
// for (int i = 0; i < PARAM.globalv.nlocal; i++)
156-
// {
157-
// for (int j = 0; j < PARAM.globalv.nlocal; j++)
158-
// {
159-
// if (std::abs(mm[i * PARAM.globalv.nlocal + j]) > 1.0e-5)
160-
// {
161-
// std::cout << std::setw(12) << mm[i * PARAM.globalv.nlocal + j];
162-
// }
163-
// else
164-
// {
165-
// std::cout << std::setw(12) << "0";
166-
// }
167-
// }
168-
// std::cout << std::endl;
169-
// }
170-
// return;
171-
//}
153+
// std::cout << "\n PRINT " << name << std::endl;
154+
// std::cout << std::setprecision(6) << std::endl;
155+
// for (int i = 0; i < PARAM.globalv.nlocal; i++)
156+
// {
157+
// for (int j = 0; j < PARAM.globalv.nlocal; j++)
158+
// {
159+
// if (std::abs(mm[i * PARAM.globalv.nlocal + j]) > 1.0e-5)
160+
// {
161+
// std::cout << std::setw(12) << mm[i * PARAM.globalv.nlocal + j];
162+
// }
163+
// else
164+
// {
165+
// std::cout << std::setw(12) << "0";
166+
// }
167+
// }
168+
// std::cout << std::endl;
169+
// }
170+
// return;
171+
// }
172172

173173
// be called in force_lo.cpp
174174
template <>
@@ -210,20 +210,40 @@ void Force_LCAO<double>::ftable(const bool isforce,
210210
// allocate DHloc_fixed_x, DHloc_fixed_y, DHloc_fixed_z
211211
this->allocate(ucell, gd, pv, fsr, two_center_bundle, orb);
212212

213-
const double* dSx[3] = { fsr.DSloc_x, fsr.DSloc_y, fsr.DSloc_z };
214-
const double* dSxy[6] = { fsr.DSloc_11, fsr.DSloc_12, fsr.DSloc_13, fsr.DSloc_22, fsr.DSloc_23, fsr.DSloc_33 };
213+
const double* dSx[3] = {fsr.DSloc_x, fsr.DSloc_y, fsr.DSloc_z};
214+
const double* dSxy[6] = {fsr.DSloc_11, fsr.DSloc_12, fsr.DSloc_13, fsr.DSloc_22, fsr.DSloc_23, fsr.DSloc_33};
215215
// calculate the force related to 'energy density matrix'.
216-
PulayForceStress::cal_pulay_fs(foverlap, soverlap,
216+
PulayForceStress::cal_pulay_fs(
217+
foverlap,
218+
soverlap,
217219
this->cal_edm(pelec, *psi, *dm, *kv, pv, PARAM.inp.nspin, PARAM.inp.nbands, ucell, *ra),
218-
ucell, pv, dSx, dSxy, isforce, isstress);
219-
220-
const double* dHx[3] = { fsr.DHloc_fixed_x, fsr.DHloc_fixed_y, fsr.DHloc_fixed_z };
221-
const double* dHxy[6] = { fsr.DHloc_fixed_11, fsr.DHloc_fixed_12, fsr.DHloc_fixed_13, fsr.DHloc_fixed_22, fsr.DHloc_fixed_23, fsr.DHloc_fixed_33 };
222-
//tvnl_dphi
220+
ucell,
221+
pv,
222+
dSx,
223+
dSxy,
224+
isforce,
225+
isstress);
226+
227+
const double* dHx[3] = {fsr.DHloc_fixed_x, fsr.DHloc_fixed_y, fsr.DHloc_fixed_z};
228+
const double* dHxy[6] = {fsr.DHloc_fixed_11,
229+
fsr.DHloc_fixed_12,
230+
fsr.DHloc_fixed_13,
231+
fsr.DHloc_fixed_22,
232+
fsr.DHloc_fixed_23,
233+
fsr.DHloc_fixed_33};
234+
// tvnl_dphi
223235
PulayForceStress::cal_pulay_fs(ftvnl_dphi, stvnl_dphi, *dm, ucell, pv, dHx, dHxy, isforce, isstress);
224236

225237
// vl_dphi
226-
PulayForceStress::cal_pulay_fs(fvl_dphi, svl_dphi, *dm, ucell, pelec->pot, gint, isforce, isstress, false/*reset dm to gint*/);
238+
PulayForceStress::cal_pulay_fs(fvl_dphi,
239+
svl_dphi,
240+
*dm,
241+
ucell,
242+
pelec->pot,
243+
gint,
244+
isforce,
245+
isstress,
246+
false /*reset dm to gint*/);
227247

228248
#ifdef __DEEPKS
229249
if (PARAM.inp.deepks_scf)
@@ -237,21 +257,21 @@ void Force_LCAO<double>::ftable(const bool isforce,
237257

238258
GlobalC::ld.cal_gedm(ucell.nat);
239259

240-
const int nks=1;
241-
DeePKS_domain::cal_f_delta_gamma(dm_gamma,
242-
ucell,
243-
orb,
244-
gd,
245-
*this->ParaV,
246-
GlobalC::ld.lmaxd,
247-
nks,
248-
kv->kvec_d,
249-
GlobalC::ld.nlm_save,
250-
GlobalC::ld.gedm,
251-
GlobalC::ld.inl_index,
252-
GlobalC::ld.F_delta,
253-
isstress,
254-
svnl_dalpha);
260+
const int nks = 1;
261+
DeePKS_domain::cal_f_delta<double>(dm_gamma,
262+
ucell,
263+
orb,
264+
gd,
265+
*this->ParaV,
266+
GlobalC::ld.lmaxd,
267+
nks,
268+
kv->kvec_d,
269+
GlobalC::ld.phialpha,
270+
GlobalC::ld.gedm,
271+
GlobalC::ld.inl_index,
272+
GlobalC::ld.F_delta,
273+
isstress,
274+
svnl_dalpha);
255275

256276
#ifdef __MPI
257277
Parallel_Reduce::reduce_all(GlobalC::ld.F_delta.c, GlobalC::ld.F_delta.nr * GlobalC::ld.F_delta.nc);
@@ -265,15 +285,15 @@ void Force_LCAO<double>::ftable(const bool isforce,
265285
if (PARAM.inp.deepks_out_unittest)
266286
{
267287
const int nks = 1; // 1 for gamma-only
268-
LCAO_deepks_io::print_dm(nks, PARAM.globalv.nlocal, this->ParaV->nrow, dm_gamma);
288+
LCAO_deepks_io::print_dm(nks, PARAM.globalv.nlocal, this->ParaV->nrow, dm_gamma);
269289

270290
GlobalC::ld.check_projected_dm();
271291

272292
GlobalC::ld.check_descriptor(ucell, PARAM.globalv.global_out_dir);
273293

274294
GlobalC::ld.check_gedm();
275295

276-
GlobalC::ld.cal_e_delta_band(dm_gamma,nks);
296+
GlobalC::ld.cal_e_delta_band(dm_gamma, nks);
277297

278298
std::ofstream ofs("E_delta_bands.dat");
279299
ofs << std::setprecision(10) << GlobalC::ld.e_delta_band;

0 commit comments

Comments
 (0)