Branch data 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 : : !!***
|
