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
1818template <>
1919void 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
174174template <>
@@ -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.psialpha ,
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