Skip to content

Commit 823c5a2

Browse files
wenfei-liwenfei-li
authored
Merge pull request #897 from wenfei-li/develop
gint : removed obsolete code; implement new gamma force & stress
2 parents 51200fe + 53a6361 commit 823c5a2

25 files changed

+708
-2222
lines changed

‎source/Makefile.Objects‎

Lines changed: 0 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -111,7 +111,6 @@ grid_bigcell.o\
111111
grid_meshball.o\
112112
gint_tools.o\
113113
gint_gamma.o\
114-
gint_gamma_common.o\
115114
gint_gamma_fvl.o\
116115
gint_gamma_vl.o\
117116
gint_gamma_rho.o\
@@ -122,7 +121,6 @@ gint_k_init.o\
122121
gint_k_vl.o\
123122
gint_k_rho.o\
124123
gint_k_fvl.o\
125-
grid_base_beta.o\
126124
ORB_control.o\
127125
ORB_read.o\
128126
ORB_atomic.o\

‎source/module_grid/Makefile.Objects‎

Lines changed: 0 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -18,7 +18,6 @@ complexmatrix.o\
1818
matrix.o\
1919

2020
OBJS_GRID=
21-
grid_base_beta.o\
2221
grid_bigcell.o\
2322
grid_meshball.o\
2423
grid_meshk.o\

‎source/src_lcao/CMakeLists.txt‎

