@@ -86,188 +86,11 @@ void Gint_k::destroy_pvpR(void)
8686 return ;
8787}
8888
89-
90-
91-
92-
9389// fold the <phi | vl |dphi(R)> * DM(R) to
9490// calculate the force.
95- void Gint_k::folding_force (
96- ModuleBase::matrix& fvl_dphi,
97- double * pvdpx,
98- double * pvdpy,
99- double * pvdpz)
100- {
101- ModuleBase::TITLE (" Gint_k" ," folding_force" );
102- ModuleBase::timer::tick (" Gint_k" ," folding_force" );
103-
104- // xiaohui modify 2013-12-17, test
105- // assert(GlobalC::GridT.lgd > 0); //mohan add 2012-06-10
106-
107- // mohan add 2014-01-20
108- const int lgd = GlobalC::GridT.lgd ;
109-
110- double ** ppx;
111- double ** ppy;
112- double ** ppz;
113-
114- if (GlobalC::GridT.lgd >0 )
115- {
116- ppx = new double *[lgd];
117- ppy = new double *[lgd];
118- ppz = new double *[lgd];
119- for (int i=0 ; i<lgd; i++)
120- {
121- ppx[i] = new double [lgd];
122- ppy[i] = new double [lgd];
123- ppz[i] = new double [lgd];
124- ModuleBase::GlobalFunc::ZEROS ( ppx[i], lgd);
125- ModuleBase::GlobalFunc::ZEROS ( ppy[i], lgd);
126- ModuleBase::GlobalFunc::ZEROS ( ppz[i], lgd);
127- }
128- }
129-
130- ModuleBase::Vector3<double > tau1, dtau;
131- for (int T1=0 ; T1<GlobalC::ucell.ntype ; ++T1)
132- {
133- Atom* atom1 = &GlobalC::ucell.atoms [T1];
134- for (int I1=0 ; I1< atom1->na ; ++I1)
135- {
136- const int iat = GlobalC::ucell.itia2iat (T1,I1);
137- if (GlobalC::GridT.in_this_processor [iat])
138- {
139- assert ( lgd > 0 );
140-
141- const int start1 = GlobalC::ucell.itiaiw2iwt (T1, I1, 0 );
142- // get the start positions of elements.
143- const int DM_start = GlobalC::GridT.nlocstartg [iat];
144- // get the coordinates of adjacent atoms.
145- tau1 = atom1->tau [I1];
146- // GlobalC::GridD.Find_atom(tau1);
147- GlobalC::GridD.Find_atom (GlobalC::ucell, tau1, T1, I1);
148- // search for the adjacent atoms.
149- int nad = 0 ;
150- for (int ad = 0 ; ad < GlobalC::GridD.getAdjacentNum ()+1 ; ++ad)
151- {
152- // get iat2
153- const int T2 = GlobalC::GridD.getType (ad);
154- const int I2 = GlobalC::GridD.getNatom (ad);
155- const int iat2 = GlobalC::ucell.itia2iat (T2, I2);
156- if (GlobalC::GridT.in_this_processor [iat2])
157- {
158- Atom* atom2 = &GlobalC::ucell.atoms [T2];
159- dtau = GlobalC::GridD.getAdjacentTau (ad) - tau1;
160- double distance = dtau.norm () * GlobalC::ucell.lat0 ;
161- double rcut = GlobalC::ORB.Phi [T1].getRcut () + GlobalC::ORB.Phi [T2].getRcut ();
162- if (distance < rcut)
163- {
164- const int start2 = GlobalC::ucell.itiaiw2iwt (T2, I2, 0 );
165- int ixxx = DM_start + GlobalC::GridT.find_R2st [iat][nad];
166- for (int iw=0 ; iw<atom1->nw ; iw++)
167- {
168- const int iw_all = start1+iw;
169- const int iw_local = GlobalC::GridT.trace_lo [iw_all];
170- // iw1_lo
171- double *vijx = ppx[iw_local];
172- double *vijy = ppy[iw_local];
173- double *vijz = ppz[iw_local];
174-
175- double *vRx = &pvdpx[ixxx]; // just fold R to normal matrix.
176- double *vRy = &pvdpy[ixxx];
177- double *vRz = &pvdpz[ixxx];
178-
179- int * iw2_lo = &GlobalC::GridT.trace_lo [start2];
180- int * iw2_end = iw2_lo + atom2->nw ;
181-
182- for (; iw2_lo<iw2_end; ++iw2_lo, ++vRx, ++vRy, ++vRz)
183- {
184- vijx[iw2_lo[0 ]] += vRx[0 ] ;
185- vijy[iw2_lo[0 ]] += vRy[0 ] ;
186- vijz[iw2_lo[0 ]] += vRz[0 ] ;
187- }
188- ixxx += atom2->nw ;
189- }
190- ++nad;
191- }// end distance<rcut
192- }
193- }// end ad
194- }
195- }// end ia
196- }// end it
197-
198-
199-
200- double * tmp = new double [GlobalV::NLOCAL*3 ];
201- for (int i=0 ; i<GlobalV::NLOCAL; ++i)
202- {
203- ModuleBase::GlobalFunc::ZEROS (tmp, 3 *GlobalV::NLOCAL);
204- const int mug = GlobalC::GridT.trace_lo [i];
205- // if the row element is on this processor
206- if (mug>=0 )
207- {
208- // GlobalV::ofs_running << " i=" << i << " mug=" << mug << std::endl;
209- for (int j=0 ; j<GlobalV::NLOCAL; ++j)
210- {
211- const int nug = GlobalC::GridT.trace_lo [j];
212- // if the col element is on this processor
213- if (nug>=0 )
214- {
215- // if(mug<nug)
216- // {
217- const int index = 3 *j;
218- tmp[index] = ppx[mug][nug];
219- tmp[index+1 ] = ppy[mug][nug];
220- tmp[index+2 ] = ppz[mug][nug];
221- // }
222- // else
223- // {
224- // tmpx[j] = 0.0;
225- // tmpy[j] = 0.0;
226- // tmpz[j] = 0.0;
227- // }
228- }
229- }
230- }
231- // collect the matrix after folding.
232- Parallel_Reduce::reduce_double_pool ( tmp, GlobalV::NLOCAL*3 );
233- for (int j=0 ; j<GlobalV::NLOCAL; j++)
234- {
235- if (!this ->LM ->ParaV ->in_this_processor (i,j))
236- {
237- continue ;
238- }
239- const int iat = GlobalC::ucell.iwt2iat [i];
240- const int index = 3 *j;
241- fvl_dphi (iat,0 ) += 2.0 *tmp[index];
242- fvl_dphi (iat,1 ) += 2.0 *tmp[index+1 ];
243- fvl_dphi (iat,2 ) += 2.0 *tmp[index+2 ];
244- }
245- }
246- delete[] tmp;
247-
248- // mohan add 2014-01-20
249- if (GlobalC::GridT.lgd > 0 )
250- {
251- // -------------------------
252- // delete the tmp matrix.
253- // -------------------------
254- for (int i=0 ; i<GlobalC::GridT.lgd ; i++)
255- {
256- delete[] ppx[i];
257- delete[] ppy[i];
258- delete[] ppz[i];
259- }
260- delete[] ppx;
261- delete[] ppy;
262- delete[] ppz;
263- }
264- ModuleBase::timer::tick (" Gint_k" ," folding_force" );
265- return ;
266- }
267-
26891// fold the <phi | vl * R_beta|dphi(R_alpha)> * DM(R) to
26992// calculate the stress.
270- void Gint_k::folding_stress (
93+ void Gint_k::folding_force (
27194 const bool isforce,
27295 const bool isstress,
27396 ModuleBase::matrix& fvl_dphi,
@@ -282,7 +105,7 @@ void Gint_k::folding_stress(
282105 double * pvdp13,
283106 double * pvdp23)
284107{
285- ModuleBase::TITLE (" Gint_k" ," folding_stress " );
108+ ModuleBase::TITLE (" Gint_k" ," folding_force " );
286109 if (!isforce&&!isstress) return ;
287110
288111 const int lgd = GlobalC::GridT.lgd ;
0 commit comments