Skip to content

Commit 85f031a

Browse files
authored
Merge branch 'develop' into neighbor2
2 parents 2709779 + 028afb0 commit 85f031a

Some content is hidden

Large Commits have some content hidden by default. Use the searchbox below for content that may be hidden.

75 files changed

+1990
-1294
lines changed

source/Makefile.Objects

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -429,6 +429,8 @@ OBJS_RELAXATION=bfgs_basic.o\
429429
relax_old.o\
430430
relax.o\
431431
line_search.o\
432+
bfgs.o\
433+
432434

433435
OBJS_SURCHEM=surchem.o\
434436
H_correction_pw.o\

source/module_elecstate/elecstate_pw.h

Lines changed: 32 additions & 18 deletions
Original file line numberDiff line numberDiff line change
@@ -17,6 +17,7 @@ class ElecStatePW : public ElecState
1717
{
1818
private:
1919
using Real = typename GetTypeReal<T>::type;
20+
2021
public:
2122
ElecStatePW(ModulePW::PW_Basis_K* wfc_basis_in,
2223
Charge* chg_in,
@@ -26,46 +27,59 @@ class ElecStatePW : public ElecState
2627
ModulePW::PW_Basis* rhodpw_in,
2728
ModulePW::PW_Basis* rhopw_in,
2829
ModulePW::PW_Basis_Big* bigpw_in);
29-
// void init(Charge* chg_in):charge(chg_in){} override;
3030

3131
~ElecStatePW();
32-
// interface for HSolver to calculate rho from Psi
32+
33+
//! interface for HSolver to calculate rho from Psi
3334
virtual void psiToRho(const psi::Psi<T, Device>& psi);
34-
// return current electronic density rho, as a input for constructing Hamiltonian
35-
// const double* getRho(int spin) const override;
36-
virtual void cal_tau(const psi::Psi<T, Device>& psi);
3735

38-
// update charge density for next scf step
39-
// void getNewRho() override;
36+
virtual void cal_tau(const psi::Psi<T, Device>& psi);
4037

4138
Real* becsum = nullptr;
42-
// init rho_data and kin_r_data
39+
40+
//! init rho_data and kin_r_data
4341
void init_rho_data();
44-
Real ** rho = nullptr, ** kin_r = nullptr; //[Device] [spin][nrxx] rho and kin_r
42+
Real** rho = nullptr;
43+
Real** kin_r = nullptr; //[Device] [spin][nrxx] rho and kin_r
4544

4645
protected:
46+
4747
ModulePW::PW_Basis* rhopw_smooth = nullptr;
48+
4849
ModulePW::PW_Basis_K* basis = nullptr;
50+
4951
UnitCell* ucell = nullptr;
52+
5053
const pseudopot_cell_vnl* ppcell = nullptr;
5154

52-
// calculate electronic charge density on grid points or density matrix in real space
53-
// the consequence charge density rho saved into rho_out, preparing for charge mixing.
55+
//! calculate electronic charge density on grid points or density matrix in real space
56+
//! the consequence charge density rho saved into rho_out, preparing for charge mixing.
5457
void updateRhoK(const psi::Psi<T, Device>& psi); // override;
55-
// sum over all pools for rho and ebands
58+
59+
//! sum over all pools for rho and ebands
5660
void parallelK();
57-
// calcualte rho for each k
61+
62+
//! calcualte rho for each k
5863
void rhoBandK(const psi::Psi<T, Device>& psi);
59-
// add to the charge density in reciprocal space the part which is due to the US augmentation.
64+
65+
//! add to the charge density in reciprocal space the part which is due to the US augmentation.
6066
void add_usrho(const psi::Psi<T, Device>& psi);
61-
// \sum_lm Q_lm(r) \sum_i <psi_i|beta_l><beta_m|psi_i> w_i
67+
68+
//! Non-local pseudopotentials
69+
//! \sum_lm Q_lm(r) \sum_i <psi_i|beta_l><beta_m|psi_i> w_i
6270
void addusdens_g(const Real* becsum, T* rhog);
6371

6472
Device * ctx = {};
73+
6574
bool init_rho = false;
75+
6676
mutable T* vkb = nullptr;
67-
Real * rho_data = nullptr, * kin_r_data = nullptr;
68-
T * wfcr = nullptr, * wfcr_another_spin = nullptr;
77+
78+
Real* rho_data = nullptr;
79+
Real* kin_r_data = nullptr;
80+
T* wfcr = nullptr;
81+
T* wfcr_another_spin = nullptr;
82+
6983
private:
7084
using meta_op = hamilt::meta_pw_op<Real, Device>;
7185
using elecstate_pw_op = elecstate::elecstate_pw_op<Real, Device>;
@@ -85,4 +99,4 @@ class ElecStatePW : public ElecState
8599

86100
} // namespace elecstate
87101

