Skip to content

Commit 30db3a4

Browse files
ESROAMERESROAMER
authored
Force calculation for hybrid gauge rt-tddft (deepmodeling#7474)
* hybrid gauge tddft force * fix no MPI compile * pass td_stype * add warning quit
1 parent 5a88cc6 commit 30db3a4

17 files changed

Lines changed: 804 additions & 10 deletions

‎source/source_esolver/esolver_ks_lcao.cpp‎

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -246,7 +246,8 @@ void ESolver_KS_LCAO<TK, TR>::cal_force(UnitCell& ucell, ModuleBase::matrix& for
246246
two_center_bundle_, orb_, force, this->scs,
247247
this->locpp, this->sf, this->kv,
248248
this->pw_rho, this->solvent, this->dftu, this->deepks,
249-
this->exx_nao, &ucell.symm);
249+
this->exx_nao, &ucell.symm, PARAM.inp.td_stype,
250+
static_cast<hamilt::Hamilt<TK>*>(this->p_hamilt));
250251

251252
// delete RA after cal_force
252253
this->RA.delete_grid();

‎source/source_esolver/esolver_ks_lcao_tddft.cpp‎

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -106,6 +106,7 @@ void ESolver_KS_LCAO_TDDFT<TR, Device>::runner(UnitCell& ucell, const int istep)
106106
//----------------------------------------------------------------
107107
this->before_scf(ucell, istep); // From ESolver_KS_LCAO
108108
td_p->initialize_phase_hybrid(ucell, dynamic_cast<hamilt::HamiltLCAO<std::complex<double>, TR>*>(this->p_hamilt)->getHR());
109+
td_p->calculate_grad_overlap(this->pv, ucell, this->gd, this->orb_.cutoffs(), this->two_center_bundle_.overlap_orb.get());
109110
// Initialize the moving spatial gauge
110111
if (use_td_moving_gauge && this->td_mg_ == nullptr)
111112
{
@@ -392,7 +393,7 @@ void ESolver_KS_LCAO_TDDFT<TR, Device>::iter_finish(UnitCell& ucell,
392393

393394
// Calculate energy-density matrix for RT-TDDFT
394395
if (conv_esolver && estep == estep_max - 1 && istep >= (PARAM.inp.init_wfc == "file" ? 0 : 1)
395-
&& PARAM.inp.td_edm == 0)
396+
&& PARAM.inp.td_edm == 0 && PARAM.inp.td_stype != 2)
396397
{
397398
if (use_tensor && use_lapack)
398399
{

‎source/source_estate/module_pot/H_TDDFT_pw.cpp‎

Lines changed: 4 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -107,6 +107,7 @@ void H_TDDFT_pw::cal_fixed_v(double* vl_pseudo)
107107
// time evolve
108108
H_TDDFT_pw::istep++;
109109
H_TDDFT_pw::istep_int = istep;
110+
global_vext_time = {0.0, 0.0, 0.0};
110111

111112
// judgement to skip vext
112113
if (!PARAM.inp.td_vext || istep > tend || istep < tstart)
@@ -122,7 +123,7 @@ void H_TDDFT_pw::cal_fixed_v(double* vl_pseudo)
122123
trigo_count = 0;
123124
heavi_count = 0;
124125

125-
global_vext_time = {0.0, 0.0, 0.0};
126+
126127

127128
for (auto direc: PARAM.inp.td_vext_dire)
128129
{
@@ -270,6 +271,7 @@ void H_TDDFT_pw::update_At()
270271
At = At + At_laststep / 2.0;
271272
At_laststep.set(0.0, 0.0, 0.0);
272273
Et.set(0.0, 0.0, 0.0);
274+
global_vext_time = { 0.0, 0.0, 0.0 };
273275

274276
// judgement to skip vext
275277
if (!PARAM.inp.td_vext || istep > tend || istep < tstart)
@@ -341,6 +343,7 @@ void H_TDDFT_pw::update_At()
341343
count++;
342344
}
343345
At = At + At_laststep / 2.0;
346+
if(stype==2)global_vext_time = { Et[0],Et[1],Et[2] };
344347

345348
ModuleBase::timer::end("H_TDDFT_pw", "update_At");
346349
return;

‎source/source_io/module_hs/cal_r_overlap_R.cpp‎

Lines changed: 40 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -516,7 +516,47 @@ ModuleBase::Vector3<double> cal_r_overlap_R::get_psi_r_psi(const ModuleBase::Vec
516516

517517
return temp_prp;
518518
}
519+
ModuleBase::Vector3<double> cal_r_overlap_R::get_psi_r_gradpsi(const ModuleBase::Vector3<double>& R1,
520+
const int& T1,
521+
const int& L1,
522+
const int& m1,
523+
const int& N1,
524+
const ModuleBase::Vector3<double>& R2,
525+
const int& T2,
526+
const int& L2,
527+
const int& m2,
528+
const int& N2,
529+
const ModuleBase::Vector3<double>& Efield,
530+
const ModuleBase::Vector3<double>& dR)
531+
{
532+
ModuleBase::Vector3<double> origin_point(0.0, 0.0, 0.0);
533+
double factor = sqrt(ModuleBase::FOUR_PI / 3.0);
534+
const ModuleBase::Vector3<double>& distance = R2 - R1;
535+
536+
ModuleBase::Vector3<double> grad_o = center2_orb11[T1][T2][L1][N1][L2].at(N2).cal_grad_overlap(origin_point, distance, m1, m2);
537+
538+
ModuleBase::Vector3<double> grad_rx = -1 * factor * center2_orb21_r[T1][T2][L1][N1][L2].at(N2).cal_grad_overlap(origin_point,
539+
distance,
540+
m1,
541+
1,
542+
m2); // m = 1
543+
544+
ModuleBase::Vector3<double> grad_ry = -1 * factor * center2_orb21_r[T1][T2][L1][N1][L2].at(N2).cal_grad_overlap(origin_point,
545+
distance,
546+
m1,
547+
2,
548+
m2); // m = -1
519549

550+
ModuleBase::Vector3<double> grad_rz = factor * center2_orb21_r[T1][T2][L1][N1][L2].at(N2).cal_grad_overlap(origin_point,
551+
distance,
552+
m1,
553+
0,
554+
m2); // m = 0
555+
556+
ModuleBase::Vector3<double> temp_prp = Efield[0] * grad_rx + Efield[1] * grad_ry + Efield[2] * grad_rz + (Efield*(R1-dR)) * grad_o;
557+
558+
return temp_prp;
559+
}
520560
void cal_r_overlap_R::get_psi_r_beta(const UnitCell& ucell,
521561
std::vector<std::vector<double>>& nlm,
522562
const ModuleBase::Vector3<double>& R1,

‎source/source_io/module_hs/cal_r_overlap_R.h‎

Lines changed: 14 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -45,6 +45,20 @@ class cal_r_overlap_R
4545
const int& m2,
4646
const int& N2
4747
);
48+
ModuleBase::Vector3<double> get_psi_r_gradpsi(
49+
const ModuleBase::Vector3<double>& R1,
50+
const int& T1,
51+
const int& L1,
52+
const int& m1,
53+
const int& N1,
54+
const ModuleBase::Vector3<double>& R2,
55+
const int& T2,
56+
const int& L2,
57+
const int& m2,
58+
const int& N2,
59+
const ModuleBase::Vector3<double>& Efield,
60+
const ModuleBase::Vector3<double>& dR
61+
);
4862
void get_psi_r_beta(
4963
const UnitCell& ucell,
5064
std::vector<std::vector<double>>& nlm,

‎source/source_lcao/FORCE_STRESS.cpp‎

Lines changed: 27 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -24,7 +24,9 @@
2424
#include "source_lcao/module_operator_lcao/nonlocal.h"
2525
#include "source_lcao/module_operator_lcao/ekinetic.h"
2626
#include "source_lcao/module_operator_lcao/overlap.h"
27+
#include "source_lcao/module_operator_lcao/td_pot_hybrid.h"
2728
#include "source_lcao/pulay_fs.h"
29+
#include "source_lcao/module_rt/force_rt_overlap.h"
2830

2931

3032
// mohan add 2025-11-04
@@ -85,7 +87,9 @@ void Force_Stress_LCAO<T>::getForceStress(UnitCell& ucell,
8587
Plus_U &dftu, // mohan add 2025-11-07
8688
Setup_DeePKS<T>& deepks,
8789
Exx_NAO<T> &exx_nao,
88-
ModuleSymmetry::Symmetry* symm)
90+
ModuleSymmetry::Symmetry* symm,
91+
const int td_stype,
92+
hamilt::Hamilt<T>* p_hamilt)
8993
{
9094
ModuleBase::TITLE("Force_Stress_LCAO", "getForceStress");
9195
ModuleBase::timer::start("Force_Stress_LCAO", "getForceStress");
@@ -113,6 +117,7 @@ void Force_Stress_LCAO<T>::getForceStress(UnitCell& ucell,
113117
ModuleBase::matrix fcc;
114118
ModuleBase::matrix fscc;
115119
ModuleBase::matrix fvnl_dalpha; // deepks
120+
ModuleBase::matrix fpothybrid;
116121

117122
fvl_dphi.create(nat, 3); // must do it now, update it later, noted by zhengdy
118123

@@ -127,6 +132,7 @@ void Force_Stress_LCAO<T>::getForceStress(UnitCell& ucell,
127132
fcc.create(nat, 3); // force due to core correction
128133
fscc.create(nat, 3); // force due to self-consistent field
129134
fvnl_dalpha.create(nat, 3); // deepks
135+
fpothybrid.create(nat, 3); // pulay force for hybrid gauge rt-tddft
130136

131137
// calculate basic terms in Force, same method with PW base
132138
this->calForcePwPart(ucell, fvl_dvl, fewalds, fcc, fscc, pelec->f_en.etxc,
@@ -199,14 +205,26 @@ void Force_Stress_LCAO<T>::getForceStress(UnitCell& ucell,
199205
hamilt::Overlap<hamilt::OperatorLCAO<T, double>> tmp_overlap(
200206
nullptr, kv.kvec_d, nullptr, nullptr, &ucell, orb.cutoffs(), &gd,
201207
two_center_bundle.overlap_orb.get());
202-
tmp_overlap.cal_force_stress(isforce, isstress, edmR, foverlap, soverlap);
208+
if(td_stype != 2)
209+
{
210+
tmp_overlap.cal_force_stress(isforce, isstress, edmR, foverlap, soverlap);
211+
}
203212

204213
// Calculate nonlocal force/stress (uses DM)
205214
hamilt::Nonlocal<hamilt::OperatorLCAO<T, double>> tmp_nonlocal(
206215
nullptr, kv.kvec_d, nullptr, &ucell, orb.cutoffs(), &gd,
207216
two_center_bundle.overlap_orb_beta.get());
208217
tmp_nonlocal.cal_force_stress(isforce, isstress, dmR, fvnl_dbeta, svnl_dbeta);
209218

219+
if(td_stype == 2)
220+
{
221+
hamilt::TD_pot_hybrid<hamilt::OperatorLCAO<T, double>> tmp_hybrid(
222+
nullptr, &kv, nullptr, nullptr, orb, &ucell, orb.cutoffs(), &gd, nullptr);
223+
tmp_hybrid.cal_force_stress(isforce, dmR, fpothybrid);
224+
225+
cal_foverlap_rt(foverlap, dmat, p_hamilt, kv, pv, ucell);
226+
}
227+
210228
// Switch back to spin channel 0
211229
if (PARAM.inp.nspin == 2)
212230
{
@@ -501,6 +519,7 @@ void Force_Stress_LCAO<T>::getForceStress(UnitCell& ucell,
501519
//---------------------------------
502520
// sum all parts of force!
503521
//---------------------------------
522+
ModuleBase::Vector3<double> net_force = {0.0, 0.0, 0.0};
504523
for (int i = 0; i < 3; i++)
505524
{
506525
double sum = 0.0;
@@ -511,7 +530,8 @@ void Force_Stress_LCAO<T>::getForceStress(UnitCell& ucell,
511530
+ fvl_dvl(iat, i) // derivative of local potential force (pw)
512531
+ fewalds(iat, i) // ewald force (pw)
513532
+ fcc(iat, i) // nonlinear core correction force (pw)
514-
+ fscc(iat, i); // self consistent corretion force (pw)
533+
+ fscc(iat, i) // self consistent corretion force (pw)
534+
+ fpothybrid(iat, i); // pulay force for hybrid gauge rt-tddft
515535

516536
// Force contribution from DFT+U, Quxin add on 20201029
517537
if (PARAM.inp.dft_plus_u)
@@ -564,7 +584,7 @@ void Force_Stress_LCAO<T>::getForceStress(UnitCell& ucell,
564584
// sum total force for correction
565585
sum += fcs(iat, i);
566586
}
567-
587+
net_force[i]=sum;
568588
if (!(PARAM.inp.gate_flag || PARAM.inp.efield_flag))
569589
{
570590
for (int iat = 0; iat < nat; ++iat)
@@ -678,6 +698,9 @@ void Force_Stress_LCAO<T>::getForceStress(UnitCell& ucell,
678698

679699
// this->printforce_total(ry, istestf, fcs);
680700
ModuleIO::print_force(GlobalV::ofs_running, ucell, "TOTAL-FORCE (eV/Angstrom)", fcs, false);
701+
net_force*= ModuleBase::Ry_to_eV / ModuleBase::BOHR_TO_A;
702+
ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, "Net force vector (eV/Ang)", net_force.x, net_force.y, net_force.z);
703+
ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, "Total drift (ev/Ang)", net_force.norm());
681704
if (istestf)
682705
{
683706
GlobalV::ofs_running << "\n FORCE INVALID TABLE." << std::endl;

‎source/source_lcao/FORCE_STRESS.h‎

Lines changed: 3 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -53,7 +53,9 @@ class Force_Stress_LCAO
5353
Plus_U &dftu, // mohan add 2025-11-07
5454
Setup_DeePKS<T> &deepks,
5555
Exx_NAO<T> &exx_nao,
56-
ModuleSymmetry::Symmetry* symm);
56+
ModuleSymmetry::Symmetry* symm,
57+
const int td_stype = 0,
58+
hamilt::Hamilt<T>* p_hamilt = nullptr);
5759

5860
private:
5961
int nat;

‎source/source_lcao/module_operator_lcao/td_pot_hybrid.cpp‎

Lines changed: 7 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -23,13 +23,17 @@ hamilt::TD_pot_hybrid<hamilt::OperatorLCAO<TK, TR>>::TD_pot_hybrid(
2323
{
2424
this->cal_type = calculation_type::lcao_tddft_periodic;
2525
this->ucell = ucell_in;
26+
this->gridD = GridD_in;
2627
#ifdef __DEBUG
2728
assert(this->ucell != nullptr);
2829
assert(this->hsk != nullptr);
2930
#endif
3031
this->init_td();
3132
// initialize HR to allocate sparse Ekinetic matrix memory
32-
this->initialize_HR(GridD_in);
33+
if(hR_in != nullptr)
34+
{
35+
this->initialize_HR(GridD_in);
36+
}
3337
}
3438

3539
// destructor
@@ -289,6 +293,7 @@ template <typename TK, typename TR>
289293
void hamilt::TD_pot_hybrid<hamilt::OperatorLCAO<TK, TR>>::contributeHk(int ik) {
290294
return;
291295
}
292-
296+
#include "td_pot_hybrid_force.hpp"
297+
template class hamilt::TD_pot_hybrid<hamilt::OperatorLCAO<double, double>>;
293298
template class hamilt::TD_pot_hybrid<hamilt::OperatorLCAO<std::complex<double>, double>>;
294299
template class hamilt::TD_pot_hybrid<hamilt::OperatorLCAO<std::complex<double>, std::complex<double>>>;

‎source/source_lcao/module_operator_lcao/td_pot_hybrid.h‎

Lines changed: 14 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -66,6 +66,9 @@ class TD_pot_hybrid<OperatorLCAO<TK, TR>> : public OperatorLCAO<TK, TR>
6666

6767
virtual void set_HR_fixed(void*) override;
6868

69+
void cal_force_stress(const bool cal_force,
70+
const HContainer<TR>* dmR,
71+
ModuleBase::matrix& force);
6972

7073
private:
7174
const UnitCell* ucell = nullptr;
@@ -120,6 +123,17 @@ class TD_pot_hybrid<OperatorLCAO<TK, TR>> : public OperatorLCAO<TK, TR>
120123

121124
/// @brief exact the nearest neighbor atoms from all adjacent atoms
122125
std::vector<AdjacentAtomInfo> adjs_all;
126+
127+
void cal_force_IJR(const int& iat1,
128+
const int& iat2,
129+
const Parallel_Orbitals* paraV,
130+
const ModuleBase::Vector3<double>& dtau,
131+
const ModuleBase::Vector3<double>& dR,
132+
TR* dmR_pointer,
133+
double* force1,
134+
double* force2);
135+
136+
const Grid_Driver* gridD = nullptr;
123137
};
124138

125139
} // namespace hamilt

0 commit comments

Comments
 (0)