Skip to content

Commit ed9ecbe

Browse files
wenfei-liwenfei-li
wenfei-li
and
wenfei-li
authored
Feature : out_mat_t, which prints the kinetic energy matrix (#1874)
* Feature : out_mat_t, which prints the kinetic energy matrix * modify creation of matrix directory * prepare for out_mat_dh --------- Co-authored-by: wenfei-li <[email protected]>
1 parent 74a6ff0 commit ed9ecbe

File tree

20 files changed

+400
-21
lines changed

20 files changed

+400
-21
lines changed

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

Lines changed: 13 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -23,7 +23,7 @@
2323
[relax_nmax](#relax_nmax) | [relax_method](#relax_method) | [relax_cg_thr](#relax_cg_thr) | [relax_bfgs_w1](#relax_bfgs_w1) | [relax_bfgs_w2](#relax_bfgs_w2) | [relax_bfgs_rmax](#relax_bfgs_rmax) | [relax_bfgs_rmin](#relax_bfgs_rmin) | [relax_bfgs_init](#relax_bfgs_init) | [cal_force](#cal_force) | [force_thr](#force_thr) | [force_thr_ev](#force_thr_ev) | [force_thr_ev2](#force_thr_ev2) | [cal_stress](#cal_stress) | [stress_thr](#stress_thr) | [press1, press2, press3](#press1-press2-press3) | [fixed_axes](#fixed_axes) | [cell_factor](#cell_factor) | [fixed_ibrav](#fixed_ibrav) | [relax_new](#relax_new) | [relax_scale_force](#relax_scale_force)
2424
- [Variables related to output information](#variables-related-to-output-information)
2525

26-
[out_force](#out_force) | [out_mul](#out_mul) | [out_freq_elec](#out_freq_elec) | [out_freq_ion](#out_freq_ion) | [out_chg](#out_chg) | [out_pot](#out_pot) | [out_dm](#out_dm) | [out_wfc_pw](#out_wfc_pw) | [out_wfc_r](#out_wfc_r) | [out_wfc_lcao](#out_wfc_lcao) | [out_dos](#out_dos) | [out_band](#out_band) | [out_proj_band](#out_proj_band) | [out_stru](#out_stru) | [out_bandgap](#out_bandgap) | [out_level](#out_level) | [out_alllog](#out_alllog) | [out_mat_hs](#out_mat_hs) | [out_mat_r](#out_mat_r) | [out_mat_hs2](#out_mat_hs2) | [out_hs2_interval](#out_hs2_interval) [out_element_info](#out_element_info) | [restart_save](#restart_save) | [restart_load](#restart_load) | [dft_plus_dmft](#dft_plus_dmft) | [rpa](#rpa)
26+
[out_force](#out_force) | [out_mul](#out_mul) | [out_freq_elec](#out_freq_elec) | [out_freq_ion](#out_freq_ion) | [out_chg](#out_chg) | [out_pot](#out_pot) | [out_dm](#out_dm) | [out_wfc_pw](#out_wfc_pw) | [out_wfc_r](#out_wfc_r) | [out_wfc_lcao](#out_wfc_lcao) | [out_dos](#out_dos) | [out_band](#out_band) | [out_proj_band](#out_proj_band) | [out_stru](#out_stru) | [out_bandgap](#out_bandgap) | [out_level](#out_level) | [out_alllog](#out_alllog) | [out_mat_hs](#out_mat_hs) | [out_mat_r](#out_mat_r) | [out_mat_hs2](#out_mat_hs2) | [out_mat_t](#out_mat_t) | [out_mat_dh](#out_mat_dh) | [out_hs2_interval](#out_hs2_interval) | [out_element_info](#out_element_info) | [restart_save](#restart_save) | [restart_load](#restart_load) | [dft_plus_dmft](#dft_plus_dmft) | [rpa](#rpa)
2727
- [Density of states](#density-of-states)
2828

2929
[dos_edelta_ev](#dos_edelta_ev) | [dos_sigma](#dos_sigma) | [dos_scale](#dos_scale) | [dos_emin_ev](#dos_emin_ev) | [dos_emax_ev](#dos_emax_ev) | [dos_nche](#dos_nche)
@@ -996,13 +996,23 @@ These variables are used to control the output of properties.
996996
- The array ROW_INDEX is of length m + 1 and encodes the index in V and COL_INDEX where the given row starts. This is equivalent to ROW_INDEX[j] encoding the total number of nonzeros above row j. The last element is NNZ , i.e., the fictitious index in V immediately after the last valid index NNZ - 1.
997997

998998
> Note: This functionality is not available for gamma_only calculations. If you want to use it in gamma_only calculations, you should turn off gamma_only, and explicitly specifies that gamma point is the only k point in the KPT file.
999-
>
999+
1000+
- **Default**: 0
1001+
1002+
### out_mat_t
1003+
- **Type**: Boolean
1004+
- **Description**: For LCAO calculations, if out_mat_t is set to 1, ABACUS will generate files containing the kinetic energy matrix T(R). The format will be the same as the Hamiltonian matrix H(R) and overlap matrix S(R) as mentioned in [out_mat_hs2](#out_mat_hs2). The name of the files will be `data-TR-sparse_SPIN0.csr` and so on. Also controled by [out_hs2_interval](#out_hs2_interval).
1005+
- **Default**: 0
1006+
1007+
### out_mat_dh
1008+
- **Type**: Boolean
1009+
- **Description**: For LCAO calculations, if out_mat_dh is set to 1, ABACUS will generate files containing the derivatives of the Hamiltonian matrix. The format will be the same as the Hamiltonian matrix H(R) and overlap matrix S(R) as mentioned in [out_mat_hs2](#out_mat_hs2). The name of the files will be `data-dHRx-sparse_SPIN0.csr` and so on. Also controled by [out_hs2_interval](#out_hs2_interval).
10001010
- **Default**: 0
10011011

10021012
### out_hs2_interval
10031013

10041014
- **Type**: Integer
1005-
- **Description**: Only relevant for printing H(R) and S(R) matrices during MD. It controls the interval at which to print. Check input parameter [out_mat_hs2](#out_mat_hs2) for more information.
1015+
- **Description**: Only relevant for printing H(R), S(R), T(R), dH(R) matrices during MD. It controls the interval at which to print. Check input parameter [out_mat_hs2](#out_mat_hs2) for more information.
10061016
- **Default**: 1
10071017

10081018
### out_element_info

‎source/module_base/global_file.cpp‎

Lines changed: 2 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -19,8 +19,7 @@ namespace ModuleBase
1919
void ModuleBase::Global_File::make_dir_out(
2020
const std::string &suffix,
2121
const std::string &calculation,
22-
const bool &out_hs,
23-
const bool &out_r,
22+
const bool &out_dir,
2423
const int rank,
2524
const bool &restart,
2625
const bool out_alllog)
@@ -124,7 +123,7 @@ void ModuleBase::Global_File::make_dir_out(
124123
}
125124

126125
// make dir for HS matrix output in md calculation
127-
if((out_hs || out_r) && calculation == "md")
126+
if((out_dir) && calculation == "md")
128127
{
129128
int make_dir_matrix = 0;
130129
std::string command1 = "test -d " + GlobalV::global_matrix_dir + " || mkdir " + GlobalV::global_matrix_dir;

‎source/module_base/global_file.h‎

Lines changed: 1 addition & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -23,8 +23,7 @@ namespace Global_File
2323
// called in input.cpp, after reading parameters.
2424
void make_dir_out(const std::string &suffix,
2525
const std::string &calculation,
26-
const bool &out_hs,
27-
const bool &out_r,
26+
const bool &out_dir,
2827
const int rank,
2928
const bool &restart,
3029
const bool out_alllog = false);

‎source/module_esolver/esolver_ks_lcao.cpp‎

Lines changed: 16 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -906,6 +906,22 @@ void ESolver_KS_LCAO::afterscf(const int istep)
906906
} // LiuXh add 2019-07-15
907907
}
908908

909+
if (hsolver::HSolverLCAO::out_mat_t)
910+
{
911+
if( !(GlobalV::CALCULATION=="md" && (istep%hsolver::HSolverLCAO::out_hsR_interval!=0)) )
912+
{
913+
ModuleIO::output_T_R(istep, this->UHM); // LiuXh add 2019-07-15
914+
} // LiuXh add 2019-07-15
915+
}
916+
917+
if (hsolver::HSolverLCAO::out_mat_dh)
918+
{
919+
if( !(GlobalV::CALCULATION=="md" && (istep%hsolver::HSolverLCAO::out_hsR_interval!=0)) )
920+
{
921+
//ModuleIO::output_DH_R(istep, this->pelec->pot->get_effective_v(), this->UHM); // LiuXh add 2019-07-15
922+
} // LiuXh add 2019-07-15
923+
}
924+
909925
// add by jingan for out r_R matrix 2019.8.14
910926
if(INPUT.out_mat_r)
911927
{

‎source/module_esolver/esolver_ks_lcao_elec.cpp‎

Lines changed: 1 addition & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -450,8 +450,7 @@ namespace ModuleESolver
450450
this->LM.allocate_HS_R(this->orb_con.ParaV.nnr);
451451
this->LM.zeros_HSR('S');
452452
this->UHM.genH.calculate_S_no(this->LM.SlocR.data());
453-
ModuleIO::output_SR("SR.csr",this->UHM);
454-
453+
ModuleIO::output_S_R(this->UHM,"SR.csr");
455454
}
456455
}
457456

‎source/module_esolver/esolver_ks_lcao_tddft.cpp‎

Lines changed: 8 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -414,6 +414,14 @@ void ESolver_KS_LCAO_TDDFT::afterscf(const int istep)
414414
}
415415
}
416416

417+
if (hsolver::HSolverLCAO::out_mat_t)
418+
{
419+
if (!(GlobalV::CALCULATION == "md" && (istep % hsolver::HSolverLCAO::out_hsR_interval != 0)))
420+
{
421+
ModuleIO::output_T_R(istep, this->UHM); // LiuXh add 2019-07-15
422+
}
423+
}
424+
417425
// add by jingan for out r_R matrix 2019.8.14
418426
if (INPUT.out_mat_r)
419427
{

‎source/module_hamilt_lcao/hamilt_lcaodft/LCAO_hamilt.cpp‎

Lines changed: 123 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -342,13 +342,131 @@ void LCAO_Hamilt::calculate_HSR_sparse(const int &current_spin, const double &sp
342342
clear_zero_elements(current_spin, sparse_threshold);
343343
}
344344

345+
void LCAO_Hamilt::calculate_STN_R_sparse_for_T(const double &sparse_threshold)
346+
{
347+
ModuleBase::TITLE("LCAO_Hamilt","calculate_STN_R_sparse_for_T");
348+
349+
int index = 0;
350+
ModuleBase::Vector3<double> dtau, tau1, tau2;
351+
ModuleBase::Vector3<double> dtau1, dtau2, tau0;
352+
353+
double temp_value_double;
354+
std::complex<double> temp_value_complex;
355+
356+
for(int T1 = 0; T1 < GlobalC::ucell.ntype; ++T1)
357+
{
358+
Atom* atom1 = &GlobalC::ucell.atoms[T1];
359+
for(int I1 = 0; I1 < atom1->na; ++I1)
360+
{
361+
tau1 = atom1->tau[I1];
362+
GlobalC::GridD.Find_atom(GlobalC::ucell, tau1, T1, I1);
363+
Atom* atom1 = &GlobalC::ucell.atoms[T1];
364+
const int start = GlobalC::ucell.itiaiw2iwt(T1,I1,0);
365+
366+
for(int ad = 0; ad < GlobalC::GridD.getAdjacentNum()+1; ++ad)
367+
{
368+
const int T2 = GlobalC::GridD.getType(ad);
369+
const int I2 = GlobalC::GridD.getNatom(ad);
370+
Atom* atom2 = &GlobalC::ucell.atoms[T2];
371+
372+
tau2 = GlobalC::GridD.getAdjacentTau(ad);
373+
dtau = tau2 - tau1;
374+
double distance = dtau.norm() * GlobalC::ucell.lat0;
375+
double rcut = GlobalC::ORB.Phi[T1].getRcut() + GlobalC::ORB.Phi[T2].getRcut();
376+
377+
bool adj = false;
378+
379+
if(distance < rcut) adj = true;
380+
else if(distance >= rcut)
381+
{
382+
for(int ad0 = 0; ad0 < GlobalC::GridD.getAdjacentNum()+1; ++ad0)
383+
{
384+
const int T0 = GlobalC::GridD.getType(ad0);
385+
386+
tau0 = GlobalC::GridD.getAdjacentTau(ad0);
387+
dtau1 = tau0 - tau1;
388+
dtau2 = tau0 - tau2;
389+
390+
double distance1 = dtau1.norm() * GlobalC::ucell.lat0;
391+
double distance2 = dtau2.norm() * GlobalC::ucell.lat0;
392+
393+
double rcut1 = GlobalC::ORB.Phi[T1].getRcut() + GlobalC::ucell.infoNL.Beta[T0].get_rcut_max();
394+
double rcut2 = GlobalC::ORB.Phi[T2].getRcut() + GlobalC::ucell.infoNL.Beta[T0].get_rcut_max();
395+
396+
if( distance1 < rcut1 && distance2 < rcut2 )
397+
{
398+
adj = true;
399+
break;
400+
}
401+
}
402+
}
403+
404+
if(adj)
405+
{
406+
const int start2 = GlobalC::ucell.itiaiw2iwt(T2,I2,0);
407+
408+
Abfs::Vector3_Order<int> dR(GlobalC::GridD.getBox(ad).x, GlobalC::GridD.getBox(ad).y, GlobalC::GridD.getBox(ad).z);
409+
410+
for(int ii=0; ii<atom1->nw*GlobalV::NPOL; ii++)
411+
{
412+
const int iw1_all = start + ii;
413+
const int mu = this->LM->ParaV->trace_loc_row[iw1_all];
414+
415+
if(mu<0)continue;
416+
417+
for(int jj=0; jj<atom2->nw*GlobalV::NPOL; jj++)
418+
{
419+
int iw2_all = start2 + jj;
420+
const int nu = this->LM->ParaV->trace_loc_col[iw2_all];
421+
422+
if(nu<0)continue;
423+
424+
if(GlobalV::NSPIN!=4)
425+
{
426+
temp_value_double = this->LM->Hloc_fixedR[index];
427+
if (std::abs(temp_value_double) > sparse_threshold)
428+
{
429+
this->LM->TR_sparse[dR][iw1_all][iw2_all] = temp_value_double;
430+
}
431+
}
432+
else
433+
{
434+
temp_value_complex = this->LM->Hloc_fixedR_soc[index];
435+
if(std::abs(temp_value_complex) > sparse_threshold)
436+
{
437+
this->LM->TR_soc_sparse[dR][iw1_all][iw2_all] = temp_value_complex;
438+
}
439+
}
440+
441+
++index;
442+
}
443+
}
444+
}
445+
}
446+
}
447+
}
448+
449+
return;
450+
}
451+
345452
void LCAO_Hamilt::calculate_SR_sparse(const double &sparse_threshold)
346453
{
347454
ModuleBase::TITLE("LCAO_Hamilt","calculate_SR_sparse");
348455
set_R_range_sparse();
349456
calculate_STN_R_sparse_for_S(sparse_threshold);
350457
}
351458

459+
void LCAO_Hamilt::calculate_TR_sparse(const double &sparse_threshold)
460+
{
461+
ModuleBase::TITLE("LCAO_Hamilt","calculate_TR_sparse");
462+
463+
//need to rebuild T(R)
464+
this->LM->zeros_HSR('T');
465+
this->genH.build_ST_new('T', 0, GlobalC::ucell, this->LM->Hloc_fixedR.data());
466+
set_R_range_sparse();
467+
calculate_STN_R_sparse_for_T(sparse_threshold);
468+
}
469+
352470
void LCAO_Hamilt::calculat_HR_dftu_sparse(const int &current_spin, const double &sparse_threshold)
353471
{
354472
ModuleBase::TITLE("LCAO_Hamilt","calculat_HR_dftu_sparse");
@@ -747,3 +865,8 @@ void LCAO_Hamilt::destroy_all_HSR_sparse(void)
747865
{
748866
this->LM->destroy_HS_R_sparse();
749867
}
868+
869+
void LCAO_Hamilt::destroy_TR_sparse(void)
870+
{
871+
this->LM->destroy_T_R_sparse();
872+
}

‎source/module_hamilt_lcao/hamilt_lcaodft/LCAO_hamilt.h‎

Lines changed: 3 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -24,6 +24,7 @@ class LCAO_Hamilt
2424
void set_R_range_sparse();
2525
void calculate_STN_R_sparse(const int &current_spin, const double &sparse_threshold);
2626
void calculate_STN_R_sparse_for_S(const double &sparse_threshold);
27+
void calculate_STN_R_sparse_for_T(const double &sparse_threshold);
2728
void calculat_HR_dftu_sparse(const int &current_spin, const double &sparse_threshold);
2829
void calculat_HR_dftu_soc_sparse(const int &current_spin, const double &sparse_threshold);
2930
void calculate_HR_exx_sparse(const int &current_spin, const double &sparse_threshold);
@@ -35,6 +36,8 @@ class LCAO_Hamilt
3536
void calculate_SR_sparse(const double &sparse_threshold);
3637
void clear_zero_elements(const int &current_spin, const double &sparse_threshold);
3738
void destroy_all_HSR_sparse(void);
39+
void calculate_TR_sparse(const double &sparse_threshold);
40+
void destroy_TR_sparse(void);
3841

3942
// used for gamma only algorithms.
4043
Gint_Gamma GG;

‎source/module_hamilt_lcao/hamilt_lcaodft/LCAO_matrix.cpp‎

Lines changed: 17 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -952,3 +952,20 @@ void LCAO_Matrix::destroy_HS_R_sparse(void)
952952

953953
return;
954954
}
955+
956+
void LCAO_Matrix::destroy_T_R_sparse(void)
957+
{
958+
ModuleBase::TITLE("LCAO_Matrix","destroy_T_R_sparse");
959+
960+
if (GlobalV::NSPIN != 4)
961+
{
962+
std::map<Abfs::Vector3_Order<int>, std::map<size_t, std::map<size_t, double>>> empty_TR_sparse;
963+
TR_sparse.swap(empty_TR_sparse);
964+
}
965+
else
966+
{
967+
std::map<Abfs::Vector3_Order<int>, std::map<size_t, std::map<size_t, std::complex<double>>>> empty_TR_soc_sparse;
968+
TR_soc_sparse.swap(empty_TR_soc_sparse);
969+
}
970+
return;
971+
}

‎source/module_hamilt_lcao/hamilt_lcaodft/LCAO_matrix.h‎

Lines changed: 3 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -97,10 +97,12 @@ class LCAO_Matrix
9797
// For HR_sparse[2], when nspin=1, only 0 is valid, when nspin=2, 0 means spin up, 1 means spin down
9898
std::map<Abfs::Vector3_Order<int>, std::map<size_t, std::map<size_t, double>>> HR_sparse[2];
9999
std::map<Abfs::Vector3_Order<int>, std::map<size_t, std::map<size_t, double>>> SR_sparse;
100+
std::map<Abfs::Vector3_Order<int>, std::map<size_t, std::map<size_t, double>>> TR_sparse;
100101

101102
// For nspin = 4
102103
std::map<Abfs::Vector3_Order<int>, std::map<size_t, std::map<size_t, std::complex<double>>>> HR_soc_sparse;
103104
std::map<Abfs::Vector3_Order<int>, std::map<size_t, std::map<size_t, std::complex<double>>>> SR_soc_sparse;
105+
std::map<Abfs::Vector3_Order<int>, std::map<size_t, std::map<size_t, std::complex<double>>>> TR_soc_sparse;
104106

105107
// Record all R direct coordinate information, even if HR or SR is a zero matrix
106108
std::set<Abfs::Vector3_Order<int>> all_R_coor;
@@ -206,6 +208,7 @@ class LCAO_Matrix
206208

207209
// jingan add 2021-6-4, modify 2021-12-2
208210
void destroy_HS_R_sparse(void);
211+
void destroy_T_R_sparse(void);
209212

210213
};
211214

0 commit comments

Comments
 (0)