33#include " ../module_base/timer.h"
44
55void 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-
118- 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" );
13240 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- if (GlobalV::CAL_STRESS)
178- {
179- ModuleBase::GlobalFunc::ZEROS (this ->UHM ->LM ->DHloc_fixed_11 , this ->ParaV ->nloc );
180- ModuleBase::GlobalFunc::ZEROS (this ->UHM ->LM ->DHloc_fixed_12 , this ->ParaV ->nloc );
181- ModuleBase::GlobalFunc::ZEROS (this ->UHM ->LM ->DHloc_fixed_13 , this ->ParaV ->nloc );
182- ModuleBase::GlobalFunc::ZEROS (this ->UHM ->LM ->DHloc_fixed_22 , this ->ParaV ->nloc );
183- ModuleBase::GlobalFunc::ZEROS (this ->UHM ->LM ->DHloc_fixed_23 , this ->ParaV ->nloc );
184- ModuleBase::GlobalFunc::ZEROS (this ->UHM ->LM ->DHloc_fixed_33 , this ->ParaV ->nloc );
185- }
186- for (int ir=0 ; ir<GlobalC::pw.nrxx ; ++ir)
187- {
188- GlobalC::pot.vr_eff1 [ir] = GlobalC::pot.vr_eff (GlobalV::CURRENT_SPIN, ir);
189- }
190-
191- // should not be set zero if VNA is used.
192- // ModuleBase::GlobalFunc::ZEROS(this->UHM->LM->DHloc_fixed_x,this->ParaV->nloc);
193- // ModuleBase::GlobalFunc::ZEROS(this->UHM->LM->DHloc_fixed_y,this->ParaV->nloc);
194- // ModuleBase::GlobalFunc::ZEROS(this->UHM->LM->DHloc_fixed_z,this->ParaV->nloc);
195- this ->UHM ->GG .cal_force (GlobalC::pot.vr_eff1 );
196-
197- for (int i=0 ; i<GlobalV::NLOCAL; i++)
198- {
199- const int iat = GlobalC::ucell.iwt2iat [i];
200- for (int j=0 ; j<GlobalV::NLOCAL; j++)
201- {
202- // const int iat2 = GlobalC::ucell.iwt2iat[j];
203- const int mu = this ->ParaV ->trace_loc_row [j];
204- const int nu = this ->ParaV ->trace_loc_col [i];
205- if (mu >= 0 && nu >= 0 )
206- {
207- const int index = mu * this ->ParaV ->ncol + nu;
208- // contribution from deriv of AO's in T+VNL term
209-
210- double dm2d2 = 2.0 * dm2d[is](nu, mu);
211-
212- if (isforce)
213- {
214- fvl_dphi (iat,0 ) -= dm2d2 * ( this ->UHM ->LM ->DHloc_fixed_x [index] + tmpDHx[index] );
215- fvl_dphi (iat,1 ) -= dm2d2 * ( this ->UHM ->LM ->DHloc_fixed_y [index] + tmpDHy[index] );
216- fvl_dphi (iat,2 ) -= dm2d2 * ( this ->UHM ->LM ->DHloc_fixed_z [index] + tmpDHz[index] );
217- }
218- if (isstress)
219- {
220- svl_dphi (0 ,0 ) += dm2d2 * this ->UHM ->LM ->DHloc_fixed_11 [index];
221- svl_dphi (0 ,1 ) += dm2d2 * this ->UHM ->LM ->DHloc_fixed_12 [index];
222- svl_dphi (0 ,2 ) += dm2d2 * this ->UHM ->LM ->DHloc_fixed_13 [index];
223- svl_dphi (1 ,1 ) += dm2d2 * this ->UHM ->LM ->DHloc_fixed_22 [index];
224- svl_dphi (1 ,2 ) += dm2d2 * this ->UHM ->LM ->DHloc_fixed_23 [index];
225- svl_dphi (2 ,2 ) += dm2d2 * this ->UHM ->LM ->DHloc_fixed_33 [index];
226- }
227- // std::cout << std::setw(5) << iat << std::setw(5) << iat2
228- // << std::setw(5) << mu << std::setw(5) << nu
229- // << std::setw(15) << this->UHM->LM->DHloc_fixed_z[index] << std::endl;
230- }
231- }
232- }
233-
234- // std::cout << "fvl_dphi:" << std::endl;
235- // for(int iat=0; iat<GlobalC::ucell.nat; ++iat)
236- // {
237- // std::cout << std::setw(5) << iat << std::setw(15) << fvl_dphi[iat][0]
238- // << std::setw(15) << fvl_dphi[iat][1]
239- // << std::setw(15) << fvl_dphi[iat][2] << std::endl;
240- // }
241-
242- } // end spin
243- // test mohan tmp
244- // test_gamma(this->UHM->LM->DHloc_fixed_x,"this->UHM->LM->DHloc_fixed_x");
245-
246- if (isstress)
247- {
248- StressTools::stress_fill (1.0 , GlobalC::ucell.omega , svl_dphi);
249- }
250-
251- delete[] tmpDHx;
252- delete[] tmpDHy;
253- delete[] tmpDHz;
254-
255- ModuleBase::timer::tick (" Force_LCAO_gamma" , " cal_fvl_dphi" );
256- return ;
257- }
41+ }
0 commit comments