Line data Source code
1 : !!****m* ABINIT/m_opernld_ylm_allwf
2 : !! NAME
3 : !! m_opernld_ylm_allwf
4 : !!
5 : !! FUNCTION
6 : !!
7 : !!
8 : !! COPYRIGHT
9 : !! Copyright (C) 2008-2026 ABINIT group (MT)
10 : !! This file is distributed under the terms of the
11 : !! GNU General Public License, see ~abinit/COPYING
12 : !! or http://www.gnu.org/copyleft/gpl.txt .
13 : !!
14 : !! SOURCE
15 :
16 : #if defined HAVE_CONFIG_H
17 : #include "config.h"
18 : #endif
19 :
20 : #include "abi_common.h"
21 :
22 : module m_opernld_ylm_allwf
23 :
24 : use defs_basis
25 : use m_errors
26 : use m_abicore
27 : use m_xmpi
28 : use m_geometry, only : strconv
29 : use defs_abitypes, only : MPI_type
30 :
31 : implicit none
32 :
33 : private
34 : !!***
35 :
36 : public :: opernld_ylm_allwf
37 : !!***
38 :
39 : contains
40 : !!***
41 :
42 : !!****f* ABINIT/opernld_ylm_allwf
43 : !! NAME
44 : !! opernld_ylm_allwf
45 : !!
46 : !! FUNCTION
47 : !! * Operate with the non-local part of the hamiltonian,
48 : !! in order to get contributions to energy/forces/stress/dyn.matrix/elst tens.
49 : !! from projected scalars
50 : !! * Operate with the non-local projectors and the overlap matrix Sij
51 : !! in order to get contributions to <c|S|c>
52 : !! from projected scalars
53 : !!
54 : !! INPUTS
55 : !! choice=chooses possible output
56 : !! cplex=1 if <p_lmn|c> scalars are real (equivalent to istwfk>1)
57 : !! 2 if <p_lmn|c> scalars are complex
58 : !! cplex_fac=1 if gxfac scalars are real, 2 if gxfac scalars are complex
59 : !! d2gxdt(cplex,nd2gxdt,nlmn,nincat,nspinor)=2nd gradients of projected scalars
60 : !! dgxdt(cplex,ndgxdt,nlmn,nincat,nspinor)=gradients of projected scalars
61 : !! dgxdtfac(cplex_fac,ndgxdtfac,nlmn,nincat,nspinor)=gradients of reduced projected scalars
62 : !! related to Vnl (NL operator)
63 : !! dgxdtfac_sij(cplex_fac,ndgxdtfac,nlmn,nincat,nspinor)=gradients of reduced projected scalars
64 : !! related to Sij (overlap)
65 : !! gx(cplex,nlmn,nincat,nspinor)= projected scalars
66 : !! gxfac(cplex_fac,nlmn,nincat,nspinor)= reduced projected scalars related to Vnl (NL operator)
67 : !! gxfac_sij(cplex,nlmn,nincat,nspinor)= reduced projected scalars related to Sij (overlap)
68 : !! ia3=gives the absolute number of the first atom in the subset presently treated
69 : !! natom=number of atoms in cell
70 : !! nd2gxdt=second dimension of d2gxdt
71 : !! ndgxdt=second dimension of dgxdt
72 : !! ndgxdtfac=second dimension of dgxdtfac
73 : !! nincat=number of atoms in the subset here treated
74 : !! nlmn=number of (l,m,n) numbers for current type of atom
75 : !! nnlout=dimension of enlout
76 : !! nspinor=number of spinorial components of the wavefunctions (on current proc)
77 : !! paw_opt= define the nonlocal operator concerned with:
78 : !! paw_opt=0 : Norm-conserving Vnl (use of Kleinman-Bylander ener.)
79 : !! paw_opt=1 : PAW nonlocal part of H (use of Dij coeffs)
80 : !! paw_opt=2 : PAW: (Vnl-lambda.Sij) (Sij=overlap matrix)
81 : !! paw_opt=3 : PAW overlap matrix (Sij)
82 : !! paw_opt=4 : both PAW nonlocal part of H (Dij) and overlap matrix (Sij)
83 : !!
84 : !! OUTPUT
85 : !! (see side effects)
86 : !!
87 : !! SIDE EFFECTS
88 : !! --If (paw_opt==0, 1 or 2)
89 : !! enlout(nnlout)= contribution to the non-local part of the following properties:
90 : !! if choice=1 : enlout(1) -> the energy
91 : !! if choice=2 : enlout(3*natom) -> 1st deriv. of energy wrt atm. pos (forces)
92 : !! if choice=3 : enlout(6) -> 1st deriv. of energy wrt strain (stresses)
93 : !! if choice=4 : enlout(6*natom) -> 2nd deriv. of energy wrt 2 atm. pos (dyn. mat.)
94 : !! if choice=23: enlout(6+3*natom) -> 1st deriv. of energy wrt atm. pos (forces) and
95 : !! 1st deriv. of energy wrt strain (stresses)
96 : !! if choice=24: enlout(9*natom) -> 1st deriv. of energy wrt atm. pos (forces) and
97 : !! 2nd deriv. of energy wrt 2 atm. pos (dyn. mat.)
98 : !! if choice=5 : enlout(3) -> 1st deriv. of energy wrt k
99 : !! if choice=53: enlout(3) -> 1st deriv. (twist) of energy wrt k
100 : !! if choice=54: enlout(18*natom) -> 2nd deriv. of energy wrt atm. pos and right k (Born eff. charge)
101 : !! if choice=55: enlout(36) -> 2nd deriv. of energy wrt strain and right k (piezoelastic tensor)
102 : !! if choice=6 : enlout(36+18*natom) -> 2nd deriv. of energy wrt 2 strains (elast. tensor) and
103 : !! 2nd deriv. of energy wrt to atm. pos and strain (internal strain)
104 : !! if choice=8 : enlout(6) -> 2nd deriv. of energy wrt 2 k
105 : !! if choice=81: enlout(18) -> 2nd deriv. of energy wrt k and right k
106 : !! --If (paw_opt==3)
107 : !! if choice=1 : enlout(1) -> contribution to <c|S|c> (note: not including <c|c>)
108 : !! if choice=2 : enlout(3*natom) -> contribution to <c|dS/d_atm.pos|c>
109 : !! if choice=53: enlout(3) -> 1st deriv. (twist) of energy wrt k
110 : !! if choice=54: enlout(18*natom) -> 2nd deriv. of energy wrt atm. pos and right k (Born eff. charge)
111 : !! if choice=55: enlout(36) -> 2nd deriv. of energy wrt strain and right k (piezoelastic tensor)
112 : !! if choice=8 : enlout(6) -> 2nd deriv. of energy wrt 2 k
113 : !! if choice=81: enlout(18) -> 2nd deriv. of energy wrt k and right k
114 : !! --If (paw_opt==4)
115 : !! not available
116 : !!
117 : !! NOTES
118 : !! Operate for one type of atom, and within this given type of atom,
119 : !! for a subset of at most nincat atoms.
120 : !!
121 : !! SOURCE
122 :
123 520 : subroutine opernld_ylm_allwf(choice,cplex,cplex_fac,ddkk,&
124 520 : & dgxdt,dgxdtfac,dgxdtfac_sij,d2gxdt,&
125 520 : & enlk,enlout,fnlk,gx,gxfac,gxfac_sij,natom,ndat,nd2gxdt,ndgxdt,&
126 520 : & ndgxdtfac,indlmn,ntypat,lmnmax,nprojs,nnlout,nspinor,paw_opt,&
127 520 : & strnlk,nattyp,gpu_option)
128 :
129 : ! Arguments ------------------------------------
130 : ! scalars
131 : integer,intent(in) :: choice,paw_opt,ntypat,ndgxdtfac,nd2gxdt,ndgxdt
132 : integer,intent(in) :: cplex,cplex_fac,natom,ndat,nnlout,nspinor,nprojs,lmnmax,gpu_option
133 :
134 : ! arrays
135 : integer,intent(in) :: indlmn(6,lmnmax,ntypat),nattyp(ntypat)
136 : real(dp),intent(in) :: d2gxdt(cplex,nd2gxdt*nprojs,ndat*nspinor)
137 : real(dp),intent(in) :: dgxdt(cplex,ndgxdt*nprojs,ndat*nspinor)
138 : real(dp),intent(in) :: dgxdtfac(cplex,ndgxdtfac*nprojs,ndat*nspinor)
139 : real(dp),intent(in) :: dgxdtfac_sij(cplex_fac,ndgxdtfac*nprojs,ndat*nspinor)
140 : real(dp),intent(in) :: gx(cplex,nprojs,ndat*nspinor)
141 : real(dp),intent(in) :: gxfac(cplex,nprojs,ndat*nspinor)
142 : real(dp),intent(in) :: gxfac_sij(cplex_fac,nprojs,ndat*nspinor)
143 : real(dp),intent(inout) :: enlout(nnlout*ndat)
144 : real(dp),intent(inout) :: enlk(ndat),ddkk(6,ndat)
145 : real(dp),intent(inout) :: fnlk(3*natom,ndat),strnlk(6,ndat)
146 :
147 : ! locals
148 : integer,parameter :: alpha(6)=(/1,2,3,3,3,2/),beta(6)=(/1,2,3,2,1,1/)
149 : integer,parameter :: gamma(3,3)=reshape((/1,6,5,6,2,4,5,4,3/),(/3,3/))
150 : integer :: force_shift, shift, nattyp_i
151 : integer :: itypat, ilmn, ia, ispinor, idat, igrad, ii, nlmn, iend, ibeg, iatm, iashift
152 : integer :: mua, mub, nu, mu, mua1, mua2, muu, mut, mushift, nushift
153 : real(dp) :: esum,esumi
154 520 : real(dp) :: d2gx(cplex)
155 :
156 : ABI_UNUSED((/gpu_option/))
157 :
158 18568 : enlout=zero
159 520 : if(paw_opt < 3) then
160 520 : if(choice==1.or.choice==3.or.choice==23.or.choice==6) then
161 : shift=0
162 : iatm=0
163 856 : esum=zero
164 856 : do itypat=1, ntypat
165 9732 : nlmn=count(indlmn(3,:,itypat)>0)
166 548 : ibeg = shift+1
167 548 : iend = shift+nattyp(itypat)*nlmn
168 548 : nattyp_i = nattyp(itypat)
169 :
170 : #ifdef HAVE_OPENMP_OFFLOAD
171 : !$OMP TARGET TEAMS DISTRIBUTE &
172 : !$OMP& MAP(to:enlk,gxfac,gx) &
173 : !$OMP& PRIVATE(idat,esum) &
174 : !$OMP& IF(gpu_option==ABI_GPU_OPENMP)
175 : #endif
176 3348 : do idat=1,ndat
177 : esum=zero
178 : !$OMP PARALLEL DO COLLAPSE(4) REDUCTION(+:esum) PRIVATE(ia,ilmn,ii,ispinor)
179 5600 : do ispinor=1,nspinor
180 16736 : do ia=1,nattyp_i
181 174384 : do ilmn=1,nlmn
182 452440 : do ii=1,cplex
183 : esum=esum +gxfac(ii,shift+(ia-1)*nlmn+ilmn,ispinor+nspinor*(idat-1)) &
184 441304 : & *gx (ii,shift+(ia-1)*nlmn+ilmn,ispinor+nspinor*(idat-1))
185 : end do
186 : end do
187 : end do
188 : end do
189 3348 : enlk(idat) = enlk(idat) + esum
190 : end do
191 :
192 548 : shift = shift + nattyp(itypat)*nlmn
193 856 : iatm = iatm+nattyp(itypat)
194 : end do
195 308 : if (choice==1) then
196 : #ifdef HAVE_OPENMP_OFFLOAD
197 : !$OMP TARGET PARALLEL DO MAP(to:enlout,enlk) PRIVATE(idat) &
198 : !$OMP& IF(gpu_option==ABI_GPU_OPENMP)
199 : #endif
200 864 : do idat=1,ndat
201 864 : enlout(idat)=enlk(idat)
202 : end do
203 : end if
204 : end if ! choice=1/3/23
205 :
206 : ! ======== Accumulate the stress tensor contributions ==========
207 520 : if (choice==3.or.choice==23) then
208 : shift=0
209 : iatm=0
210 568 : do itypat=1, ntypat
211 6084 : nlmn=count(indlmn(3,:,itypat)>0)
212 356 : ibeg = shift+1
213 356 : iend = shift+nattyp(itypat)*nlmn
214 356 : nattyp_i = nattyp(itypat)
215 :
216 : #ifdef HAVE_OPENMP_OFFLOAD
217 : !$OMP TARGET TEAMS DISTRIBUTE COLLAPSE(2) &
218 : !$OMP& MAP(to:enlout,gxfac,dgxdt) &
219 : !$OMP& PRIVATE(idat,igrad,esum) &
220 : !$OMP& IF(gpu_option==ABI_GPU_OPENMP)
221 : #endif
222 1620 : do idat=1,ndat
223 9204 : do igrad=1,6
224 : esum=zero
225 : !$OMP PARALLEL DO COLLAPSE(4) REDUCTION(+:esum) PRIVATE(ispinor,ia,ilmn,ii)
226 15168 : do ispinor=1,nspinor
227 40512 : do ia=1,nattyp_i
228 : !Following loops are a [D][Z]DOT
229 387360 : do ilmn=1,nlmn
230 1000464 : do ii=1,cplex
231 : esum=esum +gxfac(ii,shift+(ia-1)*nlmn+ilmn,ispinor+nspinor*(idat-1)) &
232 975120 : & *dgxdt(ii,ndgxdt*shift + (ia-1)*nlmn*ndgxdt + (ilmn-1)*ndgxdt + igrad,ispinor+nspinor*(idat-1))
233 : end do
234 : end do
235 : end do
236 : end do
237 8848 : enlout((idat-1)*nnlout+igrad) = enlout((idat-1)*nnlout+igrad) + two*esum
238 : end do
239 : end do
240 :
241 356 : shift = shift + nattyp(itypat)*nlmn
242 568 : iatm = iatm+nattyp(itypat)
243 : end do
244 : end if
245 :
246 : ! ============ Accumulate the forces contributions =============
247 520 : if (choice==2.or.choice==23) then
248 212 : shift=0
249 212 : iatm=0
250 212 : force_shift=0; if(choice==23) force_shift=6
251 568 : do itypat=1, ntypat
252 6084 : nlmn=count(indlmn(3,:,itypat)>0)
253 356 : nattyp_i = nattyp(itypat)
254 :
255 : #ifdef HAVE_OPENMP_OFFLOAD
256 : !$OMP TARGET TEAMS DISTRIBUTE COLLAPSE(3) &
257 : !$OMP& MAP(to:enlout,gxfac,dgxdt) &
258 : !$OMP& PRIVATE(idat,igrad,ia,esum) &
259 : !$OMP& IF(gpu_option==ABI_GPU_OPENMP)
260 : #endif
261 1620 : do idat=1,ndat
262 5844 : do ia=1,nattyp_i
263 18160 : do igrad=1,3
264 : !Following loops are a [D][Z]DOT
265 : esum=zero
266 : !$OMP PARALLEL DO COLLAPSE(3) REDUCTION(+:esum) PRIVATE(ispinor,ilmn,ii)
267 25344 : do ispinor=1,nspinor
268 202560 : do ilmn=1,nlmn
269 500232 : do ii=1,cplex
270 : esum=esum +gxfac(ii,shift+(ia-1)*nlmn+ilmn,ispinor+nspinor*(idat-1)) &
271 487560 : & *dgxdt(ii,ndgxdt*shift + (ia-1)*nlmn*ndgxdt + (ilmn-1)*ndgxdt + (igrad+force_shift) ,ispinor+nspinor*(idat-1))
272 : end do
273 : end do
274 : end do
275 : enlout((idat-1)*nnlout + force_shift + (iatm+ia-1)*3 + igrad)= &
276 16896 : & enlout((idat-1)*nnlout + force_shift + (iatm+ia-1)*3 + igrad) + two*esum
277 : end do
278 : end do
279 : end do
280 :
281 356 : shift = shift + nattyp(itypat)*nlmn
282 568 : iatm = iatm+nattyp(itypat)
283 : end do
284 : end if
285 :
286 : ! ====== Accumulate the dynamical matrix contributions =========
287 520 : if (choice==4) then
288 : shift=0
289 : iatm=0
290 0 : do itypat=1, ntypat
291 0 : nlmn=count(indlmn(3,:,itypat)>0)
292 0 : nattyp_i = nattyp(itypat)
293 :
294 : #ifdef HAVE_OPENMP_OFFLOAD
295 : !$OMP TARGET TEAMS DISTRIBUTE COLLAPSE(3) &
296 : !$OMP& MAP(to:enlout,gxfac,dgxdt,dgxdtfac,d2gxdt) &
297 : !$OMP& PRIVATE(idat,ia,esum,mu,mua,mub) &
298 : !$OMP& IF(gpu_option==ABI_GPU_OPENMP)
299 : #endif
300 0 : do idat=1,ndat
301 0 : do ia=1,nattyp_i
302 0 : do mu=1,6
303 0 : mua=alpha(mu);mub=beta(mu)
304 0 : esum=zero
305 : !$OMP PARALLEL DO COLLAPSE(3) REDUCTION(+:esum) PRIVATE(ispinor,ilmn,ii)
306 0 : do ispinor=1,nspinor
307 0 : do ilmn=1,nlmn
308 0 : do ii=1,cplex
309 : esum=esum&
310 : & +gxfac(ii,shift+(ia-1)*nlmn+ilmn,ispinor+nspinor*(idat-1))&
311 : & *d2gxdt(ii, nd2gxdt*shift+(ia-1)*nlmn*nd2gxdt+(ilmn-1)*nd2gxdt+mu, ispinor+nspinor*(idat-1))&
312 : & +dgxdtfac(ii, ndgxdt*shift+(ia-1)*nlmn*ndgxdt+(ilmn-1)*ndgxdt+mub, ispinor+nspinor*(idat-1))&
313 0 : & *dgxdt(ii, ndgxdt*shift+(ia-1)*nlmn*ndgxdt+(ilmn-1)*ndgxdt+mua, ispinor+nspinor*(idat-1))
314 : end do
315 : end do
316 : end do
317 : enlout((idat-1)*nnlout + (iatm+ia-1)*6 + mu)= &
318 0 : & enlout((idat-1)*nnlout + (iatm+ia-1)*6 + mu) + two*esum
319 : end do
320 : end do
321 : end do
322 :
323 0 : shift = shift + nattyp(itypat)*nlmn
324 0 : iatm = iatm+nattyp(itypat)
325 : end do
326 : end if
327 :
328 : ! ======= Accumulate the elastic tensor contributions ==========
329 520 : if (choice==6) then
330 : shift=0; iatm=0
331 0 : do itypat=1, ntypat
332 0 : nlmn=count(indlmn(3,:,itypat)>0)
333 0 : nattyp_i = nattyp(itypat)
334 :
335 : #ifdef HAVE_OPENMP_OFFLOAD
336 : !$OMP TARGET TEAMS DISTRIBUTE COLLAPSE(3) &
337 : !$OMP& MAP(to:fnlk,gxfac,dgxdt) &
338 : !$OMP& PRIVATE(idat,ia,esum,mu) &
339 : !$OMP& IF(gpu_option==ABI_GPU_OPENMP)
340 : #endif
341 0 : do idat=1,ndat
342 0 : do ia=1,nattyp_i
343 0 : do mu=1,3
344 : esum=zero
345 : !$OMP PARALLEL DO COLLAPSE(3) REDUCTION(+:esum) PRIVATE(ispinor,ilmn,ii)
346 0 : do ispinor=1,nspinor
347 0 : do ilmn=1,nlmn
348 0 : do ii=1,cplex
349 : esum=esum+gxfac(ii, shift+(ia-1)*nlmn+ilmn, ispinor+nspinor*(idat-1)) &
350 0 : & *dgxdt(ii, ndgxdt*shift+(ia-1)*nlmn*ndgxdt+(ilmn-1)*ndgxdt+mu+6, ispinor+nspinor*(idat-1))
351 : end do
352 : end do
353 : end do
354 0 : fnlk((iatm+ia-1)*3+mu,idat)=fnlk((iatm+ia-1)*3+mu,idat)+two*esum
355 : end do
356 : end do
357 : end do
358 :
359 0 : shift = shift + nattyp(itypat)*nlmn
360 0 : iatm = iatm+nattyp(itypat)
361 : end do
362 :
363 :
364 : shift=0; iatm=0
365 0 : do itypat=1, ntypat
366 0 : nlmn=count(indlmn(3,:,itypat)>0)
367 0 : nattyp_i = nattyp(itypat)
368 :
369 : #ifdef HAVE_OPENMP_OFFLOAD
370 : !$OMP TARGET TEAMS DISTRIBUTE COLLAPSE(2) &
371 : !$OMP& MAP(to:strnlk,gxfac,dgxdt) &
372 : !$OMP& PRIVATE(idat,esum,mu) &
373 : !$OMP& IF(gpu_option==ABI_GPU_OPENMP)
374 : #endif
375 0 : do idat=1,ndat
376 0 : do mu=1,6
377 : esum=zero
378 : !$OMP PARALLEL DO COLLAPSE(4) REDUCTION(+:esum) PRIVATE(ispinor,ilmn,ii,ia)
379 0 : do ispinor=1,nspinor
380 0 : do ia=1,nattyp_i
381 0 : do ilmn=1,nlmn
382 0 : do ii=1,cplex
383 : esum=esum+gxfac(ii, shift+(ia-1)*nlmn+ilmn, ispinor+nspinor*(idat-1)) &
384 0 : & *dgxdt(ii, ndgxdt*shift+(ia-1)*nlmn*ndgxdt+(ilmn-1)*ndgxdt+mu, ispinor+nspinor*(idat-1))
385 : end do
386 : end do
387 : end do
388 : end do
389 0 : strnlk(mu,idat)=strnlk(mu,idat)+two*esum
390 : end do
391 : end do
392 :
393 0 : shift = shift + nattyp(itypat)*nlmn
394 0 : iatm = iatm+nattyp(itypat)
395 : end do
396 :
397 : shift=0; iatm=0
398 0 : do itypat=1, ntypat
399 0 : nlmn=count(indlmn(3,:,itypat)>0)
400 0 : nattyp_i = nattyp(itypat)
401 :
402 : #ifdef HAVE_OPENMP_OFFLOAD
403 : !$OMP TARGET TEAMS DISTRIBUTE COLLAPSE(3) &
404 : !$OMP& MAP(to:enlout,gxfac,dgxdt,dgxdtfac,d2gxdt) &
405 : !$OMP& PRIVATE(idat,ia,esum,mu,mua,mub,nu,mushift,nushift) &
406 : !$OMP& IF(gpu_option==ABI_GPU_OPENMP)
407 : #endif
408 0 : do idat=1,ndat
409 0 : do mub=1,6
410 0 : do mua=1,6
411 0 : mushift=6*(mub-1);nushift=(3*natom+6)*(mub-1)
412 0 : mu=mushift+mua;nu=nushift+mua
413 0 : esum=zero
414 : !$OMP PARALLEL DO COLLAPSE(4) REDUCTION(+:esum) PRIVATE(ispinor,ia,ilmn,ii)
415 0 : do ispinor=1,nspinor
416 0 : do ia=1,nattyp_i
417 0 : do ilmn=1,nlmn
418 0 : do ii=1,cplex
419 : esum=esum+gxfac(ii,shift+(ia-1)*nlmn+ilmn,ispinor+nspinor*(idat-1))&
420 : & *d2gxdt(ii,nd2gxdt*shift+(ia-1)*nlmn*nd2gxdt+(ilmn-1)*nd2gxdt+mu,ispinor+nspinor*(idat-1))&
421 : & +dgxdtfac(ii,ndgxdtfac*shift+(ia-1)*nlmn*ndgxdtfac+(ilmn-1)*ndgxdtfac+mua,ispinor+nspinor*(idat-1))&
422 0 : & *dgxdt(ii,ndgxdt*shift+(ia-1)*nlmn*ndgxdt+(ilmn-1)*ndgxdt+mub,ispinor+nspinor*(idat-1))
423 : end do
424 : end do
425 : end do
426 : end do
427 0 : enlout((idat-1)*nnlout + nu)=enlout((idat-1)*nnlout + nu)+two*esum
428 : end do
429 : end do
430 : end do
431 :
432 0 : shift = shift + nattyp(itypat)*nlmn
433 0 : iatm = iatm+nattyp(itypat)
434 : end do
435 :
436 : shift=0; iatm=0
437 0 : do itypat=1, ntypat
438 0 : nlmn=count(indlmn(3,:,itypat)>0)
439 0 : nattyp_i = nattyp(itypat)
440 :
441 : #ifdef HAVE_OPENMP_OFFLOAD
442 : !$OMP TARGET TEAMS DISTRIBUTE COLLAPSE(4) &
443 : !$OMP& MAP(to:enlout,gxfac,dgxdt,dgxdtfac,d2gxdt) &
444 : !$OMP& PRIVATE(idat,ia,esum,mu,mua,mub,nu,mushift,nushift) &
445 : !$OMP& IF(gpu_option==ABI_GPU_OPENMP)
446 : #endif
447 0 : do idat=1,ndat
448 0 : do mub=1,6
449 0 : do ia=1,nattyp_i
450 0 : do mua=1,3
451 0 : mushift=36+3*(mub-1);nushift=6+(iatm+ia-1)*3+(3*natom+6)*(mub-1)
452 0 : mu=mushift+mua;nu=nushift+mua
453 0 : esum=zero
454 : !$OMP PARALLEL DO COLLAPSE(3) REDUCTION(+:esum) PRIVATE(ispinor,ilmn,ii)
455 0 : do ispinor=1,nspinor
456 0 : do ilmn=1,nlmn
457 0 : do ii=1,cplex
458 : esum=esum+(gxfac(ii,shift+(ia-1)*nlmn+ilmn,ispinor+nspinor*(idat-1))&
459 : & *d2gxdt(ii,nd2gxdt*shift+(ia-1)*nlmn*nd2gxdt+(ilmn-1)*nd2gxdt+mu,ispinor+nspinor*(idat-1))&
460 : & +dgxdtfac(ii,ndgxdtfac*shift+(ia-1)*nlmn*ndgxdtfac+(ilmn-1)*ndgxdtfac+mub,ispinor+nspinor*(idat-1))&
461 0 : & *dgxdt(ii,ndgxdt*shift+(ia-1)*nlmn*ndgxdt+(ilmn-1)*ndgxdt+mua+6,ispinor+nspinor*(idat-1)))
462 : end do
463 : end do
464 : end do
465 0 : enlout((idat-1)*nnlout + nu)=enlout((idat-1)*nnlout + nu)+two*esum
466 : end do
467 : end do
468 : end do
469 : end do
470 :
471 0 : shift = shift + nattyp(itypat)*nlmn
472 0 : iatm = iatm+nattyp(itypat)
473 : end do
474 : end if
475 :
476 : end if ! paw_opt < 3
477 :
478 520 : if(paw_opt==3) then
479 :
480 : ! ====== Accumulate contribution to <c|d2S/d_atm_pos d_left_k|c> =========
481 0 : if (choice==54) then
482 : shift=0
483 : iatm=0
484 0 : do itypat=1, ntypat
485 0 : nlmn=count(indlmn(3,:,itypat)>0)
486 0 : nattyp_i = nattyp(itypat)
487 :
488 : #ifdef HAVE_OPENMP_OFFLOAD
489 : !$OMP TARGET TEAMS DISTRIBUTE COLLAPSE(2) &
490 : !$OMP& MAP(to:enlout,gxfac_sij,dgxdt,dgxdtfac_sij,d2gxdt) &
491 : !$OMP& PRIVATE(idat,igrad,ia,esum,esumi,mu,nu,mua,mub,iashift) &
492 : !$OMP& IF(gpu_option==ABI_GPU_OPENMP)
493 : #endif
494 0 : do idat=1,ndat
495 0 : do ia=1,nattyp_i
496 0 : iashift=18*(ia+iatm-1)
497 0 : if(cplex==2) then
498 0 : do mua=1,3 ! atm. pos
499 0 : do mub=1,3 ! k
500 0 : mu=(mua-1)*3+mub
501 0 : nu=(mua-1)*6+mub*2-1
502 0 : esum=zero; esumi=zero
503 : !$OMP PARALLEL DO COLLAPSE(2) REDUCTION(+:esum,esumi) PRIVATE(ispinor,ilmn)
504 0 : do ispinor=1,nspinor
505 0 : do ilmn=1,nlmn
506 : esum=esum &
507 : & +dgxdt(1, ndgxdt*shift+(ia-1)*nlmn*ndgxdt+(ilmn-1)*ndgxdt+mua, ispinor+nspinor*(idat-1))&
508 : & *dgxdtfac_sij(1,ndgxdt*shift+(ia-1)*nlmn*ndgxdt+(ilmn-1)*ndgxdt+3+mub, ispinor+nspinor*(idat-1)) &
509 : & +dgxdt(2, ndgxdt*shift+(ia-1)*nlmn*ndgxdt+(ilmn-1)*ndgxdt+mua,ispinor+nspinor*(idat-1))&
510 : & *dgxdtfac_sij(2, ndgxdt*shift+(ia-1)*nlmn*ndgxdt+(ilmn-1)*ndgxdt+3+mub, ispinor+nspinor*(idat-1)) &
511 : & +gxfac_sij(1, shift+(ia-1)*nlmn+ilmn, ispinor+nspinor*(idat-1))&
512 : & *d2gxdt(1,nd2gxdt*shift+(ia-1)*nlmn*nd2gxdt+(ilmn-1)*nd2gxdt+mu,ispinor+nspinor*(idat-1)) &
513 : & +gxfac_sij(2, shift+(ia-1)*nlmn+ilmn, ispinor+nspinor*(idat-1))&
514 0 : & *d2gxdt(2,nd2gxdt*shift+(ia-1)*nlmn*nd2gxdt+(ilmn-1)*nd2gxdt+mu,ispinor+nspinor*(idat-1))
515 :
516 : esumi=esumi &
517 : & +dgxdt(1, ndgxdt*shift+(ia-1)*nlmn*ndgxdt+(ilmn-1)*ndgxdt+mua, ispinor+nspinor*(idat-1))&
518 : & *dgxdtfac_sij(2, ndgxdt*shift+(ia-1)*nlmn*ndgxdt+(ilmn-1)*ndgxdt+3+mub, ispinor+nspinor*(idat-1)) &
519 : & -dgxdt(2, ndgxdt*shift+(ia-1)*nlmn*ndgxdt+(ilmn-1)*ndgxdt+mua, ispinor+nspinor*(idat-1))&
520 : & *dgxdtfac_sij(1, ndgxdt*shift+(ia-1)*nlmn*ndgxdt+(ilmn-1)*ndgxdt+3+mub, ispinor+nspinor*(idat-1)) &
521 : & +gxfac_sij(1, shift+(ia-1)*nlmn+ilmn, ispinor+nspinor*(idat-1))&
522 : & *d2gxdt(2, nd2gxdt*shift+(ia-1)*nlmn*nd2gxdt+(ilmn-1)*nd2gxdt+mu, ispinor+nspinor*(idat-1)) &
523 : & -gxfac_sij(2, shift+(ia-1)*nlmn+ilmn, ispinor+nspinor*(idat-1))&
524 0 : & *d2gxdt(1, nd2gxdt*shift+(ia-1)*nlmn*nd2gxdt+(ilmn-1)*nd2gxdt+mu, ispinor+nspinor*(idat-1))
525 : end do
526 : end do
527 0 : enlout(nnlout*(idat-1)+iashift+nu)=enlout(nnlout*(idat-1)+iashift+nu)+esum
528 0 : enlout(nnlout*(idat-1)+iashift+nu+1)=enlout(nnlout*(idat-1)+iashift+nu+1)+esumi
529 : end do
530 : end do
531 : ! If cplex=1, dgxdt, d2gxdt and dgxdtfac_sij are real for atm. pos, pure imaginary for k
532 : else
533 0 : do mua=1,3 ! atm. pos
534 0 : do mub=1,3 ! k
535 0 : mu=(mua-1)*3+mub
536 0 : nu=(mua-1)*6+mub*2-1
537 0 : esumi=zero
538 : !$OMP PARALLEL DO COLLAPSE(2) REDUCTION(+:esumi) PRIVATE(ispinor,ilmn)
539 0 : do ispinor=1,nspinor
540 0 : do ilmn=1,nlmn
541 : esumi=esumi &
542 : & +dgxdt(1, ndgxdt*shift+(ia-1)*nlmn*ndgxdt+(ilmn-1)*ndgxdt+mua, ispinor+nspinor*(idat-1))&
543 : & *dgxdtfac_sij(1, ndgxdt*shift+(ia-1)*nlmn*ndgxdt+(ilmn-1)*ndgxdt+3+mub, ispinor+nspinor*(idat-1)) &
544 : & +gxfac_sij(1, shift+(ia-1)*nlmn+ilmn, ispinor+nspinor*(idat-1))&
545 0 : & *d2gxdt(1, nd2gxdt*shift+(ia-1)*nlmn*nd2gxdt+(ilmn-1)*nd2gxdt+mu, ispinor+nspinor*(idat-1))
546 : end do
547 : end do
548 0 : enlout(nnlout*(idat-1)+iashift+nu+1)=enlout(nnlout*(idat-1)+iashift+nu+1)+esumi
549 : end do
550 : end do
551 : end if
552 : end do ! ia
553 : end do ! idat
554 :
555 0 : shift = shift + nattyp(itypat)*nlmn
556 0 : iatm = iatm+nattyp(itypat)
557 : end do ! itypat
558 : end if
559 :
560 : ! ====== Accumulate contribution to <c|d2S/d_dstrain d_right_k|c> =========
561 0 : if (choice==55) then
562 0 : if(.not. (cplex==2.and.cplex_fac==2)) ABI_BUG("cplex==1 not supported")
563 : shift=0
564 : iatm=0
565 0 : do itypat=1, ntypat
566 0 : nlmn=count(indlmn(3,:,itypat)>0)
567 0 : nattyp_i = nattyp(itypat)
568 :
569 0 : if(cplex==2.and.cplex_fac==2) then
570 : #ifdef HAVE_OPENMP_OFFLOAD
571 : !$OMP TARGET TEAMS DISTRIBUTE &
572 : !$OMP& MAP(to:enlout,gxfac_sij,dgxdt,dgxdtfac_sij,d2gxdt) &
573 : !$OMP& PRIVATE(idat,igrad,ia,esum,esumi,mu,nu,mua,mub,mua1,mua2,muu,mut) &
574 : !$OMP& IF(gpu_option==ABI_GPU_OPENMP)
575 : #endif
576 0 : do idat=1,ndat
577 : ! If cplex=1, dgxdt is real for strain, pure imaginary for k;
578 : ! If cplex_fac=1, dgxdtfac is pure imaginary for k;
579 : ! First compute 2nd-derivative contribution
580 0 : do mua=1,6 ! strain (lambda,nu)
581 0 : do mub=1,3 ! k (mu)
582 0 : mu=(mua-1)*6+mub*2-1
583 0 : esum=zero; esumi=zero
584 0 : mua1=alpha(mua) ! (nu)
585 0 : mua2=beta(mua) ! (lambda)
586 0 : muu=3*(gamma(mua1,mub)-1)+mua2
587 0 : mut=3*(gamma(mua2,mub)-1)+mua1
588 : !$OMP PARALLEL DO COLLAPSE(3) REDUCTION(+:esum,esumi) PRIVATE(ispinor,ia,ilmn,d2gx)
589 0 : do ispinor=1,nspinor
590 0 : do ia=1,nattyp_i
591 0 : do ilmn=1,nlmn
592 : d2gx(1:cplex)=half*(&
593 : & d2gxdt(1:cplex, nd2gxdt*shift+(ia-1)*nlmn*nd2gxdt+(ilmn-1)*nd2gxdt+muu, ispinor+nspinor*(idat-1)) &
594 : & +d2gxdt(1:cplex, nd2gxdt*shift+(ia-1)*nlmn*nd2gxdt+(ilmn-1)*nd2gxdt+mut, ispinor+nspinor*(idat-1))&
595 0 : & )
596 :
597 : esum=esum &
598 : & +dgxdt(1, ndgxdt*shift+(ia-1)*nlmn*ndgxdt+(ilmn-1)*ndgxdt+mua, ispinor+nspinor*(idat-1))&
599 : & *dgxdtfac_sij(1, ndgxdt*shift+(ia-1)*nlmn*ndgxdt+(ilmn-1)*ndgxdt+6+mub, ispinor+nspinor*(idat-1)) &
600 : & +dgxdt(2, ndgxdt*shift+(ia-1)*nlmn*ndgxdt+(ilmn-1)*ndgxdt+mua, ispinor+nspinor*(idat-1))&
601 : & *dgxdtfac_sij(2, ndgxdt*shift+(ia-1)*nlmn*ndgxdt+(ilmn-1)*ndgxdt+6+mub, ispinor+nspinor*(idat-1)) &
602 : & +gxfac_sij(1, shift+(ia-1)*nlmn+ilmn, ispinor+nspinor*(idat-1))*d2gx(1)&
603 0 : & +gxfac_sij(2, shift+(ia-1)*nlmn+ilmn, ispinor+nspinor*(idat-1))*d2gx(2)
604 :
605 : esumi=esumi &
606 : & +dgxdt(1, ndgxdt*shift+(ia-1)*nlmn*ndgxdt+(ilmn-1)*ndgxdt+mua, ispinor+nspinor*(idat-1))&
607 : & *dgxdtfac_sij(2, ndgxdt*shift+(ia-1)*nlmn*ndgxdt+(ilmn-1)*ndgxdt+6+mub, ispinor+nspinor*(idat-1)) &
608 : & -dgxdt(2, ndgxdt*shift+(ia-1)*nlmn*ndgxdt+(ilmn-1)*ndgxdt+mua, ispinor+nspinor*(idat-1))&
609 : & *dgxdtfac_sij(1, ndgxdt*shift+(ia-1)*nlmn*ndgxdt+(ilmn-1)*ndgxdt+6+mub, ispinor+nspinor*(idat-1)) &
610 : & +gxfac_sij(1, shift+(ia-1)*nlmn+ilmn, ispinor+nspinor*(idat-1))*d2gx(2)&
611 0 : & -gxfac_sij(2, shift+(ia-1)*nlmn+ilmn, ispinor+nspinor*(idat-1))*d2gx(1)
612 : end do
613 : end do
614 : end do
615 0 : enlout(nnlout*(idat-1)+mu) = enlout(nnlout*(idat-1)+mu)+esum
616 0 : enlout(nnlout*(idat-1)+mu+1) = enlout(nnlout*(idat-1)+mu+1)+esumi
617 : end do
618 : end do
619 : ! Then store 1st-derivative contribution
620 0 : do nu=1,3
621 0 : mu=nu*2-1
622 0 : esum=zero; esumi=zero
623 : !$OMP PARALLEL DO COLLAPSE(3) REDUCTION(+:esum,esumi) PRIVATE(ispinor,ia,ilmn)
624 0 : do ispinor=1,nspinor
625 0 : do ia=1,nattyp_i
626 0 : do ilmn=1,nlmn
627 : esum=esum&
628 : & +gxfac_sij(1, shift+(ia-1)*nlmn+ilmn, ispinor+nspinor*(idat-1))&
629 : & *dgxdt(1, ndgxdt*shift+(ia-1)*nlmn*ndgxdt+(ilmn-1)*ndgxdt+6+nu, ispinor+nspinor*(idat-1)) &
630 : & +gxfac_sij(2, shift+(ia-1)*nlmn+ilmn, ispinor+nspinor*(idat-1))&
631 0 : & *dgxdt(2, ndgxdt*shift+(ia-1)*nlmn*ndgxdt+(ilmn-1)*ndgxdt+6+nu, ispinor+nspinor*(idat-1))
632 :
633 : esumi=esumi&
634 : & +gxfac_sij(1, shift+(ia-1)*nlmn+ilmn, ispinor+nspinor*(idat-1))&
635 : & *dgxdt(2, ndgxdt*shift+(ia-1)*nlmn*ndgxdt+(ilmn-1)*ndgxdt+6+nu, ispinor+nspinor*(idat-1)) &
636 : & -gxfac_sij(2, shift+(ia-1)*nlmn+ilmn, ispinor+nspinor*(idat-1))&
637 0 : & *dgxdt(1, ndgxdt*shift+(ia-1)*nlmn*ndgxdt+(ilmn-1)*ndgxdt+6+nu, ispinor+nspinor*(idat-1))
638 : end do
639 : end do
640 : end do
641 0 : ddkk(mu,idat) = ddkk(mu,idat)+esum
642 0 : ddkk(mu+1,idat)= ddkk(mu+1,idat)+esumi
643 : end do
644 : ! If cplex=1, dgxdt, d2gxdt and dgxdtfac_sij are real for atm. pos, pure imaginary for k
645 : ! ABI_BUG("Not implemented")
646 : ! do ilmn=1,nlmn
647 : ! mu=1
648 : ! do mua=1,6 ! strain (lambda,nu)
649 : ! mua1=alpha(mua) ! (nu)
650 : ! mua2=beta(mua) ! (lambda)
651 : ! do mub=1,3 ! k (mu)
652 : ! muu=3*(gamma(mua1,mub)-1)+mua2
653 : ! mut=3*(gamma(mua2,mub)-1)+mua1
654 : ! d2gx(1)=half*(d2gxdt(1,muu,ilmn,ia,idat)+d2gxdt(1,mut,ilmn,ia,idat))
655 : ! enljj(mu+1)=enljj(mu+1) &
656 : ! +dgxdt(1,mua,ilmn,ia,idat)*dgxdtfac_sij(1,6+mub,ilmn,ia,idat) &
657 : ! +gxfac_sij(1,ilmn,ia,idat)*d2gx(1)
658 : ! mu=mu+2
659 : ! end do
660 : ! end do
661 : ! ! Then store 1st-derivative contribution
662 : ! mu=1
663 : ! do nu=1,3
664 : ! enlj(mu+1)=enlj(mu+1)+gxfac_sij(1,ilmn,ia,idat)*dgxdt(1,6+nu,ilmn,ia,idat)
665 : ! mu=mu+2
666 : ! end do
667 : ! end do
668 : end do ! idat
669 :
670 0 : shift = shift + nattyp(itypat)*nlmn
671 0 : iatm = iatm+nattyp(itypat)
672 : end if
673 : end do ! itypat
674 : end if
675 :
676 : end if ! paw_opt == 3
677 :
678 520 : end subroutine opernld_ylm_allwf
679 : !!***
680 :
681 : end module m_opernld_ylm_allwf
682 : !!***
|