Lines changed: 0 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -12,7 +12,6 @@ list(APPEND objects
1212
FORCE_gamma_edm.cpp
1313
FORCE_gamma_tvnl.cpp
1414
FORCE_gamma_vl.cpp
15-
FORCE_gamma_vl_new.cpp
1615
FORCE_k.cpp
1716
LCAO_diago.cpp
1817
LCAO_evolve.cpp
@@ -31,10 +30,8 @@ list(APPEND objects
3130
dftu_relax.cpp
3231
dftu_yukawa.cpp
3332
gint_gamma.cpp
34-
gint_gamma_common.cpp
3533
gint_gamma_env.cpp
3634
gint_gamma_fvl.cpp
37-
gint_gamma_fvl_new.cpp
3835
gint_gamma_mull.cpp
3936
gint_gamma_rho.cpp
4037
gint_gamma_vl.cpp
@@ -45,7 +42,6 @@ list(APPEND objects
4542
gint_k_vl.cpp
4643
gint_tools.cpp
4744
global_fp.cpp
48-
grid_base_beta.cpp
4945
grid_bigcell.cpp
5046
grid_meshball.cpp
5147
grid_meshcell.cpp

‎source/src_lcao/FORCE_gamma.cpp‎

Lines changed: 3 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -50,10 +50,9 @@ void Force_LCAO_gamma::ftable_gamma (
5050

5151
this->cal_ftvnl_dphi(loc.dm_gamma, isforce, isstress, ftvnl_dphi, stvnl_dphi);
5252
this->calFvnlDbeta(loc.dm_gamma, isforce, isstress, fvnl_dbeta, svnl_dbeta, GlobalV::vnl_method);
53-
this->cal_fvl_dphi(loc.dm_gamma, isforce, isstress, fvl_dphi, svl_dphi);
54-
55-
//this->cal_fvl_dphi_new(loc.DM, isforce, isstress, fvl_dphi, svl_dphi);
56-
53+
54+
this->cal_fvl_dphi(loc.DM, isforce, isstress, fvl_dphi, svl_dphi);
55+
5756
//caoyu add for DeePKS
5857
#ifdef __DEEPKS
5958
if (GlobalV::deepks_scf)

‎source/src_lcao/FORCE_gamma.h‎

Lines changed: 0 additions & 14 deletions
Original file line numberDiff line numberDiff line change
@@ -115,20 +115,6 @@ class Force_LCAO_gamma
115115
// forces related to local pseudopotentials
116116
//-------------------------------------------
117117
void cal_fvl_dphi(
118-
ModuleBase::matrix& dm2d,
119-
const bool isforce,
120-
const bool isstress,
121-
ModuleBase::matrix& fvl_dphi,
122-
ModuleBase::matrix& svl_dphi);
123-
124-
void cal_fvl_dphi(
125-
const std::vector<ModuleBase::matrix> &dm2d,
126-
const bool isforce,
127-
const bool isstress,
128-
ModuleBase::matrix& fvl_dphi,
129-
ModuleBase::matrix& svl_dphi);
130-
131-
void cal_fvl_dphi_new(
132118
double*** DM_in,
133119
const bool isforce,
134120
const bool isstress,

‎source/src_lcao/FORCE_gamma_vl.cpp‎

Lines changed: 7 additions & 215 deletions
Original file line numberDiff line numberDiff line change
@@ -3,247 +3,39 @@
33
#include "../module_base/timer.h"
44

55
void Force_LCAO_gamma::cal_fvl_dphi(
6-
ModuleBase::matrix& dm2d,
6+
double*** DM_in,
77
const bool isforce,
88
const bool isstress,
99
ModuleBase::matrix& fvl_dphi,
1010
ModuleBase::matrix& svl_dphi)
1111
{
1212
ModuleBase::TITLE("Force_LCAO_gamma","cal_fvl_dphi");
1313
ModuleBase::timer::tick("Force_LCAO_gamma","cal_fvl_dphi");
14-
15-
if(isforce)
16-
{
17-
ModuleBase::GlobalFunc::ZEROS (this->UHM->LM->DHloc_fixed_x, this->ParaV->nloc);
18-
ModuleBase::GlobalFunc::ZEROS (this->UHM->LM->DHloc_fixed_y, this->ParaV->nloc);
19-
ModuleBase::GlobalFunc::ZEROS (this->UHM->LM->DHloc_fixed_z, this->ParaV->nloc);
20-
}
21-
if(isstress)
22-
{
23-
ModuleBase::GlobalFunc::ZEROS (this->UHM->LM->DHloc_fixed_11, this->ParaV->nloc);
24-
ModuleBase::GlobalFunc::ZEROS (this->UHM->LM->DHloc_fixed_12, this->ParaV->nloc);
25-
ModuleBase::GlobalFunc::ZEROS (this->UHM->LM->DHloc_fixed_13, this->ParaV->nloc);
26-
ModuleBase::GlobalFunc::ZEROS (this->UHM->LM->DHloc_fixed_22, this->ParaV->nloc);
27-
ModuleBase::GlobalFunc::ZEROS (this->UHM->LM->DHloc_fixed_23, this->ParaV->nloc);
28-
ModuleBase::GlobalFunc::ZEROS (this->UHM->LM->DHloc_fixed_33, this->ParaV->nloc);
29-
}
30-
31-
32-
double* tmpDHx = new double[this->ParaV->nloc];
33-
double* tmpDHy = new double[this->ParaV->nloc];
34-
double* tmpDHz = new double[this->ParaV->nloc];
35-
ModuleBase::GlobalFunc::ZEROS( tmpDHx, this->ParaV->nloc );
36-
ModuleBase::GlobalFunc::ZEROS( tmpDHy, this->ParaV->nloc );
37-
ModuleBase::GlobalFunc::ZEROS( tmpDHz, this->ParaV->nloc );
38-
for(int i=0; i<this->ParaV->nloc; ++i)
39-
{
40-
tmpDHx[i] = this->UHM->LM->DHloc_fixed_x[i];
41-
tmpDHy[i] = this->UHM->LM->DHloc_fixed_y[i];
42-
tmpDHz[i] = this->UHM->LM->DHloc_fixed_z[i];
43-
}
44-
45-
//calculate dVL
46-
//calculate <dphi | VL | phi>
47-
48-
// mohan add 2021, needs reconstruction!!!
4914
int istep = 1;
5015
GlobalC::pot.init_pot(istep, GlobalC::pw.strucFac);
51-
16+
fvl_dphi.zero_out();
17+
svl_dphi.zero_out();
5218
for(int is=0; is<GlobalV::NSPIN; ++is)
5319
{
5420
GlobalV::CURRENT_SPIN = is;
55-
56-
ModuleBase::GlobalFunc::ZEROS (this->UHM->LM->DHloc_fixed_x, this->ParaV->nloc);
57-
ModuleBase::GlobalFunc::ZEROS (this->UHM->LM->DHloc_fixed_y, this->ParaV->nloc);
58-
ModuleBase::GlobalFunc::ZEROS (this->UHM->LM->DHloc_fixed_z, this->ParaV->nloc);
59-
6021
for(int ir=0; ir<GlobalC::pw.nrxx; ++ir)
6122
{
6223
GlobalC::pot.vr_eff1[ir] = GlobalC::pot.vr_eff(GlobalV::CURRENT_SPIN, ir);
6324
}
6425

65-
this->UHM->GG.cal_force(GlobalC::pot.vr_eff1);
66-
67-
for(int i=0; i<GlobalV::NLOCAL; i++)
68-
{
69-
const int iat = GlobalC::ucell.iwt2iat[i];
70-
for(int j=0; j<GlobalV::NLOCAL; j++)
71-
{
72-
//const int iat2 = GlobalC::ucell.iwt2iat[j];
73-
const int mu = this->ParaV->trace_loc_row[j];
74-
const int nu = this->ParaV->trace_loc_col[i];
75-
if (mu >= 0 && nu >= 0 )
76-
{
77-
const int index = mu * this->ParaV->ncol + nu;
78-
//contribution from deriv of AO's in T+VNL term
79-
80-
double dm2d2 = 2.0 * dm2d(is, index);
26+
this->UHM->GG.cal_force(DM_in[is], GlobalC::pot.vr_eff1, fvl_dphi, svl_dphi, isforce, isstress);
27+
}
8128

82-
if(isforce)
83-
{
84-
fvl_dphi(iat,0) -= dm2d2 * ( this->UHM->LM->DHloc_fixed_x[index] + tmpDHx[index] );
85-
fvl_dphi(iat,1) -= dm2d2 * ( this->UHM->LM->DHloc_fixed_y[index] + tmpDHy[index] );
86-
fvl_dphi(iat,2) -= dm2d2 * ( this->UHM->LM->DHloc_fixed_z[index] + tmpDHz[index] );
87-
}
88-
if(isstress)
89-
{
90-
svl_dphi(0,0) += dm2d2 * this->UHM->LM->DHloc_fixed_11[index];
91-
svl_dphi(0,1) += dm2d2 * this->UHM->LM->DHloc_fixed_12[index];
92-
svl_dphi(0,2) += dm2d2 * this->UHM->LM->DHloc_fixed_13[index];
93-
svl_dphi(1,1) += dm2d2 * this->UHM->LM->DHloc_fixed_22[index];
94-
svl_dphi(1,2) += dm2d2 * this->UHM->LM->DHloc_fixed_23[index];
95-
svl_dphi(2,2) += dm2d2 * this->UHM->LM->DHloc_fixed_33[index];
96-
}
97-
}
98-
}
99-
}
100-
} // end spin
10129
if(isstress)
10230
{
10331
for(int i=0;i<3;i++)
10432
{
10533
for(int j=0;j<3;j++)
10634
{
10735
if(i<j) svl_dphi(j,i) = svl_dphi(i,j);
108-
svl_dphi(i,j) /= GlobalC::ucell.omega;
36+
svl_dphi(i,j) /= -GlobalC::ucell.omega;
10937
}
11038
}
11139
}
112-
113-
delete[] tmpDHx;
114-
delete[] tmpDHy;
115-
delete[] tmpDHz;
116-
117-
11840
ModuleBase::timer::tick("Force_LCAO_gamma","cal_fvl_dphi");
119-
return;
120-
}
121-
122-
123-
124-
void Force_LCAO_gamma::cal_fvl_dphi(
125-
const std::vector<ModuleBase::matrix> &dm2d,
126-
const bool isforce,
127-
const bool isstress,
128-
ModuleBase::matrix& fvl_dphi,
129-
ModuleBase::matrix& svl_dphi)
130-
{
131-
ModuleBase::TITLE("Force_LCAO_gamma","cal_fvl_dphi");
132-
ModuleBase::timer::tick("Force_LCAO_gamma","cal_fvl_dphi");
133-
134-
ModuleBase::GlobalFunc::ZEROS (this->UHM->LM->DHloc_fixed_x, this->ParaV->nloc);
135-
ModuleBase::GlobalFunc::ZEROS (this->UHM->LM->DHloc_fixed_y, this->ParaV->nloc);
136-
ModuleBase::GlobalFunc::ZEROS (this->UHM->LM->DHloc_fixed_z, this->ParaV->nloc);
137-
if(GlobalV::CAL_STRESS)
138-
{
139-
ModuleBase::GlobalFunc::ZEROS (this->UHM->LM->DHloc_fixed_11, this->ParaV->nloc);
140-
ModuleBase::GlobalFunc::ZEROS (this->UHM->LM->DHloc_fixed_12, this->ParaV->nloc);
141-
ModuleBase::GlobalFunc::ZEROS (this->UHM->LM->DHloc_fixed_13, this->ParaV->nloc);
142-
ModuleBase::GlobalFunc::ZEROS (this->UHM->LM->DHloc_fixed_22, this->ParaV->nloc);
143-
ModuleBase::GlobalFunc::ZEROS (this->UHM->LM->DHloc_fixed_23, this->ParaV->nloc);
144-
ModuleBase::GlobalFunc::ZEROS (this->UHM->LM->DHloc_fixed_33, this->ParaV->nloc);
145-
}
146-
147-
148-
double* tmpDHx = new double[this->ParaV->nloc];
149-
double* tmpDHy = new double[this->ParaV->nloc];
150-
double* tmpDHz = new double[this->ParaV->nloc];
151-
ModuleBase::GlobalFunc::ZEROS( tmpDHx, this->ParaV->nloc );
152-
ModuleBase::GlobalFunc::ZEROS( tmpDHy, this->ParaV->nloc );
153-
ModuleBase::GlobalFunc::ZEROS( tmpDHz, this->ParaV->nloc );
154-
for(int i=0; i<this->ParaV->nloc; ++i)
155-
{
156-
tmpDHx[i] = this->UHM->LM->DHloc_fixed_x[i];
157-
tmpDHy[i] = this->UHM->LM->DHloc_fixed_y[i];
158-
tmpDHz[i] = this->UHM->LM->DHloc_fixed_z[i];
159-
//std::cout << " this->UHM->LM->DHloc_fixed_x=" << this->UHM->LM->DHloc_fixed_x[i] << std::endl;
160-
//std::cout << " this->UHM->LM->DHloc_fixed_y=" << this->UHM->LM->DHloc_fixed_y[i] << std::endl;
161-
//std::cout << " this->UHM->LM->DHloc_fixed_z=" << this->UHM->LM->DHloc_fixed_z[i] << std::endl;
162-
}
163-
164-
//calculate dVL
165-
//calculate <dphi | VL | phi>
166-
167-
int istep = 1;
168-
GlobalC::pot.init_pot(istep, GlobalC::pw.strucFac);
169-
170-
for(int is=0; is<GlobalV::NSPIN; ++is)
171-
{
172-
GlobalV::CURRENT_SPIN = is;
173-
174-
ModuleBase::GlobalFunc::ZEROS (this->UHM->LM->DHloc_fixed_x, this->ParaV->nloc);
175-
ModuleBase::GlobalFunc::ZEROS (this->UHM->LM->DHloc_fixed_y, this->ParaV->nloc);
176-
ModuleBase::GlobalFunc::ZEROS (this->UHM->LM->DHloc_fixed_z, this->ParaV->nloc);
177-
178-
for(int ir=0; ir<GlobalC::pw.nrxx; ++ir)
179-
{
180-
GlobalC::pot.vr_eff1[ir] = GlobalC::pot.vr_eff(GlobalV::CURRENT_SPIN, ir);
181-
}
182-
183-
//should not be set zero if VNA is used.
184-
//ModuleBase::GlobalFunc::ZEROS(this->UHM->LM->DHloc_fixed_x,this->ParaV->nloc);
185-
//ModuleBase::GlobalFunc::ZEROS(this->UHM->LM->DHloc_fixed_y,this->ParaV->nloc);
186-
//ModuleBase::GlobalFunc::ZEROS(this->UHM->LM->DHloc_fixed_z,this->ParaV->nloc);
187-
this->UHM->GG.cal_force(GlobalC::pot.vr_eff1);
188-
189-
for(int i=0; i<GlobalV::NLOCAL; i++)
190-
{
191-
const int iat = GlobalC::ucell.iwt2iat[i];
192-
for(int j=0; j<GlobalV::NLOCAL; j++)
193-
{
194-
//const int iat2 = GlobalC::ucell.iwt2iat[j];
195-
const int mu = this->ParaV->trace_loc_row[j];
196-
const int nu = this->ParaV->trace_loc_col[i];
197-
if (mu >= 0 && nu >= 0 )
198-
{
199-
const int index = mu * this->ParaV->ncol + nu;
200-
//contribution from deriv of AO's in T+VNL term
201-
202-
double dm2d2 = 2.0 * dm2d[is](nu, mu);
203-
204-
if(isforce)
205-
{
206-
fvl_dphi(iat,0) -= dm2d2 * ( this->UHM->LM->DHloc_fixed_x[index] + tmpDHx[index] );
207-
fvl_dphi(iat,1) -= dm2d2 * ( this->UHM->LM->DHloc_fixed_y[index] + tmpDHy[index] );
208-
fvl_dphi(iat,2) -= dm2d2 * ( this->UHM->LM->DHloc_fixed_z[index] + tmpDHz[index] );
209-
}
210-
if(isstress)
211-
{
212-
svl_dphi(0,0) += dm2d2 * this->UHM->LM->DHloc_fixed_11[index];
213-
svl_dphi(0,1) += dm2d2 * this->UHM->LM->DHloc_fixed_12[index];
214-
svl_dphi(0,2) += dm2d2 * this->UHM->LM->DHloc_fixed_13[index];
215-
svl_dphi(1,1) += dm2d2 * this->UHM->LM->DHloc_fixed_22[index];
216-
svl_dphi(1,2) += dm2d2 * this->UHM->LM->DHloc_fixed_23[index];
217-
svl_dphi(2,2) += dm2d2 * this->UHM->LM->DHloc_fixed_33[index];
218-
}
219-
//std::cout << std::setw(5) << iat << std::setw(5) << iat2
220-
//<< std::setw(5) << mu << std::setw(5) << nu
221-
//<< std::setw(15) << this->UHM->LM->DHloc_fixed_z[index] << std::endl;
222-
}
223-
}
224-
}
225-
226-
//std::cout << "fvl_dphi:" << std::endl;
227-
//for(int iat=0; iat<GlobalC::ucell.nat; ++iat)
228-
//{
229-
//std::cout << std::setw(5) << iat << std::setw(15) << fvl_dphi[iat][0]
230-
//<< std::setw(15) << fvl_dphi[iat][1]
231-
//<< std::setw(15) << fvl_dphi[iat][2] << std::endl;
232-
//}
233-
234-
} // end spin
235-
// test mohan tmp
236-
//test_gamma(this->UHM->LM->DHloc_fixed_x,"this->UHM->LM->DHloc_fixed_x");
237-
238-
if(isstress)
239-
{
240-
StressTools::stress_fill(1.0, GlobalC::ucell.omega, svl_dphi);
241-
}
242-
243-
delete[] tmpDHx;
244-
delete[] tmpDHy;
245-
delete[] tmpDHz;
246-
247-
ModuleBase::timer::tick("Force_LCAO_gamma", "cal_fvl_dphi");
248-
return;
249-
}
41+
}

‎source/src_lcao/FORCE_gamma_vl_new.cpp‎

Lines changed: 0 additions & 28 deletions
This file was deleted.

‎source/src_lcao/LCAO_hamilt.cpp‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -37,7 +37,7 @@ void LCAO_Hamilt::set_lcao_matrices(void)
3737
{
3838
// mohan add 2012-03-29
3939
// calculate the grid integration of 'Vl' matrix for gamma algorithms.
40-
this->GG.prepare(GlobalC::ucell.latvec, GlobalC::ucell.lat0, this->LM);
40+
this->GG.prepare(GlobalC::ucell.latvec, GlobalC::ucell.lat0);
4141

4242
// calulate the 'S', 'T' and 'Vnl' matrix for gamma algorithms.
4343
this->calculate_STNR_gamma();
@@ -87,7 +87,7 @@ void LCAO_Hamilt::calculate_Hgamma( const int &ik , vector<ModuleBase::matrix> d
8787
// calculate the 'Vl' matrix using gamma-algorithms.
8888
if(GlobalV::VL_IN_H)
8989
{
90-
this->GG.cal_vlocal(GlobalC::pot.vr_eff1);
90+
this->GG.cal_vlocal(GlobalC::pot.vr_eff1, *this->LM);
9191

9292
#ifdef __MPI //liyuanbo 2022/2/23
9393
// Peize Lin add 2016-12-03

0 commit comments

Comments
 (0)