Branch data Line data Source code
1 : : !!****m* ABINIT/m_paw_nhat
2 : : !! NAME
3 : : !! m_paw_nhat
4 : : !!
5 : : !! FUNCTION
6 : : !! This module contains several routines related to the PAW compensation
7 : : !! charge density (i.e. n^hat(r)).
8 : : !!
9 : : !! COPYRIGHT
10 : : !! Copyright (C) 2018-2026 ABINIT group (FJ, MT, MG, TRangel)
11 : : !! This file is distributed under the terms of the
12 : : !! GNU General Public License, see ~abinit/COPYING
13 : : !! or http://www.gnu.org/copyleft/gpl.txt .
14 : : !!
15 : : !! SOURCE
16 : :
17 : : #if defined HAVE_CONFIG_H
18 : : #include "config.h"
19 : : #endif
20 : :
21 : : #include "abi_common.h"
22 : :
23 : : MODULE m_paw_nhat
24 : :
25 : : use defs_basis
26 : : use m_abicore
27 : : use m_errors
28 : : use m_xmpi
29 : : use m_xomp
30 : : use m_gputk
31 : : use m_abi_linalg
32 : : use, intrinsic :: iso_c_binding, only: c_size_t,c_loc
33 : :
34 : : use defs_abitypes, only : MPI_type
35 : : use m_time, only : timab
36 : : use m_pawang, only : pawang_type
37 : : use m_pawtab, only : pawtab_type
38 : : use m_pawfgrtab, only : pawfgrtab_type
39 : : use m_pawrhoij, only : pawrhoij_type
40 : : use m_pawcprj, only : pawcprj_type
41 : : use m_paw_finegrid, only : pawgylm,pawrfgd_fft,pawrfgd_wvl,pawexpiqr
42 : : use m_paral_atom, only : get_my_atmtab, free_my_atmtab
43 : : use m_distribfft, only : distribfft_type
44 : : use m_geometry, only : xred2xcart
45 : : use m_cgtools, only : mean_fftr
46 : : use m_mpinfo, only : set_mpi_enreg_fft,unset_mpi_enreg_fft,initmpi_seq
47 : : use m_fft, only : zerosym, fourwf, fourdp
48 : : use m_paw_lmn, only : klmn2ijlmn
49 : :
50 : : implicit none
51 : :
52 : : private
53 : :
54 : : !public procedures.
55 : : public :: pawmknhat ! Compute compensation charge density on the real space (fine) grid
56 : : public :: pawmknhat_psipsi ! Compute compensation charge density associated to the product of two WF
57 : : public :: pawmknhat_psipsi_ndat ! Batched variant, assumes compact nhat12 storage
58 : : public :: pawnhatfr ! Compute frozen part of 1st-order compensation charge density nhat^(1) (DFPT)
59 : : public :: pawdijhat_ndat ! Compute compensation charge contribution
60 : : public :: pawsushat ! Compute contrib. to the product of two WF from compensation charge density
61 : : public :: nhatgrid ! Determine points of the (fine) grid that are located around atoms - PW version
62 : : public :: wvl_nhatgrid ! Determine points of the (fine) grid that are located around atoms - WVL version
63 : :
64 : : CONTAINS !========================================================================================
65 : : !!***
66 : :
67 : : !----------------------------------------------------------------------
68 : :
69 : : !!****f* m_paw_nhat/pawmknhat
70 : : !! NAME
71 : : !! pawmknhat
72 : : !!
73 : : !! FUNCTION
74 : : !! PAW only:
75 : : !! Compute compensation charge density (and derivatives) on the fine FFT grid
76 : : !! Can also compute first-order compensation charge density (RF calculations)
77 : : !!
78 : : !! INPUTS
79 : : !! cplex: if 1, real space 1-order functions on FFT grid are REAL, if 2, COMPLEX
80 : : !! distribfft<type(distribfft_type)>=--optional-- contains infos related to FFT parallelism
81 : : !! gprimd(3,3)=dimensional primitive translations for reciprocal space
82 : : !! ider= 0: nhat(r) is computed
83 : : !! 1: cartesian derivatives of nhat(r) are computed
84 : : !! 2: nhat(r) and derivatives are computed
85 : : !! idir=direction of atomic displacement (in case of atomic displ. perturb.)
86 : : !! ipert=index of perturbation; must be 0 for ground-state calculations
87 : : !! izero=if 1, unbalanced components of nhat(g) have to be set to zero
88 : : !! me_g0=--optional-- 1 if the current process treat the g=0 plane-wave (only needed when comm_fft is present)
89 : : !! mpi_atmtab(:)=--optional-- indexes of the atoms treated by current proc
90 : : !! comm_atom=--optional-- MPI communicator over atoms
91 : : !! comm_fft=--optional-- MPI communicator over FFT components
92 : : !! my_natom=number of atoms treated by current processor
93 : : !! natom=total number of atoms in cell
94 : : !! nfft=number of point on the rectangular fft grid
95 : : !! ngfft(18)=contain all needed information about 3D FFT, see ~abinit/doc/variables/vargs.htm#ngfft
96 : : !! nhatgrdim= -PAW only- 0 if pawgrnhat array is not used ; 1 otherwise
97 : : !! ntypat=number of types of atoms in unit cell.
98 : : !! paral_kgb=--optional-- 1 if "band-FFT" parallelism is activated (only needed when comm_fft is present)
99 : : !! pawang <type(pawang_type)>=paw angular mesh and related data
100 : : !! pawfgrtab(my_natom) <type(pawfgrtab_type)>=atomic data given on fine rectangular grid
101 : : !! pawrhoij(my_natom) <type(pawrhoij_type)>= paw rhoij occupancies and related data
102 : : !! (1st-order occupancies if ipert>0)
103 : : !! pawrhoij0(my_natom) <type(pawrhoij_type)>= GS paw rhoij occupancies and related data (used only if ipert>0)
104 : : !! set equat to pawrhoij for GS calculations
105 : : !! pawtab(ntypat) <type(pawtab_type)>=paw tabulated starting data
106 : : !! qphon(3)=wavevector of the phonon (RF only)
107 : : !! rprimd(3,3)=dimensional primitive translations for real space
108 : : !! ucvol=volume of the unit cell
109 : : !! xred(3,natom)= reduced atomic coordinates
110 : : !!
111 : : !! OUTPUT
112 : : !! === if ider=0 or 2
113 : : !! compch_fft=compensation charge inside spheres computed over fine fft grid
114 : : !! pawnhat(nfft,ispden)=nhat on fine rectangular grid
115 : : !! === if ider=1 or 2
116 : : !! pawgrnhat(nfft,ispden,3)=derivatives of nhat on fine rectangular grid (and derivatives)
117 : : !!
118 : : !! SOURCE
119 : :
120 : 33421 : subroutine pawmknhat(compch_fft,cplex,ider,idir,ipert,izero,gprimd,&
121 : 33421 : & my_natom,natom,nfft,ngfft,nhatgrdim,nspden,ntypat,pawang,pawfgrtab,&
122 : 66842 : & pawgrnhat,pawnhat,pawrhoij,pawrhoij0,pawtab,qphon,rprimd,ucvol,usewvl,xred,&
123 [ + + + + : 33421 : & mpi_atmtab,comm_atom,comm_fft,mpi_comm_wvl,me_g0,paral_kgb,distribfft,gpu_thread_limit) ! optional arguments
- + ]
124 : :
125 : : !Arguments ---------------------------------------------
126 : : !scalars
127 : : integer,intent(in) :: cplex,ider,idir,ipert,izero,my_natom,natom,nfft
128 : : integer,intent(in) :: usewvl
129 : : integer,intent(in) :: nhatgrdim,nspden,ntypat
130 : : integer,optional,intent(in) :: me_g0,comm_atom,comm_fft,mpi_comm_wvl,paral_kgb,gpu_thread_limit
131 : : real(dp),intent(in) :: ucvol
132 : : real(dp),intent(inout) :: compch_fft
133 : : type(distribfft_type),optional,intent(in),target :: distribfft
134 : : type(pawang_type),intent(in) :: pawang
135 : : !arrays
136 : : integer,intent(in) :: ngfft(18)
137 : : integer,optional,target,intent(in) :: mpi_atmtab(:)
138 : : real(dp),intent(in) :: gprimd(3,3),qphon(3),rprimd(3,3),xred(3,natom)
139 : : real(dp),intent(out) :: pawgrnhat(cplex*nfft,nspden,3*nhatgrdim)
140 : : real(dp),intent(inout) :: pawnhat(cplex*nfft,nspden) !vz_i
141 : : type(pawfgrtab_type),intent(inout) :: pawfgrtab(my_natom)
142 : : type(pawrhoij_type),intent(in) :: pawrhoij(my_natom),pawrhoij0(my_natom)
143 : : type(pawtab_type),intent(in) :: pawtab(ntypat)
144 : :
145 : : !Local variables ---------------------------------------
146 : : !scalars
147 : : integer :: cplex_rhoij,iatom,iatom_tot,ic,ierr,ii,ils,ilslm,iq0,irhoij,ispden,itypat
148 : : integer :: jc,jrhoij,kc,klm,klmn,lmax,lmin,lm_size,mfgd,mm,mpi_comm_sphgrid
149 : : integer :: my_comm_atom,my_comm_fft,nfgd,nfftot,option,optgr0,optgr1,optgr2,paral_kgb_fft
150 : : logical :: compute_grad,compute_nhat,my_atmtab_allocated,need_frozen,paral_atom,qeq0
151 : : logical :: compute_phonons,has_phase
152 : : type(distribfft_type),pointer :: my_distribfft
153 : 33421 : type(mpi_type) :: mpi_enreg_fft
154 : : !arrays
155 : 33421 : integer,pointer :: my_atmtab(:)
156 : 66842 : real(dp) :: ro(cplex),ro_ql(cplex),tmp_compch_fft(nspden),tsec(2)
157 : 33421 : real(dp),allocatable :: pawgrnhat_atm(:,:),pawnhat_atm(:),work(:,:)
158 : :
159 : : ! *************************************************************************
160 : :
161 : : DBG_ENTER("COLL")
162 : :
163 : 33421 : compute_nhat=(ider==0.or.ider==2)
164 : 33421 : compute_grad=(ider==1.or.ider==2)
165 [ + + ]: 6579 : compute_phonons=(ipert>0.and.ipert<=natom)
166 : :
167 : : !Compatibility tests
168 : 33421 : qeq0=(qphon(1)**2+qphon(2)**2+qphon(3)**2<1.d-15)
169 [ + + ]: 33421 : if (present(comm_fft)) then
170 [ - + ]: 32309 : if ((.not.present(paral_kgb)).or.(.not.present(me_g0))) then
171 : 0 : ABI_BUG('Need paral_kgb and me_g0 with comm_fft !')
172 : : end if
173 : : end if
174 [ + + - + ]: 33421 : if(ider>0.and.nhatgrdim==0) then
175 : 0 : ABI_BUG('Gradients of nhat required but not allocated!')
176 : : end if
177 [ + + ]: 33421 : if (my_natom>0) then
178 [ + + - + ]: 30149 : if(nspden>1.and.nspden/=pawrhoij(1)%nspden) then
179 : 0 : ABI_BUG('Wrong values for nspden and pawrhoij%nspden!')
180 : : end if
181 [ + + - + ]: 30149 : if(nspden>1.and.nspden/=pawfgrtab(1)%nspden) then
182 : 0 : ABI_BUG('Wrong values for nspden and pawfgrtab%nspden!')
183 : : end if
184 [ - + ]: 30149 : if(pawrhoij(1)%qphase<cplex) then
185 : 0 : ABI_BUG('Must have pawrhoij()%qphase >= cplex!')
186 : : end if
187 [ + + ]: 30149 : if (compute_phonons.and.(.not.qeq0)) then
188 [ - + ]: 611 : if (pawfgrtab(1)%rfgd_allocated==0) then
189 : 0 : ABI_BUG('pawfgrtab()%rfgd array must be allocated!')
190 : : end if
191 [ - + ]: 611 : if (compute_grad.and.(.not.compute_nhat)) then
192 : 0 : ABI_BUG('When q<>0, nhat gradients need nhat!')
193 : : end if
194 : : end if
195 : : end if
196 : :
197 : : !nhat1 does not have to be computed for ddk or d2dk
198 [ + - - + ]: 33421 : if (ipert==natom+1.or.ipert==natom+10) return
199 : :
200 : : !Set up parallelism over atoms
201 [ + + + + ]: 33421 : paral_atom=(present(comm_atom).and.(my_natom/=natom))
202 [ + + ]: 33421 : nullify(my_atmtab);if (present(mpi_atmtab)) my_atmtab => mpi_atmtab
203 [ + + ]: 33421 : my_comm_atom=xmpi_comm_self;if (present(comm_atom)) my_comm_atom=comm_atom
204 : : call get_my_atmtab(my_comm_atom,my_atmtab,my_atmtab_allocated,paral_atom,natom,&
205 : 33421 : & my_natom_ref=my_natom)
206 : :
207 : : !Initialisations
208 [ - + ]: 33421 : if ((.not.compute_nhat).and.(.not.compute_grad)) return
209 [ + + + + ]: 105734 : mfgd=zero;if (my_natom>0) mfgd=maxval(pawfgrtab(1:my_natom)%nfgd)
210 [ + + ]: 33421 : if (compute_nhat) then
211 [ + + + - : 99462 : ABI_MALLOC(pawnhat_atm,(cplex*mfgd))
+ - + + -
+ - + ]
212 [ + + + + ]: 1401353832 : pawnhat=zero
213 : : end if
214 [ + + ]: 33421 : if (compute_grad) then
215 [ - + + - : 2334 : ABI_MALLOC(pawgrnhat_atm,(cplex*mfgd,3))
+ - + - -
+ - + ]
216 [ + + + + : 87201547 : pawgrnhat=zero
+ + ]
217 : : end if
218 : :
219 : : !mpi communicators for spherical grid:
220 : 33421 : mpi_comm_sphgrid=xmpi_comm_self !no communicators passed
221 [ + + + - ]: 33421 : if(present(comm_fft) .and. usewvl==0) mpi_comm_sphgrid=comm_fft
222 [ + + - + ]: 33421 : if(present(mpi_comm_wvl) .and. usewvl==1) mpi_comm_sphgrid=mpi_comm_wvl
223 : :
224 : : !------------------------------------------------------------------------
225 : : !----- Loop over atoms
226 : : !------------------------------------------------------------------------
227 : :
228 [ + + ]: 105734 : do iatom=1,my_natom
229 [ + + ]: 72313 : iatom_tot=iatom;if (paral_atom) iatom_tot=my_atmtab(iatom)
230 : :
231 : 72313 : itypat=pawrhoij(iatom)%itypat
232 : 72313 : lm_size=pawfgrtab(iatom)%l_size**2
233 [ + + + + : 72313 : need_frozen=((compute_nhat).and.(ipert==iatom_tot.or.ipert==natom+3.or.ipert==natom+4))
+ + + + ]
234 : 72313 : nfgd=pawfgrtab(iatom)%nfgd
235 : 72313 : cplex_rhoij=pawrhoij(iatom)%cplex_rhoij
236 : 72313 : iq0=cplex_rhoij*pawrhoij(iatom)%lmn2_size
237 : :
238 : : ! Eventually compute g_l(r).Y_lm(r) factors for the current atom (if not already done)
239 : : if (((compute_nhat).and.(pawfgrtab(iatom)%gylm_allocated==0)).or.&
240 [ + + + + : 72313 : & ((compute_grad).and.(pawfgrtab(iatom)%gylmgr_allocated==0)).or.&
+ + + - +
+ - + ]
241 : : & ((compute_grad.and.need_frozen).and.(pawfgrtab(iatom)%gylmgr2_allocated==0))) then
242 : 18 : optgr0=0;optgr1=0;optgr2=0
243 [ + - + - ]: 18 : if ((compute_nhat).and.(pawfgrtab(iatom)%gylm_allocated==0)) then
244 [ + - ]: 18 : if (allocated(pawfgrtab(iatom)%gylm)) then
245 : 18 : ABI_FREE(pawfgrtab(iatom)%gylm)
246 : : end if
247 [ - + + - : 72 : ABI_MALLOC(pawfgrtab(iatom)%gylm,(nfgd,pawfgrtab(iatom)%l_size**2))
- + + - +
- + - - +
- + - + ]
248 : 18 : pawfgrtab(iatom)%gylm_allocated=2;optgr0=1
249 : : end if
250 [ + + + - ]: 18 : if ((compute_grad).and.(pawfgrtab(iatom)%gylmgr_allocated==0)) then
251 [ + - ]: 9 : if (allocated(pawfgrtab(iatom)%gylmgr)) then
252 : 9 : ABI_FREE(pawfgrtab(iatom)%gylmgr)
253 : : end if
254 [ - + + - : 36 : ABI_MALLOC(pawfgrtab(iatom)%gylmgr,(3,nfgd,pawfgrtab(iatom)%l_size**2))
- + + - +
- + - - +
- + - + ]
255 : 9 : pawfgrtab(iatom)%gylmgr_allocated=2;optgr1=1
256 : : end if
257 [ - + - - ]: 18 : if ((compute_grad.and.need_frozen).and.(pawfgrtab(iatom)%gylmgr2_allocated==0)) then
258 [ # # ]: 0 : if (allocated(pawfgrtab(iatom)%gylmgr2)) then
259 : 0 : ABI_FREE(pawfgrtab(iatom)%gylmgr2)
260 : : end if
261 [ # # # # : 0 : ABI_MALLOC(pawfgrtab(iatom)%gylmgr2,(6,nfgd,pawfgrtab(iatom)%l_size**2))
# # # # #
# # # # #
# # # # ]
262 : 0 : pawfgrtab(iatom)%gylmgr2_allocated=2;optgr2=1
263 : : end if
264 [ + - ]: 18 : if (optgr0+optgr1+optgr2>0) then
265 : : call pawgylm(pawfgrtab(iatom)%gylm,pawfgrtab(iatom)%gylmgr,pawfgrtab(iatom)%gylmgr2,&
266 : 18 : & lm_size,nfgd,optgr0,optgr1,optgr2,pawtab(itypat),pawfgrtab(iatom)%rfgd)
267 : : end if
268 : : end if
269 : :
270 : :
271 : : ! Eventually compute exp(-i.q.r) factors for the current atom (if not already done)
272 [ + + - + ]: 72313 : if (compute_phonons.and.(.not.qeq0).and.pawfgrtab(iatom)%expiqr_allocated==0) then
273 [ # # ]: 0 : if (allocated(pawfgrtab(iatom)%expiqr)) then
274 : 0 : ABI_FREE(pawfgrtab(iatom)%expiqr)
275 : : end if
276 [ # # # # : 0 : ABI_MALLOC(pawfgrtab(iatom)%expiqr,(2,nfgd))
# # # # #
# # # #
# ]
277 : : call pawexpiqr(pawfgrtab(iatom)%expiqr,gprimd,nfgd,qphon,&
278 : 0 : & pawfgrtab(iatom)%rfgd,xred(:,iatom_tot))
279 : 0 : pawfgrtab(iatom)%expiqr_allocated=2
280 : : end if
281 [ + + + + ]: 72313 : has_phase=(compute_phonons.and.pawfgrtab(iatom)%expiqr_allocated/=0)
282 : :
283 : : ! Eventually compute frozen part of nhat for the current atom (if not already done)
284 [ + + + - : 72313 : if ((need_frozen).and.((pawfgrtab(iatom)%nhatfr_allocated==0).or.&
+ + - + ]
285 : : & (compute_grad.and.pawfgrtab(iatom)%nhatfrgr_allocated==0))) then
286 [ # # ]: 0 : if (allocated(pawfgrtab(iatom)%nhatfr)) then
287 : 0 : ABI_FREE(pawfgrtab(iatom)%nhatfr)
288 : : end if
289 [ # # # # : 0 : ABI_MALLOC(pawfgrtab(iatom)%nhatfr,(nfgd,pawfgrtab(iatom)%nspden))
# # # # #
# # # # #
# # # # ]
290 : 0 : option=0;pawfgrtab(iatom)%nhatfr_allocated=2
291 [ # # ]: 0 : if (compute_grad) then
292 : 0 : option=1
293 [ # # ]: 0 : if (allocated(pawfgrtab(iatom)%nhatfrgr)) then
294 : 0 : ABI_FREE(pawfgrtab(iatom)%nhatfrgr)
295 : : end if
296 [ # # # # : 0 : ABI_MALLOC(pawfgrtab(iatom)%nhatfrgr,(3,nfgd,pawfgrtab(iatom)%nspden))
# # # # #
# # # # #
# # # # ]
297 : 0 : pawfgrtab(iatom)%nhatfrgr_allocated=2
298 : : end if
299 : : call pawnhatfr(option,idir,ipert,1,natom,nspden,ntypat,pawang,pawfgrtab(iatom),&
300 : 0 : & pawrhoij0(iatom),pawtab,rprimd)
301 : : end if
302 : :
303 : : ! ------------------------------------------------------------------------
304 : : ! ----- Loop over density components
305 : : ! ------------------------------------------------------------------------
306 : :
307 [ + + ]: 161077 : do ispden=1,nspden
308 : :
309 [ + + + + ]: 155432441 : if (compute_nhat) pawnhat_atm(1:cplex*nfgd)=zero
310 [ + + + + : 33724215 : if (compute_grad) pawgrnhat_atm(1:cplex*nfgd,1:3)=zero
+ + ]
311 : :
312 : : ! ------------------------------------------------------------------------
313 : : ! ----- Loop over ij channels (basis components)
314 : : ! ------------------------------------------------------------------------
315 : 88764 : jrhoij=1
316 [ + + ]: 3731195 : do irhoij=1,pawrhoij(iatom)%nrhoijsel
317 : 3642431 : klmn=pawrhoij(iatom)%rhoijselect(irhoij)
318 : 3642431 : klm =pawtab(itypat)%indklmn(1,klmn)
319 : 3642431 : lmin=pawtab(itypat)%indklmn(3,klmn)
320 : 3642431 : lmax=pawtab(itypat)%indklmn(4,klmn)
321 : :
322 : : ! Retrieve rhoij
323 [ + + ]: 3642431 : if (pawrhoij(iatom)%nspden/=2) then
324 : 2739525 : ro(1)=pawrhoij(iatom)%rhoijp(jrhoij,ispden)
325 [ + + ]: 2739525 : if (cplex==2) ro(2)=pawrhoij(iatom)%rhoijp(iq0+jrhoij,ispden)
326 : : else
327 [ + + ]: 902906 : if (ispden==1) then
328 : 451453 : ro(1)=pawrhoij(iatom)%rhoijp(jrhoij,1)+pawrhoij(iatom)%rhoijp(jrhoij,2)
329 [ + + ]: 451453 : if (cplex==2) ro(2)=pawrhoij(iatom)%rhoijp(iq0+jrhoij,1)+pawrhoij(iatom)%rhoijp(iq0+jrhoij,2)
330 [ + - ]: 451453 : else if (ispden==2) then
331 : 451453 : ro(1)=pawrhoij(iatom)%rhoijp(jrhoij,1)
332 [ + + ]: 451453 : if (cplex==2) ro(2)=pawrhoij(iatom)%rhoijp(iq0+jrhoij,1)
333 : : end if
334 : : end if
335 [ + + ]: 7319873 : ro(1:cplex)=pawtab(itypat)%dltij(klmn)*ro(1:cplex)
336 : :
337 [ + + ]: 3642431 : if (compute_nhat) then
338 [ + + ]: 3591731 : if (cplex==1) then
339 [ - + + + ]: 10552113 : do ils=lmin,lmax,2
340 [ + + ]: 39562846 : do mm=-ils,ils
341 : 29010733 : ilslm=ils*ils+ils+mm+1
342 [ + + ]: 36006126 : if (pawang%gntselect(ilslm,klm)>0) then
343 : 7074416 : ro_ql(1)=ro(1)*pawtab(itypat)%qijl(ilslm,klmn)
344 : : !$OMP PARALLEL DO PRIVATE(ic)
345 [ + + ]: 26738324291 : do ic=1,nfgd
346 : 26738324291 : pawnhat_atm(ic)=pawnhat_atm(ic)+ro_ql(1)*pawfgrtab(iatom)%gylm(ic,ilslm)
347 : : end do
348 : : end if
349 : : end do
350 : : end do
351 : : else
352 [ - + + + ]: 94625 : do ils=lmin,lmax,2
353 [ + + ]: 296175 : do mm=-ils,ils
354 : 201550 : ilslm=ils*ils+ils+mm+1
355 [ + + ]: 261164 : if (pawang%gntselect(ilslm,klm)>0) then
356 [ + + ]: 164316 : ro_ql(1:2)=ro(1:2)*pawtab(itypat)%qijl(ilslm,klmn)
357 : : !$OMP PARALLEL DO PRIVATE(ic,jc)
358 [ + + ]: 109296806 : do ic=1,nfgd
359 : 109242034 : jc=2*ic-1
360 [ + + ]: 327780874 : pawnhat_atm(jc:jc+1)=pawnhat_atm(jc:jc+1)+ro_ql(1:2)*pawfgrtab(iatom)%gylm(ic,ilslm)
361 : : end do
362 : : end if
363 : : end do
364 : : end do
365 : : end if
366 : : end if
367 : :
368 [ + + ]: 3642431 : if (compute_grad) then
369 [ + - ]: 119809 : if (cplex==1) then
370 [ - + + + ]: 383252 : do ils=lmin,lmax,2
371 [ + + ]: 1534923 : do mm=-ils,ils
372 : 1151671 : ilslm=ils*ils+ils+mm+1
373 [ + + ]: 1415114 : if (pawang%gntselect(ilslm,klm)>0) then
374 : 309093 : ro_ql(1)=ro(1)*pawtab(itypat)%qijl(ilslm,klmn)
375 [ + + ]: 1332987074 : do ic=1,nfgd
376 : 1332677981 : pawgrnhat_atm(ic,1)=pawgrnhat_atm(ic,1)+ro_ql(1)*pawfgrtab(iatom)%gylmgr(1,ic,ilslm)
377 : 1332677981 : pawgrnhat_atm(ic,2)=pawgrnhat_atm(ic,2)+ro_ql(1)*pawfgrtab(iatom)%gylmgr(2,ic,ilslm)
378 : 1332987074 : pawgrnhat_atm(ic,3)=pawgrnhat_atm(ic,3)+ro_ql(1)*pawfgrtab(iatom)%gylmgr(3,ic,ilslm)
379 : : end do
380 : : end if
381 : : end do
382 : : end do
383 : : else
384 [ # # # # ]: 0 : do ils=lmin,lmax,2
385 [ # # ]: 0 : do mm=-ils,ils
386 : 0 : ilslm=ils*ils+ils+mm+1
387 [ # # ]: 0 : if (pawang%gntselect(ilslm,klm)>0) then
388 [ # # ]: 0 : ro_ql(1:2)=ro(1:2)*pawtab(itypat)%qijl(ilslm,klmn)
389 [ # # ]: 0 : do ic=1,nfgd
390 : 0 : jc=2*ic-1
391 : : pawgrnhat_atm(jc:jc+1,1)=pawgrnhat_atm(jc:jc+1,1) &
392 [ # # ]: 0 : & +ro_ql(1:2)*pawfgrtab(iatom)%gylmgr(1,ic,ilslm)
393 : : pawgrnhat_atm(jc:jc+1,2)=pawgrnhat_atm(jc:jc+1,2) &
394 [ # # ]: 0 : & +ro_ql(1:2)*pawfgrtab(iatom)%gylmgr(2,ic,ilslm)
395 : : pawgrnhat_atm(jc:jc+1,3)=pawgrnhat_atm(jc:jc+1,3) &
396 [ # # ]: 0 : & +ro_ql(1:2)*pawfgrtab(iatom)%gylmgr(3,ic,ilslm)
397 : : end do
398 : : end if
399 : : end do
400 : : end do
401 : : end if
402 : : end if
403 : :
404 : : ! ------------------------------------------------------------------------
405 : : ! ----- End loop over ij channels
406 : : ! ------------------------------------------------------------------------
407 : 3731195 : jrhoij=jrhoij+cplex_rhoij
408 : : end do
409 : :
410 : : ! If RF calculation, add frozen part of 1st-order compensation density
411 [ + + ]: 88764 : if (need_frozen) then
412 [ + + ]: 6571 : if (cplex==1) then
413 [ + + ]: 5036419 : do ic=1,nfgd
414 : 5036419 : pawnhat_atm(ic)=pawnhat_atm(ic)+pawfgrtab(iatom)%nhatfr(ic,ispden)
415 : : end do
416 : : else
417 [ + + ]: 848260 : do ic=1,nfgd
418 : 847760 : jc=2*ic-1
419 : 848260 : pawnhat_atm(jc)=pawnhat_atm(jc)+pawfgrtab(iatom)%nhatfr(ic,ispden)
420 : : end do
421 : : end if
422 [ + + ]: 6571 : if (compute_grad) then
423 [ - + ]: 4 : if (cplex==1) then
424 [ # # ]: 0 : do ic=1,nfgd
425 : 0 : pawgrnhat_atm(ic,1)=pawgrnhat_atm(ic,1)+pawfgrtab(iatom)%nhatfrgr(1,ic,ispden)
426 : 0 : pawgrnhat_atm(ic,2)=pawgrnhat_atm(ic,2)+pawfgrtab(iatom)%nhatfrgr(2,ic,ispden)
427 : 0 : pawgrnhat_atm(ic,3)=pawgrnhat_atm(ic,3)+pawfgrtab(iatom)%nhatfrgr(3,ic,ispden)
428 : : end do
429 : : else
430 [ + + ]: 14784 : do ic=1,nfgd
431 : 14780 : jc=2*ic-1
432 : 14780 : pawgrnhat_atm(jc,1)=pawgrnhat_atm(jc,1)+pawfgrtab(iatom)%nhatfrgr(1,ic,ispden)
433 : 14780 : pawgrnhat_atm(jc,2)=pawgrnhat_atm(jc,2)+pawfgrtab(iatom)%nhatfrgr(2,ic,ispden)
434 : 14784 : pawgrnhat_atm(jc,3)=pawgrnhat_atm(jc,3)+pawfgrtab(iatom)%nhatfrgr(3,ic,ispden)
435 : : end do
436 : : end if
437 : : end if
438 : : end if
439 : :
440 : : ! If needed, multiply eventually by exp(-i.q.r) phase
441 [ + + ]: 88764 : if (has_phase) then
442 [ - + ]: 944 : if (cplex==1) then
443 [ # # ]: 0 : do ic=1,nfgd
444 : 0 : pawnhat_atm(ic)=pawnhat_atm(ic)*pawfgrtab(iatom)%expiqr(1,ic)
445 : : end do
446 : : else
447 [ + + ]: 1515196 : do ic=1,nfgd
448 : 1514252 : jc=2*ic-1
449 : 1514252 : ro_ql(1)= pawfgrtab(iatom)%expiqr(1,ic)
450 : 1514252 : ro_ql(2)=-pawfgrtab(iatom)%expiqr(2,ic)
451 [ + + ]: 4542756 : ro(1:2)=pawnhat_atm(jc:jc+1)
452 : 1514252 : pawnhat_atm(jc )=ro(1)*ro_ql(1)-ro(2)*ro_ql(2)
453 : 1515196 : pawnhat_atm(jc+1)=ro(2)*ro_ql(1)+ro(1)*ro_ql(2)
454 : : end do
455 : : end if
456 [ + + ]: 944 : if (compute_grad) then
457 [ - + ]: 4 : if (cplex==1) then
458 [ # # ]: 0 : do ic=1,nfgd
459 [ # # ]: 0 : pawgrnhat_atm(ic,1:3)=pawgrnhat_atm(ic,1:3)*pawfgrtab(iatom)%expiqr(1,ic)
460 : : end do
461 : : else
462 [ + + ]: 14784 : do ic=1,nfgd
463 : 14780 : jc=2*ic-1
464 : : ! dn^hat(r)/dr_i * exp(-i.q.r)
465 : 14780 : ro_ql(1)= pawfgrtab(iatom)%expiqr(1,ic)
466 : 14780 : ro_ql(2)=-pawfgrtab(iatom)%expiqr(2,ic)
467 [ + + ]: 59120 : do ii=1,3
468 [ + + ]: 133020 : ro(1:2)=pawgrnhat_atm(jc:jc+1,ii)
469 : 44340 : pawgrnhat_atm(jc ,ii)=ro(1)*ro_ql(1)-ro(2)*ro_ql(2)
470 : 59120 : pawgrnhat_atm(jc+1,ii)=ro(2)*ro_ql(1)+ro(1)*ro_ql(2)
471 : : end do
472 : : ! -i.q_i * [n^hat(r).exp(-i.q.r)]
473 [ + + ]: 44340 : ro(1:2)=pawnhat_atm(jc:jc+1)
474 [ + + ]: 59124 : do ii=1,3
475 : 44340 : pawgrnhat_atm(jc ,ii)=pawgrnhat_atm(jc ,ii)+qphon(ii)*ro(2)
476 : 59120 : pawgrnhat_atm(jc+1,ii)=pawgrnhat_atm(jc+1,ii)-qphon(ii)*ro(1)
477 : : end do
478 : : end do
479 : : end if
480 : : end if
481 : : end if
482 : :
483 : : ! Add the contribution of the atom to the compensation charge
484 : : ! LB-2025-12-11 : if the PAW sphere overlaps with itself (true if it is larger than the unit cell, rare case but possible...),
485 : : ! then several values of ic can give the same kc in the following loops, so the iterations are not independent,
486 : : ! and cannot be parallelized (for example with OpenMP directives)
487 [ + + ]: 88764 : if (compute_nhat) then
488 [ + + ]: 87700 : if (cplex==1) then
489 : : ! Not possible to parallelize here (see comment above)
490 [ + + ]: 152401929 : do ic=1,nfgd
491 : 152315173 : kc=pawfgrtab(iatom)%ifftsph(ic)
492 : 152401929 : pawnhat(kc,ispden)=pawnhat(kc,ispden)+pawnhat_atm(ic)
493 : : end do
494 : : else
495 : : ! Not possible to parallelize here (see comment above)
496 [ + + ]: 1515196 : do ic=1,nfgd
497 : 1514252 : jc=2*ic-1;kc=2*pawfgrtab(iatom)%ifftsph(ic)-1
498 [ + + ]: 4543700 : pawnhat(kc:kc+1,ispden)=pawnhat(kc:kc+1,ispden)+pawnhat_atm(jc:jc+1)
499 : : end do
500 : : end if
501 : : end if
502 [ + + ]: 161077 : if (compute_grad) then
503 [ + + ]: 2698 : if (cplex==1) then
504 : : ! Not possible to parallelize here (see comment above)
505 [ + + ]: 11182253 : do ic=1,nfgd
506 : 11179559 : kc=pawfgrtab(iatom)%ifftsph(ic)
507 [ + + ]: 44720930 : pawgrnhat(kc,ispden,1:3)=pawgrnhat(kc,ispden,1:3)+pawgrnhat_atm(ic,1:3)
508 : : end do
509 : : else
510 : : ! Not possible to parallelize here (see comment above)
511 [ + + ]: 14784 : do ic=1,nfgd
512 : 14780 : jc=2*ic-1;kc=2*pawfgrtab(iatom)%ifftsph(ic)-1
513 [ + + ]: 59124 : do ii=1,3
514 [ + + ]: 147800 : pawgrnhat(kc:kc+1,ispden,ii)=pawgrnhat(kc:kc+1,ispden,ii)+pawgrnhat_atm(jc:jc+1,ii)
515 : : end do
516 : : end do
517 : : end if
518 : : end if
519 : : ! ------------------------------------------------------------------------
520 : : ! ----- End loop over density components
521 : : ! ------------------------------------------------------------------------
522 : : end do
523 : :
524 [ + + ]: 72313 : if (pawfgrtab(iatom)%gylm_allocated==2) then
525 [ - + ]: 18 : ABI_FREE(pawfgrtab(iatom)%gylm)
526 [ - + ]: 18 : ABI_MALLOC(pawfgrtab(iatom)%gylm,(0,0))
527 : 18 : pawfgrtab(iatom)%gylm_allocated=0
528 : : end if
529 [ + + ]: 72313 : if (pawfgrtab(iatom)%gylmgr_allocated==2) then
530 [ - + ]: 9 : ABI_FREE(pawfgrtab(iatom)%gylmgr)
531 [ - + ]: 9 : ABI_MALLOC(pawfgrtab(iatom)%gylmgr,(0,0,0))
532 : 9 : pawfgrtab(iatom)%gylmgr_allocated=0
533 : : end if
534 [ - + ]: 72313 : if (pawfgrtab(iatom)%gylmgr2_allocated==2) then
535 [ # # ]: 0 : ABI_FREE(pawfgrtab(iatom)%gylmgr2)
536 [ # # ]: 0 : ABI_MALLOC(pawfgrtab(iatom)%gylmgr2,(0,0,0))
537 : 0 : pawfgrtab(iatom)%gylmgr2_allocated=0
538 : : end if
539 [ - + ]: 72313 : if (pawfgrtab(iatom)%nhatfr_allocated==2) then
540 [ # # ]: 0 : ABI_FREE(pawfgrtab(iatom)%nhatfr)
541 [ # # ]: 0 : ABI_MALLOC(pawfgrtab(iatom)%nhatfr,(0,0))
542 : 0 : pawfgrtab(iatom)%nhatfr_allocated=0
543 : : end if
544 [ - + ]: 72313 : if (pawfgrtab(iatom)%nhatfrgr_allocated==2) then
545 [ # # ]: 0 : ABI_FREE(pawfgrtab(iatom)%nhatfrgr)
546 [ # # ]: 0 : ABI_MALLOC(pawfgrtab(iatom)%nhatfrgr,(0,0,0))
547 : 0 : pawfgrtab(iatom)%nhatfrgr_allocated=0
548 : : end if
549 [ - + ]: 105734 : if (pawfgrtab(iatom)%expiqr_allocated==2) then
550 [ # # ]: 0 : ABI_FREE(pawfgrtab(iatom)%expiqr)
551 [ # # ]: 0 : ABI_MALLOC(pawfgrtab(iatom)%expiqr,(0,0))
552 : 0 : pawfgrtab(iatom)%expiqr_allocated=0
553 : : end if
554 : :
555 : : ! ------------------------------------------------------------------------
556 : : ! ----- End loop over atoms
557 : : ! ------------------------------------------------------------------------
558 : : end do
559 : :
560 : : !----- Free some memory
561 [ + + ]: 33421 : if (compute_nhat) then
562 [ - + ]: 33154 : ABI_FREE(pawnhat_atm)
563 : : end if
564 [ + + ]: 33421 : if (compute_grad) then
565 [ - + ]: 778 : ABI_FREE(pawgrnhat_atm)
566 : : end if
567 : :
568 : : !----- Reduction in case of parallelism
569 [ + + ]: 33421 : if (paral_atom) then
570 : 11200 : call timab(48,1,tsec)
571 [ + - ]: 11200 : if (compute_nhat) then
572 : 11200 : call xmpi_sum(pawnhat,my_comm_atom,ierr)
573 : : end if
574 [ + + ]: 11200 : if (compute_grad) then
575 : 8 : call xmpi_sum(pawgrnhat,my_comm_atom,ierr)
576 : : end if
577 : 11200 : call timab(48,2,tsec)
578 : : end if
579 : :
580 : : !----- Avoid unbalanced g-components numerical errors
581 [ - + - - : 33421 : if (izero==1.and.compute_nhat.and.usewvl==0) then
- - ]
582 : : ! Create fake mpi_enreg to wrap fourdp
583 [ # # ]: 0 : if (present(distribfft)) then
584 : : my_distribfft => distribfft
585 : : else
586 [ # # ]: 0 : ABI_MALLOC(my_distribfft,)
587 : 0 : call my_distribfft%init_seq('f',ngfft(2),ngfft(3),'fourdp')
588 : : end if
589 : 0 : call initmpi_seq(mpi_enreg_fft)
590 [ # # # # : 0 : ABI_FREE(mpi_enreg_fft%distribfft)
# # # # #
# # # # #
# # # # #
# # # # #
# # ]
591 [ # # ]: 0 : if (present(comm_fft)) then
592 : 0 : call set_mpi_enreg_fft(mpi_enreg_fft,comm_fft,my_distribfft,me_g0,paral_kgb)
593 : 0 : my_comm_fft=comm_fft;paral_kgb_fft=paral_kgb
594 : : else
595 : 0 : my_comm_fft=xmpi_comm_self;paral_kgb_fft=0;
596 : 0 : mpi_enreg_fft%distribfft => my_distribfft
597 : : end if
598 : : ! do FFT
599 [ # # # # : 0 : ABI_MALLOC(work,(2,nfft))
# # # # #
# # # ]
600 [ # # ]: 0 : do ispden=1,min(2,nspden)
601 : 0 : call fourdp(cplex,work,pawnhat(:,ispden),-1,mpi_enreg_fft,nfft,1,ngfft,0)
602 : 0 : call zerosym(work,2,ngfft(1),ngfft(2),ngfft(3),comm_fft=my_comm_fft,distribfft=my_distribfft)
603 : 0 : call fourdp(cplex,work,pawnhat(:,ispden),+1,mpi_enreg_fft,nfft,1,ngfft,0)
604 : : end do
605 [ # # ]: 0 : ABI_FREE(work)
606 : : ! Destroy fake mpi_enreg
607 : 0 : call unset_mpi_enreg_fft(mpi_enreg_fft)
608 [ # # ]: 0 : if (.not.present(distribfft)) then
609 : 0 : call my_distribfft%free()
610 [ # # # # : 0 : ABI_FREE(my_distribfft)
# # # # #
# # # # #
# # # # #
# # # #
# ]
611 : : end if
612 : : end if
613 : :
614 : : !----- Computation of compensation charge over real space grid
615 [ + + + + ]: 33421 : if (compute_nhat.and.ipert==0) then
616 [ + + ]: 106300 : nfftot=PRODUCT(ngfft(1:3))
617 : : call mean_fftr(pawnhat,tmp_compch_fft,nfft,nfftot,1,&
618 : 26575 : & mpi_comm_sphgrid=mpi_comm_sphgrid,gpu_thread_limit=gpu_thread_limit)
619 : 26575 : compch_fft = tmp_compch_fft(1)
620 : 26575 : compch_fft=compch_fft*ucvol
621 : : end if
622 : :
623 : : !Destroy atom table used for parallelism
624 : 33421 : call free_my_atmtab(my_atmtab,my_atmtab_allocated)
625 : :
626 : : DBG_EXIT("COLL")
627 : :
628 [ + + - + : 66842 : end subroutine pawmknhat
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + ]
629 : : !!***
630 : :
631 : : !----------------------------------------------------------------------
632 : :
633 : : !!****f* m_paw_nhat/pawmknhat_psipsi_ndat
634 : : !! NAME
635 : : !! pawmknhat_psipsi_ndat
636 : : !!
637 : : !! FUNCTION
638 : : !! PAW only:
639 : : !! Compute on the fine FFT grid the compensation charge density (and derivatives) associated
640 : : !! to the product of two wavefunctions n_{12}(r) = \Psi_1* \Psi_2.
641 : : !! Based on pawmknhat_psipsi, but:
642 : : !! - support batching
643 : : !! - execution on GPU
644 : : !!
645 : : !! INPUTS
646 : : !! cprj1(natom,nspinor), cprj2(natom,nspinor) <type(pawcprj_type)>=
647 : : !! projected input wave functions <Proj_i|Cnk> with all NL projectors corresponding to
648 : : !! the \Psi_1 and \Psi_2, respectively.
649 : : !! distribfft<type(distribfft_type)>=--optional-- contains infos related to FFT parallelism
650 : : !! ider= 0: nhat(r) is computed
651 : : !! 1: cartesian derivatives of nhat(r) are computed
652 : : !! 2: nhat(r) and derivatives are computed
653 : : !! 3: nhat(r) and gradients of nhat wrt atomic coordinates are computed
654 : : !! Note: ider>0 not compatible with ipert>0
655 : : !! izero=if 1, unbalanced components of nhat(g) have to be set to zero
656 : : !! mpi_atmtab(:)=--optional-- indexes of the atoms treated by current proc
657 : : !! comm_atom=--optional-- MPI communicator over atoms
658 : : !! my_natom=number of atoms treated by current processor
659 : : !! natom=total number of atoms in cell
660 : : !! nfft=number of point on the rectangular fft grid
661 : : !! ngfft(18)=contain all needed information about 3D FFT, see ~abinit/doc/variables/vargs.htm#ngfft
662 : : !! nhat12_grdim= 0 if grnhat12 array is not used ; 1 otherwise
663 : : !! ntypat=number of types of atoms in unit cell.
664 : : !! pawang <type(pawang_type)>=paw angular mesh and related data
665 : : !! pawfgrtab(my_natom) <type(pawfgrtab_type)>=atomic data given on fine rectangular grid
666 : : !! pawtab(ntypat) <type(pawtab_type)>=paw tabulated starting data
667 : : !! ndat1=first batching size (from number of band)
668 : : !! ndat2=second batching size (from number of occupated states)
669 : : !! gpu_option=if equal to ABI_GPU_OPENMP, run on GPU
670 : : !!
671 : : !! OUTPUT
672 : : !! === if ider=0 or 2
673 : : !! nhat12(2,nfgd_max,nspinor**2,ndat2,ndat1,natom)=nhat on fine rectangular grid*exp(iqr),
674 : : !! === if ider=1 or 2 (not implemented)
675 : : !! grnhat12(nfft,nspinor**2,3)=gradient of (nhat*exp(iqr)) on fine rectangular grid (derivative versus r)
676 : : !! === if ider=3
677 : : !! grnhat_12(2,nfgd_max,nspinor**2,3,natom*(ider/3),ndat2,ndat1)=derivatives of nhat on fine rectangular grid versus R*exp(iqr).
678 : : !!
679 : : !! NOTES
680 : : !! izero=1 is not supported by this routine: it would require a full-grid FFT smoothing of nhat12,
681 : : !! which is incompatible with the compact (nfgd_max-sized) storage used here.
682 : : !!
683 : : !! nhat12 and grnhat_12 are expected to be sized after each atom's PAW augmentation sphere (nfgd_max points)
684 : : !! and indexed by the local in-sphere point index.
685 : : !! (see pawfgrtab(iatom)%ifftsph for the mapping to the full FFT grid).
686 : : !!
687 : : !! Only cases with ider in {0,3} are supported, grnhat12 isn't touched but kept for sticking with pawmknhat_psipsi prototype
688 : : !!
689 : : !!
690 : : !! SOURCE
691 : :
692 : 19356 : subroutine pawmknhat_psipsi_ndat(cprj1,cprj2,ider,izero,my_natom,natom,nfft,ngfft,nhat12_grdim,&
693 [ - + - + ]: 19356 : & nspinor,ntypat,ndat1,ndat2,pawang,pawfgrtab,grnhat12,nhat12,nattyp,pawtab, &
694 [ - + - - : 19356 : & gprimd,grnhat_12,qphon,xred,atindx,mpi_atmtab,comm_atom,gpu_option) ! optional arguments
- - + - +
- - + ]
695 : :
696 : : !Arguments ---------------------------------------------
697 : : !scalars
698 : : integer,intent(in) :: ider,izero,my_natom,natom,nfft,nhat12_grdim,ntypat,nspinor,ndat1,ndat2
699 : : integer,optional,intent(in) :: gpu_option
700 : : integer,optional,intent(in) :: comm_atom
701 : : type(pawang_type),intent(in),target :: pawang
702 : : !arrays
703 : : integer,intent(in) :: ngfft(18),nattyp(ntypat)
704 : : integer,optional,intent(in) ::atindx(natom)
705 : : integer,optional,target,intent(in) :: mpi_atmtab(:)
706 : : real(dp),optional, intent(in) ::gprimd(3,3),qphon(3),xred(3,natom)
707 : : real(dp),intent(out) :: grnhat12(:,:,:,:,:,:)
708 : : real(dp),optional,target,intent(out) :: grnhat_12(:,:,:,:,:,:,:)
709 : : real(dp),target,intent(out) :: nhat12(:,:,:,:,:,:)
710 : : type(pawfgrtab_type),intent(inout),target :: pawfgrtab(my_natom)
711 : : type(pawtab_type),intent(in),target :: pawtab(ntypat)
712 : : type(pawcprj_type),intent(in) :: cprj1(natom,nspinor*ndat1),cprj2(natom,nspinor*ndat2)
713 : :
714 : : !Local variables ---------------------------------------
715 : : !scalars
716 : : complex(dp), parameter :: cminusone = (-1._dp,0._dp)
717 : : integer :: iatm,iatom,iatom_tot,ic,ierr,ils,ilslm,isp1,isp2,isploop,itypat,klm,klmn,idat1,idat2,ia,nfgd_max
718 : : integer :: lmax,lmin,lm_size,mm,my_comm_atom,optgr0,optgr1
719 : : integer :: ilmn,jlmn,lmn_size,lmn2_size,gpu_option_,nprojs,shift,nlmn,nfgd
720 : : logical :: compute_grad,compute_grad1,compute_nhat,my_atmtab_allocated,paral_atom,qeq0,compute_phonon,order
721 : : real(dp) :: wgt,wgt1,wgt2,wgt3
722 : : integer :: gemm_n
723 : : #ifdef HAVE_OPENMP_OFFLOAD
724 : : integer :: gemm_batch
725 : : #endif
726 : : !arrays
727 : : integer,parameter :: spinor_idxs(2,4)=RESHAPE((/1,1,2,2,1,2,2,1/),(/2,4/))
728 : 9678 : integer,pointer :: my_atmtab(:)
729 : : real(dp) :: rdum(1),tsec(2),ro(2),ro_ql(2)
730 : 9678 : real(dp),allocatable :: qijl(:,:),projs1(:,:,:),projs2(:,:,:),gnt_scal(:,:)
731 : 9678 : real(dp),allocatable, target :: cpf_re(:,:,:,:),cpf_im(:,:,:,:)
732 : 9678 : real(dp),allocatable, target :: gemm_re(:,:,:,:),gemm_im(:,:,:,:)
733 : 9678 : real(dp),allocatable, target :: gemm_gr_re(:,:,:,:,:),gemm_gr_im(:,:,:,:,:)
734 : 9678 : real(dp),allocatable, target :: atom_wgylm(:,:,:),atom_wgylmgr(:,:,:,:)
735 : 9678 : real(dp), ABI_CONTIGUOUS pointer :: nhat12_atm(:,:,:,:,:,:)
736 : 9678 : real(dp), ABI_CONTIGUOUS pointer :: atom_expiqr(:,:,:),atom_gylm(:,:,:),atom_dltij(:),atom_gylmgr(:,:,:,:)
737 : 9678 : integer, ABI_CONTIGUOUS pointer :: atom_nfgd(:),atom_indklmn(:,:)
738 : :
739 : : ! *************************************************************************
740 : :
741 : : DBG_ENTER("COLL")
742 : :
743 : : ABI_UNUSED(nfft)
744 : : ABI_UNUSED(ngfft)
745 : : !Compatibility tests
746 : : if (ider>0.and.nhat12_grdim==0) then
747 : : ! ABI_BUG('Gradients of nhat required but not allocated !')
748 : : end if
749 : 9678 : if (nspinor==2) then
750 : 0 : ABI_BUG('nspinor==2 not coded!')
751 : : end if
752 [ + - ]: 9678 : gpu_option_=ABI_GPU_DISABLED; if (present(gpu_option)) gpu_option_=gpu_option
753 [ + - ]: 9678 : if(gpu_option_/=ABI_GPU_OPENMP) gpu_option_=ABI_GPU_DISABLED ! Only OpenMP variant supported
754 [ - + ]: 9678 : if(ider==1 .or. ider==2) then
755 : 0 : ABI_BUG('ider=={1,2} not coded with GPU!')
756 : : end if
757 [ - + - - : 9678 : if (izero==1.and.(ider==0.or.ider==2.or.ider==3)) then
- - ]
758 : : ! nhat12 is stored compactly (nfgd_max-sized, per atom): the full-grid FFT/zerosym
759 : : ! smoothing of unbalanced g-components cannot be performed on this compact storage.
760 : 0 : ABI_BUG('izero=1 not supported by pawmknhat_psipsi_ndat (compact nhat12 storage)!')
761 : : end if
762 : :
763 : 9678 : compute_phonon=.false.;qeq0=.false.
764 [ + - + - ]: 9678 : if (present(gprimd).and.present(qphon).and.present(xred)) compute_phonon=.true.
765 : 9678 : if (compute_phonon) qeq0=(qphon(1)**2+qphon(2)**2+qphon(3)**2<1.d-15)
766 : 9678 : if (present(atindx)) order=.true.
767 : : !Set up parallelism over atoms
768 [ - + - - ]: 9678 : paral_atom=(present(comm_atom).and.(my_natom/=natom))
769 [ - + ]: 9678 : nullify(my_atmtab);if (present(mpi_atmtab)) my_atmtab => mpi_atmtab
770 [ - + ]: 9678 : my_comm_atom=xmpi_comm_self;if (present(comm_atom)) my_comm_atom=comm_atom
771 : 9678 : call get_my_atmtab(my_comm_atom,my_atmtab,my_atmtab_allocated,paral_atom,natom,my_natom_ref=my_natom)
772 : :
773 : : !Initialisations
774 [ + + + - ]: 9678 : compute_nhat=(ider==0.or.ider==2.or.ider==3)
775 : 9678 : compute_grad=(ider==1.or.ider==2)
776 : 9678 : compute_grad1=(ider==3)
777 [ + - ]: 9678 : if ((.not.compute_nhat).and.(.not.compute_grad)) return
778 : :
779 [ - + - - : 9678 : if (compute_grad) grnhat12=zero
- - - - -
- - - -
- ]
780 [ + + ]: 9678 : if (compute_grad1) then
781 : : select case(gpu_option_)
782 : : case (ABI_GPU_DISABLED)
783 [ + + ]: 5764 : do idat1=1,ndat1
784 [ + + + + : 200627845 : grnhat_12(:,:,:,:,:,:,idat1) = zero
+ + + + +
+ + + ]
785 : : end do
786 : : case (ABI_GPU_OPENMP)
787 : : !FIXME grnhat_12 assumed to be mapped on GPU
788 [ # # ]: 0 : do idat1=1,ndat1
789 [ # # # # : 0 : call gpu_set_to_zero(grnhat_12(:,:,:,:,:,:,idat1),int(2,c_size_t)*size(grnhat_12,2)*nspinor**2*3*natom*ndat2)
# # # # #
# # # # #
# # # # #
# ]
790 : : end do
791 : : case default
792 [ + - - ]: 2882 : ABI_BUG("Unsupported GPU option")
793 : : end select
794 : : end if
795 : :
796 [ - + + - : 38712 : ABI_MALLOC(gnt_scal,(size(pawang%gntselect,1),size(pawang%gntselect,2)))
- + + - +
- + - - +
- + ]
797 [ + + + + ]: 977478 : gnt_scal=0
798 [ + + ]: 106458 : do klm=1,size(pawang%gntselect,2)
799 [ + + ]: 977478 : do ilslm=1,size(pawang%gntselect,1)
800 [ + + ]: 967800 : if(pawang%gntselect(ilslm,klm)>0) gnt_scal(ilslm,klm)=1
801 : : end do
802 : : end do
803 : :
804 : : nprojs=0
805 [ + + ]: 29034 : do iatom = 1,my_natom
806 : 29034 : nprojs = nprojs + cprj1(iatom, 1)%nlmn
807 : : end do
808 [ - + + - : 38712 : ABI_MALLOC(projs1,(2,nprojs,nspinor*ndat1))
- + + - +
- + - - +
- + ]
809 [ - + + - : 38712 : ABI_MALLOC(projs2,(2,nprojs,nspinor*ndat2))
- + + - +
- + - - +
- + ]
810 : : !$OMP PARALLEL DO PRIVATE(shift,idat2,iatom,nlmn)
811 [ + + ]: 19356 : do idat1=1, ndat1*nspinor
812 : : shift = 0
813 [ + + ]: 38712 : do iatom = 1,my_natom
814 : 19356 : nlmn = cprj1(iatom, idat1)%nlmn
815 [ + + + + ]: 471048 : projs1(:, shift+1:shift+nlmn, idat1) = cprj1(iatom, idat1)%cp(:, 1:nlmn)
816 : 29034 : shift = shift + nlmn
817 : : end do
818 : : end do
819 : : !$OMP PARALLEL DO PRIVATE(shift,idat2,iatom,nlmn)
820 [ + + ]: 48051 : do idat2=1, ndat2*nspinor
821 : : shift = 0
822 [ + + ]: 124797 : do iatom = 1,my_natom
823 : 76746 : nlmn = cprj2(iatom, idat2)%nlmn
824 [ + + + + ]: 1878384 : projs2(:, shift+1:shift+nlmn, idat2) = cprj2(iatom, idat2)%cp(:, 1:nlmn)
825 : 115119 : shift = shift + nlmn
826 : : end do
827 : : end do
828 : : #ifdef HAVE_OPENMP_OFFLOAD
829 : : !$OMP TARGET ENTER DATA MAP(to:projs1,projs2,nattyp) IF(gpu_option_==ABI_GPU_OPENMP)
830 : : #endif
831 : : !------------------------------------------------------------------------
832 : : !----- Loop over atoms types
833 : : !------------------------------------------------------------------------
834 : : shift = 0; iatm=0
835 [ + + ]: 20070 : do itypat=1,ntypat
836 : 10392 : atom_dltij => pawtab(itypat)%dltij
837 : 10392 : atom_indklmn => pawtab(itypat)%indklmn
838 : 10392 : lm_size = pawtab(itypat)%l_size**2
839 : 10392 : lmn_size = pawtab(itypat)%lmn_size
840 : 10392 : lmn2_size = pawtab(itypat)%lmn2_size
841 [ - + + - : 41568 : ABI_MALLOC(qijl,(lm_size,lmn2_size))
- + + - +
- + - - +
- + ]
842 [ + + + + ]: 3498756 : qijl=zero
843 [ + - + - : 3509148 : qijl=pawtab(itypat)%qijl
- + - - -
- - - + +
+ + ]
844 : 10392 : nlmn = cprj1(iatm+1, 1)%nlmn
845 : :
846 : 10392 : nhat12_atm => nhat12(:,:,:,:,:,iatm+1:iatm+nattyp(itypat))
847 : :
848 [ + - ]: 10392 : if (compute_nhat) then
849 [ + - ]: 10392 : if(gpu_option_==ABI_GPU_DISABLED) then
850 [ + + ]: 29748 : do ia=1,nattyp(itypat)
851 [ + + + + : 192121446 : nhat12_atm(:,:,:,:,:,ia) = zero
+ + + + +
+ ]
852 : : end do
853 [ # # ]: 0 : else if(gpu_option_==ABI_GPU_OPENMP) then
854 [ # # ]: 0 : do ia=1,nattyp(itypat)
855 : 0 : call gpu_set_to_zero(nhat12_atm(:,:,:,:,:,ia),int(2,c_size_t)*size(nhat12_atm,2)*(nspinor**2)*ndat2*ndat1)
856 : : end do
857 : : end if
858 : : end if
859 : :
860 : : ! Separate real and imaginary parts to perform a batched GEMM later.
861 [ - + + - : 62352 : ABI_MALLOC(cpf_re,(lmn2_size,ndat2,ndat1,nattyp(itypat)))
- + + - -
+ + - - +
+ - + - +
- - + -
+ ]
862 [ - + + - : 51960 : ABI_MALLOC(cpf_im,(lmn2_size,ndat2,ndat1,nattyp(itypat)))
- + + - -
+ + - - +
+ - - + -
+ ]
863 : : #ifdef HAVE_OPENMP_OFFLOAD
864 : : !$OMP TARGET ENTER DATA MAP(alloc:cpf_re,cpf_im) IF(gpu_option_==ABI_GPU_OPENMP)
865 : : #endif
866 : :
867 : : !------------------------------------------------------------------------
868 : : !----- Loop over atoms (init)
869 : : !------------------------------------------------------------------------
870 : 10392 : nfgd_max=1
871 [ + + ]: 29748 : do ia=1,nattyp(itypat)
872 : 19356 : iatom=iatm+ia
873 : 19356 : iatom_tot=iatom
874 : :
875 : 19356 : nfgd_max = MAX(pawfgrtab(iatom)%nfgd,nfgd_max)
876 : :
877 : : ! Eventually compute g_l(r).Y_lm(r) factors for the current atom (if not already done)
878 [ + - - + : 19356 : if (((compute_nhat).and.(pawfgrtab(iatom)%gylm_allocated==0)).or.&
- - - - ]
879 : : & (((compute_grad).or.(compute_grad1)).and.(pawfgrtab(iatom)%gylmgr_allocated==0))) then
880 : 19356 : optgr0=0; optgr1=0
881 [ + - + - ]: 19356 : if ((compute_nhat).and.(pawfgrtab(iatom)%gylm_allocated==0)) then
882 [ + - ]: 19356 : if (allocated(pawfgrtab(iatom)%gylm)) then
883 : 19356 : ABI_FREE(pawfgrtab(iatom)%gylm)
884 : : end if
885 [ - + + - : 77424 : ABI_MALLOC(pawfgrtab(iatom)%gylm,(pawfgrtab(iatom)%nfgd,pawfgrtab(iatom)%l_size**2))
- + + - +
- + - - +
- + - + ]
886 : 19356 : pawfgrtab(iatom)%gylm_allocated=2;optgr0=1
887 : : end if
888 [ + + + - ]: 19356 : if (((compute_grad).or.(compute_grad1)).and.(pawfgrtab(iatom)%gylmgr_allocated==0)) then
889 [ + - ]: 5764 : if (allocated(pawfgrtab(iatom)%gylmgr)) then
890 : 5764 : ABI_FREE(pawfgrtab(iatom)%gylmgr)
891 : : end if
892 [ - + + - : 23056 : ABI_MALLOC(pawfgrtab(iatom)%gylmgr,(3,pawfgrtab(iatom)%nfgd,pawfgrtab(iatom)%l_size**2))
- + + - +
- + - - +
- + - + ]
893 : 5764 : pawfgrtab(iatom)%gylmgr_allocated=2;optgr1=1
894 : : end if
895 [ + - ]: 19356 : if (optgr0+optgr1>0) then
896 : : call pawgylm(pawfgrtab(iatom)%gylm,pawfgrtab(iatom)%gylmgr,rdum,&
897 : : & lm_size,pawfgrtab(iatom)%nfgd,optgr0,optgr1,0,pawtab(itypat),&
898 : 19356 : & pawfgrtab(iatom)%rfgd)
899 : : end if
900 : : end if
901 [ + + + - ]: 29748 : if (compute_phonon.and.(.not.qeq0).and.(pawfgrtab(iatom)%expiqr_allocated==0)) then
902 [ + - ]: 15642 : if (allocated(pawfgrtab(iatom)%expiqr)) then
903 : 15642 : ABI_FREE(pawfgrtab(iatom)%expiqr)
904 : : end if
905 [ - + + - : 46926 : ABI_MALLOC(pawfgrtab(iatom)%expiqr,(2,pawfgrtab(iatom)%nfgd))
+ - + - -
+ - + -
+ ]
906 : : call pawexpiqr(pawfgrtab(iatom)%expiqr,gprimd,pawfgrtab(iatom)%nfgd,qphon,&
907 [ + - ]: 15642 : & pawfgrtab(iatom)%rfgd,xred(:,iatom_tot))
908 : 15642 : pawfgrtab(iatom)%expiqr_allocated=2
909 : : end if
910 : : end do
911 : :
912 [ + - ]: 10392 : if (compute_nhat) then
913 : : ! GEMM output buffers: C(ic,idat2,idat1) per atom, batched over ia.
914 [ - + + - : 62352 : ABI_MALLOC(gemm_re,(nfgd_max,ndat2,ndat1,nattyp(itypat)))
- + + - -
+ + - - +
+ - + - +
- - + - +
- + ]
915 [ - + + - : 51960 : ABI_MALLOC(gemm_im,(nfgd_max,ndat2,ndat1,nattyp(itypat)))
- + + - -
+ + - - +
+ - - + -
+ - + ]
916 : : #ifdef HAVE_OPENMP_OFFLOAD
917 : : !$OMP TARGET ENTER DATA MAP(alloc:gemm_re,gemm_im) IF(gpu_option_==ABI_GPU_OPENMP)
918 : : #endif
919 : : end if
920 [ + + ]: 10392 : if (compute_grad1) then
921 : : ! GEMM output buffers for the gradient path: C(dir,ic,idat2,idat1) per
922 : : ! atom, batched over ia.
923 [ - + + - : 21576 : ABI_MALLOC(gemm_gr_re,(3,nfgd_max,ndat2,ndat1,nattyp(itypat)))
- + + - -
+ + - - +
+ - + - +
- - + - +
- + ]
924 [ - + + - : 17980 : ABI_MALLOC(gemm_gr_im,(3,nfgd_max,ndat2,ndat1,nattyp(itypat)))
- + + - -
+ + - - +
+ - - + -
+ - + ]
925 : : #ifdef HAVE_OPENMP_OFFLOAD
926 : : !$OMP TARGET ENTER DATA MAP(alloc:gemm_gr_re,gemm_gr_im) IF(gpu_option_==ABI_GPU_OPENMP)
927 : : #endif
928 : : end if
929 : :
930 [ - + + - : 31176 : ABI_MALLOC(atom_nfgd, (nfgd_max))
+ - + - -
+ - + ]
931 [ - + + - : 51960 : ABI_MALLOC(atom_gylm, ( nfgd_max,lm_size,nattyp(itypat)))
- + + - -
+ + - + -
+ - - + -
+ ]
932 [ + - ]: 10392 : if (compute_nhat) then
933 [ - + + - : 51960 : ABI_MALLOC(atom_wgylm,(nfgd_max,lmn2_size,nattyp(itypat)))
- + + - -
+ + - + -
+ - - + -
+ - + ]
934 : : end if
935 [ + + ]: 10392 : if(compute_phonon.and.(.not.qeq0)) then
936 [ - + + - : 31284 : ABI_MALLOC(atom_expiqr, (2,nfgd_max,nattyp(itypat)))
- + + - +
- + - - +
- + ]
937 : : end if
938 [ + + ]: 10392 : if(compute_grad1) then
939 [ - + + - : 17980 : ABI_MALLOC(atom_gylmgr, (3,nfgd_max,lm_size,nattyp(itypat)))
- + + - -
+ + - + -
+ - - + -
+ ]
940 [ - + + - : 17980 : ABI_MALLOC(atom_wgylmgr,(3,nfgd_max,lmn2_size,nattyp(itypat)))
- + + - -
+ + - + -
+ - - + -
+ - + ]
941 : : end if
942 : :
943 [ + + ]: 29748 : do ia=1,nattyp(itypat)
944 : 19356 : iatom=iatm+ia
945 : 19356 : nfgd = pawfgrtab(iatom)%nfgd
946 : :
947 : 19356 : atom_nfgd(ia) = pawfgrtab(iatom)%nfgd
948 [ + + + + ]: 136667140 : atom_gylm(1:nfgd,1:lm_size,ia) = pawfgrtab(iatom)%gylm(1:nfgd,1:lm_size)
949 [ + + ]: 19356 : if(compute_phonon.and.(.not.qeq0)) then
950 [ + + + + ]: 36007884 : atom_expiqr(1:2,1:nfgd,ia) = pawfgrtab(iatom)%expiqr(1:2,1:nfgd)
951 : : end if
952 [ + + ]: 29748 : if(compute_grad1) then
953 [ + + + + : 170666792 : atom_gylmgr(1:3,1:nfgd,1:lm_size,ia)= pawfgrtab(iatom)%gylmgr(1:3,1:nfgd,1:lm_size)
+ + ]
954 : : end if
955 : : end do
956 : :
957 : : #ifdef HAVE_OPENMP_OFFLOAD
958 : : !$OMP TARGET ENTER DATA MAP(to:atom_gylm,atom_indklmn,atom_nfgd,atom_dltij,qijl,gnt_scal) IF(gpu_option_==ABI_GPU_OPENMP)
959 : : !$OMP TARGET ENTER DATA MAP(to:atom_expiqr) IF(gpu_option_==ABI_GPU_OPENMP .and. compute_phonon .and. (.not.qeq0))
960 : : !$OMP TARGET ENTER DATA MAP(to:atom_gylmgr) IF(gpu_option_==ABI_GPU_OPENMP .and. compute_grad1)
961 : : !$OMP TARGET ENTER DATA MAP(alloc:atom_wgylm) IF(gpu_option_==ABI_GPU_OPENMP .and. compute_nhat)
962 : : !$OMP TARGET ENTER DATA MAP(alloc:atom_wgylmgr) IF(gpu_option_==ABI_GPU_OPENMP .and. compute_grad1)
963 : : #endif
964 : :
965 : : ! Precompute the weight factors once per atom here, in order to simplify next compute kernels
966 : : ! for making them more GPU-friendly
967 [ + - ]: 10392 : if (compute_nhat) then
968 [ + - ]: 10392 : if(gpu_option_==ABI_GPU_DISABLED) then
969 : : !$OMP PARALLEL DO PRIVATE(ic,klmn,klm,lmin,lmax,ils,mm,ilslm,wgt)
970 [ + + ]: 29748 : do ia=1,nattyp(itypat)
971 [ + + ]: 15470624 : do ic=1,atom_nfgd(ia)
972 [ + + ]: 561066326 : do klmn=1,lmn2_size
973 : 545606094 : klm =atom_indklmn(1,klmn)
974 : 545606094 : lmin=atom_indklmn(3,klmn)
975 : 545606094 : lmax=atom_indklmn(4,klmn)
976 : 545606094 : wgt=zero
977 [ - + + + ]: 1408938030 : do ils=lmin,lmax,2
978 [ + + ]: 3906288582 : do mm=-ils,ils
979 : 2497350552 : ilslm=ils*ils+ils+mm+1
980 [ + + ]: 3360682488 : if (pawang%gntselect(ilslm,klm)>0) then
981 : 772553124 : wgt=wgt+qijl(ilslm,klmn)*atom_gylm(ic,ilslm,ia)
982 : : end if
983 : : end do
984 : : end do
985 : 561046970 : atom_wgylm(ic,klmn,ia)=atom_dltij(klmn)*half*wgt
986 : : end do
987 : : end do
988 : : end do
989 : : else if(gpu_option_==ABI_GPU_OPENMP) then
990 : : #ifdef HAVE_OPENMP_OFFLOAD
991 : : ! Init to zero to avoid random uninitialised values
992 : : call gpu_set_to_zero(atom_wgylm, int(nfgd_max,c_size_t)*lmn2_size*nattyp(itypat))
993 : : !$OMP TARGET TEAMS DISTRIBUTE MAP(to:atom_gylm,atom_indklmn,atom_dltij,qijl,gnt_scal,atom_nfgd,nattyp,atom_wgylm)&
994 : : !$OMP& PRIVATE(ia)
995 : : do ia=1,nattyp(itypat)
996 : : !$OMP PARALLEL DO COLLAPSE(2) PRIVATE(ic,klmn,klm,lmin,lmax,ils,mm,ilslm,wgt)
997 : : do ic=1,atom_nfgd(ia)
998 : : do klmn=1,lmn2_size
999 : : klm =atom_indklmn(1,klmn)
1000 : : lmin=atom_indklmn(3,klmn)
1001 : : lmax=atom_indklmn(4,klmn)
1002 : : wgt=zero
1003 : : do ils=lmin,lmax,2
1004 : : do mm=-ils,ils
1005 : : ilslm=ils*ils+ils+mm+1
1006 : : wgt=wgt+qijl(ilslm,klmn)*atom_gylm(ic,ilslm,ia)*gnt_scal(ilslm,klm)
1007 : : end do
1008 : : end do
1009 : : atom_wgylm(ic,klmn,ia)=atom_dltij(klmn)*half*wgt
1010 : : end do
1011 : : end do
1012 : : end do
1013 : : #endif
1014 : : end if
1015 : : end if
1016 : :
1017 [ + + ]: 10392 : if (compute_grad1) then
1018 [ + - ]: 3596 : if(gpu_option_==ABI_GPU_DISABLED) then
1019 : : !$OMP PARALLEL DO PRIVATE(ic,klmn,klm,lmin,lmax,ils,mm,ilslm,wgt1,wgt2,wgt3)
1020 [ + + ]: 9360 : do ia=1,nattyp(itypat)
1021 [ + + ]: 5025172 : do ic=1,atom_nfgd(ia)
1022 [ + + ]: 175325366 : do klmn=1,lmn2_size
1023 : 170303790 : klm =atom_indklmn(1,klmn)
1024 : 170303790 : lmin=atom_indklmn(3,klmn)
1025 : 170303790 : lmax=atom_indklmn(4,klmn)
1026 : 170303790 : wgt1=zero; wgt2=zero; wgt3=zero
1027 [ - + + + ]: 439407078 : do ils=lmin,lmax,2
1028 [ + + ]: 1216622070 : do mm=-ils,ils
1029 : 777214992 : ilslm=ils*ils+ils+mm+1
1030 [ + + ]: 1046318280 : if (pawang%gntselect(ilslm,klm)>0) then
1031 : 240874860 : wgt1=wgt1+qijl(ilslm,klmn)*atom_gylmgr(1,ic,ilslm,ia)
1032 : 240874860 : wgt2=wgt2+qijl(ilslm,klmn)*atom_gylmgr(2,ic,ilslm,ia)
1033 : 240874860 : wgt3=wgt3+qijl(ilslm,klmn)*atom_gylmgr(3,ic,ilslm,ia)
1034 : : end if
1035 : : end do
1036 : : end do
1037 : 170303790 : atom_wgylmgr(1,ic,klmn,ia)=atom_dltij(klmn)*half*wgt1
1038 : 170303790 : atom_wgylmgr(2,ic,klmn,ia)=atom_dltij(klmn)*half*wgt2
1039 : 175319602 : atom_wgylmgr(3,ic,klmn,ia)=atom_dltij(klmn)*half*wgt3
1040 : : end do
1041 : : end do
1042 : : end do
1043 [ # # ]: 0 : else if(gpu_option_==ABI_GPU_OPENMP) then
1044 : : ! Init to zero to avoid random uninitialised values
1045 : 0 : call gpu_set_to_zero(atom_wgylmgr, int(3,c_size_t)*nfgd_max*lmn2_size*nattyp(itypat))
1046 : : #ifdef HAVE_OPENMP_OFFLOAD
1047 : : !$OMP TARGET TEAMS DISTRIBUTE MAP(to:atom_gylmgr,atom_indklmn,atom_dltij,qijl,gnt_scal,atom_nfgd,nattyp,atom_wgylmgr)&
1048 : : !$OMP& PRIVATE(ia)
1049 : : do ia=1,nattyp(itypat)
1050 : : !$OMP PARALLEL DO COLLAPSE(2) PRIVATE(ic,klmn,klm,lmin,lmax,ils,mm,ilslm,wgt1,wgt2,wgt3)
1051 : : do ic=1,atom_nfgd(ia)
1052 : : do klmn=1,lmn2_size
1053 : : klm =atom_indklmn(1,klmn)
1054 : : lmin=atom_indklmn(3,klmn)
1055 : : lmax=atom_indklmn(4,klmn)
1056 : : wgt1=zero; wgt2=zero; wgt3=zero
1057 : : do ils=lmin,lmax,2
1058 : : do mm=-ils,ils
1059 : : ilslm=ils*ils+ils+mm+1
1060 : : wgt1=wgt1+qijl(ilslm,klmn)*atom_gylmgr(1,ic,ilslm,ia)*gnt_scal(ilslm,klm)
1061 : : wgt2=wgt2+qijl(ilslm,klmn)*atom_gylmgr(2,ic,ilslm,ia)*gnt_scal(ilslm,klm)
1062 : : wgt3=wgt3+qijl(ilslm,klmn)*atom_gylmgr(3,ic,ilslm,ia)*gnt_scal(ilslm,klm)
1063 : : end do
1064 : : end do
1065 : : atom_wgylmgr(1,ic,klmn,ia)=atom_dltij(klmn)*half*wgt1
1066 : : atom_wgylmgr(2,ic,klmn,ia)=atom_dltij(klmn)*half*wgt2
1067 : : atom_wgylmgr(3,ic,klmn,ia)=atom_dltij(klmn)*half*wgt3
1068 : : end do
1069 : : end do
1070 : : end do
1071 : : #endif
1072 : : end if
1073 : : end if
1074 : :
1075 [ + + ]: 20784 : do isploop=1,nspinor**2 ! Loop over density components of the compensation charge.
1076 : : ! TODO Here we might take advantage of symmetry relations between the four components if nspinor==2
1077 : 10392 : isp1=spinor_idxs(1,isploop)
1078 : 10392 : isp2=spinor_idxs(2,isploop)
1079 : :
1080 [ + + ]: 29748 : do ia=1,nattyp(itypat)
1081 : 19356 : iatom=iatm+ia
1082 [ + - ]: 19356 : if(gpu_option_==ABI_GPU_DISABLED) then
1083 : : !$OMP PARALLEL DO PRIVATE(idat1,idat2,ilmn,jlmn,klmn)
1084 [ + + ]: 38712 : do idat1=1,ndat1
1085 [ + + ]: 115458 : do idat2=1,ndat2
1086 [ + + ]: 2785137 : do klmn=1,lmn2_size ! Loop over ij channels of this atom type.
1087 : 2689035 : ilmn=atom_indklmn(7,klmn)
1088 : 2689035 : jlmn=atom_indklmn(8,klmn)
1089 : : cpf_re(klmn,idat2,idat1,ia) = &
1090 : : & (projs1(1,shift+ilmn,isp1+(idat1-1)*nspinor) * projs2(1,shift+jlmn,isp2+(idat2-1)*nspinor)&
1091 : : & +projs1(2,shift+ilmn,isp1+(idat1-1)*nspinor) * projs2(2,shift+jlmn,isp2+(idat2-1)*nspinor)&
1092 : : & +projs1(1,shift+jlmn,isp1+(idat1-1)*nspinor) * projs2(1,shift+ilmn,isp2+(idat2-1)*nspinor)&
1093 : 2689035 : & +projs1(2,shift+jlmn,isp1+(idat1-1)*nspinor) * projs2(2,shift+ilmn,isp2+(idat2-1)*nspinor))
1094 : :
1095 : : cpf_im(klmn,idat2,idat1,ia) = &
1096 : : & (projs1(1,shift+ilmn,isp1+(idat1-1)*nspinor) * projs2(2,shift+jlmn,isp2+(idat2-1)*nspinor)&
1097 : : & -projs1(2,shift+ilmn,isp1+(idat1-1)*nspinor) * projs2(1,shift+jlmn,isp2+(idat2-1)*nspinor)&
1098 : : & +projs1(1,shift+jlmn,isp1+(idat1-1)*nspinor) * projs2(2,shift+ilmn,isp2+(idat2-1)*nspinor)&
1099 : 2765781 : & -projs1(2,shift+jlmn,isp1+(idat1-1)*nspinor) * projs2(1,shift+ilmn,isp2+(idat2-1)*nspinor))
1100 : : end do
1101 : : end do
1102 : : end do
1103 : : else if(gpu_option_==ABI_GPU_OPENMP) then
1104 : : #ifdef HAVE_OPENMP_OFFLOAD
1105 : : !$OMP TARGET TEAMS DISTRIBUTE COLLAPSE(2) &
1106 : : !$OMP& PRIVATE(idat1,idat2) MAP(to:cpf_re,cpf_im,projs1,projs2,atom_indklmn,atom_nfgd)
1107 : : do idat1=1,ndat1
1108 : : do idat2=1,ndat2
1109 : : !$OMP PARALLEL DO PRIVATE(ilmn,jlmn,klmn)
1110 : : do klmn=1,lmn2_size ! Loop over ij channels of this atom type.
1111 : : ilmn=atom_indklmn(7,klmn)
1112 : : jlmn=atom_indklmn(8,klmn)
1113 : : cpf_re(klmn,idat2,idat1,ia) = &
1114 : : & (projs1(1,shift+ilmn,isp1+(idat1-1)*nspinor) * projs2(1,shift+jlmn,isp2+(idat2-1)*nspinor)&
1115 : : & +projs1(2,shift+ilmn,isp1+(idat1-1)*nspinor) * projs2(2,shift+jlmn,isp2+(idat2-1)*nspinor)&
1116 : : & +projs1(1,shift+jlmn,isp1+(idat1-1)*nspinor) * projs2(1,shift+ilmn,isp2+(idat2-1)*nspinor)&
1117 : : & +projs1(2,shift+jlmn,isp1+(idat1-1)*nspinor) * projs2(2,shift+ilmn,isp2+(idat2-1)*nspinor))
1118 : :
1119 : : cpf_im(klmn,idat2,idat1,ia) = &
1120 : : & (projs1(1,shift+ilmn,isp1+(idat1-1)*nspinor) * projs2(2,shift+jlmn,isp2+(idat2-1)*nspinor)&
1121 : : & -projs1(2,shift+ilmn,isp1+(idat1-1)*nspinor) * projs2(1,shift+jlmn,isp2+(idat2-1)*nspinor)&
1122 : : & +projs1(1,shift+jlmn,isp1+(idat1-1)*nspinor) * projs2(2,shift+ilmn,isp2+(idat2-1)*nspinor)&
1123 : : & -projs1(2,shift+jlmn,isp1+(idat1-1)*nspinor) * projs2(1,shift+ilmn,isp2+(idat2-1)*nspinor))
1124 : : end do
1125 : : end do
1126 : : end do
1127 : : #endif
1128 : : end if
1129 : 29748 : shift = shift + nlmn
1130 : : end do ! ia
1131 : :
1132 [ + - ]: 10392 : if (compute_nhat) then
1133 [ + - ]: 10392 : if(gpu_option_==ABI_GPU_DISABLED) then
1134 : 10392 : gemm_n = ndat2*ndat1
1135 [ + + ]: 29748 : do ia=1,nattyp(itypat)
1136 : : call abi_xgemm('n','n',nfgd_max,gemm_n,lmn2_size,cone,&
1137 : : & atom_wgylm(:,:,ia),nfgd_max,&
1138 : : & cpf_re(:,:,:,ia),lmn2_size,czero,&
1139 : 19356 : & gemm_re(:,:,:,ia),nfgd_max,x_cplx=1)
1140 : : call abi_xgemm('n','n',nfgd_max,gemm_n,lmn2_size,cone,&
1141 : : & atom_wgylm(:,:,ia),nfgd_max,&
1142 : : & cpf_im(:,:,:,ia),lmn2_size,czero,&
1143 : 29748 : & gemm_im(:,:,:,ia),nfgd_max,x_cplx=1)
1144 : : end do
1145 : : !$OMP PARALLEL DO COLLAPSE(2) PRIVATE(idat1,idat2,ic)
1146 [ + + ]: 29748 : do ia=1,nattyp(itypat)
1147 [ + + ]: 49104 : do idat1=1,ndat1
1148 [ + + ]: 115458 : do idat2=1,ndat2
1149 [ + + ]: 60913706 : do ic=1,atom_nfgd(ia)
1150 : 60817604 : nhat12_atm(1,ic,isploop,idat2,idat1,ia)=nhat12_atm(1,ic,isploop,idat2,idat1,ia)+gemm_re(ic,idat2,idat1,ia)
1151 : 60894350 : nhat12_atm(2,ic,isploop,idat2,idat1,ia)=nhat12_atm(2,ic,isploop,idat2,idat1,ia)+gemm_im(ic,idat2,idat1,ia)
1152 : : end do
1153 : : end do
1154 : : end do
1155 : : end do
1156 : : else if(gpu_option_==ABI_GPU_OPENMP) then
1157 : : #ifdef HAVE_OPENMP_OFFLOAD
1158 : : gemm_n = ndat2*ndat1
1159 : : gemm_batch = nattyp(itypat)
1160 : : !$OMP TARGET DATA USE_DEVICE_ADDR(atom_wgylm,cpf_re,cpf_im,gemm_re,gemm_im)
1161 : : call abi_gpu_xgemm_strided(1,'n','n',nfgd_max,gemm_n,lmn2_size,cone,&
1162 : : & c_loc(atom_wgylm),nfgd_max,nfgd_max*lmn2_size,&
1163 : : & c_loc(cpf_re),lmn2_size,lmn2_size*gemm_n,czero,&
1164 : : & c_loc(gemm_re),nfgd_max,nfgd_max*gemm_n,gemm_batch)
1165 : : call abi_gpu_xgemm_strided(1,'n','n',nfgd_max,gemm_n,lmn2_size,cone,&
1166 : : & c_loc(atom_wgylm),nfgd_max,nfgd_max*lmn2_size,&
1167 : : & c_loc(cpf_im),lmn2_size,lmn2_size*gemm_n,czero,&
1168 : : & c_loc(gemm_im),nfgd_max,nfgd_max*gemm_n,gemm_batch)
1169 : : !$OMP END TARGET DATA
1170 : :
1171 : : !$OMP TARGET TEAMS DISTRIBUTE COLLAPSE(3) &
1172 : : !$OMP& MAP(to:nhat12_atm,gemm_re,gemm_im,atom_nfgd,nattyp)&
1173 : : !$OMP& PRIVATE(idat1,idat2,ia)
1174 : : do ia=1,nattyp(itypat)
1175 : : do idat1=1,ndat1
1176 : : do idat2=1,ndat2
1177 : : !$OMP PARALLEL DO PRIVATE(ic)
1178 : : do ic=1,atom_nfgd(ia)
1179 : : nhat12_atm(1,ic,isploop,idat2,idat1,ia)=nhat12_atm(1,ic,isploop,idat2,idat1,ia)+gemm_re(ic,idat2,idat1,ia)
1180 : : nhat12_atm(2,ic,isploop,idat2,idat1,ia)=nhat12_atm(2,ic,isploop,idat2,idat1,ia)+gemm_im(ic,idat2,idat1,ia)
1181 : : end do
1182 : : end do
1183 : : end do
1184 : : end do ! ia
1185 : : #endif
1186 : : end if
1187 : : end if ! compute_nhat
1188 : :
1189 [ + + ]: 10392 : if (compute_grad1) then
1190 [ + - ]: 3596 : if(gpu_option_==ABI_GPU_DISABLED) then
1191 : 3596 : gemm_n = ndat2*ndat1
1192 [ + + ]: 9360 : do ia=1,nattyp(itypat)
1193 : : call abi_xgemm('n','n',3*nfgd_max,gemm_n,lmn2_size,cone,&
1194 : : & atom_wgylmgr(:,:,:,ia),3*nfgd_max,&
1195 : : & cpf_re(:,:,:,ia),lmn2_size,czero,&
1196 : 5764 : & gemm_gr_re(:,:,:,:,ia),3*nfgd_max,x_cplx=1)
1197 : : call abi_xgemm('n','n',3*nfgd_max,gemm_n,lmn2_size,cone,&
1198 : : & atom_wgylmgr(:,:,:,ia),3*nfgd_max,&
1199 : : & cpf_im(:,:,:,ia),lmn2_size,czero,&
1200 : 9360 : & gemm_gr_im(:,:,:,:,ia),3*nfgd_max,x_cplx=1)
1201 : : end do
1202 : : !$OMP PARALLEL DO COLLAPSE(2) PRIVATE(iatom,idat1,idat2,ic)
1203 [ + + ]: 9360 : do ia=1,nattyp(itypat)
1204 [ + + ]: 15124 : do idat1=1,ndat1
1205 [ + + ]: 33906 : do idat2=1,ndat2
1206 : 22378 : iatom=iatm+ia
1207 [ + + ]: 19145490 : do ic=1,atom_nfgd(ia)
1208 : 19117348 : grnhat_12(1,ic,isploop,1,iatom,idat2,idat1)=grnhat_12(1,ic,isploop,1,iatom,idat2,idat1)+gemm_gr_re(1,ic,idat2,idat1,ia)
1209 : 19117348 : grnhat_12(1,ic,isploop,2,iatom,idat2,idat1)=grnhat_12(1,ic,isploop,2,iatom,idat2,idat1)+gemm_gr_re(2,ic,idat2,idat1,ia)
1210 : 19117348 : grnhat_12(1,ic,isploop,3,iatom,idat2,idat1)=grnhat_12(1,ic,isploop,3,iatom,idat2,idat1)+gemm_gr_re(3,ic,idat2,idat1,ia)
1211 : 19117348 : grnhat_12(2,ic,isploop,1,iatom,idat2,idat1)=grnhat_12(2,ic,isploop,1,iatom,idat2,idat1)+gemm_gr_im(1,ic,idat2,idat1,ia)
1212 : 19117348 : grnhat_12(2,ic,isploop,2,iatom,idat2,idat1)=grnhat_12(2,ic,isploop,2,iatom,idat2,idat1)+gemm_gr_im(2,ic,idat2,idat1,ia)
1213 : 19139726 : grnhat_12(2,ic,isploop,3,iatom,idat2,idat1)=grnhat_12(2,ic,isploop,3,iatom,idat2,idat1)+gemm_gr_im(3,ic,idat2,idat1,ia)
1214 : : end do
1215 : : end do
1216 : : end do
1217 : : end do
1218 : : else if(gpu_option_==ABI_GPU_OPENMP) then
1219 : : #ifdef HAVE_OPENMP_OFFLOAD
1220 : : gemm_n = ndat2*ndat1
1221 : : gemm_batch = nattyp(itypat)
1222 : : !$OMP TARGET DATA USE_DEVICE_ADDR(atom_wgylmgr,cpf_re,cpf_im,gemm_gr_re,gemm_gr_im)
1223 : : call abi_gpu_xgemm_strided(1,'n','n',3*nfgd_max,gemm_n,lmn2_size,cone,&
1224 : : & c_loc(atom_wgylmgr),3*nfgd_max,3*nfgd_max*lmn2_size,&
1225 : : & c_loc(cpf_re),lmn2_size,lmn2_size*gemm_n,czero,&
1226 : : & c_loc(gemm_gr_re),3*nfgd_max,3*nfgd_max*gemm_n,gemm_batch)
1227 : : call abi_gpu_xgemm_strided(1,'n','n',3*nfgd_max,gemm_n,lmn2_size,cone,&
1228 : : & c_loc(atom_wgylmgr),3*nfgd_max,3*nfgd_max*lmn2_size,&
1229 : : & c_loc(cpf_im),lmn2_size,lmn2_size*gemm_n,czero,&
1230 : : & c_loc(gemm_gr_im),3*nfgd_max,3*nfgd_max*gemm_n,gemm_batch)
1231 : : !$OMP END TARGET DATA
1232 : : ! Scatter the dense (nfgd_max-wide) GEMM output back into
1233 : : ! grnhat_12, honoring the real per-atom sphere size atom_nfgd(ia).
1234 : : ! grnhat_12 is itself sized on nfgd_max (per-atom sphere), so the local
1235 : : ! sphere point index ic is used directly (no mapping to the global FFT grid).
1236 : : !$OMP TARGET TEAMS DISTRIBUTE COLLAPSE(3) &
1237 : : !$OMP& MAP(to:grnhat_12,gemm_gr_re,gemm_gr_im,atom_nfgd,nattyp)&
1238 : : !$OMP& PRIVATE(idat1,idat2,ia,iatom)
1239 : : do ia=1,nattyp(itypat)
1240 : : do idat1=1,ndat1
1241 : : do idat2=1,ndat2
1242 : : !$OMP PARALLEL DO PRIVATE(ic,iatom)
1243 : : do ic=1,atom_nfgd(ia)
1244 : : iatom=iatm+ia
1245 : : grnhat_12(1,ic,isploop,1,iatom,idat2,idat1)=grnhat_12(1,ic,isploop,1,iatom,idat2,idat1)+gemm_gr_re(1,ic,idat2,idat1,ia)
1246 : : grnhat_12(1,ic,isploop,2,iatom,idat2,idat1)=grnhat_12(1,ic,isploop,2,iatom,idat2,idat1)+gemm_gr_re(2,ic,idat2,idat1,ia)
1247 : : grnhat_12(1,ic,isploop,3,iatom,idat2,idat1)=grnhat_12(1,ic,isploop,3,iatom,idat2,idat1)+gemm_gr_re(3,ic,idat2,idat1,ia)
1248 : : grnhat_12(2,ic,isploop,1,iatom,idat2,idat1)=grnhat_12(2,ic,isploop,1,iatom,idat2,idat1)+gemm_gr_im(1,ic,idat2,idat1,ia)
1249 : : grnhat_12(2,ic,isploop,2,iatom,idat2,idat1)=grnhat_12(2,ic,isploop,2,iatom,idat2,idat1)+gemm_gr_im(2,ic,idat2,idat1,ia)
1250 : : grnhat_12(2,ic,isploop,3,iatom,idat2,idat1)=grnhat_12(2,ic,isploop,3,iatom,idat2,idat1)+gemm_gr_im(3,ic,idat2,idat1,ia)
1251 : : end do
1252 : : end do
1253 : : end do
1254 : : end do ! ia
1255 : : #endif
1256 : : end if
1257 : : end if ! compute_grad1
1258 : :
1259 [ + - ]: 10392 : if (compute_nhat) then
1260 : : ! If needed, multiply eventually by exp(-i.q.r) phase
1261 [ + + + - ]: 10392 : if(compute_phonon.and.(.not.qeq0).and.pawfgrtab(iatom)%expiqr_allocated/=0) then
1262 [ + - ]: 7821 : if(gpu_option_==ABI_GPU_DISABLED) then
1263 : : !$OMP PARALLEL DO COLLAPSE(2) PRIVATE(ro,ro_ql,ic)
1264 [ + + ]: 23463 : do ia=1,nattyp(itypat)
1265 [ + + ]: 39105 : do idat1=1,ndat1
1266 [ + + ]: 93852 : do idat2=1,ndat2
1267 [ + + ]: 48067866 : do ic=1,atom_nfgd(ia)
1268 : 47989656 : iatom=iatm+ia
1269 [ + + ]: 143968968 : ro(1:2)=nhat12_atm(1:2,ic,isploop,idat2,idat1,ia)
1270 : 47989656 : nhat12_atm(1,ic,isploop,idat2,idat1,ia)=ro(1)*atom_expiqr(1,ic,ia)-ro(2)*atom_expiqr(2,ic,ia)
1271 : 48052224 : nhat12_atm(2,ic,isploop,idat2,idat1,ia)=ro(2)*atom_expiqr(1,ic,ia)+ro(1)*atom_expiqr(2,ic,ia)
1272 : : end do
1273 : : end do
1274 : : end do
1275 : : end do ! ia
1276 : : else if(gpu_option_==ABI_GPU_OPENMP) then
1277 : : #ifdef HAVE_OPENMP_OFFLOAD
1278 : : !$OMP TARGET TEAMS DISTRIBUTE COLLAPSE(2) &
1279 : : !$OMP& MAP(to:atom_expiqr,atom_nfgd,nhat12_atm,nattyp)
1280 : : do ia=1,nattyp(itypat)
1281 : : do idat1=1,ndat1
1282 : : !$OMP PARALLEL DO COLLAPSE(2) PRIVATE(ic,ro)
1283 : : do idat2=1,ndat2
1284 : : do ic=1,atom_nfgd(ia)
1285 : : ro(1)=nhat12_atm(1,ic,isploop,idat2,idat1,ia)
1286 : : ro(2)=nhat12_atm(2,ic,isploop,idat2,idat1,ia)
1287 : : nhat12_atm(1,ic,isploop,idat2,idat1,ia)=ro(1)*atom_expiqr(1,ic,ia)-ro(2)*atom_expiqr(2,ic,ia)
1288 : : nhat12_atm(2,ic,isploop,idat2,idat1,ia)=ro(2)*atom_expiqr(1,ic,ia)+ro(1)*atom_expiqr(2,ic,ia)
1289 : : end do
1290 : : end do
1291 : : end do
1292 : : end do ! ia
1293 : : #endif
1294 : : end if
1295 : : end if
1296 : : end if
1297 [ + + ]: 20784 : if (compute_grad1) then
1298 [ + + + - ]: 3596 : if(compute_phonon.and.(.not.qeq0).and.pawfgrtab(iatom)%expiqr_allocated/=0) then
1299 [ + - ]: 1764 : if(gpu_option_==ABI_GPU_DISABLED) then
1300 : : !$OMP PARALLEL DO COLLAPSE(2) PRIVATE(ia,idat1,idat2,ro,ro_ql,ic)
1301 [ + + ]: 5292 : do ia=1,nattyp(itypat)
1302 [ + + ]: 8820 : do idat1=1,ndat1
1303 [ + + ]: 21168 : do idat2=1,ndat2
1304 [ + + ]: 10841544 : do ic=1,atom_nfgd(ia)
1305 : 10823904 : iatom=iatm+ia
1306 : 10823904 : ro_ql(1)= atom_expiqr(1,ic,ia)
1307 : 10823904 : ro_ql(2)= atom_expiqr(2,ic,ia)
1308 : 10823904 : ro(1)=grnhat_12(1,ic,isploop,1,iatom,idat2,idat1)
1309 : 10823904 : ro(2)=grnhat_12(2,ic,isploop,1,iatom,idat2,idat1)
1310 : 10823904 : grnhat_12(1,ic,isploop,1,iatom,idat2,idat1)=ro(1)*ro_ql(1)-ro(2)*ro_ql(2)
1311 : 10823904 : grnhat_12(2,ic,isploop,1,iatom,idat2,idat1)=ro(2)*ro_ql(1)+ro(1)*ro_ql(2)
1312 : 10823904 : ro(1)=grnhat_12(1,ic,isploop,2,iatom,idat2,idat1)
1313 : 10823904 : ro(2)=grnhat_12(2,ic,isploop,2,iatom,idat2,idat1)
1314 : 10823904 : grnhat_12(1,ic,isploop,2,iatom,idat2,idat1)=ro(1)*ro_ql(1)-ro(2)*ro_ql(2)
1315 : 10823904 : grnhat_12(2,ic,isploop,2,iatom,idat2,idat1)=ro(2)*ro_ql(1)+ro(1)*ro_ql(2)
1316 : 10823904 : ro(1)=grnhat_12(1,ic,isploop,3,iatom,idat2,idat1)
1317 : 10823904 : ro(2)=grnhat_12(2,ic,isploop,3,iatom,idat2,idat1)
1318 : 10823904 : grnhat_12(1,ic,isploop,3,iatom,idat2,idat1)=ro(1)*ro_ql(1)-ro(2)*ro_ql(2)
1319 : 10838016 : grnhat_12(2,ic,isploop,3,iatom,idat2,idat1)=ro(2)*ro_ql(1)+ro(1)*ro_ql(2)
1320 : : end do
1321 : : end do
1322 : : end do
1323 : : end do ! ia
1324 : : else if(gpu_option_==ABI_GPU_OPENMP) then
1325 : : #ifdef HAVE_OPENMP_OFFLOAD
1326 : : !$OMP TARGET TEAMS DISTRIBUTE COLLAPSE(3) &
1327 : : !$OMP& MAP(to:atom_expiqr,atom_nfgd,nhat12_atm,grnhat_12,nattyp) PRIVATE(idat1,idat2)
1328 : : do ia=1,nattyp(itypat)
1329 : : do idat1=1,ndat1
1330 : : do idat2=1,ndat2
1331 : : !$OMP PARALLEL DO PRIVATE(iatom,ic,ro,ro_ql)
1332 : : do ic=1,atom_nfgd(ia)
1333 : : iatom=iatm+ia
1334 : : ro_ql(1)= atom_expiqr(1,ic,ia)
1335 : : ro_ql(2)= atom_expiqr(2,ic,ia)
1336 : : ro(1)=grnhat_12(1,ic,isploop,1,iatom,idat2,idat1)
1337 : : ro(2)=grnhat_12(2,ic,isploop,1,iatom,idat2,idat1)
1338 : : grnhat_12(1,ic,isploop,1,iatom,idat2,idat1)=ro(1)*ro_ql(1)-ro(2)*ro_ql(2)
1339 : : grnhat_12(2,ic,isploop,1,iatom,idat2,idat1)=ro(2)*ro_ql(1)+ro(1)*ro_ql(2)
1340 : : ro(1)=grnhat_12(1,ic,isploop,2,iatom,idat2,idat1)
1341 : : ro(2)=grnhat_12(2,ic,isploop,2,iatom,idat2,idat1)
1342 : : grnhat_12(1,ic,isploop,2,iatom,idat2,idat1)=ro(1)*ro_ql(1)-ro(2)*ro_ql(2)
1343 : : grnhat_12(2,ic,isploop,2,iatom,idat2,idat1)=ro(2)*ro_ql(1)+ro(1)*ro_ql(2)
1344 : : ro(1)=grnhat_12(1,ic,isploop,3,iatom,idat2,idat1)
1345 : : ro(2)=grnhat_12(2,ic,isploop,3,iatom,idat2,idat1)
1346 : : grnhat_12(1,ic,isploop,3,iatom,idat2,idat1)=ro(1)*ro_ql(1)-ro(2)*ro_ql(2)
1347 : : grnhat_12(2,ic,isploop,3,iatom,idat2,idat1)=ro(2)*ro_ql(1)+ro(1)*ro_ql(2)
1348 : : end do
1349 : : end do
1350 : : end do
1351 : : end do ! ia
1352 : : #endif
1353 : : end if
1354 : : end if
1355 : : end if
1356 : :
1357 : : end do ! isploop (density components of the compensation charge)
1358 : :
1359 [ + + ]: 29748 : do ia=1,nattyp(itypat)
1360 : 19356 : iatom=iatm+ia
1361 [ + - ]: 19356 : if (pawfgrtab(iatom)%gylm_allocated==2) then
1362 [ - + ]: 19356 : ABI_FREE(pawfgrtab(iatom)%gylm)
1363 [ - + ]: 19356 : ABI_MALLOC(pawfgrtab(iatom)%gylm,(0,0))
1364 : 19356 : pawfgrtab(iatom)%gylm_allocated=0
1365 : : end if
1366 [ + + ]: 19356 : if (pawfgrtab(iatom)%gylmgr_allocated==2) then
1367 [ - + ]: 5764 : ABI_FREE(pawfgrtab(iatom)%gylmgr)
1368 [ - + ]: 5764 : ABI_MALLOC(pawfgrtab(iatom)%gylmgr,(0,0,0))
1369 : 5764 : pawfgrtab(iatom)%gylmgr_allocated=0
1370 : : end if
1371 [ + + ]: 29748 : if (pawfgrtab(iatom)%expiqr_allocated==2) then
1372 [ - + ]: 15642 : ABI_FREE(pawfgrtab(iatom)%expiqr)
1373 [ - + ]: 15642 : ABI_MALLOC(pawfgrtab(iatom)%expiqr,(0,0))
1374 : 15642 : pawfgrtab(iatom)%expiqr_allocated=0
1375 : : end if
1376 : : end do ! ia
1377 : :
1378 : 10392 : iatm=iatm+nattyp(itypat)
1379 : : #ifdef HAVE_OPENMP_OFFLOAD
1380 : : !$OMP TARGET EXIT DATA MAP(delete:atom_nfgd,atom_indklmn,atom_gylm,atom_dltij,qijl,gnt_scal) IF(gpu_option_==ABI_GPU_OPENMP)
1381 : : !$OMP TARGET EXIT DATA MAP(delete:atom_expiqr) IF(gpu_option_==ABI_GPU_OPENMP .and. compute_phonon .and. (.not.qeq0))
1382 : : !$OMP TARGET EXIT DATA MAP(delete:atom_gylmgr) IF(gpu_option_==ABI_GPU_OPENMP .and. compute_grad1)
1383 : : !$OMP TARGET EXIT DATA MAP(delete:atom_wgylm) IF(gpu_option_==ABI_GPU_OPENMP .and. compute_nhat)
1384 : : !$OMP TARGET EXIT DATA MAP(delete:atom_wgylmgr) IF(gpu_option_==ABI_GPU_OPENMP .and. compute_grad1)
1385 : : #endif
1386 [ - + ]: 10392 : ABI_FREE(atom_nfgd)
1387 [ - + ]: 10392 : ABI_FREE(atom_gylm)
1388 [ + + ]: 10392 : if (compute_grad1) then
1389 [ - + ]: 3596 : ABI_FREE(atom_gylmgr)
1390 : : end if
1391 [ + + ]: 10392 : if (compute_phonon.and.(.not.qeq0)) then
1392 [ - + ]: 7821 : ABI_FREE(atom_expiqr)
1393 : : end if
1394 [ - + ]: 10392 : ABI_FREE(qijl)
1395 [ + - ]: 10392 : if (compute_nhat) then
1396 [ - + ]: 10392 : ABI_FREE(atom_wgylm)
1397 : : end if
1398 [ + + ]: 10392 : if (compute_grad1) then
1399 [ - + ]: 3596 : ABI_FREE(atom_wgylmgr)
1400 : : end if
1401 : : #ifdef HAVE_OPENMP_OFFLOAD
1402 : : !$OMP TARGET EXIT DATA MAP(delete:cpf_re,cpf_im) IF(gpu_option_==ABI_GPU_OPENMP)
1403 : : !$OMP TARGET EXIT DATA MAP(delete:gemm_re,gemm_im) IF(compute_nhat.and.gpu_option_==ABI_GPU_OPENMP)
1404 : : !$OMP TARGET EXIT DATA MAP(delete:gemm_gr_re,gemm_gr_im) IF(compute_grad1.and.gpu_option_==ABI_GPU_OPENMP)
1405 : : #endif
1406 [ - + ]: 10392 : ABI_FREE(cpf_re)
1407 [ - + ]: 10392 : ABI_FREE(cpf_im)
1408 [ + - ]: 10392 : if (compute_nhat) then
1409 [ - + ]: 10392 : ABI_FREE(gemm_re)
1410 [ - + ]: 10392 : ABI_FREE(gemm_im)
1411 : : end if
1412 [ + + ]: 10392 : if (compute_grad1) then
1413 [ - + ]: 3596 : ABI_FREE(gemm_gr_re)
1414 [ - + ]: 3596 : ABI_FREE(gemm_gr_im)
1415 : : end if
1416 : 20070 : nullify(nhat12_atm)
1417 : : end do ! itypat
1418 : :
1419 : : #ifdef HAVE_OPENMP_OFFLOAD
1420 : : !$OMP TARGET EXIT DATA MAP(delete:projs1,projs2,nattyp) IF(gpu_option_==ABI_GPU_OPENMP)
1421 : : #endif
1422 [ - + ]: 9678 : ABI_FREE(projs1)
1423 [ - + ]: 9678 : ABI_FREE(projs2)
1424 [ - + ]: 9678 : ABI_FREE(gnt_scal)
1425 : :
1426 [ + + ]: 9678 : if (compute_grad1) then
1427 : : select case (gpu_option_)
1428 : : case (ABI_GPU_DISABLED)
1429 [ + + + + : 200627845 : grnhat_12=-grnhat_12
+ + + + +
+ + + +
+ ]
1430 : : case (ABI_GPU_OPENMP)
1431 [ # # ]: 0 : call abi_xscal(size(grnhat_12),cminusone,grnhat_12,1,x_cplx=1,gpu_option=gpu_option_)
1432 : : case default
1433 [ + - - ]: 2882 : ABI_BUG("Unsupported GPU option")
1434 : : end select
1435 : : end if
1436 : :
1437 : : !----- Reduction in case of parallelism -----!
1438 [ - + ]: 9678 : if (paral_atom)then
1439 : 0 : call timab(48,1,tsec)
1440 [ # # ]: 0 : if (compute_nhat) then
1441 : 0 : call xmpi_sum(nhat12,my_comm_atom,ierr)
1442 : : end if
1443 [ # # ]: 0 : if (compute_grad) then
1444 : 0 : call xmpi_sum(grnhat12,my_comm_atom,ierr)
1445 : : end if
1446 [ # # ]: 0 : if (compute_grad1) then
1447 : 0 : call xmpi_sum(grnhat_12,my_comm_atom,ierr)
1448 : : end if
1449 : 0 : call timab(48,2,tsec)
1450 : : end if
1451 : :
1452 : : !Destroy atom table used for parallelism
1453 : 9678 : call free_my_atmtab(my_atmtab,my_atmtab_allocated)
1454 : :
1455 : : DBG_EXIT("COLL")
1456 : :
1457 [ - + - + : 19356 : end subroutine pawmknhat_psipsi_ndat
- + - + -
+ - + -
+ ]
1458 : : !!***
1459 : :
1460 : 0 : subroutine pawmknhat_psipsi(cprj1,cprj2,ider,izero,my_natom,natom,nfft,ngfft,nhat12_grdim,&
1461 : 0 : & nspinor,ntypat,ndat1,ndat2,pawang,pawfgrtab,grnhat12,nhat12,pawtab, &
1462 [ # # # # : 0 : & gprimd,grnhat_12,qphon,xred,atindx,mpi_atmtab,comm_atom,comm_fft,me_g0,paral_kgb,&
# # # # #
# # # ]
1463 : : & distribfft) ! optional arguments
1464 : :
1465 : : !Arguments ---------------------------------------------
1466 : : !scalars
1467 : : integer,intent(in) :: ider,izero,my_natom,natom,nfft,nhat12_grdim,ntypat,nspinor,ndat1,ndat2
1468 : : integer,optional,intent(in) :: me_g0,comm_fft,paral_kgb
1469 : : integer,optional,intent(in) :: comm_atom
1470 : : type(distribfft_type),optional,intent(in),target :: distribfft
1471 : : type(pawang_type),intent(in) :: pawang
1472 : : !arrays
1473 : : integer,intent(in) :: ngfft(18)
1474 : : integer,optional,intent(in) ::atindx(natom)
1475 : : integer,optional,target,intent(in) :: mpi_atmtab(:)
1476 : : real(dp),optional, intent(in) ::gprimd(3,3),qphon(3),xred(3,natom)
1477 : : real(dp),intent(out) :: grnhat12(2,nfft,nspinor**2,3*nhat12_grdim,ndat2,ndat1)
1478 : : real(dp),optional,intent(out) :: grnhat_12(:,:,:,:,:,:,:)
1479 : : real(dp),intent(out) :: nhat12(2,nfft,nspinor**2,ndat2,ndat1)
1480 : : type(pawfgrtab_type),intent(inout) :: pawfgrtab(my_natom)
1481 : : type(pawtab_type),intent(in) :: pawtab(ntypat)
1482 : : type(pawcprj_type),intent(in) :: cprj1(natom,nspinor*ndat1),cprj2(natom,nspinor*ndat2)
1483 : :
1484 : : !Local variables ---------------------------------------
1485 : : !scalars
1486 : : integer :: iatm,iatom,iatom_tot,ic,ierr,ils,ilslm,isp1,isp2,isploop,itypat,jc,klm,klmn,idat1,idat2
1487 : : integer :: lmax,lmin,lm_size,mm,my_comm_atom,my_comm_fft,optgr0,optgr1,paral_kgb_fft
1488 : : integer :: cplex,ilmn,jlmn,lmn_size,lmn2_size
1489 : : real(dp) :: re_p,im_p
1490 : : logical :: compute_grad,compute_grad1,compute_nhat,my_atmtab_allocated,paral_atom,qeq0,compute_phonon,order
1491 : : type(distribfft_type),pointer :: my_distribfft
1492 : 0 : type(mpi_type) :: mpi_enreg_fft
1493 : : !arrays
1494 : : integer,parameter :: spinor_idxs(2,4)=RESHAPE((/1,1,2,2,1,2,2,1/),(/2,4/))
1495 : 0 : integer,pointer :: my_atmtab(:)
1496 : : real(dp) :: rdum(1),cpf(2),cpf_ql(2),tsec(2),ro(2),ro_ql(2)
1497 : 0 : real(dp),allocatable :: work(:,:), qijl(:,:), nhat12_atm(:,:,:,:,:)
1498 : :
1499 : : ! *************************************************************************
1500 : :
1501 : : DBG_ENTER("COLL")
1502 : :
1503 : : !Compatibility tests
1504 : 0 : if (present(comm_fft)) then
1505 [ # # ]: 0 : if ((.not.present(paral_kgb)).or.(.not.present(me_g0))) then
1506 : 0 : ABI_BUG('Need paral_kgb and me_g0 with comm_fft!')
1507 : : end if
1508 [ # # ]: 0 : if (present(paral_kgb)) then
1509 [ # # ]: 0 : if (paral_kgb/=0) then
1510 : 0 : ABI_BUG('paral_kgb/=0 not coded!')
1511 : : end if
1512 : : end if
1513 : : end if
1514 : : if (ider>0.and.nhat12_grdim==0) then
1515 : : ! ABI_BUG('Gradients of nhat required but not allocated !')
1516 : : end if
1517 [ # # ]: 0 : if (nspinor==2) then
1518 : 0 : ABI_BUG('nspinor==2 not coded!')
1519 : : end if
1520 : :
1521 : 0 : compute_phonon=.false.;qeq0=.false.
1522 [ # # # # ]: 0 : if (present(gprimd).and.present(qphon).and.present(xred)) compute_phonon=.true.
1523 : 0 : if (compute_phonon) qeq0=(qphon(1)**2+qphon(2)**2+qphon(3)**2<1.d-15)
1524 : : if (present(atindx)) order=.true.
1525 : : !Set up parallelism over atoms
1526 [ # # # # ]: 0 : paral_atom=(present(comm_atom).and.(my_natom/=natom))
1527 [ # # ]: 0 : nullify(my_atmtab);if (present(mpi_atmtab)) my_atmtab => mpi_atmtab
1528 [ # # ]: 0 : my_comm_atom=xmpi_comm_self;if (present(comm_atom)) my_comm_atom=comm_atom
1529 : 0 : call get_my_atmtab(my_comm_atom,my_atmtab,my_atmtab_allocated,paral_atom,natom,my_natom_ref=my_natom)
1530 : :
1531 : : !Initialisations
1532 [ # # # # ]: 0 : compute_nhat=(ider==0.or.ider==2.or.ider==3)
1533 : 0 : compute_grad=(ider==1.or.ider==2)
1534 : 0 : compute_grad1=(ider==3)
1535 [ # # ]: 0 : if ((.not.compute_nhat).and.(.not.compute_grad)) return
1536 : :
1537 [ # # # # : 0 : if (compute_nhat) nhat12=zero
# # # # #
# # # ]
1538 [ # # # # : 0 : if (compute_grad) grnhat12=zero
# # # # #
# # # #
# ]
1539 [ # # # # : 0 : if (compute_grad1) grnhat_12=zero
# # # # #
# # # # #
# # ]
1540 : :
1541 : : if (compute_grad) then
1542 : : ! ABI_BUG('compute_grad not tested!')
1543 : : end if
1544 : :
1545 : : !------------------------------------------------------------------------
1546 : : !----- Loop over atoms
1547 : : !------------------------------------------------------------------------
1548 [ # # ]: 0 : do iatom=1,my_natom
1549 [ # # ]: 0 : iatom_tot=iatom;if (paral_atom) iatom_tot=my_atmtab(iatom)
1550 : 0 : iatm=iatom_tot
1551 : 0 : if (order) iatm=atindx(iatom_tot)
1552 : 0 : itypat = pawfgrtab(iatom)%itypat
1553 : 0 : lm_size = pawfgrtab(iatom)%l_size**2
1554 : 0 : lmn_size = pawtab(itypat)%lmn_size
1555 : 0 : lmn2_size = pawtab(itypat)%lmn2_size
1556 [ # # # # : 0 : ABI_MALLOC(qijl,(lm_size,lmn2_size))
# # # # #
# # # # #
# # ]
1557 [ # # # # ]: 0 : qijl=zero
1558 [ # # # # : 0 : qijl=pawtab(itypat)%qijl
# # # # #
# # # # #
# # ]
1559 [ # # # # : 0 : ABI_MALLOC(nhat12_atm, (2,nfft,nspinor**2,ndat2,ndat1))
# # # # #
# # # # #
# # # # #
# # # #
# ]
1560 [ # # # # : 0 : if (compute_nhat) nhat12_atm=zero
# # # # #
# # # ]
1561 : :
1562 : : ! Eventually compute g_l(r).Y_lm(r) factors for the current atom (if not already done)
1563 [ # # # # : 0 : if (((compute_nhat).and.(pawfgrtab(iatom)%gylm_allocated==0)).or.&
# # ]
1564 : : & (((compute_grad).or.(compute_grad1)).and.(pawfgrtab(iatom)%gylmgr_allocated==0))) then
1565 : 0 : optgr0=0; optgr1=0
1566 [ # # # # ]: 0 : if ((compute_nhat).and.(pawfgrtab(iatom)%gylm_allocated==0)) then
1567 [ # # ]: 0 : if (allocated(pawfgrtab(iatom)%gylm)) then
1568 : 0 : ABI_FREE(pawfgrtab(iatom)%gylm)
1569 : : end if
1570 [ # # # # : 0 : ABI_MALLOC(pawfgrtab(iatom)%gylm,(pawfgrtab(iatom)%nfgd,pawfgrtab(iatom)%l_size**2))
# # # # #
# # # # #
# # # # ]
1571 : 0 : pawfgrtab(iatom)%gylm_allocated=2;optgr0=1
1572 : : end if
1573 [ # # # # ]: 0 : if (((compute_grad).or.(compute_grad1)).and.(pawfgrtab(iatom)%gylmgr_allocated==0)) then
1574 [ # # ]: 0 : if (allocated(pawfgrtab(iatom)%gylmgr)) then
1575 : 0 : ABI_FREE(pawfgrtab(iatom)%gylmgr)
1576 : : end if
1577 [ # # # # : 0 : ABI_MALLOC(pawfgrtab(iatom)%gylmgr,(3,pawfgrtab(iatom)%nfgd,pawfgrtab(iatom)%l_size**2))
# # # # #
# # # # #
# # # # ]
1578 : 0 : pawfgrtab(iatom)%gylmgr_allocated=2;optgr1=1
1579 : : end if
1580 [ # # ]: 0 : if (optgr0+optgr1>0) then
1581 : : call pawgylm(pawfgrtab(iatom)%gylm,pawfgrtab(iatom)%gylmgr,rdum,&
1582 : : & lm_size,pawfgrtab(iatom)%nfgd,optgr0,optgr1,0,pawtab(itypat),&
1583 : 0 : & pawfgrtab(iatom)%rfgd)
1584 : : end if
1585 : :
1586 : : end if
1587 [ # # # # ]: 0 : if (compute_phonon.and.(.not.qeq0).and.(pawfgrtab(iatom)%expiqr_allocated==0)) then
1588 [ # # ]: 0 : if (allocated(pawfgrtab(iatom)%expiqr)) then
1589 : 0 : ABI_FREE(pawfgrtab(iatom)%expiqr)
1590 : : end if
1591 [ # # # # : 0 : ABI_MALLOC(pawfgrtab(iatom)%expiqr,(2,pawfgrtab(iatom)%nfgd))
# # # # #
# # # #
# ]
1592 : : call pawexpiqr(pawfgrtab(iatom)%expiqr,gprimd,pawfgrtab(iatom)%nfgd,qphon,&
1593 [ # # ]: 0 : & pawfgrtab(iatom)%rfgd,xred(:,iatom_tot))
1594 : 0 : pawfgrtab(iatom)%expiqr_allocated=2
1595 : : end if
1596 : :
1597 : : !$OMP PARALLEL DO COLLAPSE(2) PRIVATE(ilslm,ils,mm,ic,jc,cpf_ql) &
1598 : : !$OMP& PRIVATE(isp1,isp2,klm,lmin,lmax,ilmn,jlmn,re_p,im_p,cpf,ro,ro_ql)
1599 [ # # ]: 0 : do idat1=1,ndat1
1600 [ # # ]: 0 : do idat2=1,ndat2
1601 [ # # ]: 0 : do isploop=1,nspinor**2 ! Loop over density components of the compensation charge.
1602 : : ! TODO Here we might take advantage of symmetry relations between the four components if nspinor==2
1603 : 0 : isp1=spinor_idxs(1,isploop)
1604 : 0 : isp2=spinor_idxs(2,isploop)
1605 : :
1606 [ # # ]: 0 : do klmn=1,lmn2_size ! Loop over ij channels of this atom type.
1607 : 0 : klm =pawtab(itypat)%indklmn(1,klmn)
1608 : 0 : lmin=pawtab(itypat)%indklmn(3,klmn) ! abs(il-jl)
1609 : 0 : lmax=pawtab(itypat)%indklmn(4,klmn) ! il+jl
1610 : 0 : ilmn=pawtab(itypat)%indklmn(7,klmn)
1611 : 0 : jlmn=pawtab(itypat)%indklmn(8,klmn)
1612 : : ! call klmn2ijlmn(klmn,lmn_size,ilmn,jlmn) ! This mapping should be stored in pawtab_type
1613 : :
1614 : : ! Retrieve the factor due to the PAW projections.
1615 : : re_p = cprj1(iatm,isp1+(idat1-1)*nspinor)%cp(1,ilmn) * cprj2(iatm,isp2+(idat2-1)*nspinor)%cp(1,jlmn) &
1616 : : & +cprj1(iatm,isp1+(idat1-1)*nspinor)%cp(2,ilmn) * cprj2(iatm,isp2+(idat2-1)*nspinor)%cp(2,jlmn) &
1617 : : & +cprj1(iatm,isp1+(idat1-1)*nspinor)%cp(1,jlmn) * cprj2(iatm,isp2+(idat2-1)*nspinor)%cp(1,ilmn) &
1618 : 0 : & +cprj1(iatm,isp1+(idat1-1)*nspinor)%cp(2,jlmn) * cprj2(iatm,isp2+(idat2-1)*nspinor)%cp(2,ilmn)
1619 : :
1620 : : im_p = cprj1(iatm,isp1+(idat1-1)*nspinor)%cp(1,ilmn) * cprj2(iatm,isp2+(idat2-1)*nspinor)%cp(2,jlmn) &
1621 : : & -cprj1(iatm,isp1+(idat1-1)*nspinor)%cp(2,ilmn) * cprj2(iatm,isp2+(idat2-1)*nspinor)%cp(1,jlmn) &
1622 : : & +cprj1(iatm,isp1+(idat1-1)*nspinor)%cp(1,jlmn) * cprj2(iatm,isp2+(idat2-1)*nspinor)%cp(2,ilmn) &
1623 : 0 : & -cprj1(iatm,isp1+(idat1-1)*nspinor)%cp(2,jlmn) * cprj2(iatm,isp2+(idat2-1)*nspinor)%cp(1,ilmn)
1624 : :
1625 : 0 : cpf(1)=re_p*pawtab(itypat)%dltij(klmn)*half
1626 : 0 : cpf(2)=im_p*pawtab(itypat)%dltij(klmn)*half
1627 : :
1628 [ # # ]: 0 : if (compute_nhat) then
1629 [ # # # # ]: 0 : do ils=lmin,lmax,2 ! Sum over (L,M)
1630 [ # # ]: 0 : do mm=-ils,ils
1631 : 0 : ilslm=ils*ils+ils+mm+1
1632 [ # # ]: 0 : if (pawang%gntselect(ilslm,klm)>0) then
1633 : 0 : cpf_ql(1)=cpf(1)*qijl(ilslm,klmn)
1634 : 0 : cpf_ql(2)=cpf(2)*qijl(ilslm,klmn)
1635 : : !!$OMP PARALLEL DO PRIVATE(ic,jc)
1636 [ # # ]: 0 : do ic=1,pawfgrtab(iatom)%nfgd
1637 : 0 : jc=pawfgrtab(iatom)%ifftsph(ic)
1638 : 0 : nhat12_atm(1,jc,isploop,idat2,idat1)=nhat12_atm(1,jc,isploop,idat2,idat1)+cpf_ql(1)*pawfgrtab(iatom)%gylm(ic,ilslm)
1639 : 0 : nhat12_atm(2,jc,isploop,idat2,idat1)=nhat12_atm(2,jc,isploop,idat2,idat1)+cpf_ql(2)*pawfgrtab(iatom)%gylm(ic,ilslm)
1640 : : end do
1641 : : end if
1642 : : end do
1643 : : end do
1644 : : end if ! compute_nhat
1645 : :
1646 [ # # ]: 0 : if (compute_grad) then
1647 [ # # # # ]: 0 : do ils=lmin,lmax,2 ! Sum over (L,M)
1648 [ # # ]: 0 : do mm=-ils,ils
1649 : 0 : ilslm=ils*ils+ils+mm+1
1650 [ # # ]: 0 : if (pawang%gntselect(ilslm,klm)>0) then
1651 : 0 : cpf_ql(1)=cpf(1)*qijl(ilslm,klmn)
1652 : 0 : cpf_ql(2)=cpf(2)*qijl(ilslm,klmn)
1653 [ # # ]: 0 : do ic=1,pawfgrtab(iatom)%nfgd
1654 : 0 : jc=pawfgrtab(iatom)%ifftsph(ic)
1655 : 0 : grnhat12(1,jc,isploop,1,idat2,idat1)=grnhat12(1,jc,isploop,1,idat2,idat1)+cpf_ql(1)*pawfgrtab(iatom)%gylmgr(1,ic,ilslm)
1656 : 0 : grnhat12(1,jc,isploop,2,idat2,idat1)=grnhat12(1,jc,isploop,2,idat2,idat1)+cpf_ql(1)*pawfgrtab(iatom)%gylmgr(2,ic,ilslm)
1657 : 0 : grnhat12(1,jc,isploop,3,idat2,idat1)=grnhat12(1,jc,isploop,3,idat2,idat1)+cpf_ql(1)*pawfgrtab(iatom)%gylmgr(3,ic,ilslm)
1658 : :
1659 : 0 : grnhat12(2,jc,isploop,1,idat2,idat1)=grnhat12(2,jc,isploop,1,idat2,idat1)+cpf_ql(2)*pawfgrtab(iatom)%gylmgr(1,ic,ilslm)
1660 : 0 : grnhat12(2,jc,isploop,2,idat2,idat1)=grnhat12(2,jc,isploop,2,idat2,idat1)+cpf_ql(2)*pawfgrtab(iatom)%gylmgr(2,ic,ilslm)
1661 : 0 : grnhat12(2,jc,isploop,3,idat2,idat1)=grnhat12(2,jc,isploop,3,idat2,idat1)+cpf_ql(2)*pawfgrtab(iatom)%gylmgr(3,ic,ilslm)
1662 : : end do
1663 : : end if
1664 : : end do
1665 : : end do
1666 : : end if ! compute_grad
1667 [ # # ]: 0 : if (compute_grad1) then
1668 [ # # # # ]: 0 : do ils=lmin,lmax,2 ! Sum over (L,M)
1669 [ # # ]: 0 : do mm=-ils,ils
1670 : 0 : ilslm=ils*ils+ils+mm+1
1671 [ # # ]: 0 : if (pawang%gntselect(ilslm,klm)>0) then
1672 : 0 : cpf_ql(1)=cpf(1)*qijl(ilslm,klmn)
1673 : 0 : cpf_ql(2)=cpf(2)*qijl(ilslm,klmn)
1674 [ # # ]: 0 : do ic=1,pawfgrtab(iatom)%nfgd
1675 : 0 : grnhat_12(1,ic,isploop,1,iatom,idat2,idat1)=grnhat_12(1,ic,isploop,1,iatom,idat2,idat1)+cpf_ql(1)*pawfgrtab(iatom)%gylmgr(1,ic,ilslm)
1676 : 0 : grnhat_12(1,ic,isploop,2,iatom,idat2,idat1)=grnhat_12(1,ic,isploop,2,iatom,idat2,idat1)+cpf_ql(1)*pawfgrtab(iatom)%gylmgr(2,ic,ilslm)
1677 : 0 : grnhat_12(1,ic,isploop,3,iatom,idat2,idat1)=grnhat_12(1,ic,isploop,3,iatom,idat2,idat1)+cpf_ql(1)*pawfgrtab(iatom)%gylmgr(3,ic,ilslm)
1678 : :
1679 : 0 : grnhat_12(2,ic,isploop,1,iatom,idat2,idat1)=grnhat_12(2,ic,isploop,1,iatom,idat2,idat1)+cpf_ql(2)*pawfgrtab(iatom)%gylmgr(1,ic,ilslm)
1680 : 0 : grnhat_12(2,ic,isploop,2,iatom,idat2,idat1)=grnhat_12(2,ic,isploop,2,iatom,idat2,idat1)+cpf_ql(2)*pawfgrtab(iatom)%gylmgr(2,ic,ilslm)
1681 : 0 : grnhat_12(2,ic,isploop,3,iatom,idat2,idat1)=grnhat_12(2,ic,isploop,3,iatom,idat2,idat1)+cpf_ql(2)*pawfgrtab(iatom)%gylmgr(3,ic,ilslm)
1682 : : end do
1683 : : end if
1684 : : end do
1685 : : end do
1686 : : end if ! compute_grad1
1687 : : end do ! klmn (ij channels)
1688 : : ! If needed, multiply eventually by exp(-i.q.r) phase
1689 [ # # ]: 0 : if (compute_nhat) then
1690 [ # # # # ]: 0 : if(compute_phonon.and.(.not.qeq0).and.pawfgrtab(iatom)%expiqr_allocated/=0) then
1691 : : !$OMP PARALLEL DO PRIVATE(ro,ro_ql,ic,jc)
1692 [ # # ]: 0 : do ic=1,pawfgrtab(iatom)%nfgd
1693 : 0 : jc=pawfgrtab(iatom)%ifftsph(ic)
1694 : 0 : ro_ql(1)= pawfgrtab(iatom)%expiqr(1,ic)
1695 : 0 : ro_ql(2)= pawfgrtab(iatom)%expiqr(2,ic)
1696 [ # # ]: 0 : ro(1:2)=nhat12_atm(1:2,jc,isploop,idat2,idat1)
1697 : 0 : nhat12_atm(1,jc,isploop,idat2,idat1)=ro(1)*ro_ql(1)-ro(2)*ro_ql(2)
1698 : 0 : nhat12_atm(2,jc,isploop,idat2,idat1)=ro(2)*ro_ql(1)+ro(1)*ro_ql(2)
1699 : : end do
1700 : : end if
1701 : : end if
1702 : :
1703 [ # # ]: 0 : if (compute_grad) then
1704 [ # # # # ]: 0 : if(compute_phonon.and.(.not.qeq0).and.pawfgrtab(iatom)%expiqr_allocated/=0) then
1705 : : !$OMP PARALLEL DO PRIVATE(ro,ro_ql,ic,jc)
1706 [ # # ]: 0 : do ic=1,pawfgrtab(iatom)%nfgd
1707 : 0 : jc=pawfgrtab(iatom)%ifftsph(ic)
1708 : 0 : ro_ql(1)= pawfgrtab(iatom)%expiqr(1,ic)
1709 : 0 : ro_ql(2)= pawfgrtab(iatom)%expiqr(2,ic)
1710 : 0 : ro(1)=grnhat12(1,jc,isploop,1,idat2,idat1)-qphon(1)*nhat12_atm(2,jc,isploop,idat2,idat1)
1711 : 0 : ro(2)=grnhat12(2,jc,isploop,1,idat2,idat1)+qphon(1)*nhat12_atm(1,jc,isploop,idat2,idat1)
1712 : 0 : grnhat12(1,jc,isploop,1,idat2,idat1)=ro(1)*ro_ql(1)-ro(2)*ro_ql(2)
1713 : 0 : grnhat12(2,jc,isploop,1,idat2,idat1)=ro(2)*ro_ql(1)+ro(1)*ro_ql(2)
1714 : 0 : ro(1)=grnhat12(1,jc,isploop,2,idat2,idat1)-qphon(2)*nhat12_atm(2,jc,isploop,idat2,idat1)
1715 : 0 : ro(2)=grnhat12(2,jc,isploop,2,idat2,idat1)+qphon(2)*nhat12_atm(1,jc,isploop,idat2,idat1)
1716 : 0 : grnhat12(1,jc,isploop,2,idat2,idat1)=ro(1)*ro_ql(1)-ro(2)*ro_ql(2)
1717 : 0 : grnhat12(2,jc,isploop,2,idat2,idat1)=ro(2)*ro_ql(1)+ro(1)*ro_ql(2)
1718 : 0 : ro(1)=grnhat12(1,jc,isploop,3,idat2,idat1)-qphon(3)*nhat12_atm(2,jc,isploop,idat2,idat1)
1719 : 0 : ro(2)=grnhat12(2,jc,isploop,3,idat2,idat1)+qphon(3)*nhat12_atm(1,jc,isploop,idat2,idat1)
1720 : 0 : grnhat12(1,jc,isploop,3,idat2,idat1)=ro(1)*ro_ql(1)-ro(2)*ro_ql(2)
1721 : 0 : grnhat12(2,jc,isploop,3,idat2,idat1)=ro(2)*ro_ql(1)+ro(1)*ro_ql(2)
1722 : : end do
1723 : : end if
1724 : : end if
1725 [ # # ]: 0 : if (compute_grad1) then
1726 [ # # # # ]: 0 : if(compute_phonon.and.(.not.qeq0).and.pawfgrtab(iatom)%expiqr_allocated/=0) then
1727 : : !$OMP PARALLEL DO PRIVATE(ro,ro_ql,ic)
1728 [ # # ]: 0 : do ic=1,pawfgrtab(iatom)%nfgd
1729 : 0 : ro_ql(1)= pawfgrtab(iatom)%expiqr(1,ic)
1730 : 0 : ro_ql(2)= pawfgrtab(iatom)%expiqr(2,ic)
1731 : 0 : ro(1)=grnhat_12(1,ic,isploop,1,iatom,idat2,idat1)
1732 : 0 : ro(2)=grnhat_12(2,ic,isploop,1,iatom,idat2,idat1)
1733 : 0 : grnhat_12(1,ic,isploop,1,iatom,idat2,idat1)=ro(1)*ro_ql(1)-ro(2)*ro_ql(2)
1734 : 0 : grnhat_12(2,ic,isploop,1,iatom,idat2,idat1)=ro(2)*ro_ql(1)+ro(1)*ro_ql(2)
1735 : 0 : ro(1)=grnhat_12(1,ic,isploop,2,iatom,idat2,idat1)
1736 : 0 : ro(2)=grnhat_12(2,ic,isploop,2,iatom,idat2,idat1)
1737 : 0 : grnhat_12(1,ic,isploop,2,iatom,idat2,idat1)=ro(1)*ro_ql(1)-ro(2)*ro_ql(2)
1738 : 0 : grnhat_12(2,ic,isploop,2,iatom,idat2,idat1)=ro(2)*ro_ql(1)+ro(1)*ro_ql(2)
1739 : 0 : ro(1)=grnhat_12(1,ic,isploop,3,iatom,idat2,idat1)
1740 : 0 : ro(2)=grnhat_12(2,ic,isploop,3,iatom,idat2,idat1)
1741 : 0 : grnhat_12(1,ic,isploop,3,iatom,idat2,idat1)=ro(1)*ro_ql(1)-ro(2)*ro_ql(2)
1742 : 0 : grnhat_12(2,ic,isploop,3,iatom,idat2,idat1)=ro(2)*ro_ql(1)+ro(1)*ro_ql(2)
1743 : : end do
1744 : : end if
1745 : : end if
1746 : : end do ! isploop (density components of the compensation charge)
1747 : : end do ! idat2
1748 : : end do ! idat1
1749 : : ! accumlate nhat12 for all the atoms
1750 : : !nhat12(2,nfft,nspinor**2,ndat2)
1751 [ # # # # : 0 : if (compute_nhat) nhat12=nhat12+nhat12_atm
# # # # #
# # # ]
1752 : :
1753 [ # # ]: 0 : if (pawfgrtab(iatom)%gylm_allocated==2) then
1754 [ # # ]: 0 : ABI_FREE(pawfgrtab(iatom)%gylm)
1755 [ # # ]: 0 : ABI_MALLOC(pawfgrtab(iatom)%gylm,(0,0))
1756 : 0 : pawfgrtab(iatom)%gylm_allocated=0
1757 : : end if
1758 [ # # ]: 0 : if (pawfgrtab(iatom)%gylmgr_allocated==2) then
1759 [ # # ]: 0 : ABI_FREE(pawfgrtab(iatom)%gylmgr)
1760 [ # # ]: 0 : ABI_MALLOC(pawfgrtab(iatom)%gylmgr,(0,0,0))
1761 : 0 : pawfgrtab(iatom)%gylmgr_allocated=0
1762 : : end if
1763 [ # # ]: 0 : ABI_FREE(qijl)
1764 [ # # ]: 0 : ABI_FREE(nhat12_atm)
1765 [ # # ]: 0 : if (pawfgrtab(iatom)%expiqr_allocated==2) then
1766 [ # # ]: 0 : ABI_FREE(pawfgrtab(iatom)%expiqr)
1767 [ # # ]: 0 : ABI_MALLOC(pawfgrtab(iatom)%expiqr,(0,0))
1768 : 0 : pawfgrtab(iatom)%expiqr_allocated=0
1769 : : end if
1770 : :
1771 : : end do ! iatom
1772 : :
1773 [ # # # # : 0 : if (compute_grad1) grnhat_12=-grnhat_12
# # # # #
# # # # #
# # ]
1774 : :
1775 : : !----- Reduction in case of parallelism -----!
1776 [ # # ]: 0 : if (paral_atom)then
1777 : 0 : call timab(48,1,tsec)
1778 [ # # ]: 0 : if (compute_nhat) then
1779 : 0 : call xmpi_sum(nhat12,my_comm_atom,ierr)
1780 : : end if
1781 [ # # ]: 0 : if (compute_grad) then
1782 : 0 : call xmpi_sum(grnhat12,my_comm_atom,ierr)
1783 : : end if
1784 [ # # ]: 0 : if (compute_grad1) then
1785 : 0 : call xmpi_sum(grnhat_12,my_comm_atom,ierr)
1786 : : end if
1787 : 0 : call timab(48,2,tsec)
1788 : : end if
1789 : :
1790 : : !----- Avoid unbalanced g-components numerical errors -----!
1791 : :
1792 [ # # # # ]: 0 : if (izero==1.and.compute_nhat) then
1793 : : ! Create fake mpi_enreg to wrap fourdp
1794 [ # # ]: 0 : if (present(distribfft)) then
1795 : : my_distribfft => distribfft
1796 : : else
1797 [ # # ]: 0 : ABI_MALLOC(my_distribfft,)
1798 : 0 : call my_distribfft%init_seq('f',ngfft(2),ngfft(3),'fourdp')
1799 : : end if
1800 : 0 : call initmpi_seq(mpi_enreg_fft)
1801 [ # # # # : 0 : ABI_FREE(mpi_enreg_fft%distribfft)
# # # # #
# # # # #
# # # # #
# # # # #
# # ]
1802 [ # # ]: 0 : if (present(comm_fft)) then
1803 : 0 : call set_mpi_enreg_fft(mpi_enreg_fft,comm_fft,my_distribfft,me_g0,paral_kgb)
1804 : 0 : my_comm_fft=comm_fft;paral_kgb_fft=paral_kgb
1805 : : else
1806 : 0 : my_comm_fft=xmpi_comm_self;paral_kgb_fft=0;
1807 : 0 : mpi_enreg_fft%distribfft => my_distribfft
1808 : : end if
1809 : : ! Do FFT
1810 [ # # # # : 0 : ABI_MALLOC(work,(2,nfft))
# # # # #
# # # ]
1811 : 0 : cplex=2
1812 [ # # ]: 0 : do idat1=1,ndat1
1813 [ # # ]: 0 : do idat2=1,ndat2
1814 [ # # ]: 0 : do isp1=1,MIN(2,nspinor**2)
1815 : 0 : call fourdp(cplex,work,nhat12(:,:,isp1,idat2,idat1),-1,mpi_enreg_fft,nfft,1,ngfft,0)
1816 : 0 : call zerosym(work,cplex,ngfft(1),ngfft(2),ngfft(3),comm_fft=my_comm_fft,distribfft=my_distribfft)
1817 : 0 : call fourdp(cplex,work,nhat12(:,:,isp1,idat2,idat1),+1,mpi_enreg_fft,nfft,1,ngfft,0)
1818 : : end do
1819 : : end do ! idat2
1820 : : end do ! idat1
1821 [ # # ]: 0 : ABI_FREE(work)
1822 : : ! Destroy fake mpi_enreg
1823 : 0 : call unset_mpi_enreg_fft(mpi_enreg_fft)
1824 [ # # ]: 0 : if (.not.present(distribfft)) then
1825 : 0 : call my_distribfft%free()
1826 [ # # # # : 0 : ABI_FREE(my_distribfft)
# # # # #
# # # # #
# # # # #
# # # #
# ]
1827 : : end if
1828 : : end if
1829 : :
1830 : : !Destroy atom table used for parallelism
1831 : 0 : call free_my_atmtab(my_atmtab,my_atmtab_allocated)
1832 : :
1833 : : DBG_EXIT("COLL")
1834 : :
1835 [ # # # # : 0 : end subroutine pawmknhat_psipsi
# # # # #
# # # # #
# # # # #
# # # #
# ]
1836 : : !!***
1837 : :
1838 : : !----------------------------------------------------------------------
1839 : :
1840 : : !!****f* m_paw_nhat/pawnhatfr
1841 : : !!
1842 : : !! NAME
1843 : : !! pawnhatfr
1844 : : !!
1845 : : !! FUNCTION
1846 : : !! PAW: Compute frozen part of 1st-order compensation charge density nhat^(1)
1847 : : !! nhatfr(r)=Sum_ij,lm[rhoij_ij.q_ij^l.(g_l(r).Y_lm(r))^(1)]
1848 : : !! Depends on q wave vector but not on first-order wave-function.
1849 : : !!
1850 : : !! INPUTS
1851 : : !! ider=0: computes frozen part of compensation density
1852 : : !! 1: computes frozen part of compensation density and cartesian gradients
1853 : : !! idir=direction of atomic displacement (in case of phonons perturb.)
1854 : : !! ipert=nindex of perturbation
1855 : : !! mpi_atmtab(:)=--optional-- indexes of the atoms treated by current proc
1856 : : !! comm_atom=--optional-- MPI communicator over atoms
1857 : : !! my_natom=number of atoms treated by current processor
1858 : : !! natom=total number of atoms in cell
1859 : : !! nspden=number of spin-density components
1860 : : !! ntypat=number of types of atoms
1861 : : !! pawang <type(pawang_type)>=paw angular mesh and related data
1862 : : !! pawfgrtab(my_natom) <type(pawfgrtab_type)>=atomic data given on fine rectangular grid
1863 : : !! pawrhoij(my_natom) <type(pawrhoij_type)>= Ground-State paw rhoij occupancies and related data
1864 : : !! pawtab(ntypat) <type(pawtab_type)>=paw tabulated starting data
1865 : : !! rprimd(3,3)=dimensional primitive translations for real space
1866 : : !!
1867 : : !! OUTPUT
1868 : : !! pawfgrtab(iatom)%nhatfr(nfgd,nspden)
1869 : : !! frozen part of charge compensation density (inside PAW spheres)
1870 : : !! =Sum_ij,lm[rhoij_ij.q_ij^l.(g_l(r).Y_lm(r))^(1)]
1871 : : !! === If ider==1
1872 : : !! pawfgrtab(iatom)%nhatfrgr(3,nfgd,nspden)
1873 : : !! gradients of frozen part of charge compensation density (inside PAW spheres)
1874 : : !! =Sum_ij,lm[rhoij_ij.q_ij^l . d/dr((g_l(r).Y_lm(r))^(1))]
1875 : : !!
1876 : : !! SOURCE
1877 : :
1878 : 3492 : subroutine pawnhatfr(ider,idir,ipert,my_natom,natom,nspden,ntypat,&
1879 : 3492 : & pawang,pawfgrtab,pawrhoij,pawtab,rprimd, &
1880 [ + + + + : 3492 : & mpi_atmtab,comm_atom) ! optional arguments (parallelism)
- + ]
1881 : :
1882 : : !Arguments ------------------------------------
1883 : : !scalars
1884 : : integer,intent(in) :: ider,idir,ipert,my_natom,natom,nspden,ntypat
1885 : : integer,optional,intent(in) :: comm_atom
1886 : : type(pawang_type),intent(in) :: pawang
1887 : : !arrays
1888 : : integer,optional,target,intent(in) :: mpi_atmtab(:)
1889 : : real(dp),intent(in) :: rprimd(3,3)
1890 : : type(pawfgrtab_type),intent(inout) :: pawfgrtab(my_natom)
1891 : : type(pawrhoij_type),intent(in) :: pawrhoij(my_natom)
1892 : : type(pawtab_type),intent(in) :: pawtab(ntypat)
1893 : :
1894 : : !Local variables-------------------------------
1895 : : !scalars
1896 : : integer :: iatom,iatom_tot,ic,ils,ilslm,irhoij,isel,ispden,istr,itypat,jrhoij
1897 : : integer :: klm,klmn,lm_size,lmn2_size,lm0,lmax,lmin,mua,mub,mm,mu,my_comm_atom,nfgd,nu,optgr0,optgr1,optgr2
1898 : : logical :: my_atmtab_allocated,my_pert,paral_atom
1899 : : real(dp) :: contrib,ro
1900 : : !arrays
1901 : : integer,parameter :: voigt(3,3)=reshape((/1,6,5,6,2,4,5,4,3/),(/3,3/))
1902 : : integer,parameter :: alpha(9)=(/1,2,3,3,3,2,2,1,1/),beta(9)=(/1,2,3,2,1,1,3,3,2/)
1903 : 3492 : integer,pointer :: my_atmtab(:)
1904 : 3492 : real(dp),allocatable :: nhatfr_tmp(:,:),nhatfrgr_tmp(:,:,:)
1905 : :
1906 : : ! *************************************************************************
1907 : :
1908 : : DBG_ENTER("COLL")
1909 : :
1910 : : !Only relevant for atomic displacement and strain perturbation
1911 [ + + + + ]: 856 : if (ipert>natom.and.ipert/=natom+3.and.ipert/=natom+4) return
1912 : :
1913 : : !Set up parallelism over atoms
1914 [ + + + + ]: 3332 : paral_atom=(present(comm_atom).and.(my_natom/=natom))
1915 [ + + ]: 3332 : nullify(my_atmtab);if (present(mpi_atmtab)) my_atmtab => mpi_atmtab
1916 [ + + ]: 3332 : my_comm_atom=xmpi_comm_self;if (present(comm_atom)) my_comm_atom=comm_atom
1917 : 3332 : call get_my_atmtab(my_comm_atom,my_atmtab,my_atmtab_allocated,paral_atom,natom,my_natom_ref=my_natom)
1918 : :
1919 : : !Compatibility tests
1920 [ + + ]: 3332 : if (my_natom>0) then
1921 [ + - - + : 3308 : if ((pawfgrtab(1)%gylm_allocated==0.or.pawfgrtab(1)%gylmgr_allocated==0).and. &
- - ]
1922 : : & pawfgrtab(1)%rfgd_allocated==0) then
1923 : 0 : ABI_BUG('pawnhatfr: pawfgrtab()%rfgd array must be allocated!')
1924 : : end if
1925 [ - + ]: 3308 : if (pawrhoij(1)%qphase/=1) then
1926 : 0 : ABI_BUG('pawnhatfr: not supposed to be called with qphase=2!')
1927 : : end if
1928 : : end if
1929 : :
1930 [ + + + + : 3332 : my_pert = (ipert<=natom).or.ipert==natom+3.or.ipert==natom+4
+ - ]
1931 : :
1932 : : !Get correct index of strain pertubation
1933 [ + + ]: 3332 : if (ipert==natom+3) istr = idir
1934 [ + + ]: 3332 : if (ipert==natom+4) istr = idir + 3
1935 : :
1936 : : !Loops over atoms
1937 [ + + ]: 10188 : do iatom=1,my_natom
1938 : 6856 : iatom_tot=iatom;if (paral_atom) iatom_tot=my_atmtab(iatom)
1939 : :
1940 : : ! Eventually allocate frozen nhat points
1941 [ - + ]: 6856 : if (my_pert) then
1942 [ + + ]: 6856 : if (pawfgrtab(iatom)%nhatfr_allocated==0) then
1943 [ + + ]: 2710 : if (allocated(pawfgrtab(iatom)%nhatfr)) then
1944 : 391 : ABI_FREE(pawfgrtab(iatom)%nhatfr)
1945 : : end if
1946 [ - + + - : 10840 : ABI_MALLOC(pawfgrtab(iatom)%nhatfr,(pawfgrtab(iatom)%nfgd,nspden))
- + + - +
- + - - +
- + - + ]
1947 : 2710 : pawfgrtab(iatom)%nhatfr_allocated=1
1948 : : end if
1949 [ + + + + ]: 6856 : if (ider==1.and.pawfgrtab(iatom)%nhatfrgr_allocated==0) then
1950 [ + - ]: 1 : if (allocated(pawfgrtab(iatom)%nhatfrgr)) then
1951 : 1 : ABI_FREE(pawfgrtab(iatom)%nhatfrgr)
1952 : : end if
1953 [ - + + - : 4 : ABI_MALLOC(pawfgrtab(iatom)%nhatfrgr,(3,pawfgrtab(iatom)%nfgd,nspden))
- + + - +
- + - - +
- + - + ]
1954 : 1 : pawfgrtab(iatom)%nhatfrgr_allocated=1
1955 : : end if
1956 : : end if
1957 : :
1958 : : ! Select if frozen part of nhat exists for the current perturbation
1959 [ - + ]: 6856 : if ((.not.my_pert).or.(pawfgrtab(iatom)%nhatfr_allocated==0)) cycle
1960 : :
1961 : : ! Some atom-dependent quantities
1962 : 6856 : itypat=pawfgrtab(iatom)%itypat
1963 : 6856 : lm_size=pawfgrtab(iatom)%l_size**2
1964 : 6856 : lmn2_size=pawtab(itypat)%lmn2_size
1965 : :
1966 : : ! Eventually compute g_l(r).Y_lm(r) factors for the current atom (if not already done)
1967 : 6856 : nfgd=pawfgrtab(iatom)%nfgd
1968 [ + - + + : 6856 : if ((pawfgrtab(iatom)%gylmgr_allocated==0).or. &
- + ]
1969 : : & (pawfgrtab(iatom)%gylmgr2_allocated==0.and.ider==1)) then
1970 : 0 : optgr0=0;optgr1=0;optgr2=0
1971 [ # # # # ]: 0 : if(ipert==natom+3.or.ipert==natom+4)then
1972 [ # # ]: 0 : if (pawfgrtab(iatom)%gylm_allocated==0) then
1973 [ # # ]: 0 : if (allocated(pawfgrtab(iatom)%gylm)) then
1974 : 0 : ABI_FREE(pawfgrtab(iatom)%gylm)
1975 : : end if
1976 [ # # # # : 0 : ABI_MALLOC(pawfgrtab(iatom)%gylm,(nfgd,lm_size))
# # # # #
# # # # #
# # # # ]
1977 : 0 : pawfgrtab(iatom)%gylm_allocated=2;optgr0=1
1978 : : end if
1979 : : end if
1980 [ # # ]: 0 : if (pawfgrtab(iatom)%gylmgr_allocated==0) then
1981 [ # # ]: 0 : if (allocated(pawfgrtab(iatom)%gylmgr)) then
1982 : 0 : ABI_FREE(pawfgrtab(iatom)%gylmgr)
1983 : : end if
1984 [ # # # # : 0 : ABI_MALLOC(pawfgrtab(iatom)%gylmgr,(3,nfgd,lm_size))
# # # # #
# # # # #
# # # # ]
1985 : 0 : pawfgrtab(iatom)%gylmgr_allocated=2;optgr1=1
1986 : : end if
1987 [ # # # # ]: 0 : if (ider==1.and.pawfgrtab(iatom)%gylmgr2_allocated==0) then
1988 [ # # ]: 0 : if (allocated(pawfgrtab(iatom)%gylmgr2)) then
1989 : 0 : ABI_FREE(pawfgrtab(iatom)%gylmgr2)
1990 : : end if
1991 [ # # # # : 0 : ABI_MALLOC(pawfgrtab(iatom)%gylmgr2,(6,nfgd,lm_size))
# # # # #
# # # # #
# # # # ]
1992 : 0 : pawfgrtab(iatom)%gylmgr2_allocated=2;optgr2=1
1993 : : end if
1994 [ # # ]: 0 : if (optgr0+optgr1+optgr2>0) then
1995 : : call pawgylm(pawfgrtab(iatom)%gylm,pawfgrtab(iatom)%gylmgr,pawfgrtab(iatom)%gylmgr2,&
1996 : 0 : & lm_size,nfgd,optgr0,optgr1,optgr2,pawtab(itypat),pawfgrtab(iatom)%rfgd)
1997 : : end if
1998 : : end if
1999 : :
2000 : :
2001 : : ! ============ Phonons ====================================
2002 [ + + ]: 6856 : if (ipert<=natom) then
2003 : :
2004 : : ! Loop over spin components
2005 [ + + ]: 11450 : do ispden=1,nspden
2006 : :
2007 [ - + + - : 18102 : ABI_MALLOC(nhatfr_tmp,(3,nfgd))
+ - + - -
+ - + ]
2008 [ + + + + ]: 20623774 : nhatfr_tmp=zero
2009 [ + + ]: 6034 : if (ider==1) then
2010 [ - + + - : 16 : ABI_MALLOC(nhatfrgr_tmp,(3,nfgd,3))
+ - + - +
- - + - +
- + ]
2011 [ + + + + : 177376 : nhatfrgr_tmp=zero
+ + ]
2012 : : end if
2013 : :
2014 : 6034 : jrhoij=1
2015 [ + + ]: 135062 : do irhoij=1,pawrhoij(iatom)%nrhoijsel
2016 : 129028 : klmn=pawrhoij(iatom)%rhoijselect(irhoij)
2017 : 129028 : klm =pawtab(itypat)%indklmn(1,klmn)
2018 : 129028 : lmin=pawtab(itypat)%indklmn(3,klmn)
2019 : 129028 : lmax=pawtab(itypat)%indklmn(4,klmn)
2020 : :
2021 [ + + ]: 129028 : if (nspden/=2) then
2022 : 88632 : ro=pawrhoij(iatom)%rhoijp(jrhoij,ispden)
2023 : : else
2024 [ + + ]: 40396 : if (ispden==1) then
2025 : 20198 : ro=pawrhoij(iatom)%rhoijp(jrhoij,1)+pawrhoij(iatom)%rhoijp(jrhoij,2)
2026 [ + - ]: 20198 : else if (ispden==2) then
2027 : 20198 : ro=pawrhoij(iatom)%rhoijp(jrhoij,1)
2028 : : end if
2029 : : end if
2030 : 129028 : ro=pawtab(itypat)%dltij(klmn)*ro
2031 : :
2032 [ - + + + ]: 336078 : do ils=lmin,lmax,2
2033 : 207050 : lm0=ils**2+ils+1
2034 [ + + ]: 928960 : do mm=-ils,ils
2035 : 592882 : ilslm=lm0+mm;isel=pawang%gntselect(lm0+mm,klm)
2036 [ + + ]: 799932 : if (isel>0) then
2037 [ + + ]: 235708531 : do ic=1,nfgd
2038 [ + + ]: 942164356 : do mu=1,3
2039 : : contrib=-ro*pawtab(itypat)%qijl(ilslm,klmn)&
2040 : 706455825 : & *pawfgrtab(iatom)%gylmgr(mu,ic,ilslm)
2041 : 941941100 : nhatfr_tmp(mu,ic)=nhatfr_tmp(mu,ic)+contrib
2042 : : end do
2043 : : end do
2044 [ + + ]: 223256 : if (ider==1) then
2045 [ + + ]: 1330560 : do ic=1,nfgd
2046 [ + + ]: 5321160 : do nu=1,3
2047 [ + + ]: 17292600 : do mu=1,3
2048 : : contrib=-ro*pawtab(itypat)%qijl(ilslm,klmn) &
2049 : 11971800 : & *pawfgrtab(iatom)%gylmgr2(voigt(mu,nu),ic,ilslm)
2050 : 15962400 : nhatfrgr_tmp(mu,ic,nu)=nhatfrgr_tmp(mu,ic,nu)+contrib
2051 : : end do
2052 : : end do
2053 : : end do
2054 : : end if
2055 : : end if
2056 : : end do
2057 : : end do
2058 : 135062 : jrhoij=jrhoij+pawrhoij(iatom)%cplex_rhoij
2059 : : end do
2060 : :
2061 : : ! Convert from cartesian to reduced coordinates
2062 [ + + ]: 5160469 : do ic=1,nfgd
2063 : : pawfgrtab(iatom)%nhatfr(ic,ispden)= &
2064 : : & rprimd(1,idir)*nhatfr_tmp(1,ic) &
2065 : : & +rprimd(2,idir)*nhatfr_tmp(2,ic) &
2066 : 5160469 : & +rprimd(3,idir)*nhatfr_tmp(3,ic)
2067 : : end do
2068 [ + + ]: 6034 : if (ider==1) then
2069 [ + + ]: 16 : do nu=1,3
2070 [ + + ]: 44356 : do ic=1,nfgd
2071 : : pawfgrtab(iatom)%nhatfrgr(nu,ic,ispden)= &
2072 : : & rprimd(1,idir)*nhatfrgr_tmp(1,ic,nu) &
2073 : : & +rprimd(2,idir)*nhatfrgr_tmp(2,ic,nu) &
2074 : 44352 : & +rprimd(3,idir)*nhatfrgr_tmp(3,ic,nu)
2075 : : end do
2076 : : end do
2077 : : end if
2078 [ - + ]: 6034 : ABI_FREE(nhatfr_tmp)
2079 [ + + ]: 11450 : if (ider==1) then
2080 [ - + ]: 4 : ABI_FREE(nhatfrgr_tmp)
2081 : : end if
2082 : : ! End loop over spin components
2083 : : end do ! ispden
2084 : :
2085 : :
2086 : : ! ============ Elastic tensor ===============================
2087 [ + + - + ]: 1440 : else if (ipert==natom+3.or.ipert==natom+4) then
2088 : : ! Loop over spin components
2089 [ + + + + ]: 1389678 : pawfgrtab(iatom)%nhatfr(:,:) = zero
2090 [ + + ]: 2952 : do ispden=1,nspden
2091 : 1512 : jrhoij=1
2092 [ + + ]: 27894 : do irhoij=1,pawrhoij(iatom)%nrhoijsel
2093 : 24942 : klmn=pawrhoij(iatom)%rhoijselect(irhoij)
2094 : 24942 : klm =pawtab(itypat)%indklmn(1,klmn)
2095 : 24942 : lmin=pawtab(itypat)%indklmn(3,klmn)
2096 : 24942 : lmax=pawtab(itypat)%indklmn(4,klmn)
2097 [ + + ]: 24942 : if (nspden/=2) then
2098 : 21486 : ro=pawrhoij(iatom)%rhoijp(jrhoij,ispden)
2099 : : else
2100 [ + + ]: 3456 : if (ispden==1) then
2101 : 1728 : ro=pawrhoij(iatom)%rhoijp(jrhoij,1)+pawrhoij(iatom)%rhoijp(jrhoij,2)
2102 [ + - ]: 1728 : else if (ispden==2) then
2103 : 1728 : ro=pawrhoij(iatom)%rhoijp(jrhoij,1)
2104 : : end if
2105 : : end if
2106 : 24942 : ro=pawtab(itypat)%dltij(klmn)*ro
2107 [ - + + + ]: 64770 : do ils=lmin,lmax,2
2108 : 39828 : lm0=ils**2+ils+1
2109 [ + + ]: 176190 : do mm=-ils,ils
2110 : 111420 : ilslm=lm0+mm;isel=pawang%gntselect(lm0+mm,klm)
2111 [ + + ]: 151248 : if (isel>0) then
2112 : : ! Sum{[Q_ij_q^LM^(1)]}
2113 [ + + ]: 56473602 : do ic=1,nfgd
2114 : 56425314 : mua=alpha(istr);mub=beta(istr)
2115 : : pawfgrtab(iatom)%nhatfr(ic,ispden) = pawfgrtab(iatom)%nhatfr(ic,ispden)+&
2116 : : & ro*pawtab(itypat)%qijl(ilslm,klmn)*half*(&
2117 : : & pawfgrtab(iatom)%gylmgr(mua,ic,ilslm)*pawfgrtab(iatom)%rfgd(mub,ic)&
2118 : 56473602 : & +pawfgrtab(iatom)%gylmgr(mub,ic,ilslm)*pawfgrtab(iatom)%rfgd(mua,ic))
2119 : : end do
2120 : : ! Add volume contribution
2121 [ + + ]: 48288 : if(istr<=3)then
2122 [ + + ]: 28236801 : do ic=1,nfgd
2123 : : pawfgrtab(iatom)%nhatfr(ic,ispden) = pawfgrtab(iatom)%nhatfr(ic,ispden)+&
2124 : 28236801 : & ro*pawtab(itypat)%qijl(ilslm,klmn)*pawfgrtab(iatom)%gylm(ic,ilslm)
2125 : : end do
2126 : : end if
2127 [ - + ]: 48288 : if (ider==1) then
2128 : 0 : ABI_ERROR("nhatgr not implemented for strain perturbationxs")
2129 : : ! do ic=1,nfgd
2130 : : ! do nu=1,6
2131 : : ! do mu=1,6
2132 : : ! contrib=-ro*pawtab(itypat)%qijl(ilslm,klmn) &
2133 : : !& *pawfgrtab(iatom)%gylmgr2(voigt(mu,nu),ic,ilslm)
2134 : : ! nhatfrgr_tmp(mu,ic,nu)=nhatfrgr_tmp(mu,ic,nu)+contrib
2135 : : ! end do
2136 : : ! end do
2137 : : ! end do
2138 : : end if
2139 : : end if
2140 : : end do
2141 : : end do
2142 : 26454 : jrhoij=jrhoij+pawrhoij(iatom)%cplex_rhoij
2143 : : end do
2144 : : end do ! ispden
2145 : : end if
2146 : :
2147 : : ! Eventually free temporary space for g_l(r).Y_lm(r) gradients and exp(-i.q.r)
2148 [ - + ]: 6856 : if (pawfgrtab(iatom)%gylmgr_allocated==2) then
2149 [ # # ]: 0 : ABI_FREE(pawfgrtab(iatom)%gylmgr)
2150 [ # # ]: 0 : ABI_MALLOC(pawfgrtab(iatom)%gylmgr,(0,0,0))
2151 : 0 : pawfgrtab(iatom)%gylmgr_allocated=0
2152 : : end if
2153 [ + - ]: 10188 : if (pawfgrtab(iatom)%gylmgr2_allocated==2) then
2154 [ # # ]: 0 : ABI_FREE(pawfgrtab(iatom)%gylmgr2)
2155 [ # # ]: 0 : ABI_MALLOC(pawfgrtab(iatom)%gylmgr2,(0,0,0))
2156 : 0 : pawfgrtab(iatom)%gylmgr2_allocated=0
2157 : : end if
2158 : :
2159 : : ! End loop on atoms
2160 : : end do
2161 : :
2162 : : !Destroy atom table used for parallelism
2163 : 3332 : call free_my_atmtab(my_atmtab,my_atmtab_allocated)
2164 : :
2165 : : DBG_EXIT("COLL")
2166 : :
2167 [ + + - + ]: 6824 : end subroutine pawnhatfr
2168 : : !!***
2169 : :
2170 : : !----------------------------------------------------------------------
2171 : :
2172 : : !!****f* m_pawdij/pawdijhat_ndat
2173 : : !! NAME
2174 : : !! pawdijhat_ndat
2175 : : !!
2176 : : !! FUNCTION
2177 : : !! Compute the "hat" contribution to the PAW pseudopotential strength Dij,
2178 : : !! i.e. the compensation charge contribution (for one atom only):
2179 : : !! D_ij^hat=Intg_R [ V(r). Sum_L(Qij^L(r)). dr]
2180 : : !!
2181 : : !! INPUTS
2182 : : !! cplex_dij=2 if dij is COMPLEX (as in the spin-orbit case), 1 if dij is REAL
2183 : : !! qphase=2 if dij contains a exp(-i.q.r) phase (as in the q<>0 RF case), 1 if not
2184 : : !! gprimd(3,3)=dimensional primitive translations for reciprocal space
2185 : : !! iatom=absolute index of current atom (between 1 and natom)
2186 : : !! natom=total number of atoms
2187 : : !! ndij= number of spin components
2188 : : !! ngrid=number of points of the real space grid (FFT, WVL, ...) treated by current proc
2189 : : !! ngridtot=total number of points of the real space grid (FFT, WVL, ...)
2190 : : !! For the FFT grid, thi should be equal to ngfft1*ngfft2*ngfft3
2191 : : !! nspden=number of spin density components
2192 : : !! nsppol=number of independent spin WF components
2193 : : !! pawang <type(pawang_type)>=paw angular mesh and related data
2194 : : !! pawfgrtab<type(pawfgrtab_type)>=atomic data given on fine rectangular grid for current atom
2195 : : !! pawtab(ntypat) <type(pawtab_type)>=paw tabulated starting data, for current atom
2196 : : !! Pot(qphase*ngrid,nspden)=potential on real space grid
2197 : : !! qphon(3)=(RF calculations only) - wavevector of the phonon
2198 : : !! ucvol=unit cell volume
2199 : : !! xred(3,my_natom)= reduced atomic coordinates
2200 : : !!
2201 : : !! OUTPUT
2202 : : !! dijhat(cplex_dij*qphase*lmn2_size,ndij)= D_ij^hat terms
2203 : : !! When Dij is complex (cplex_dij=2):
2204 : : !! dij(2*i-1,:) contains the real part, dij(2*i,:) contains the imaginary part
2205 : : !! When a exp(-i.q.r) phase is included (qphase=2):
2206 : : !! dij(1:cplex_dij*lmn2_size,:)
2207 : : !! contains the real part of the phase, i.e. D_ij*cos(q.r)
2208 : : !! dij(cplex_dij*lmn2_size+1:2*cplex_dij*lmn2_size,:)
2209 : : !! contains the imaginary part of the phase, i.e. D_ij*sin(q.r)
2210 : : !!
2211 : : !! SOURCE
2212 : :
2213 [ - + ]: 10392 : subroutine pawdijhat_ndat(dijhat,cplex_dij,qphase,gprimd,iatm,&
2214 : : & natom,ndij,ngrid,ngridtot,nspden,nsppol,ndat,nattyp,&
2215 : 10392 : & pawang,pawfgrtab,pawtab,Pot,qphon,ucvol,xred,&
2216 : : & gpu_option) ! Optional argument
2217 : :
2218 : : !Arguments ---------------------------------------------
2219 : : !scalars
2220 : : integer,intent(in) :: cplex_dij,iatm,natom,ndij,nattyp
2221 : : integer,intent(in) :: ngrid,ngridtot,nspden,nsppol,ndat,qphase
2222 : : integer,intent(in),optional :: gpu_option
2223 : : real(dp),intent(in) :: ucvol
2224 : : type(pawang_type),intent(in),target :: pawang
2225 : : !arrays
2226 : : real(dp),intent(in) :: gprimd(3,3),Pot(qphase*ngrid,nspden,ndat),qphon(3),xred(3,natom)
2227 : : real(dp),intent(out),target :: dijhat(:,:,:)
2228 : : type(pawtab_type),intent(in),target :: pawtab
2229 : : type(pawfgrtab_type),intent(inout),target :: pawfgrtab(natom)
2230 : :
2231 : : !Local variables ---------------------------------------
2232 : : !scalars
2233 : : integer :: ic,idij,idijend,icpq,ispden,iatom,idat,ia,ilslm,klm
2234 : : integer :: jc,klmn,klmn1,nfgd_max,iatom_tot
2235 : : integer :: l_size_full,lm_size,lmn2_size,nfgd,nsploop,optgr0,gpu_option_
2236 : : logical :: has_qphase,qne0
2237 : : real(dp) :: vi,vr,scal_r
2238 : : complex(dp) :: gemm_alpha
2239 : : character(len=500) :: msg
2240 : : !arrays
2241 : : real(dp) :: rdum1(1),rdum2(2)
2242 : 10392 : real(dp),allocatable,target :: dijhat_idij(:,:,:,:),prod(:,:,:,:),atom_potg(:,:,:,:),gnt_scal(:,:)
2243 : 10392 : real(dp),allocatable,target :: atom_expiqr(:,:,:),atom_gylm(:,:,:),atom_qijl(:,:)
2244 : 10392 : integer,allocatable,target :: atom_ifftsph(:,:),atom_nfgd(:)
2245 : :
2246 : : ! *************************************************************************
2247 : :
2248 : : !Useful data
2249 : 10392 : lm_size = pawtab%lcut_size**2
2250 : 10392 : lmn2_size = pawtab%lmn2_size
2251 : 10392 : qne0=(qphon(1)**2+qphon(2)**2+qphon(3)**2>=1.d-15)
2252 [ + + - + ]: 10392 : has_qphase=(qne0.and.qphase==2)
2253 : 10392 : scal_r=ucvol/dble(ngridtot)
2254 [ + - ]: 10392 : gpu_option_=ABI_GPU_DISABLED; if (present(gpu_option)) gpu_option_=gpu_option
2255 : :
2256 [ - + + - : 41568 : ABI_MALLOC(gnt_scal,(size(pawang%gntselect,1),size(pawang%gntselect,2)))
- + + - +
- + - - +
+ - ]
2257 [ + + + + ]: 1049592 : gnt_scal=0
2258 [ + + ]: 114312 : do klm=1,size(pawang%gntselect,2)
2259 [ + + ]: 1049592 : do ilslm=1,size(pawang%gntselect,1)
2260 [ + + ]: 1039200 : if(pawang%gntselect(ilslm,klm)>0) gnt_scal(ilslm,klm)=1
2261 : : end do
2262 : : end do
2263 : :
2264 : 10392 : l_size_full = size(pawtab%qijl,1)
2265 [ - + + - : 41568 : ABI_MALLOC(atom_qijl,(l_size_full,lmn2_size))
- + + - +
- + - - +
- + ]
2266 [ + - + - : 3509148 : atom_qijl = pawtab%qijl
- + - - -
- - - + +
+ + ]
2267 : : !Init memory
2268 : : #ifdef HAVE_OPENMP_OFFLOAD
2269 : : !$OMP TARGET ENTER DATA MAP(alloc:dijhat) IF(gpu_option_==ABI_GPU_OPENMP)
2270 : : #endif
2271 [ + - ]: 10392 : if(gpu_option_==ABI_GPU_DISABLED) then
2272 [ + + + + : 5484564 : dijhat=zero
+ + ]
2273 [ # # ]: 0 : else if(gpu_option_==ABI_GPU_OPENMP) then
2274 [ # # # # : 0 : call gpu_set_to_zero(dijhat,int(cplex_dij,c_size_t)*qphase*lmn2_size*ndij*ndat*nattyp)
# # # # #
# # # ]
2275 : : end if
2276 : :
2277 : : !------------------------------------------------------------------------
2278 : : !----- Loop over atoms (init)
2279 : : !------------------------------------------------------------------------
2280 : 10392 : nfgd_max=1
2281 [ + + ]: 29748 : do ia=1,nattyp
2282 : 19356 : iatom=iatm+ia
2283 : 19356 : iatom_tot=iatom
2284 : :
2285 : 19356 : nfgd_max=MAX(pawfgrtab(iatom)%nfgd,nfgd_max)
2286 : 19356 : nfgd=pawfgrtab(iatom)%nfgd
2287 : :
2288 : : !Eventually compute g_l(r).Y_lm(r) factors for the current atom (if not already done)
2289 [ + - ]: 19356 : if (pawfgrtab(iatom)%gylm_allocated==0) then
2290 [ + - ]: 19356 : if (allocated(pawfgrtab(iatom)%gylm)) then
2291 : 19356 : ABI_FREE(pawfgrtab(iatom)%gylm)
2292 : : end if
2293 [ - + + - : 77424 : ABI_MALLOC(pawfgrtab(iatom)%gylm,(nfgd,lm_size))
- + + - +
- + - - +
- + - + ]
2294 : 19356 : pawfgrtab(iatom)%gylm_allocated=2;optgr0=1
2295 : 19356 : call pawgylm(pawfgrtab(iatom)%gylm,rdum1,rdum2,lm_size,nfgd,optgr0,0,0,pawtab,pawfgrtab(iatom)%rfgd)
2296 : : end if
2297 : :
2298 : : !Eventually compute exp(i.q.r) factors for the current atom (if not already done)
2299 [ + + + - ]: 29748 : if (has_qphase.and.pawfgrtab(iatom)%expiqr_allocated==0) then
2300 [ - + ]: 15642 : if (pawfgrtab(iatom)%rfgd_allocated==0) then
2301 : 0 : msg='pawfgrtab()%rfgd array must be allocated !'
2302 : 0 : ABI_BUG(msg)
2303 : : end if
2304 [ + - ]: 15642 : if (allocated(pawfgrtab(iatom)%expiqr)) then
2305 : 15642 : ABI_FREE(pawfgrtab(iatom)%expiqr)
2306 : : end if
2307 [ - + + - : 46926 : ABI_MALLOC(pawfgrtab(iatom)%expiqr,(2,nfgd))
+ - + - -
+ - + -
+ ]
2308 : 15642 : call pawexpiqr(pawfgrtab(iatom)%expiqr,gprimd,nfgd,qphon,pawfgrtab(iatom)%rfgd,xred(:,iatom))
2309 : 15642 : pawfgrtab(iatom)%expiqr_allocated=2
2310 : : end if
2311 : : end do ! ia
2312 : :
2313 [ - + + - : 31176 : ABI_MALLOC(atom_nfgd, (nattyp))
+ - + - -
+ - + ]
2314 [ - + + - : 51960 : ABI_MALLOC(atom_gylm, ( nfgd_max,lm_size,nattyp))
- + + - -
+ + - + -
+ - - + -
+ ]
2315 [ - + + - : 41568 : ABI_MALLOC(atom_ifftsph,(nfgd_max,nattyp))
- + + - +
- + - - +
- + ]
2316 [ + + ]: 10392 : if(has_qphase) then
2317 [ - + + - : 33855 : ABI_MALLOC(atom_expiqr, (2,nfgd_max,nattyp))
- + + - +
- + - - +
- + ]
2318 : : end if
2319 : :
2320 : : !Zero-pad atom_gylm beyond each atom's actual nfgd as GEMMs covers
2321 : : !nfgd max value among atoms.
2322 [ + + + + : 136677532 : atom_gylm=zero
+ + ]
2323 [ + + ]: 29748 : do ia=1,nattyp
2324 : 19356 : iatom=iatm+ia
2325 : 19356 : nfgd=pawfgrtab(iatom)%nfgd
2326 : :
2327 : 19356 : atom_nfgd(ia) = pawfgrtab(iatom)%nfgd
2328 [ + + + + ]: 136667140 : atom_gylm(1:nfgd,1:lm_size,ia) = pawfgrtab(iatom)%gylm(1:nfgd,1:lm_size)
2329 [ + + ]: 15460232 : atom_ifftsph(1:nfgd,ia) = pawfgrtab(iatom)%ifftsph(1:nfgd)
2330 [ + + ]: 29748 : if(has_qphase) then
2331 [ + + + + ]: 36007884 : atom_expiqr(1:2,1:nfgd,ia) = pawfgrtab(iatom)%expiqr(1:2,1:nfgd)
2332 : : end if
2333 : : end do ! ia
2334 : :
2335 : : !prod/dijhat_idij/atom_potg store the "phase" (qphase) index as the outermost
2336 : : !(4th) array dimension so that, for a fixed phase component, the (idat,ia)
2337 : : ![or (idat,ia,ic)] sub-block is contiguous and can be passed directly to
2338 : : !GEMM as a single dense (or strided-batched) matrix, with no repacking.
2339 [ - + + - : 62352 : ABI_MALLOC(prod,(lm_size,ndat,nattyp,qphase))
- + + - -
+ + - - +
+ - + - +
- - + -
+ ]
2340 [ - + + - : 62352 : ABI_MALLOC(dijhat_idij,(lmn2_size,ndat,nattyp,qphase))
- + + - -
+ + - - +
+ - + - +
- - + -
+ ]
2341 [ - + + - : 62352 : ABI_MALLOC(atom_potg,(nfgd_max,ndat,nattyp,qphase))
- + + - -
+ + - - +
+ - + - +
- - + -
+ ]
2342 : :
2343 : : #ifdef HAVE_OPENMP_OFFLOAD
2344 : : !$OMP TARGET ENTER DATA MAP(alloc:prod,dijhat_idij,atom_potg) IF(gpu_option_==ABI_GPU_OPENMP)
2345 : : !$OMP TARGET ENTER DATA MAP(to:atom_gylm,atom_ifftsph,atom_qijl,atom_nfgd) IF(gpu_option_==ABI_GPU_OPENMP)
2346 : : !$OMP TARGET ENTER DATA MAP(to:atom_expiqr) IF(gpu_option_==ABI_GPU_OPENMP .and. has_qphase)
2347 : : #endif
2348 : : !----------------------------------------------------------
2349 : : !Loop over spin components
2350 : : !----------------------------------------------------------
2351 [ - + ]: 10392 : nsploop=nsppol;if (ndij==4) nsploop=4
2352 [ + + ]: 20784 : do idij=1,nsploop
2353 [ - + - - : 20784 : if (idij<=nsppol.or.(nspden==4.and.idij<=3)) then
- - ]
2354 : :
2355 : 10392 : idijend=idij+idij/3
2356 [ + + ]: 20784 : do ispden=idij,idijend
2357 : :
2358 : : ! ------------------------------------------------------
2359 : : ! Gather V(r) on the atom's fine-grid points (and apply the
2360 : : ! exp(-i.q.r) phase rotation when relevant), zero-padded up
2361 : : ! to nfgd_max so the GEMMs below can use a common K=nfgd_max.
2362 : : ! ------------------------------------------------------
2363 : : ! Note for non-collinear magnetism:
2364 : : ! We compute Int[V^(alpha,beta)(r).g_l(r).Y_lm(r)]
2365 : : ! Remember: if nspden=4, V is stored as : V^11, V^22, V^12, i.V^21
2366 : : ! Note: has_qphase implies qphase==2 (has_qphase=qne0.and.qphase==2),
2367 : : ! so the qphase==1 & has_qphase combination can never occur.
2368 : :
2369 [ + - ]: 10392 : if(gpu_option_==ABI_GPU_DISABLED) then
2370 [ + + + + : 121858588 : atom_potg=zero
+ + + + ]
2371 [ # # ]: 0 : else if(gpu_option_==ABI_GPU_OPENMP) then
2372 : 0 : call gpu_set_to_zero(atom_potg,int(nfgd_max,c_size_t)*ndat*nattyp*qphase)
2373 : : end if
2374 : :
2375 [ - + ]: 10392 : if (qphase==1) then
2376 : : ! ===== Standard case, real Pot/dijhat ============================
2377 [ # # ]: 0 : if(gpu_option_==ABI_GPU_DISABLED) then
2378 : : !$OMP PARALLEL DO COLLAPSE(2) PRIVATE(ic,idat)
2379 [ # # ]: 0 : do ia=1,nattyp
2380 [ # # ]: 0 : do idat=1,ndat
2381 [ # # ]: 0 : do ic=1,atom_nfgd(ia)
2382 : 0 : atom_potg(ic,idat,ia,1)=Pot(atom_ifftsph(ic,ia),ispden,idat)
2383 : : end do
2384 : : end do
2385 : : end do ! ia
2386 : : else if(gpu_option_==ABI_GPU_OPENMP) then
2387 : : #ifdef HAVE_OPENMP_OFFLOAD
2388 : : !$OMP TARGET TEAMS DISTRIBUTE COLLAPSE(2) &
2389 : : !$OMP& PRIVATE(ia,idat) MAP(to:atom_potg,Pot,atom_ifftsph,atom_nfgd)
2390 : : do ia=1,nattyp
2391 : : do idat=1,ndat
2392 : : !$OMP PARALLEL DO PRIVATE(ic)
2393 : : do ic=1,atom_nfgd(ia)
2394 : : atom_potg(ic,idat,ia,1)=Pot(atom_ifftsph(ic,ia),ispden,idat)
2395 : : end do
2396 : : end do
2397 : : end do ! ia
2398 : : #endif
2399 : : end if
2400 [ + + ]: 10392 : else if (.not.has_qphase) then
2401 : : ! ===== qphase==2, no exp(iqr) phase ================================
2402 [ + - ]: 2571 : if(gpu_option_==ABI_GPU_DISABLED) then
2403 : : !$OMP PARALLEL DO COLLAPSE(2) PRIVATE(ic,idat,jc)
2404 [ + + ]: 6285 : do ia=1,nattyp
2405 [ + + ]: 20463 : do idat=1,ndat
2406 [ + + ]: 12845840 : do ic=1,atom_nfgd(ia)
2407 : 12827948 : jc=2*atom_ifftsph(ic,ia)
2408 : 12827948 : atom_potg(ic,idat,ia,1)=Pot(jc-1,ispden,idat)
2409 : 12842126 : atom_potg(ic,idat,ia,2)=Pot(jc ,ispden,idat)
2410 : : end do
2411 : : end do
2412 : : end do ! ia
2413 : : else if(gpu_option_==ABI_GPU_OPENMP) then
2414 : : #ifdef HAVE_OPENMP_OFFLOAD
2415 : : !$OMP TARGET TEAMS DISTRIBUTE COLLAPSE(2) &
2416 : : !$OMP& PRIVATE(ia,idat) MAP(to:atom_potg,Pot,atom_ifftsph,atom_nfgd)
2417 : : do ia=1,nattyp
2418 : : do idat=1,ndat
2419 : : !$OMP PARALLEL DO PRIVATE(ic,jc)
2420 : : do ic=1,atom_nfgd(ia)
2421 : : jc=2*atom_ifftsph(ic,ia)
2422 : : atom_potg(ic,idat,ia,1)=Pot(jc-1,ispden,idat)
2423 : : atom_potg(ic,idat,ia,2)=Pot(jc ,ispden,idat)
2424 : : end do
2425 : : end do
2426 : : end do ! ia
2427 : : #endif
2428 : : end if
2429 : : else
2430 : : ! ===== qphase==2, including exp(iqr) phase (DFPT only) =============
2431 : : ! The exp(-i.q.r) rotation is folded directly into the gather, so
2432 : : ! it disappears from the later GEMM-based reduction entirely.
2433 [ + - ]: 7821 : if(gpu_option_==ABI_GPU_DISABLED) then
2434 : : !$OMP PARALLEL DO COLLAPSE(2) PRIVATE(ic,idat,jc,vr,vi)
2435 [ + + ]: 23463 : do ia=1,nattyp
2436 [ + + ]: 86031 : do idat=1,ndat
2437 [ + + ]: 48067866 : do ic=1,atom_nfgd(ia)
2438 : 47989656 : jc=2*atom_ifftsph(ic,ia)
2439 : 47989656 : vr=Pot(jc-1,ispden,idat);vi=Pot(jc,ispden,idat)
2440 : 47989656 : atom_potg(ic,idat,ia,1)=vr*atom_expiqr(1,ic,ia)-vi*atom_expiqr(2,ic,ia)
2441 : 48052224 : atom_potg(ic,idat,ia,2)=vr*atom_expiqr(2,ic,ia)+vi*atom_expiqr(1,ic,ia)
2442 : : end do
2443 : : end do
2444 : : end do ! ia
2445 : : else if(gpu_option_==ABI_GPU_OPENMP) then
2446 : : #ifdef HAVE_OPENMP_OFFLOAD
2447 : : !$OMP TARGET TEAMS DISTRIBUTE COLLAPSE(2) &
2448 : : !$OMP& PRIVATE(ia,idat) MAP(to:atom_potg,Pot,atom_ifftsph,atom_expiqr,atom_nfgd)
2449 : : do ia=1,nattyp
2450 : : do idat=1,ndat
2451 : : !$OMP PARALLEL DO PRIVATE(ic,jc,vr,vi)
2452 : : do ic=1,atom_nfgd(ia)
2453 : : jc=2*atom_ifftsph(ic,ia)
2454 : : vr=Pot(jc-1,ispden,idat);vi=Pot(jc,ispden,idat)
2455 : : atom_potg(ic,idat,ia,1)=vr*atom_expiqr(1,ic,ia)-vi*atom_expiqr(2,ic,ia)
2456 : : atom_potg(ic,idat,ia,2)=vr*atom_expiqr(2,ic,ia)+vi*atom_expiqr(1,ic,ia)
2457 : : end do
2458 : : end do
2459 : : end do ! ia
2460 : : #endif
2461 : : end if
2462 : : end if
2463 : :
2464 : : ! ----------------------------------------------------------
2465 : : ! Compute Int[V(r).g_l(r).Y_lm(r)] = gylm^T . atom_potg
2466 : : ! as one GEMM per atom (CPU) or one batched-strided GEMM
2467 : : ! covering all atoms of the type (GPU).
2468 : : ! The ucvol/ngridtot scaling factor is folded in as GEMM alpha.
2469 : : ! ----------------------------------------------------------
2470 : 10392 : gemm_alpha=dcmplx(scal_r,zero)
2471 [ + - ]: 10392 : if(gpu_option_==ABI_GPU_DISABLED) then
2472 [ + + ]: 31176 : do icpq=1,qphase
2473 [ + + ]: 69888 : do ia=1,nattyp
2474 : : call abi_xgemm('T','N',lm_size,ndat,nfgd_max,gemm_alpha,&
2475 : : & atom_gylm(:,:,ia),nfgd_max,&
2476 : : & atom_potg(:,:,ia,icpq),nfgd_max,czero,&
2477 : 59496 : & prod(:,:,ia,icpq),lm_size,x_cplx=1)
2478 : : end do
2479 : : end do
2480 : : else if(gpu_option_==ABI_GPU_OPENMP) then
2481 : : #ifdef HAVE_OPENMP_OFFLOAD
2482 : : !$OMP TARGET DATA USE_DEVICE_ADDR(atom_gylm,atom_potg,prod)
2483 : : do icpq=1,qphase
2484 : : call abi_gpu_xgemm_strided(1,'T','N',lm_size,ndat,nfgd_max,gemm_alpha,&
2485 : : & c_loc(atom_gylm),nfgd_max,nfgd_max*lm_size,&
2486 : : & c_loc(atom_potg(1,1,1,icpq)),nfgd_max,nfgd_max*ndat,czero,&
2487 : : & c_loc(prod(1,1,1,icpq)),lm_size,lm_size*ndat,nattyp)
2488 : : end do
2489 : : !$OMP END TARGET DATA
2490 : : #endif
2491 : : end if
2492 : :
2493 : : ! ----------------------------------------------------------
2494 : : ! Compute Sum_(i,j)_LM { q_ij^L Int[V(r).g_l(r).Y_lm(r)] } = qijl^T . prod
2495 : : ! This reduction does not depend on the atom index (qijl is a
2496 : : ! per-atom-type quantity), so the whole (ndat,nattyp) batch is
2497 : : ! handled by a SINGLE dense GEMM call (no atom loop, no per-klmn
2498 : : ! Gaunt-selection loop needed: qijl is already zero wherever the
2499 : : ! Gaunt selection rule forbids a given (ilslm,klmn) pair).
2500 : : ! ----------------------------------------------------------
2501 [ + + ]: 31176 : do icpq=1,qphase
2502 : : call abi_xgemm('T','N',lmn2_size,ndat*nattyp,lm_size,cone,&
2503 : : & atom_qijl(:,:),l_size_full,&
2504 : : & prod(:,:,:,icpq),lm_size,czero,&
2505 : : & dijhat_idij(:,:,:,icpq),lmn2_size,&
2506 : 31176 : & x_cplx=1,gpu_option=gpu_option_)
2507 : : end do
2508 : :
2509 : : ! ----------------------------------------------------------
2510 : : ! Deduce some part of Dij according to symmetries
2511 : : ! ----------------------------------------------------------
2512 : :
2513 : : !if ispden=1 => real part of D^11_ij
2514 : : !if ispden=2 => real part of D^22_ij
2515 : : !if ispden=3 => real part of D^12_ij
2516 : : !if ispden=4 => imaginary part of D^12_ij
2517 : : !icpq=1 => real part of the exp(-i.q.r) phase (or the only part if qphase==1)
2518 : : !icpq=2 => imaginary part of the exp(-i.q.r) phase (qphase==2 only)
2519 [ + + ]: 41568 : do icpq=1,qphase
2520 [ + - ]: 31176 : if(gpu_option_==ABI_GPU_DISABLED) then
2521 : : !$OMP PARALLEL DO COLLAPSE(2) PRIVATE(idat,klmn,klmn1)
2522 [ + + ]: 59496 : do ia=1,nattyp
2523 [ + + ]: 212988 : do idat=1,ndat
2524 [ + + ]: 5570274 : do klmn=1,lmn2_size
2525 : 5378070 : klmn1=max(1,ispden-2)+(klmn-1+(icpq-1)*lmn2_size)*cplex_dij
2526 : 5531562 : dijhat(klmn1,idij+(idat-1)*ndij,ia)=dijhat_idij(klmn,idat,ia,icpq)
2527 : : end do
2528 : : end do
2529 : : end do ! ia
2530 : : else if(gpu_option_==ABI_GPU_OPENMP) then
2531 : : #ifdef HAVE_OPENMP_OFFLOAD
2532 : : !$OMP TARGET TEAMS DISTRIBUTE COLLAPSE(2) &
2533 : : !$OMP& MAP(to:dijhat,dijhat_idij) PRIVATE(idat,ia)
2534 : : do ia=1,nattyp
2535 : : do idat=1,ndat
2536 : : !$OMP PARALLEL DO PRIVATE(klmn,klmn1)
2537 : : do klmn=1,lmn2_size
2538 : : klmn1=max(1,ispden-2)+(klmn-1+(icpq-1)*lmn2_size)*cplex_dij
2539 : : dijhat(klmn1,idij+(idat-1)*ndij,ia)=dijhat_idij(klmn,idat,ia,icpq)
2540 : : end do
2541 : : end do
2542 : : end do ! ia
2543 : : #endif
2544 : : end if
2545 : : end do ! icpq
2546 : :
2547 : : end do !ispden
2548 : :
2549 : : !Non-collinear: D_ij(:,4)=D^21_ij=D^12_ij^*
2550 [ # # # # ]: 0 : else if (nspden==4.and.idij==4) then
2551 [ # # ]: 0 : do ia=1,nattyp
2552 [ # # ]: 0 : do idat=1,ndat
2553 [ # # ]: 0 : dijhat(:,idij+(idat-1)*ndij,ia)=dijhat(:,idij-1+(idat-1)*ndij,ia)
2554 : : end do
2555 : : end do ! ia
2556 [ # # ]: 0 : if (cplex_dij==2) then
2557 [ # # ]: 0 : do ia=1,nattyp
2558 [ # # ]: 0 : do idat=1,ndat
2559 [ # # ]: 0 : do klmn=2,lmn2_size*cplex_dij,cplex_dij
2560 [ # # ]: 0 : dijhat(klmn,idij+(idat-1)*ndij,ia)=-dijhat(klmn,idij+(idat-1)*ndij,ia)
2561 : : end do
2562 : : end do
2563 : : end do ! ia
2564 [ # # ]: 0 : if (qphase==2) then
2565 [ # # ]: 0 : do ia=1,nattyp
2566 [ # # ]: 0 : do idat=1,ndat
2567 [ # # ]: 0 : do klmn=2+lmn2_size*cplex_dij,2*lmn2_size*cplex_dij,cplex_dij
2568 [ # # ]: 0 : dijhat(klmn,idij+(idat-1)*ndij,ia)=-dijhat(klmn,idij+(idat-1)*ndij,ia)
2569 : : end do
2570 : : end do
2571 : : end do ! ia
2572 : : end if
2573 : : end if
2574 : :
2575 : : !Antiferro: D_ij(:,2)=D^down_ij=D^up_ij
2576 [ # # # # ]: 0 : else if (nsppol==1.and.idij==2) then
2577 [ # # ]: 0 : do ia=1,nattyp
2578 [ # # ]: 0 : do idat=1,ndat
2579 [ # # ]: 0 : dijhat(:,idij+(idat-1)*ndij,ia)=dijhat(:,idij-1+(idat-1)*ndij,ia)
2580 : : end do
2581 : : end do ! ia
2582 : : end if
2583 : :
2584 : : !----------------------------------------------------------
2585 : : !End loop on spin density components
2586 : : end do
2587 : :
2588 : : #ifdef HAVE_OPENMP_OFFLOAD
2589 : : !$OMP TARGET EXIT DATA MAP(delete:prod,dijhat_idij,atom_potg) IF(gpu_option_==ABI_GPU_OPENMP)
2590 : : !$OMP TARGET EXIT DATA MAP(delete:atom_gylm,atom_ifftsph,atom_qijl,atom_nfgd) IF(gpu_option_==ABI_GPU_OPENMP)
2591 : : !$OMP TARGET EXIT DATA MAP(delete:atom_expiqr) IF(gpu_option_==ABI_GPU_OPENMP .and. has_qphase)
2592 : : !$OMP TARGET EXIT DATA MAP(from:dijhat) IF(gpu_option_==ABI_GPU_OPENMP)
2593 : : #endif
2594 [ - + ]: 10392 : ABI_FREE(atom_nfgd)
2595 [ - + ]: 10392 : ABI_FREE(atom_gylm)
2596 [ - + ]: 10392 : ABI_FREE(atom_ifftsph)
2597 [ + + ]: 10392 : if(has_qphase) then
2598 [ - + ]: 7821 : ABI_FREE(atom_expiqr)
2599 : : end if
2600 : : !Free temporary memory spaces
2601 [ - + ]: 10392 : ABI_FREE(gnt_scal)
2602 [ - + ]: 10392 : ABI_FREE(atom_qijl)
2603 [ - + ]: 10392 : ABI_FREE(prod)
2604 [ - + ]: 10392 : ABI_FREE(dijhat_idij)
2605 [ - + ]: 10392 : ABI_FREE(atom_potg)
2606 : :
2607 [ + + ]: 29748 : do ia=1,nattyp
2608 : 19356 : iatom=iatm+ia
2609 [ + - ]: 19356 : if (pawfgrtab(iatom)%gylm_allocated==2) then
2610 [ - + ]: 19356 : ABI_FREE(pawfgrtab(iatom)%gylm)
2611 [ - + ]: 19356 : ABI_MALLOC(pawfgrtab(iatom)%gylm,(0,0))
2612 : 19356 : pawfgrtab(iatom)%gylm_allocated=0
2613 : : end if
2614 [ + + ]: 29748 : if (pawfgrtab(iatom)%expiqr_allocated==2) then
2615 [ - + ]: 15642 : ABI_FREE(pawfgrtab(iatom)%expiqr)
2616 [ - + ]: 15642 : ABI_MALLOC(pawfgrtab(iatom)%expiqr,(0,0))
2617 : 15642 : pawfgrtab(iatom)%expiqr_allocated=0
2618 : : end if
2619 : : end do ! ia
2620 : :
2621 [ - + ]: 10392 : end subroutine pawdijhat_ndat
2622 : : !!***
2623 : :
2624 : : !----------------------------------------------------------------------
2625 : :
2626 : : !!****f* m_paw_nhat/pawsushat
2627 : : !! NAME
2628 : : !! pawsushat
2629 : : !!
2630 : : !! FUNCTION
2631 : : !! PAW only, for susceptibility matrix:
2632 : : !! Compute contribution to the product of two wavefunctions (exchange charge density)
2633 : : !! from hat (compensation charge) density (in reciprocal space and eventually in real space):
2634 : : !! sushat_{ij,R}(g)=Sum_{L}[Q^L_ijR(g)]
2635 : : !!
2636 : : !! INPUTS
2637 : : !! atindx(natom)=index table for atoms, inverse of atindx
2638 : : !! cprj_k(natom,nspinor*nband_k)= wave functions projected with non-local projectors:
2639 : : !! cprj_k=<p_i|Cnk> where p_i is a non-local projector.
2640 : : !! WARNING: cprj(iatom,:) ARE SORTED BY ATOM TYPE !!!
2641 : : !! distribfft<type(distribfft_type)>=--optional-- contains infos related to FFT parallelism
2642 : : !! gbound_diel(2*mgfftdiel+8,2)=G sphere boundary for small FFT sphere.
2643 : : !! gylmg_diel(npwdiel,lmax_diel**2,ntypat)= -PAW only- Fourier transform of g_l(r).Y_ml(r) shape functions
2644 : : !! iband1,iband2= indices of the bands concerned with
2645 : : !! ispinor1,ispinor2= indices of spinorial components concerned with
2646 : : !! istwf_k=input option parameter that describes the storage of wfs
2647 : : !! kg_diel(3,npwdiel)=reduced planewave coordinates for the dielectric matrix.
2648 : : !! lmax_diel=1+max. value of l angular momentum used for dielectric matrix
2649 : : !! me_g0=--optional-- 1 if the current process treat the g=0 plane-wave (only needed when comm_fft is present)
2650 : : !! mgfftdiel=maximum size of 1D FFTs, for the computation of the dielectric matrix
2651 : : !! mpi_atmtab(:)=--optional-- indexes of the atoms treated by current proc
2652 : : !! comm_atom=--optional-- MPI communicator over atoms
2653 : : !! comm_fft=--optional-- MPI communicator over FT components
2654 : : !! natom=number of atoms in cell
2655 : : !! nband=number of bands at this k point for that spin polarization
2656 : : !! ndiel4,ndiel5,ndiel6= FFT dimensions, modified to avoid cache trashing
2657 : : !! nfftdiel=number of FFT grid points for the small (diel) grid
2658 : : !! ngfftdiel(18)=contain all needed information about 3D FFT, for dielectric matrix
2659 : : !! nspinor=number of spinorial components of the wavefunctions
2660 : : !! ntypat=number of types of atoms in unit cell.
2661 : : !! optreal=0 if WF product has to be output in reciprocal space
2662 : : !! 1 if WF product has to be output in real space
2663 : : !! paral_kgb=--optional-- 1 if "band-FFT" parallelism is activated (only needed when comm_fft is present)
2664 : : !! pawang <type(pawang_type)>=paw angular mesh and related data
2665 : : !! pawtab(ntypat) <type(pawtab_type)>=paw tabulated starting data
2666 : : !! ph3d_diel(2,npwdiel,natom*usepaw)=3-dim structure factors, for each atom and plane wave, for dielectric matrix
2667 : : !! typat(natom)=type (integer) for each atom
2668 : : !!
2669 : : !! SIDE EFFECTS
2670 : : !! === if optreal=0
2671 : : !! wfprod(2,npwdiel)=PAW contrib. to product of two wavefunctions (iband1,iband2):
2672 : : !! is added (in reciprocal space)
2673 : : !! === if optreal=1
2674 : : !! wfraug(2,ndiel4,ndiel5,ndiel6)=PAW contrib. to product of two wavefunctions (iband1,iband2)
2675 : : !! is added (in real space)
2676 : : !!
2677 : : !! SOURCE
2678 : :
2679 : 0 : subroutine pawsushat(atindx,cprj_k,gbound_diel,gylmg_diel,iband1,iband2,ispinor1,ispinor2,istwf_k,kg_diel,&
2680 : : & lmax_diel,mgfftdiel,natom,nband,ndiel4,ndiel5,ndiel6,&
2681 : : & ngfftdiel,npwdiel,nspinor,ntypat,optreal,&
2682 : 0 : & pawang,pawtab,ph3d_diel,typat,wfprod,wfraug, &
2683 [ # # # # : 0 : & mpi_atmtab,comm_atom,comm_fft,me_g0,paral_kgb,distribfft) ! optional arguments (parallelism)
# # ]
2684 : :
2685 : : !Arguments ---------------------------------------------
2686 : : !scalars
2687 : : integer,intent(in) :: iband1,iband2,ispinor1,ispinor2,istwf_k,lmax_diel,mgfftdiel
2688 : : integer,intent(in) :: natom,nband,ndiel4,ndiel5,ndiel6,npwdiel,nspinor
2689 : : integer,intent(in) :: ntypat,optreal
2690 : : integer,optional,intent(in) :: me_g0,comm_atom,comm_fft,paral_kgb
2691 : : type(distribfft_type),optional,intent(in),target :: distribfft
2692 : : type(pawang_type),intent(in) :: pawang
2693 : : !arrays
2694 : : integer,intent(in) :: atindx(natom),gbound_diel(2*mgfftdiel+8,2)
2695 : : integer,intent(in) :: kg_diel(3,npwdiel),ngfftdiel(18),typat(natom)
2696 : : integer,optional,target,intent(in) :: mpi_atmtab(:)
2697 : : real(dp),intent(in) :: gylmg_diel(npwdiel,lmax_diel**2,ntypat)
2698 : : real(dp),intent(in) :: ph3d_diel(2,npwdiel,natom)
2699 : : real(dp),intent(inout) :: wfprod(2,npwdiel*(1-optreal))
2700 : : real(dp),intent(inout) :: wfraug(2,ndiel4,ndiel5,ndiel6*optreal)
2701 : : type(pawcprj_type),intent(in) :: cprj_k(natom,nspinor*nband)
2702 : : type(pawtab_type),intent(in) :: pawtab(ntypat)
2703 : :
2704 : : !Local variables ---------------------------------------
2705 : : !scalars
2706 : : integer :: cplex,iatm,iatom,iatom_tot,ibsp1,ibsp2,ierr,il,ilmn,ils,ilslm,ipw
2707 : : integer :: itypat,j0lmn,jlmn,klm,klmn,lmax,lmin,mm,my_comm_atom,my_comm_fft,my_natom,tim_fourwf
2708 : : real(dp) :: phil1,phil2,sgn,weight_dum,wf1,wf2
2709 : : logical :: my_atmtab_allocated,parity,paral_atom
2710 : : type(distribfft_type),pointer :: my_distribfft
2711 : 0 : type(mpi_type) :: mpi_enreg_fft
2712 : : !arrays
2713 : 0 : integer,pointer :: my_atmtab(:)
2714 : : real(dp) :: ro(2),ro_ql(2)
2715 : 0 : real(dp),allocatable :: dummy(:,:),wfprod_paw(:,:),wfraug_paw(:,:,:,:)
2716 : :
2717 : : ! *************************************************************************
2718 : :
2719 : : DBG_ENTER("COLL")
2720 : :
2721 : 0 : if (present(comm_fft)) then
2722 [ # # ]: 0 : if ((.not.present(paral_kgb)).or.(.not.present(me_g0))) then
2723 : 0 : ABI_BUG('Need paral_kgb and me_g0 with comm_fft !')
2724 : : end if
2725 : : end if
2726 : :
2727 : : !Set up parallelism over atoms
2728 [ # # # # ]: 0 : paral_atom=(present(comm_atom).and.(my_natom/=natom))
2729 [ # # ]: 0 : nullify(my_atmtab);if (present(mpi_atmtab)) my_atmtab => mpi_atmtab
2730 [ # # ]: 0 : my_comm_atom=xmpi_comm_self;if (present(comm_atom)) my_comm_atom=comm_atom
2731 : 0 : call get_my_atmtab(my_comm_atom,my_atmtab,my_atmtab_allocated,paral_atom,natom)
2732 [ # # ]: 0 : my_natom=natom;if (paral_atom) my_natom=size(my_atmtab)
2733 : :
2734 [ # # ]: 0 : cplex=1;if (istwf_k>1) cplex=2
2735 [ # # # # : 0 : ABI_MALLOC(wfprod_paw,(2,npwdiel))
# # # # #
# # # ]
2736 [ # # # # ]: 0 : wfprod_paw(:,:)=zero
2737 : 0 : ibsp1=(iband1-1)*nspinor+ispinor1
2738 : 0 : ibsp2=(iband2-1)*nspinor+ispinor2
2739 : :
2740 : : !------------------------------------------------------------------------
2741 : : !----- Loop over atoms
2742 : : !------------------------------------------------------------------------
2743 [ # # ]: 0 : do iatom=1,my_natom
2744 [ # # ]: 0 : iatom_tot=iatom;if (paral_atom) iatom_tot=my_atmtab(iatom)
2745 : 0 : iatm=atindx(iatom_tot)
2746 : 0 : itypat=typat(iatom_tot)
2747 : :
2748 : : ! ------------------------------------------------------------------------
2749 : : ! ----- Loop over ij channels (basis components)
2750 : : ! ------------------------------------------------------------------------
2751 [ # # ]: 0 : do jlmn=1,pawtab(itypat)%lmn_size
2752 : 0 : j0lmn=jlmn*(jlmn-1)/2
2753 [ # # ]: 0 : do ilmn=1,jlmn
2754 : 0 : klmn=j0lmn+ilmn
2755 : 0 : klm =pawtab(itypat)%indklmn(1,klmn)
2756 : 0 : lmin=pawtab(itypat)%indklmn(3,klmn)
2757 : 0 : lmax=pawtab(itypat)%indklmn(4,klmn)
2758 : :
2759 : 0 : ro(1)=cprj_k(iatm,ibsp1)%cp(1,ilmn)*cprj_k(iatm,ibsp2)%cp(1,jlmn)
2760 [ # # ]: 0 : if (cplex==2) then
2761 : 0 : ro(1)=ro(1)+cprj_k(iatm,ibsp1)%cp(2,ilmn)*cprj_k(iatm,ibsp2)%cp(2,jlmn)
2762 : : ro(2)=cprj_k(iatm,ibsp1)%cp(2,ilmn)*cprj_k(iatm,ibsp2)%cp(1,jlmn) &
2763 : 0 : & -cprj_k(iatm,ibsp1)%cp(1,ilmn)*cprj_k(iatm,ibsp2)%cp(2,jlmn)
2764 : : end if
2765 [ # # ]: 0 : ro(1:cplex)=ro(1:cplex)*pawtab(itypat)%dltij(klmn)
2766 : :
2767 [ # # # # ]: 0 : do ils=lmin,lmax,2
2768 : 0 : il=mod(ils,4);parity=(mod(il,2)==0)
2769 [ # # ]: 0 : sgn=one;if (il>1) sgn=-one
2770 : :
2771 [ # # ]: 0 : do mm=-ils,ils
2772 : 0 : ilslm=ils*ils+ils+mm+1
2773 [ # # ]: 0 : if (pawang%gntselect(ilslm,klm)>0) then
2774 : :
2775 [ # # ]: 0 : ro_ql(1:cplex)=pawtab(itypat)%qijl(ilslm,klmn)*ro(1:cplex)
2776 : :
2777 : : ! Compute: Sum_{ijR} [ cpi* cpj qij^l (-i)^l g_l(g) S_lm(g) ]
2778 : :
2779 [ # # ]: 0 : if (cplex==1) then
2780 [ # # ]: 0 : if (parity) then
2781 [ # # ]: 0 : do ipw=1,npwdiel
2782 : 0 : phil1= sgn*ph3d_diel(1,ipw,iatm) ! (i)^l.exp(i.g.R)
2783 : 0 : phil2= sgn*ph3d_diel(2,ipw,iatm)
2784 : 0 : wf1= phil1*ro_ql(1) ! cpi* cpj qij^l (-i)^l.exp(-i.g.R)
2785 : 0 : wf2=-phil2*ro_ql(1)
2786 : 0 : wfprod_paw(1,ipw)=wfprod_paw(1,ipw)+wf1*gylmg_diel(ipw,ilslm,itypat)
2787 : 0 : wfprod_paw(2,ipw)=wfprod_paw(2,ipw)+wf2*gylmg_diel(ipw,ilslm,itypat)
2788 : : end do
2789 : : else
2790 [ # # ]: 0 : do ipw=1,npwdiel
2791 : 0 : phil1=-sgn*ph3d_diel(2,ipw,iatm) ! (i)^l.exp(i.g.R)
2792 : 0 : phil2= sgn*ph3d_diel(1,ipw,iatm)
2793 : 0 : wf1= phil1*ro_ql(1) ! cpi* cpj qij^l (-i)^l.exp(-i.g.R)
2794 : 0 : wf2=-phil2*ro_ql(1)
2795 : 0 : wfprod_paw(1,ipw)=wfprod_paw(1,ipw)+wf1*gylmg_diel(ipw,ilslm,itypat)
2796 : 0 : wfprod_paw(2,ipw)=wfprod_paw(2,ipw)+wf2*gylmg_diel(ipw,ilslm,itypat)
2797 : : end do
2798 : : end if
2799 : :
2800 : : else
2801 : :
2802 [ # # ]: 0 : if (parity) then
2803 [ # # ]: 0 : do ipw=1,npwdiel
2804 : 0 : phil1= sgn*ph3d_diel(1,ipw,iatm) ! (i)^l.exp(i.g.R)
2805 : 0 : phil2= sgn*ph3d_diel(2,ipw,iatm)
2806 : 0 : wf1=phil1*ro_ql(1)+phil2*ro_ql(2) ! cpi* cpj qij^l (-i)^l.exp(-i.g.R)
2807 : 0 : wf2=phil1*ro_ql(2)-phil2*ro_ql(1)
2808 : 0 : wfprod_paw(1,ipw)=wfprod_paw(1,ipw)+wf1*gylmg_diel(ipw,ilslm,itypat)
2809 : 0 : wfprod_paw(2,ipw)=wfprod_paw(2,ipw)+wf2*gylmg_diel(ipw,ilslm,itypat)
2810 : : end do
2811 : : else
2812 [ # # ]: 0 : do ipw=1,npwdiel
2813 : 0 : phil1=-sgn*ph3d_diel(2,ipw,iatm) ! (i)^l.exp(i.g.R)
2814 : 0 : phil2= sgn*ph3d_diel(1,ipw,iatm)
2815 : 0 : wf1=phil1*ro_ql(1)+phil2*ro_ql(2) ! cpi* cpj qij^l (-i)^l.exp(-i.g.R)
2816 : 0 : wf2=phil1*ro_ql(2)-phil2*ro_ql(1)
2817 : 0 : wfprod_paw(1,ipw)=wfprod_paw(1,ipw)+wf1*gylmg_diel(ipw,ilslm,itypat)
2818 : 0 : wfprod_paw(2,ipw)=wfprod_paw(2,ipw)+wf2*gylmg_diel(ipw,ilslm,itypat)
2819 : : end do
2820 : : end if
2821 : :
2822 : : end if
2823 : : end if
2824 : : end do
2825 : : end do
2826 : :
2827 : : ! ----- End loop over ij channels
2828 : : end do
2829 : : end do
2830 : :
2831 : : ! ----- End loop over atoms
2832 : : end do
2833 : :
2834 : : !Reduction in case of parallelism over atoms
2835 [ # # ]: 0 : if (paral_atom) then
2836 : 0 : call xmpi_sum(wfprod_paw,my_comm_atom,ierr)
2837 : : end if
2838 : :
2839 [ # # ]: 0 : if (optreal==0) then
2840 : :
2841 : : ! === Output in reciprocal space
2842 [ # # # # ]: 0 : wfprod(:,:)=wfprod(:,:)+wfprod_paw(:,:)
2843 : :
2844 : : else
2845 : : ! === Output in reciprocal space
2846 : 0 : tim_fourwf=17;weight_dum=0
2847 : : ! Create fake mpi_enreg to wrap fourdp
2848 [ # # ]: 0 : if (present(distribfft)) then
2849 : : my_distribfft => distribfft
2850 : : else
2851 [ # # ]: 0 : ABI_MALLOC(my_distribfft,)
2852 : 0 : call my_distribfft%init_seq('c',ngfftdiel(2),ngfftdiel(3),'fourwf')
2853 : : end if
2854 : 0 : call initmpi_seq(mpi_enreg_fft)
2855 [ # # # # : 0 : ABI_FREE(mpi_enreg_fft%distribfft)
# # # # #
# # # # #
# # # # #
# # # # #
# # ]
2856 [ # # ]: 0 : if (present(comm_fft)) then
2857 : 0 : call set_mpi_enreg_fft(mpi_enreg_fft,comm_fft,my_distribfft,me_g0,paral_kgb)
2858 : 0 : my_comm_fft=comm_fft
2859 : 0 : mpi_enreg_fft%paral_kgb = paral_kgb
2860 : : else
2861 : 0 : my_comm_fft=xmpi_comm_self
2862 : 0 : mpi_enreg_fft%paral_kgb = 0
2863 : 0 : mpi_enreg_fft%distribfft => my_distribfft
2864 : : end if
2865 : : ! do FFT
2866 [ # # # # : 0 : ABI_MALLOC(wfraug_paw,(2,ndiel4,ndiel5,ndiel6))
# # # # #
# # # # #
# # # # #
# ]
2867 : : call fourwf(1,dummy,wfprod_paw,dummy,wfraug_paw,gbound_diel,gbound_diel,&
2868 : : & istwf_k,kg_diel,kg_diel,mgfftdiel,mpi_enreg_fft,1,ngfftdiel,1,npwdiel,&
2869 : 0 : & ndiel4,ndiel5,ndiel6,0,tim_fourwf,weight_dum,weight_dum)
2870 [ # # # # : 0 : wfraug(:,:,:,:)=wfraug(:,:,:,:)+wfraug_paw(:,:,:,:)
# # # # ]
2871 [ # # ]: 0 : ABI_FREE(wfraug_paw)
2872 : 0 : call unset_mpi_enreg_fft(mpi_enreg_fft)
2873 [ # # ]: 0 : if (.not.present(distribfft)) then
2874 : 0 : call my_distribfft%free()
2875 [ # # # # : 0 : ABI_FREE(my_distribfft)
# # # # #
# # # # #
# # # # #
# # # #
# ]
2876 : : end if
2877 : : end if
2878 : :
2879 [ # # ]: 0 : ABI_FREE(wfprod_paw)
2880 : :
2881 : : !Destroy atom table used for parallelism
2882 : 0 : call free_my_atmtab(my_atmtab,my_atmtab_allocated)
2883 : :
2884 : : DBG_EXIT("COLL")
2885 : :
2886 [ # # # # : 0 : end subroutine pawsushat
# # # # #
# # # # #
# # # # #
# # # #
# ]
2887 : : !!***
2888 : :
2889 : : !----------------------------------------------------------------------
2890 : :
2891 : : !!****f* m_paw_nhat/nhatgrid
2892 : : !! NAME
2893 : : !! nhatgrid
2894 : : !!
2895 : : !! FUNCTION
2896 : : !! Determine parts of the rectangular (fine) grid that are contained
2897 : : !! inside spheres around atoms (used to compute n_hat density).
2898 : : !! If corresponding option is selected, compute also g_l(r)*Y_lm(r)
2899 : : !! (and derivatives) on this grid (g_l=radial shape function).
2900 : : !!
2901 : : !! INPUTS
2902 : : !! atindx1(natom)=index table for atoms, inverse of atindx (see gstate.f)
2903 : : !! distribfft<type(distribfft_type)>=--optional-- contains all the information related
2904 : : !! to the FFT parallelism and plane sharing
2905 : : !! gmet(3,3)=reciprocal space metric tensor in bohr**-2
2906 : : !! mpi_atmtab(:)=--optional-- indexes of the atoms treated by current proc
2907 : : !! comm_atom=--optional-- MPI communicator over atoms
2908 : : !! comm_fft=--optional-- MPI communicator over FFT components
2909 : : !! my_natom=number of atoms treated by current processor
2910 : : !! natom=total number of atoms in cell
2911 : : !! nattyp(ntypat)= # atoms of each type.
2912 : : !! ngfft(18)=contain all needed information about 3D FFT, see ~abinit/doc/variables/vargs.htm#ngfft
2913 : : !! ntypat=number of types of atoms in unit cell
2914 : : !! optcut= option for the cut-off radius of spheres:
2915 : : !! if optcut=0, cut-off radius=pawtab%rshp=cut-off radius of compensation charge
2916 : : !! if optcut=1, cut-off radius=pawtab%rpaw=radius of PAW augmentation regions
2917 : : !! optgr0= 1 if g_l(r)*Y_lm(r) are computed
2918 : : !! optgr1= 1 if first derivatives of g_l(r)*Y_lm(r) are computed
2919 : : !! optgr2= 1 if second derivatives of g_l(r)*Y_lm(r) are computed
2920 : : !! optrad= 1 if vectors (r-r_atom) on the fine grid around atoms have to be stored
2921 : : !! pawfgrtab(natom) <type(pawfgrtab_type)>=atomic data given on fine rectangular grid
2922 : : !! pawtab(ntypat) <type(pawtab_type)>=paw tabulated starting data
2923 : : !! rprimd(3,3)=dimensional primitive translations in real space (bohr)
2924 : : !! typat(natom)=type (integer) for each atom
2925 : : !! typord=1 if the output is ordered by type of atoms, 0 otherwise
2926 : : !! ucvol=unit cell volume in bohr**3
2927 : : !! xred(3,natom)=reduced dimensionless atomic coordinates
2928 : : !!
2929 : : !! OUTPUT
2930 : : !! pawfgrtab(natom)%ifftsph(nfgd)=FFT index (fine grid) of a points in paw spheres around each atom
2931 : : !! pawfgrtab(natom)%nfgd= number of (fine grid) FFT points in paw spheres around atoms
2932 : : !! if (optgr0==1)
2933 : : !! pawfgrtab(natom)%gylm(nfgd,l_size**2)= g_l(r)*Y_lm(r) around each atom
2934 : : !! if (optgr1==1)
2935 : : !! pawfgrtab(natom)%gylmgr(3,nfgd,l_size**2)= derivatives of g_l(r)*Y_lm(r) wrt cart. coordinates
2936 : : !! if (optgr2==1)
2937 : : !! pawfgrtab(natom)%gylmgr2(6,nfgd,l_size**2)= second derivatives of g_l(r)*Y_lm(r) wrt cart. coordinates
2938 : : !! if (optrad==1)
2939 : : !! pawfgrtab(natom)%rfgd(3,nfgd)= coordinates of r-r_atom around each atom
2940 : : !!
2941 : : !! SOURCE
2942 : :
2943 : 4720 : subroutine nhatgrid(atindx1,gmet,my_natom,natom,nattyp,ngfft,ntypat,&
2944 : 2360 : & optcut,optgr0,optgr1,optgr2,optrad,pawfgrtab,pawtab,rprimd,typat,ucvol,xred, &
2945 [ + + + + : 2360 : & mpi_atmtab,comm_atom,comm_fft,distribfft,typord) ! optional arguments (parallelism)
- + ]
2946 : :
2947 : : !Arguments ---------------------------------------------
2948 : : !scalars
2949 : : integer,intent(in) :: my_natom,natom,ntypat,optcut,optgr0,optgr1,optgr2,optrad
2950 : : integer,optional,intent(in) :: comm_atom,comm_fft,typord
2951 : : real(dp),intent(in) :: ucvol
2952 : : type(distribfft_type),optional,target,intent(in) :: distribfft
2953 : : !arrays
2954 : : integer,intent(in) :: ngfft(18),typat(natom)
2955 : : integer,intent(in),target :: atindx1(natom),nattyp(ntypat)
2956 : : integer,optional,target,intent(in) :: mpi_atmtab(:)
2957 : : real(dp),intent(in) :: gmet(3,3),rprimd(3,3),xred(3,natom)
2958 : : type(pawfgrtab_type),intent(inout) :: pawfgrtab(my_natom)
2959 : : type(pawtab_type),intent(in) :: pawtab(ntypat)
2960 : :
2961 : : !Local variables ------------------------------
2962 : : !scalars
2963 : : integer :: i3,iat,iatm,iatom,iatom_,iatom_tot,itypat,lm_size,me_fft,my_comm_atom,n1,n2,n3,nfgd
2964 : : logical :: grid_found,my_atmtab_allocated,paral_atom
2965 : : real(dp) :: rcut
2966 : : character(len=500) :: msg
2967 : : !arrays
2968 : 2360 : integer,allocatable :: ifftsph_tmp(:)
2969 : 2360 : integer,pointer :: my_atindx1(:),my_atmtab(:),my_nattyp(:)
2970 : 2360 : integer, ABI_CONTIGUOUS pointer :: fftn3_distrib(:),ffti3_local(:)
2971 : : real(dp) :: tsec(2)
2972 : 2360 : real(dp),allocatable :: rfgd_tmp(:,:)
2973 : :
2974 : : ! *************************************************************************
2975 : :
2976 : : DBG_ENTER("COLL")
2977 : :
2978 : 2360 : call timab(559,1,tsec)
2979 [ + + ]: 2360 : if (my_natom==0) return
2980 : :
2981 : : !Set up parallelism over FFT
2982 : 2202 : me_fft=0
2983 [ + + ]: 2202 : if (present(comm_fft)) then
2984 : 2050 : me_fft=xmpi_comm_rank(comm_fft)
2985 : : end if
2986 : :
2987 : : !Set up parallelism over atoms
2988 [ + + + + ]: 2202 : paral_atom=(present(comm_atom).and.(my_natom/=natom))
2989 [ + + ]: 2202 : nullify(my_atmtab);if (present(mpi_atmtab)) my_atmtab => mpi_atmtab
2990 [ + + ]: 2202 : my_comm_atom=xmpi_comm_self;if (present(comm_atom)) my_comm_atom=comm_atom
2991 : : call get_my_atmtab(my_comm_atom,my_atmtab,my_atmtab_allocated,paral_atom,natom,&
2992 : 2202 : & my_natom_ref=my_natom)
2993 [ + + ]: 2202 : if (paral_atom) then
2994 [ - + + - : 1860 : ABI_MALLOC(my_atindx1,(natom))
+ - + - -
+ - + ]
2995 [ - + + - : 1860 : ABI_MALLOC(my_nattyp,(ntypat))
+ - + - -
+ - + ]
2996 [ + + + + ]: 3320 : my_atindx1(:)=0;my_nattyp(:)=0
2997 : : iat=1
2998 [ + + ]: 1452 : do itypat=1,ntypat
2999 [ + - ]: 1452 : if (my_natom>0) then
3000 [ + + ]: 1866 : do iatom=1,my_natom
3001 [ + + ]: 1866 : if(typat(my_atmtab(iatom))==itypat)then
3002 : 771 : my_nattyp(itypat)=my_nattyp(itypat)+1
3003 : 771 : my_atindx1(iat)=iatom
3004 : 771 : iat=iat+1
3005 : : end if
3006 : : end do
3007 : : end if
3008 : : end do
3009 : : else
3010 : 1582 : my_atindx1 => atindx1
3011 : 1582 : my_nattyp => nattyp
3012 : : end if
3013 : :
3014 : : !Get the distrib associated with this fft_grid
3015 : 2202 : n1=ngfft(1);n2=ngfft(2);n3=ngfft(3)
3016 [ + + ]: 2202 : if (present(distribfft)) then
3017 : 2050 : grid_found=.false.
3018 [ + + ]: 2050 : if (n2 == distribfft%n2_coarse) then
3019 [ + + ]: 131 : if (n3== size(distribfft%tab_fftdp3_distrib)) then
3020 : 124 : fftn3_distrib => distribfft%tab_fftdp3_distrib
3021 : 124 : ffti3_local => distribfft%tab_fftdp3_local
3022 : 124 : grid_found=.true.
3023 : : end if
3024 : : end if
3025 [ + - ]: 2050 : if (n2 == distribfft%n2_fine) then
3026 [ + - ]: 2050 : if (n3 == size(distribfft%tab_fftdp3dg_distrib)) then
3027 : 2050 : fftn3_distrib => distribfft%tab_fftdp3dg_distrib
3028 : 2050 : ffti3_local => distribfft%tab_fftdp3dg_local
3029 : : grid_found = .true.
3030 : : end if
3031 : : end if
3032 [ # # ]: 0 : if (.not.(grid_found)) then
3033 : 0 : msg='Unable to find an allocated distrib for this fft grid!'
3034 : 0 : ABI_BUG(msg)
3035 : : end if
3036 : : else
3037 [ - + + - : 456 : ABI_MALLOC(fftn3_distrib,(n3))
+ - + - -
+ - + ]
3038 [ - + + - : 304 : ABI_MALLOC(ffti3_local,(n3))
- + - + ]
3039 [ + + + - : 13907 : fftn3_distrib=0;ffti3_local=(/(i3,i3=1,n3)/)
+ + + + ]
3040 : : end if
3041 : :
3042 : : !Loop over types of atom
3043 : : !-------------------------------------------
3044 : 2202 : iatm=0
3045 [ + + ]: 5264 : do itypat=1,ntypat
3046 : :
3047 [ + + ]: 3062 : if (optcut==1) then
3048 : 4 : rcut=pawtab(itypat)%rpaw
3049 : : else
3050 : 3058 : rcut=pawtab(itypat)%rshp
3051 : : end if
3052 : :
3053 : : ! Loop over atoms
3054 : : ! -------------------------------------------
3055 [ + + ]: 9997 : do iat=1,my_nattyp(itypat)
3056 : 4733 : iatm=iatm+1;iatom=my_atindx1(iatm)
3057 [ + + ]: 4733 : iatom_tot=iatom;if (paral_atom) iatom_tot=my_atmtab(iatom)
3058 [ + + - + ]: 4733 : iatom_=iatom;if(present(typord)) iatom_=merge(iatm,iatom,typord==1)
3059 : 4733 : lm_size=pawfgrtab(iatom_)%l_size**2
3060 : :
3061 : : ! ------------------------------------------------------------------
3062 : : ! A-Determine FFT points and r-R vectors around the atom
3063 : : ! ------------------------------------------------------------------
3064 : :
3065 : : call pawrfgd_fft(ifftsph_tmp,gmet,n1,n2,n3,nfgd,rcut,rfgd_tmp,rprimd,ucvol,&
3066 [ - + - + ]: 4733 : & xred(:,iatom_tot),fft_distrib=fftn3_distrib,fft_index=ffti3_local,me_fft=me_fft)
3067 : :
3068 : : ! Allocate arrays defining sphere (and related data) around current atom
3069 [ + - ]: 4733 : if (allocated(pawfgrtab(iatom_)%ifftsph)) then
3070 : 4733 : ABI_FREE(pawfgrtab(iatom_)%ifftsph)
3071 : : end if
3072 [ + + + - : 14199 : ABI_MALLOC(pawfgrtab(iatom_)%ifftsph,(nfgd))
+ - + + -
+ - + -
+ ]
3073 : 4733 : pawfgrtab(iatom_)%nfgd=nfgd
3074 [ + + ]: 9483857 : pawfgrtab(iatom_)%ifftsph(1:nfgd)=ifftsph_tmp(1:nfgd)
3075 : :
3076 [ + + ]: 4733 : if (optrad==1) then
3077 [ + - ]: 2762 : if (allocated(pawfgrtab(iatom_)%rfgd)) then
3078 : 2762 : ABI_FREE(pawfgrtab(iatom_)%rfgd)
3079 : : end if
3080 [ + + + - : 8286 : ABI_MALLOC(pawfgrtab(iatom_)%rfgd,(3,nfgd))
+ - + + -
+ - + -
+ ]
3081 : 2762 : pawfgrtab(iatom_)%rfgd_allocated=1
3082 [ + + + + ]: 19807382 : pawfgrtab(iatom_)%rfgd(1:3,1:nfgd)=rfgd_tmp(1:3,1:nfgd)
3083 : : end if
3084 : :
3085 [ + + ]: 4733 : if (optgr0==1) then
3086 [ + - ]: 3394 : if (allocated(pawfgrtab(iatom_)%gylm)) then
3087 : 3394 : ABI_FREE(pawfgrtab(iatom_)%gylm)
3088 : : end if
3089 [ + + + - : 13576 : ABI_MALLOC(pawfgrtab(iatom_)%gylm,(nfgd,lm_size))
- + + - +
- + + - +
- + - + ]
3090 : 3394 : pawfgrtab(iatom_)%gylm_allocated=1
3091 : : end if
3092 : :
3093 [ + + ]: 4733 : if (optgr1==1) then
3094 [ + - ]: 2955 : if (allocated(pawfgrtab(iatom_)%gylmgr)) then
3095 : 2955 : ABI_FREE(pawfgrtab(iatom_)%gylmgr)
3096 : : end if
3097 [ + + + - : 11820 : ABI_MALLOC(pawfgrtab(iatom_)%gylmgr,(3,nfgd,lm_size))
- + + - +
- + + - +
- + - + ]
3098 : 2955 : pawfgrtab(iatom_)%gylmgr_allocated=1
3099 : : end if
3100 : :
3101 [ + + ]: 4733 : if (optgr2==1) then
3102 [ + - ]: 97 : if (allocated(pawfgrtab(iatom_)%gylmgr2)) then
3103 : 97 : ABI_FREE(pawfgrtab(iatom_)%gylmgr2)
3104 : : end if
3105 [ - + + - : 388 : ABI_MALLOC(pawfgrtab(iatom_)%gylmgr2,(6,nfgd,lm_size))
- + + - +
- + - - +
- + - + ]
3106 : 97 : pawfgrtab(iatom_)%gylmgr2_allocated=1
3107 : : end if
3108 : :
3109 : : ! ------------------------------------------------------------------
3110 : : ! B-Calculate g_l(r-R)*Y_lm(r-R) for each r around the atom R
3111 : : ! ------------------------------------------------------------------
3112 [ + + ]: 4733 : if (optgr0+optgr1+optgr2>0) then
3113 : : call pawgylm(pawfgrtab(iatom_)%gylm,pawfgrtab(iatom_)%gylmgr,pawfgrtab(iatom_)%gylmgr2,&
3114 : 4675 : & lm_size,nfgd,optgr0,optgr1,optgr2,pawtab(itypat),rfgd_tmp(:,1:nfgd))
3115 : : end if
3116 : :
3117 : : ! End loops over types/atoms
3118 : : ! -------------------------------------------
3119 [ - + ]: 4733 : ABI_FREE(ifftsph_tmp)
3120 [ - + ]: 7795 : ABI_FREE(rfgd_tmp)
3121 : : end do
3122 : : end do
3123 : :
3124 : : !Destroy atom tables used for parallelism
3125 : 2202 : call free_my_atmtab(my_atmtab,my_atmtab_allocated)
3126 [ + + ]: 2202 : if (paral_atom) then
3127 [ - + ]: 620 : ABI_FREE(my_atindx1)
3128 [ - + ]: 620 : ABI_FREE(my_nattyp)
3129 : : end if
3130 : :
3131 [ + + ]: 2202 : if (.not.present(distribfft)) then
3132 [ - + ]: 152 : ABI_FREE(fftn3_distrib)
3133 [ - + ]: 152 : ABI_FREE(ffti3_local)
3134 : : end if
3135 : :
3136 : 2202 : call timab(559,2,tsec)
3137 : :
3138 : : DBG_EXIT("COLL")
3139 : :
3140 [ - + - + ]: 4562 : end subroutine nhatgrid
3141 : : !!***
3142 : :
3143 : : !----------------------------------------------------------------------
3144 : :
3145 : : !!****f* m_paw_nhat/wvl_nhatgrid
3146 : : !! NAME
3147 : : !! wvl_nhatgrid
3148 : : !!
3149 : : !! FUNCTION
3150 : : !! Determine parts of the rectangular (fine) grid that are contained
3151 : : !! inside spheres around atoms (used to compute n_hat density).
3152 : : !! If corresponding option is selected, compute also g_l(r)*Y_lm(r)
3153 : : !! (and derivatives) on this grid (g_l=radial shape function).
3154 : : !!
3155 : : !! INPUTS
3156 : : !!
3157 : : !! OUTPUT
3158 : : !! pawfgrtab(natom)%ifftsph(nfgd)=FFT index (fine grid) of a points in paw spheres around each atom
3159 : : !! pawfgrtab(natom)%nfgd= number of (fine grid) FFT points in paw spheres around atoms
3160 : : !! if (optgr0==1)
3161 : : !! pawfgrtab(natom)%gylm(nfgd,l_size**2)= g_l(r)*Y_lm(r) around each atom
3162 : : !! if (optgr1==1)
3163 : : !! pawfgrtab(natom)%gylmgr(3,nfgd,l_size**2)= derivatives of g_l(r)*Y_lm(r) wrt cart. coordinates
3164 : : !! if (optgr2==1)
3165 : : !! pawfgrtab(natom)%gylmgr2(6,nfgd,l_size**2)= second derivatives of g_l(r)*Y_lm(r) wrt cart. coordinates
3166 : : !! if (optrad==1)
3167 : : !! pawfgrtab(natom)%rfgd(3,nfgd)= coordinates of r-r_atom around each atom
3168 : : !!
3169 : : !! NOTES
3170 : : !! PENDING: ADD PARALELLISM OVER ATOMS:
3171 : : !! COPY NHATGRID
3172 : : !!
3173 : : !! SOURCE
3174 : :
3175 : 0 : subroutine wvl_nhatgrid(atindx1,geocode,h,i3s,natom,natom_tot,&
3176 : 0 : & nattyp,ntypat,n1,n1i,n2,n2i,n3,n3pi,optcut,optgr0,optgr1,optgr2,optrad,&
3177 : 0 : & pawfgrtab,pawtab,psppar,rprimd,shift,xred)
3178 : :
3179 : : !Arguments ---------------------------------------------
3180 : : !scalars
3181 : : integer,intent(in) :: i3s,natom,natom_tot,ntypat,optcut,optgr0,optgr1,optgr2,optrad
3182 : : integer,intent(in) :: n1,n2,n3,n1i,n2i,n3pi,shift
3183 : : real(dp),intent(in) :: h(3)
3184 : : character(1),intent(in) :: geocode
3185 : : !integer,intent(in),optional :: mpi_comm_wvl
3186 : : !arrays
3187 : : integer,intent(in) :: atindx1(natom),nattyp(ntypat)
3188 : : real(dp),intent(in) :: psppar(0:4,0:6,ntypat),rprimd(3,3)
3189 : : real(dp),intent(inout) :: xred(3,natom)
3190 : : type(pawfgrtab_type),intent(inout) :: pawfgrtab(natom)
3191 : : type(pawtab_type),intent(in) :: pawtab(ntypat)
3192 : :
3193 : : !Local variables ------------------------------
3194 : : !scalars
3195 : : !buffer to be added at the end of the last dimension of an array to control bounds_check
3196 : : integer :: iat,iatm,iatom,iatom_tot,itypat,lm_size,nfgd
3197 : 0 : real(dp) :: rloc,rshp,xcart(3,natom)
3198 : : !arrays
3199 : 0 : integer,allocatable :: ifftsph_tmp(:)
3200 : : real(dp) :: hh(3) !fine grid spacing for wavelets
3201 : : real(dp) :: tsec(2)
3202 : 0 : real(dp),allocatable :: rfgd_tmp(:,:)
3203 : :
3204 : : ! *************************************************************************
3205 : :
3206 : : DBG_ENTER("COLL")
3207 : :
3208 : : #if !defined HAVE_BIGDFT
3209 : 0 : BIGDFT_NOTENABLED_ERROR()
3210 : : #endif
3211 : :
3212 : 0 : call timab(559,1,tsec)
3213 : :
3214 : : !Set up parallelism for wvl
3215 : : !for debug: use me_wvl=xmpi_comm_rank(MPI_COMM_WORLD)
3216 : : !if (present(mpi_comm_wvl)) then
3217 : : !me_wvl=xmpi_comm_rank(mpi_comm_wvl)
3218 : : !nproc_wvl=xmpi_comm_size(mpi_comm_wvl)
3219 : : !else
3220 : : !me_wvl=0;nproc_wvl=1
3221 : : !end if
3222 : : !Pending: parallelism over atoms: see nhatgrid
3223 : :
3224 [ # # ]: 0 : if (natom_tot<natom) then ! This test has to be remove when natom_tot is used
3225 : 0 : ABI_BUG(' natom_tot<natom !')
3226 : : end if
3227 : :
3228 : : !Fine grid
3229 [ # # ]: 0 : hh(:)=0.5d0*h(:)
3230 : :
3231 : : !Compute xcart from xred
3232 : 0 : call xred2xcart(natom,rprimd,xcart,xred)
3233 : :
3234 : : !Loop over types of atom
3235 : 0 : iatm=0
3236 [ # # ]: 0 : do itypat=1,ntypat
3237 : :
3238 : 0 : rloc=psppar(0,0,itypat)
3239 [ # # ]: 0 : if (optcut==1) then
3240 : 0 : rshp=pawtab(itypat)%rpaw
3241 : : else
3242 : 0 : rshp=pawtab(itypat)%rshp
3243 : : end if
3244 : :
3245 : : ! Loop over atoms
3246 [ # # ]: 0 : do iat=1,nattyp(itypat)
3247 : 0 : iatm=iatm+1;iatom=atindx1(iatm)
3248 : 0 : iatom_tot=iatom; !if (paral_atom) iatom_tot=my_atmtab(iatom)
3249 : 0 : lm_size=pawfgrtab(iatom)%l_size**2
3250 : :
3251 : : ! Determine FFT points and r-R vectors around the atom
3252 : : call pawrfgd_wvl(geocode,hh,ifftsph_tmp,i3s,n1,n1i,n2,n2i,n3,n3pi,nfgd,rshp,rloc,&
3253 [ # # # # ]: 0 : & rfgd_tmp,shift,xcart(:,iatom_tot))
3254 : :
3255 : : ! Allocate arrays defining sphere (and related data) around current atom
3256 [ # # ]: 0 : if (allocated(pawfgrtab(iatom)%ifftsph)) then
3257 : 0 : ABI_FREE(pawfgrtab(iatom)%ifftsph)
3258 : : end if
3259 [ # # # # : 0 : ABI_MALLOC(pawfgrtab(iatom)%ifftsph,(nfgd))
# # # # #
# # # #
# ]
3260 : 0 : pawfgrtab(iatom)%nfgd=nfgd
3261 [ # # ]: 0 : pawfgrtab(iatom)%ifftsph(1:nfgd)=ifftsph_tmp(1:nfgd)
3262 : :
3263 [ # # ]: 0 : if (optrad==1) then
3264 [ # # ]: 0 : if (allocated(pawfgrtab(iatom)%rfgd)) then
3265 : 0 : ABI_FREE(pawfgrtab(iatom)%rfgd)
3266 : : end if
3267 [ # # # # : 0 : ABI_MALLOC(pawfgrtab(iatom)%rfgd,(3,nfgd))
# # # # #
# # # #
# ]
3268 : 0 : pawfgrtab(iatom)%rfgd_allocated=1
3269 [ # # # # ]: 0 : pawfgrtab(iatom)%rfgd(1:3,1:nfgd)=rfgd_tmp(1:3,1:nfgd)
3270 : : end if
3271 : :
3272 [ # # ]: 0 : if (optgr0==1) then
3273 [ # # ]: 0 : if (allocated(pawfgrtab(iatom)%gylm)) then
3274 : 0 : ABI_FREE(pawfgrtab(iatom)%gylm)
3275 : : end if
3276 [ # # # # : 0 : ABI_MALLOC(pawfgrtab(iatom)%gylm,(nfgd,lm_size))
# # # # #
# # # # #
# # # # ]
3277 : 0 : pawfgrtab(iatom)%gylm_allocated=1
3278 : : end if
3279 : :
3280 [ # # ]: 0 : if (optgr1==1) then
3281 [ # # ]: 0 : if (allocated(pawfgrtab(iatom)%gylmgr)) then
3282 : 0 : ABI_FREE(pawfgrtab(iatom)%gylmgr)
3283 : : end if
3284 [ # # # # : 0 : ABI_MALLOC(pawfgrtab(iatom)%gylmgr,(3,nfgd,lm_size))
# # # # #
# # # # #
# # # # ]
3285 : 0 : pawfgrtab(iatom)%gylmgr_allocated=1
3286 : : end if
3287 : :
3288 [ # # ]: 0 : if (optgr2==1) then
3289 [ # # ]: 0 : if (allocated(pawfgrtab(iatom)%gylmgr2)) then
3290 : 0 : ABI_FREE(pawfgrtab(iatom)%gylmgr2)
3291 : : end if
3292 [ # # # # : 0 : ABI_MALLOC(pawfgrtab(iatom)%gylmgr2,(6,nfgd,lm_size))
# # # # #
# # # # #
# # # # ]
3293 : 0 : pawfgrtab(iatom)%gylmgr2_allocated=1
3294 : : end if
3295 : :
3296 : : ! Calculate g_l(r-R)*Y_lm(r-R) for each r around the atom R
3297 [ # # ]: 0 : if (optgr0+optgr1+optgr2>0) then
3298 : : call pawgylm(pawfgrtab(iatom)%gylm,pawfgrtab(iatom)%gylmgr,pawfgrtab(iatom)%gylmgr2,&
3299 : 0 : & lm_size,nfgd,optgr0,optgr1,optgr2,pawtab(itypat),rfgd_tmp(:,1:nfgd))
3300 : : end if
3301 : :
3302 : : ! End loops over types/atoms
3303 [ # # ]: 0 : ABI_FREE(ifftsph_tmp)
3304 [ # # ]: 0 : ABI_FREE(rfgd_tmp)
3305 : : end do
3306 : : end do
3307 : :
3308 : 0 : call timab(559,2,tsec)
3309 : :
3310 : : DBG_EXIT("COLL")
3311 : :
3312 [ # # # # ]: 0 : end subroutine wvl_nhatgrid
3313 : : !!***
3314 : :
3315 : : !----------------------------------------------------------------------
3316 : :
3317 : : END MODULE m_paw_nhat
3318 : : !!***
|
