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 33217 : subroutine pawmknhat(compch_fft,cplex,ider,idir,ipert,izero,gprimd,&
121 33217 : & my_natom,natom,nfft,ngfft,nhatgrdim,nspden,ntypat,pawang,pawfgrtab,&
122 66434 : & pawgrnhat,pawnhat,pawrhoij,pawrhoij0,pawtab,qphon,rprimd,ucvol,usewvl,xred,&
123 33217 : & 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 33217 : type(mpi_type) :: mpi_enreg_fft
154 : !arrays
155 33217 : integer,pointer :: my_atmtab(:)
156 66434 : real(dp) :: ro(cplex),ro_ql(cplex),tmp_compch_fft(nspden),tsec(2)
157 33217 : real(dp),allocatable :: pawgrnhat_atm(:,:),pawnhat_atm(:),work(:,:)
158 :
159 : ! *************************************************************************
160 :
161 : DBG_ENTER("COLL")
162 :
163 33217 : compute_nhat=(ider==0.or.ider==2)
164 33217 : compute_grad=(ider==1.or.ider==2)
165 6579 : compute_phonons=(ipert>0.and.ipert<=natom)
166 :
167 : !Compatibility tests
168 33217 : qeq0=(qphon(1)**2+qphon(2)**2+qphon(3)**2<1.d-15)
169 33217 : if (present(comm_fft)) then
170 32109 : 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 33217 : if(ider>0.and.nhatgrdim==0) then
175 0 : ABI_BUG('Gradients of nhat required but not allocated!')
176 : end if
177 33217 : if (my_natom>0) then
178 29945 : if(nspden>1.and.nspden/=pawrhoij(1)%nspden) then
179 0 : ABI_BUG('Wrong values for nspden and pawrhoij%nspden!')
180 : end if
181 29945 : if(nspden>1.and.nspden/=pawfgrtab(1)%nspden) then
182 0 : ABI_BUG('Wrong values for nspden and pawfgrtab%nspden!')
183 : end if
184 29945 : if(pawrhoij(1)%qphase<cplex) then
185 0 : ABI_BUG('Must have pawrhoij()%qphase >= cplex!')
186 : end if
187 29945 : 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 33217 : if (ipert==natom+1.or.ipert==natom+10) return
199 :
200 : !Set up parallelism over atoms
201 33217 : paral_atom=(present(comm_atom).and.(my_natom/=natom))
202 33217 : nullify(my_atmtab);if (present(mpi_atmtab)) my_atmtab => mpi_atmtab
203 33217 : 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 33217 : & my_natom_ref=my_natom)
206 :
207 : !Initialisations
208 33217 : if ((.not.compute_nhat).and.(.not.compute_grad)) return
209 105326 : mfgd=zero;if (my_natom>0) mfgd=maxval(pawfgrtab(1:my_natom)%nfgd)
210 33217 : if (compute_nhat) then
211 98850 : ABI_MALLOC(pawnhat_atm,(cplex*mfgd))
212 1297353424 : pawnhat=zero
213 : end if
214 33217 : 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 33217 : mpi_comm_sphgrid=xmpi_comm_self !no communicators passed
221 33217 : if(present(comm_fft) .and. usewvl==0) mpi_comm_sphgrid=comm_fft
222 33217 : if(present(mpi_comm_wvl) .and. usewvl==1) mpi_comm_sphgrid=mpi_comm_wvl
223 :
224 : !------------------------------------------------------------------------
225 : !----- Loop over atoms
226 : !------------------------------------------------------------------------
227 :
228 105326 : do iatom=1,my_natom
229 72109 : iatom_tot=iatom;if (paral_atom) iatom_tot=my_atmtab(iatom)
230 :
231 72109 : itypat=pawrhoij(iatom)%itypat
232 72109 : lm_size=pawfgrtab(iatom)%l_size**2
233 72109 : need_frozen=((compute_nhat).and.(ipert==iatom_tot.or.ipert==natom+3.or.ipert==natom+4))
234 72109 : nfgd=pawfgrtab(iatom)%nfgd
235 72109 : cplex_rhoij=pawrhoij(iatom)%cplex_rhoij
236 72109 : 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 72109 : & ((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 72109 : 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 72109 : 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 72109 : 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 160669 : do ispden=1,nspden
308 :
309 155212901 : if (compute_nhat) pawnhat_atm(1:cplex*nfgd)=zero
310 33724011 : if (compute_grad) pawgrnhat_atm(1:cplex*nfgd,1:3)=zero
311 :
312 : ! ------------------------------------------------------------------------
313 : ! ----- Loop over ij channels (basis components)
314 : ! ------------------------------------------------------------------------
315 88560 : jrhoij=1
316 3725423 : do irhoij=1,pawrhoij(iatom)%nrhoijsel
317 3636863 : klmn=pawrhoij(iatom)%rhoijselect(irhoij)
318 3636863 : klm =pawtab(itypat)%indklmn(1,klmn)
319 3636863 : lmin=pawtab(itypat)%indklmn(3,klmn)
320 3636863 : lmax=pawtab(itypat)%indklmn(4,klmn)
321 :
322 : ! Retrieve rhoij
323 3636863 : if (pawrhoij(iatom)%nspden/=2) then
324 2733957 : ro(1)=pawrhoij(iatom)%rhoijp(jrhoij,ispden)
325 2733957 : 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 7308737 : ro(1:cplex)=pawtab(itypat)%dltij(klmn)*ro(1:cplex)
336 :
337 3636863 : if (compute_nhat) then
338 3586163 : if (cplex==1) then
339 10537585 : do ils=lmin,lmax,2
340 39522646 : do mm=-ils,ils
341 28985061 : ilslm=ils*ils+ils+mm+1
342 35971494 : if (pawang%gntselect(ilslm,klm)>0) then
343 7065828 : ro_ql(1)=ro(1)*pawtab(itypat)%qijl(ilslm,klmn)
344 : !$OMP PARALLEL DO PRIVATE(ic)
345 26729128899 : do ic=1,nfgd
346 26729128899 : 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 3636863 : 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 3725423 : jrhoij=jrhoij+cplex_rhoij
408 : end do
409 :
410 : ! If RF calculation, add frozen part of 1st-order compensation density
411 88560 : 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 88560 : 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 88560 : if (compute_nhat) then
488 87496 : if (cplex==1) then
489 : ! Not possible to parallelize here (see comment above)
490 152182389 : do ic=1,nfgd
491 152095837 : kc=pawfgrtab(iatom)%ifftsph(ic)
492 152182389 : 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 160669 : 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 72109 : 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 72109 : 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 72109 : 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 72109 : 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 72109 : 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 105326 : 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 33217 : if (compute_nhat) then
562 32950 : ABI_FREE(pawnhat_atm)
563 : end if
564 33217 : if (compute_grad) then
565 778 : ABI_FREE(pawgrnhat_atm)
566 : end if
567 :
568 : !----- Reduction in case of parallelism
569 33217 : 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 33217 : 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 33217 : if (compute_nhat.and.ipert==0) then
616 105484 : nfftot=PRODUCT(ngfft(1:3))
617 : call mean_fftr(pawnhat,tmp_compch_fft,nfft,nfftot,1,&
618 26371 : & mpi_comm_sphgrid=mpi_comm_sphgrid,gpu_thread_limit=gpu_thread_limit)
619 26371 : compch_fft = tmp_compch_fft(1)
620 26371 : compch_fft=compch_fft*ucvol
621 : end if
622 :
623 : !Destroy atom table used for parallelism
624 33217 : call free_my_atmtab(my_atmtab,my_atmtab_allocated)
625 :
626 : DBG_EXIT("COLL")
627 :
628 66434 : 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 3480 : subroutine pawnhatfr(ider,idir,ipert,my_natom,natom,nspden,ntypat,&
1879 3480 : & pawang,pawfgrtab,pawrhoij,pawtab,rprimd, &
1880 3480 : & 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 3480 : integer,pointer :: my_atmtab(:)
1904 3480 : 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 844 : 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 6812 : 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 4704 : subroutine nhatgrid(atindx1,gmet,my_natom,natom,nattyp,ngfft,ntypat,&
2944 2352 : & optcut,optgr0,optgr1,optgr2,optrad,pawfgrtab,pawtab,rprimd,typat,ucvol,xred, &
2945 2352 : & 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 2352 : integer,allocatable :: ifftsph_tmp(:)
2969 2352 : integer,pointer :: my_atindx1(:),my_atmtab(:),my_nattyp(:)
2970 2352 : integer, ABI_CONTIGUOUS pointer :: fftn3_distrib(:),ffti3_local(:)
2971 : real(dp) :: tsec(2)
2972 2352 : real(dp),allocatable :: rfgd_tmp(:,:)
2973 :
2974 : ! *************************************************************************
2975 :
2976 : DBG_ENTER("COLL")
2977 :
2978 2352 : call timab(559,1,tsec)
2979 2352 : if (my_natom==0) return
2980 :
2981 : !Set up parallelism over FFT
2982 2194 : me_fft=0
2983 2194 : if (present(comm_fft)) then
2984 2046 : me_fft=xmpi_comm_rank(comm_fft)
2985 : end if
2986 :
2987 : !Set up parallelism over atoms
2988 2194 : paral_atom=(present(comm_atom).and.(my_natom/=natom))
2989 2194 : nullify(my_atmtab);if (present(mpi_atmtab)) my_atmtab => mpi_atmtab
2990 2194 : 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 2194 : & my_natom_ref=my_natom)
2993 2194 : 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 1574 : my_atindx1 => atindx1
3011 1574 : my_nattyp => nattyp
3012 : end if
3013 :
3014 : !Get the distrib associated with this fft_grid
3015 2194 : n1=ngfft(1);n2=ngfft(2);n3=ngfft(3)
3016 2194 : if (present(distribfft)) then
3017 2046 : grid_found=.false.
3018 2046 : 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 2046 : if (n2 == distribfft%n2_fine) then
3026 2046 : if (n3 == size(distribfft%tab_fftdp3dg_distrib)) then
3027 2046 : fftn3_distrib => distribfft%tab_fftdp3dg_distrib
3028 2046 : 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 444 : ABI_MALLOC(fftn3_distrib,(n3))
3038 296 : ABI_MALLOC(ffti3_local,(n3))
3039 12703 : fftn3_distrib=0;ffti3_local=(/(i3,i3=1,n3)/)
3040 : end if
3041 :
3042 : !Loop over types of atom
3043 : !-------------------------------------------
3044 2194 : iatm=0
3045 5248 : do itypat=1,ntypat
3046 :
3047 3054 : if (optcut==1) then
3048 4 : rcut=pawtab(itypat)%rpaw
3049 : else
3050 3050 : rcut=pawtab(itypat)%rshp
3051 : end if
3052 :
3053 : ! Loop over atoms
3054 : ! -------------------------------------------
3055 9973 : do iat=1,my_nattyp(itypat)
3056 4725 : iatm=iatm+1;iatom=my_atindx1(iatm)
3057 4725 : iatom_tot=iatom;if (paral_atom) iatom_tot=my_atmtab(iatom)
3058 4725 : iatom_=iatom;if(present(typord)) iatom_=merge(iatm,iatom,typord==1)
3059 4725 : 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 4725 : & 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 4725 : if (allocated(pawfgrtab(iatom_)%ifftsph)) then
3070 4725 : ABI_FREE(pawfgrtab(iatom_)%ifftsph)
3071 : end if
3072 14175 : ABI_MALLOC(pawfgrtab(iatom_)%ifftsph,(nfgd))
3073 4725 : pawfgrtab(iatom_)%nfgd=nfgd
3074 9471195 : pawfgrtab(iatom_)%ifftsph(1:nfgd)=ifftsph_tmp(1:nfgd)
3075 :
3076 4725 : 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 4725 : if (optgr0==1) then
3086 3386 : if (allocated(pawfgrtab(iatom_)%gylm)) then
3087 3386 : ABI_FREE(pawfgrtab(iatom_)%gylm)
3088 : end if
3089 13544 : ABI_MALLOC(pawfgrtab(iatom_)%gylm,(nfgd,lm_size))
3090 3386 : pawfgrtab(iatom_)%gylm_allocated=1
3091 : end if
3092 :
3093 4725 : 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 4725 : 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 4725 : if (optgr0+optgr1+optgr2>0) then
3113 : call pawgylm(pawfgrtab(iatom_)%gylm,pawfgrtab(iatom_)%gylmgr,pawfgrtab(iatom_)%gylmgr2,&
3114 4667 : & lm_size,nfgd,optgr0,optgr1,optgr2,pawtab(itypat),rfgd_tmp(:,1:nfgd))
3115 : end if
3116 :
3117 : ! End loops over types/atoms
3118 : ! -------------------------------------------
3119 4725 : ABI_FREE(ifftsph_tmp)
3120 7779 : ABI_FREE(rfgd_tmp)
3121 : end do
3122 : end do
3123 :
3124 : !Destroy atom tables used for parallelism
3125 2194 : call free_my_atmtab(my_atmtab,my_atmtab_allocated)
3126 2194 : if (paral_atom) then
3127 620 : ABI_FREE(my_atindx1)
3128 620 : ABI_FREE(my_nattyp)
3129 : end if
3130 :
3131 2194 : if (.not.present(distribfft)) then
3132 148 : ABI_FREE(fftn3_distrib)
3133 148 : ABI_FREE(ffti3_local)
3134 : end if
3135 :
3136 2194 : call timab(559,2,tsec)
3137 :
3138 : DBG_EXIT("COLL")
3139 :
3140 4546 : 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 : !!***
|