88-
#endif
102+
#endif

source/module_elecstate/test/charge_extra_test.cpp

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -87,7 +87,7 @@ Structure_Factor::Structure_Factor()
8787
Structure_Factor::~Structure_Factor()
8888
{
8989
}
90-
void Structure_Factor::setup_structure_factor(UnitCell* Ucell, const ModulePW::PW_Basis* rho_basis)
90+
void Structure_Factor::setup_structure_factor(const UnitCell* Ucell, const ModulePW::PW_Basis* rho_basis)
9191
{
9292
}
9393

source/module_esolver/esolver_gets.cpp

Lines changed: 5 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -101,14 +101,15 @@ void ESolver_GetS::runner(UnitCell& ucell, const int istep)
101101
PARAM.inp.test_atom_input);
102102

103103
Record_adj RA;
104-
RA.for_2d(this->pv, PARAM.globalv.gamma_only_local, orb_.cutoffs());
104+
RA.for_2d(ucell,this->pv, PARAM.globalv.gamma_only_local, orb_.cutoffs());
105105

106106
if (this->p_hamilt == nullptr)
107107
{
108108
if (PARAM.inp.nspin == 4)
109109
{
110110
this->p_hamilt
111-
= new hamilt::HamiltLCAO<std::complex<double>, std::complex<double>>(&this->pv,
111+
= new hamilt::HamiltLCAO<std::complex<double>, std::complex<double>>(ucell,
112+
&this->pv,
112113
this->kv,
113114
*(two_center_bundle_.overlap_orb),
114115
orb_.cutoffs());
@@ -117,7 +118,8 @@ void ESolver_GetS::runner(UnitCell& ucell, const int istep)
117118
}
118119
else
119120
{
120-
this->p_hamilt = new hamilt::HamiltLCAO<std::complex<double>, double>(&this->pv,
121+
this->p_hamilt = new hamilt::HamiltLCAO<std::complex<double>, double>(ucell,
122+
&this->pv,
121123
this->kv,
122124
*(two_center_bundle_.overlap_orb),
123125
orb_.cutoffs());

source/module_esolver/esolver_ks_lcao.cpp

Lines changed: 6 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -164,7 +164,7 @@ void ESolver_KS_LCAO<TK, TR>::before_all_runners(UnitCell& ucell, const Input_pa
164164
// 5) initialize Hamilt in LCAO
165165
// * allocate H and S matrices according to computational resources
166166
// * set the 'trace' between local H/S and global H/S
167-
LCAO_domain::divide_HS_in_frag(PARAM.globalv.gamma_only_local, pv, this->kv.get_nks(), orb_);
167+
LCAO_domain::divide_HS_in_frag(PARAM.globalv.gamma_only_local, ucell , pv, this->kv.get_nks(), orb_);
168168

169169
#ifdef __EXX
170170
// 6) initialize exx
@@ -294,6 +294,7 @@ void ESolver_KS_LCAO<TK, TR>::cal_force(UnitCell& ucell, ModuleBase::matrix& for
294294
PARAM.inp.cal_stress,
295295
PARAM.inp.test_force,
296296
PARAM.inp.test_stress,
297+
ucell,
297298
this->pv,
298299
this->pelec,
299300
this->psi,
@@ -643,7 +644,7 @@ void ESolver_KS_LCAO<TK, TR>::iter_init(UnitCell& ucell, const int istep, const
643644
GlobalC::dftu.set_dmr(dynamic_cast<elecstate::ElecStateLCAO<TK>*>(this->pelec)->get_DM());
644645
}
645646
// Calculate U and J if Yukawa potential is used
646-
GlobalC::dftu.cal_slater_UJ(this->pelec->charge->rho, this->pw_rho->nrxx);
647+
GlobalC::dftu.cal_slater_UJ(ucell,this->pelec->charge->rho, this->pw_rho->nrxx);
647648
}
648649

649650
#ifdef __DEEPKS
@@ -874,11 +875,11 @@ void ESolver_KS_LCAO<TK, TR>::iter_finish(UnitCell& ucell, const int istep, int&
874875
{
875876
const std::vector<std::vector<TK>>& tmp_dm
876877
= dynamic_cast<elecstate::ElecStateLCAO<TK>*>(this->pelec)->get_DM()->get_DMK_vector();
877-
ModuleDFTU::dftu_cal_occup_m(iter, tmp_dm, this->kv, this->p_chgmix->get_mixing_beta(), this->p_hamilt);
878+
ModuleDFTU::dftu_cal_occup_m(iter, ucell,tmp_dm, this->kv, this->p_chgmix->get_mixing_beta(), this->p_hamilt);
878879
}
879-
GlobalC::dftu.cal_energy_correction(istep);
880+
GlobalC::dftu.cal_energy_correction(ucell,istep);
880881
}
881-
GlobalC::dftu.output();
882+
GlobalC::dftu.output(ucell);
882883
}
883884

884885
// (7) for deepks, calculate delta_e

source/module_esolver/lcao_before_scf.cpp

Lines changed: 3 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -127,7 +127,7 @@ void ESolver_KS_LCAO<TK, TR>::before_scf(UnitCell& ucell, const int istep)
127127
// (2)For each atom, calculate the adjacent atoms in different cells
128128
// and allocate the space for H(R) and S(R).
129129
// If k point is used here, allocate HlocR after atom_arrange.
130-
this->RA.for_2d(this->pv, PARAM.globalv.gamma_only_local, orb_.cutoffs());
130+
this->RA.for_2d(ucell,this->pv, PARAM.globalv.gamma_only_local, orb_.cutoffs());
131131

132132
// 2. density matrix extrapolation
133133

@@ -174,7 +174,7 @@ void ESolver_KS_LCAO<TK, TR>::before_scf(UnitCell& ucell, const int istep)
174174
}
175175

176176
// prepare grid in Gint
177-
LCAO_domain::grid_prepare(this->GridT, this->GG, this->GK, orb_, *this->pw_rho, *this->pw_big);
177+
LCAO_domain::grid_prepare(this->GridT, this->GG, this->GK, ucell, orb_, *this->pw_rho, *this->pw_big);
178178

179179
// init Hamiltonian
180180
if (this->p_hamilt != nullptr)
@@ -188,6 +188,7 @@ void ESolver_KS_LCAO<TK, TR>::before_scf(UnitCell& ucell, const int istep)
188188
this->p_hamilt = new hamilt::HamiltLCAO<TK, TR>(
189189
PARAM.globalv.gamma_only_local ? &(this->GG) : nullptr,
190190
PARAM.globalv.gamma_only_local ? nullptr : &(this->GK),
191+
ucell,
191192
&this->pv,
192193
this->pelec->pot,
193194
this->kv,

source/module_esolver/lcao_others.cpp

Lines changed: 3 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -128,7 +128,7 @@ void ESolver_KS_LCAO<TK, TR>::others(UnitCell& ucell, const int istep)
128128
// (2)For each atom, calculate the adjacent atoms in different cells
129129
// and allocate the space for H(R) and S(R).
130130
// If k point is used here, allocate HlocR after atom_arrange.
131-
this->RA.for_2d(this->pv, PARAM.globalv.gamma_only_local, orb_.cutoffs());
131+
this->RA.for_2d(ucell,this->pv, PARAM.globalv.gamma_only_local, orb_.cutoffs());
132132

133133
// 2. density matrix extrapolation
134134

@@ -175,7 +175,7 @@ void ESolver_KS_LCAO<TK, TR>::others(UnitCell& ucell, const int istep)
175175
}
176176

177177
// prepare grid in Gint
178-
LCAO_domain::grid_prepare(this->GridT, this->GG, this->GK, orb_, *this->pw_rho, *this->pw_big);
178+
LCAO_domain::grid_prepare(this->GridT, this->GG, this->GK, ucell,orb_, *this->pw_rho, *this->pw_big);
179179

180180
// init Hamiltonian
181181
if (this->p_hamilt != nullptr)
@@ -189,6 +189,7 @@ void ESolver_KS_LCAO<TK, TR>::others(UnitCell& ucell, const int istep)
189189
this->p_hamilt = new hamilt::HamiltLCAO<TK, TR>(
190190
PARAM.globalv.gamma_only_local ? &(this->GG) : nullptr,
191191
PARAM.globalv.gamma_only_local ? nullptr : &(this->GK),
192+
ucell,
192193
&this->pv,
193194
this->pelec->pot,
194195
this->kv,

source/module_hamilt_general/module_surchem/cal_pseudo.cpp

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -8,7 +8,7 @@ void surchem::gauss_charge(const UnitCell& cell,
88
complex<double>* N,
99
Structure_Factor* sf)
1010
{
11-
sf->setup_structure_factor(&GlobalC::ucell, rho_basis); // call strucFac(ntype,ngmc)
11+
sf->setup_structure_factor(&cell, rho_basis); // call strucFac(ntype,ngmc)
1212
for (int it = 0; it < cell.ntype; it++)
1313
{
1414
double RCS = GetAtom.atom_RCS[cell.atoms[it].ncpp.psd];

source/module_hamilt_general/module_surchem/cal_vcav.cpp

Lines changed: 10 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -2,7 +2,10 @@
22
#include "module_hamilt_general/module_xc/xc_functional.h"
33
#include "surchem.h"
44

5-
void lapl_rho(const std::complex<double>* rhog, double* lapn, const ModulePW::PW_Basis* rho_basis)
5+
void lapl_rho(const double& tpiba2,
6+
const std::complex<double>* rhog,
7+
double* lapn,
8+
const ModulePW::PW_Basis* rho_basis)
69
{
710
std::complex<double> *gdrtmpg = new std::complex<double>[rho_basis->npw];
811
ModuleBase::GlobalFunc::ZEROS(lapn, rho_basis->nrxx);
@@ -22,9 +25,9 @@ void lapl_rho(const std::complex<double>* rhog, double* lapn, const ModulePW::PW
2225
// bring the gdr from G --> R
2326
rho_basis->recip2real(aux, aux);
2427
// remember to multily 2pi/a0, which belongs to G vectors.
25-
for (int ir = 0; ir < rho_basis->nrxx; ir++) {
26-
lapn[ir] -= aux[ir].real() * GlobalC::ucell.tpiba2;
27-
}
28+
for (int ir = 0; ir < rho_basis->nrxx; ir++)
29+
lapn[ir] -= aux[ir].real() * tpiba2;
30+
2831
}
2932

3033
delete[] gdrtmpg;
@@ -70,7 +73,7 @@ void surchem::createcavity(const UnitCell& ucell,
7073
ModuleBase::GlobalFunc::ZEROS(lapn, rho_basis->nrxx);
7174

7275
// nabla n
73-
XC_Functional::grad_rho(PS_TOTN, nablan, rho_basis, GlobalC::ucell.tpiba);
76+
XC_Functional::grad_rho(PS_TOTN, nablan, rho_basis, ucell.tpiba);
7477

7578
// |\nabla n |^2 = nablan_2
7679
for (int ir = 0; ir < rho_basis->nrxx; ir++)
@@ -79,7 +82,7 @@ void surchem::createcavity(const UnitCell& ucell,
7982
}
8083

8184
// Laplacian of n
82-
lapl_rho(PS_TOTN, lapn, rho_basis);
85+
lapl_rho(ucell.tpiba2,PS_TOTN, lapn, rho_basis);
8386

8487
//-------------------------------------------------------------
8588
// add -Lap(n)/|\nabla n| to vwork and copy \sqrt(|\nabla n|^2)
@@ -131,7 +134,7 @@ void surchem::createcavity(const UnitCell& ucell,
131134

132135
// \nabla(1 / |\nabla n|), ggn in real space
133136
ModuleBase::Vector3<double> *ggn = new ModuleBase::Vector3<double>[rho_basis->nrxx];
134-
XC_Functional::grad_rho(inv_gn, ggn, rho_basis, GlobalC::ucell.tpiba);
137+
XC_Functional::grad_rho(inv_gn, ggn, rho_basis, ucell.tpiba);
135138

136139
//-------------------------------------------------------------
137140
// add -(\nabla n . \nabla(1/ |\nabla n|)) to Vcav in real space

source/module_hamilt_general/module_surchem/cal_vel.cpp

Lines changed: 3 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -17,6 +17,7 @@ void shape_gradn(const double* PS_TOTN_real, const ModulePW::PW_Basis* rho_basis
1717
}
1818

1919
void eps_pot(const double* PS_TOTN_real,
20+
const double& tpiba,
2021
const complex<double>* phi,
2122
const ModulePW::PW_Basis* rho_basis,
2223
double* d_eps,
@@ -36,7 +37,7 @@ void eps_pot(const double* PS_TOTN_real,
3637
double *phisq = new double[rho_basis->nrxx];
3738

3839
// nabla phi
39-
XC_Functional::grad_rho(phi, nabla_phi, rho_basis, GlobalC::ucell.tpiba);
40+
XC_Functional::grad_rho(phi, nabla_phi, rho_basis, tpiba);
4041

4142
for (int ir = 0; ir < rho_basis->nrxx; ir++)
4243
{
@@ -118,7 +119,7 @@ ModuleBase::matrix surchem::cal_vel(const UnitCell& cell,
118119
this->Ael *= cell.omega / rho_basis->nxyz;
119120

120121
// the 2nd item of tmp_Vel
121-
eps_pot(PS_TOTN_real, Sol_phi, rho_basis, epsilon, epspot);
122+
eps_pot(PS_TOTN_real, cell.tpiba, Sol_phi, rho_basis, epsilon, epspot);
122123

123124
for (int i = 0; i < rho_basis->nrxx; i++)
124125
{

0 commit comments

Comments
 (0)