Line data Source code
1 : !!****m* ABINIT/m_cgprj
2 : !! NAME
3 : !! m_cgprj
4 : !!
5 : !! FUNCTION
6 : !! Routines to compute <Proj_i|Cnk> with |Cnk> expressed in reciprocal space.
7 : !!
8 : !! COPYRIGHT
9 : !! Copyright (C) 1998-2026 ABINIT group (MT)
10 : !! This file is distributed under the terms of the
11 : !! GNU General Public License, see ~abinit/COPYING
12 : !! or http://www.gnu.org/copyleft/gpl.txt .
13 : !!
14 : !! SOURCE
15 :
16 : #if defined HAVE_CONFIG_H
17 : #include "config.h"
18 : #endif
19 :
20 : #include "abi_common.h"
21 :
22 : module m_cgprj
23 :
24 : use defs_basis
25 : use m_abicore
26 : use m_errors
27 : use m_xmpi
28 :
29 : use defs_abitypes, only : MPI_type
30 : use defs_datatypes, only : pseudopotential_type
31 : use m_fstrings, only : itoa, sjoin
32 : use m_kg, only : ph1d3d, mkkpg
33 : use m_geometry, only : strconv
34 : use m_mkffnl, only : mkffnl
35 : use m_mpinfo, only : proc_distrb_cycle
36 : use m_pawcprj, only : pawcprj_type, pawcprj_alloc, pawcprj_put, pawcprj_free, &
37 : pawcprj_set_zero, pawcprj_mpi_sum, pawcprj_copy, pawcprj_lincom
38 : use m_opernla_ylm, only : opernla_ylm
39 : use m_opernla_ylm_mv, only : opernla_ylm_mv
40 : use m_opernla_gemm, only : opernla_gemm
41 : use m_gemm_nonlop_projectors
42 : use m_time, only : timab
43 : use m_io_tools, only : flush_unit
44 :
45 : use m_xg
46 : use m_xg_nonlop
47 :
48 : implicit none
49 :
50 : private
51 : !!***
52 :
53 : public :: getcprj
54 : public :: cprj_rotate
55 : public :: ctocprj
56 : public :: xg_cprj_copy
57 :
58 : integer,parameter,public :: XG_TO_CPRJ=1
59 : integer,parameter,public :: CPRJ_TO_XG=2
60 : !!***
61 :
62 : contains
63 : !!***
64 :
65 : !!****f* ABINIT/getcprj
66 : !! NAME
67 : !! getcprj
68 : !!
69 : !! FUNCTION
70 : !! Compute <Proj_i|Cnk> for one wave function |Cnk> expressed in reciprocal space.
71 : !! Compute also derivatives of <Proj_i|Cnk>.
72 : !! |Proj_i> are non-local projectors (for each atom and each l,m,n)
73 : !!
74 : !! INPUTS
75 : !! choice=chooses possible output:
76 : !! In addition to projected wave function:
77 : !! choice=1 => nothing else
78 : !! =2 => 1st gradients with respect to atomic position(s)
79 : !! =3 => 1st gradients with respect to strain(s)
80 : !! =23=> 1st gradients with respect to strain(s) and atm pos
81 : !! =4 => 2nd derivatives with respect to atomic pos.
82 : !! =24=> 1st and 2nd derivatives with respect to atomic pos.
83 : !! =5 => 1st gradients with respect to k wavevector
84 : !! =6 => 2nd derivatives with respect to strain and atm. pos.
85 : !! cpopt=1 if <Proj_i|Cnk> are already in memory; see below (side effects).
86 : !! cwavef(2,nspinor*npw_k)=input cmplx wavefunction coefficients <G|Cnk>
87 : !! ffnl(npw_k,dimffnl,lmnmax,ntypat)=nonlocal form factors to be used for the application of the nl operator
88 : !! idir=direction of the derivative, i.e. dir. of - atom to be moved in the case choice=2
89 : !! - strain component in the case choice=3
90 : !! - k point direction in the case choice=5
91 : !! Compatible only with choice=2,3,5; if idir=0, all derivatives are computed
92 : !! indlmn(6,i,ntypat)= array giving l,m,n,lm,ln,s for i=lmn
93 : !! istwf_k=option parameter that describes the storage of wfs
94 : !! kg_k(3,npw_k)=reduced planewave coordinates
95 : !! kpg(npw_k,npkg)=(k+G) components and related data
96 : !! kpoint(3)=k point in terms of recip. translations
97 : !! lmnmax=max. number of (l,m,n) components over all types of atoms
98 : !! mgfft=maximum size of 1D FFTs
99 : !! mpi_enreg=information about MPI parallelization
100 : !! natom=number of atoms in cell
101 : !! nattyp(ntypat)=number of atoms of each type
102 : !! ngfft(18)=contain all needed information about 3D FFT, see ~ABINIT/Infos/vargs.htm#ngfft
103 : !! nloalg(3)=governs the choice of the algorithm for nonlocal operator
104 : !! npw_k=number of planewaves for given k point
105 : !! nspinor=number of spinorial components of the wavefunctions (on current proc)
106 : !! ntypat=number of types of atoms in unit cell
107 : !! phkxred(2,natom)=phase factors exp(2 pi kpoint.xred)
108 : !! ph1d(2,3*(2*mgfft+1)*natom)=1-dim structure factor phase information
109 : !! ph3d(2,npw_k,natom)=3D structure factors, for each atom and plane wave
110 : !! only used if nloalg(2)>0
111 : !! ucvol= unit cell volume
112 : !! useylm=governs the way the nonlocal operator is to be applied
113 : !! is_kprime=wether provided arrays relate to K or Kprime for Hamiltonian operator (optional, for GEMM nonlop)
114 : !!
115 : !! SIDE EFFECTS
116 : !! cwaveprj(natom,nspinor) <type(pawcprj_type)>=projected input wave function <Proj_i|Cnk> with all NL projectors
117 : !! (and derivatives)
118 : !! if cpopt=1 the projected scalars have been already been computed and
119 : !! only derivatives are computed here
120 : !! if cpopt=0 the projected scalars and derivatives are computed here
121 : !!
122 : !! TODO
123 : !! Spin-orbit
124 : !!
125 : !! SOURCE
126 :
127 4264420 : subroutine getcprj(choice,cpopt,cwavef,cwaveprj,ffnl,&
128 1066105 : & idir,indlmn,istwf_k,kg_k,kpg,kpoint,lmnmax,mgfft,mpi_enreg,ndat,&
129 1066105 : & natom,nattyp,ngfft,nloalg,npw_k,nspinor,ntypat,&
130 1066105 : & phkxred,ph1d,ph3d,ucvol,useylm,is_kprime,gpu_option)
131 :
132 : !Arguments -------------------------------
133 : !scalars
134 : integer,intent(in) :: choice,cpopt,idir,istwf_k,lmnmax,ndat
135 : integer,intent(in) :: mgfft,natom,npw_k,nspinor,ntypat,useylm
136 : integer,intent(in),optional :: gpu_option
137 : real(dp),intent(in) :: ucvol
138 : logical,intent(in),optional :: is_kprime
139 : type(MPI_type),intent(in) :: mpi_enreg
140 : !arrays
141 : integer,intent(in) :: indlmn(6,lmnmax,ntypat),kg_k(3,npw_k),nattyp(ntypat)
142 : integer,intent(in) :: ngfft(18),nloalg(3)
143 : real(dp),intent(in) :: cwavef(2,npw_k*nspinor*ndat)
144 : real(dp),intent(in),target :: ffnl(:,:,:,:),kpg(:,:),ph3d(:,:,:)
145 : real(dp),intent(in) :: kpoint(3),ph1d(2,3*(2*mgfft+1)*natom),phkxred(2,natom)
146 : type(pawcprj_type),intent(inout) :: cwaveprj(natom,nspinor*ndat)
147 :
148 : !Local variables-------------------------------
149 : !scalars
150 : logical :: no_opernla_mv,no_opernla_gemm,l_is_kprime
151 : integer :: choice_,cplex,dimffnl,ia,ia1,ia2,ia3,ia4,iatm,ic,ii,ilmn,ishift,ispinor,itypat,idat,nprojs,shift,iatom,igrad
152 : integer :: jc,matblk,mincat,nd2gxdt,ndgxdt,nincat,nkpg,nkpg_,nlmn,signs,l_gpu_option
153 : !arrays
154 : real(dp) :: tsec(2)
155 1066105 : integer,allocatable :: cplex_dgxdt(:),cplex_d2gxdt(:),indlmn_typ(:,:)
156 1066105 : real(dp),allocatable :: d2gxdt(:,:,:,:,:),dgxdt(:,:,:,:,:),ffnl_typ(:,:,:)
157 1066105 : real(dp),allocatable :: gx(:,:,:,:)
158 1066105 : real(dp),allocatable :: vgx(:,:,:),vdgxdt(:,:,:)
159 1066105 : real(dp), pointer :: kpg_(:,:),ph3d_(:,:,:)
160 1066105 : real(dp), allocatable :: temp_realvec_r(:),temp_realvec_i(:)
161 : real(dp) :: d2gxdt_dum_in(1,1,1,1,1)
162 :
163 : ! *********************************************************************
164 :
165 : DBG_ENTER('COLL')
166 :
167 1066105 : call timab(1290,1,tsec)
168 :
169 : !Nothing to do in that case
170 1066105 : if (cpopt==1.and.choice==1) return
171 :
172 1066105 : l_gpu_option=ABI_GPU_DISABLED; if(present(gpu_option)) l_gpu_option = gpu_option
173 1066105 : l_is_kprime=.false.; if(present(is_kprime)) l_is_kprime = is_kprime
174 :
175 : !Not available for useylm=0
176 1066105 : if (useylm==0) then
177 0 : ABI_ERROR('Not available for useylm=0 !')
178 : end if
179 :
180 : !Error on bad choice
181 1066105 : if ((choice<1.or.choice>6).and.choice/=23.and.choice/=24) then
182 0 : ABI_BUG('Does not presently support this choice !')
183 : end if
184 :
185 : !Error on bad idir
186 1066105 : if (idir>0.and.choice/=2.and.choice/=3.and.choice/=5) then
187 0 : ABI_BUG('Does not support idir>0 for this choice')
188 : end if
189 :
190 : !Error on sizes
191 1066105 : nkpg=size(kpg,2)
192 1066105 : if (nkpg>0) then
193 : if( (choice==2.and.nkpg<3) .or. &
194 177384 : & ((choice==4.or.choice==24).and.nkpg<9) .or. &
195 : & ((choice==6.or.choice==3.or.choice==23).and.nkpg<3) ) then
196 0 : ABI_BUG('Incorrect size for kpg array !')
197 : end if
198 : end if
199 1066105 : if (size(ffnl,1)/=npw_k.or.size(ffnl,3)/=lmnmax) then
200 0 : ABI_BUG('Incorrect size for ffnl!')
201 : end if
202 1066105 : if (size(ph3d,dim=1)>0) then
203 1066105 : if (size(ph3d,2)/=npw_k) then
204 0 : ABI_BUG('Incorrect size for ph3d!')
205 : end if
206 : end if
207 :
208 1066105 : no_opernla_mv = nloalg(1)==4.or.nloalg(1)==8.or.nloalg(1)==10 ! have to be consistent with nonlop_ylm
209 : no_opernla_gemm = (.not. gemm_nonlop_use_gemm) &
210 : & .or. choice==4 .or. choice==6 .or. ndat==1 &
211 1066105 : & .or. l_gpu_option/=gemm_nonlop_gpu_option
212 :
213 : !Define dimensions of projected scalars
214 1066105 : dimffnl=size(ffnl,2)
215 1066105 : ndgxdt=0;nd2gxdt=0
216 1066105 : if (idir==0) then
217 919057 : if (choice==2) ndgxdt=3
218 919057 : if (choice==3) ndgxdt=6
219 919057 : if (choice==23) ndgxdt=9
220 919057 : if (choice==4) nd2gxdt=6
221 919057 : if (choice==24) then
222 0 : ndgxdt=3;nd2gxdt=6
223 : end if
224 919057 : if (choice==5) ndgxdt=3
225 919057 : if (choice==6) then
226 0 : ndgxdt=9;nd2gxdt=54
227 : end if
228 : else
229 147048 : ndgxdt=1
230 : end if
231 1066105 : if(cwaveprj(1,1)%ncpgr<ndgxdt+nd2gxdt) then
232 0 : ABI_BUG('Incorrect size for ncpgr')
233 : end if
234 : !Eventually re-compute (k+G) vectors (and related data)
235 1066105 : if (nkpg==0) then
236 888721 : nkpg_=0
237 888721 : if (choice==4.or.choice==24) nkpg_=9
238 888721 : if (choice==2.or.choice==3.or.choice==23) nkpg_=3
239 3554884 : ABI_MALLOC(kpg_,(npw_k,nkpg_))
240 888721 : if (nkpg_>0) then
241 16128 : call mkkpg(kg_k,kpg_,kpoint,nkpg_,npw_k)
242 : end if
243 : else
244 177384 : nkpg_=nkpg
245 177384 : kpg_ => kpg
246 : end if
247 :
248 : !Some other dims
249 2950572 : mincat=min(NLO_MINCAT,maxval(nattyp))
250 1066105 : cplex=2;if (istwf_k>1) cplex=1
251 1066105 : choice_=choice;if (cpopt==1) choice_=-choice
252 1066105 : signs=1;if (idir>0) signs=2
253 : !Eventually allocate temporary array for ph3d
254 1066105 : if (nloalg(2)<=0) then
255 0 : matblk=mincat
256 0 : ABI_MALLOC(ph3d_,(2,npw_k,matblk))
257 : else
258 1066105 : matblk=size(ph3d,3)
259 1066105 : ph3d_ => ph3d
260 : end if
261 :
262 1066105 : if(no_opernla_gemm) then
263 : !Loop over atom types
264 : ia1=1;iatm=0
265 2950572 : do itypat=1,ntypat
266 1884467 : ia2=ia1+nattyp(itypat)-1;if (ia2<ia1) cycle
267 17490577 : nlmn=count(indlmn(3,:,itypat)>0)
268 :
269 : ! Retrieve some data for this type of atom
270 5653401 : ABI_MALLOC(indlmn_typ,(6,nlmn))
271 9422335 : ABI_MALLOC(ffnl_typ,(npw_k,dimffnl,nlmn))
272 110771749 : indlmn_typ(:,1:nlmn)=indlmn(:,1:nlmn,itypat)
273 2948221065 : ffnl_typ(:,:,1:nlmn)=ffnl(:,:,1:nlmn,itypat)
274 :
275 : ! Loop on blocks of atoms inside type
276 1884467 : do ia3=ia1,ia2,mincat
277 1884467 : ia4=min(ia2,ia3+mincat-1);nincat=ia4-ia3+1
278 : ! Prepare the phase factors if they were not already computed
279 1884467 : if (nloalg(2)<=0) then
280 : call ph1d3d(ia3,ia4,kg_k,matblk,natom,npw_k,ngfft(1),ngfft(2),ngfft(3),&
281 0 : & phkxred,ph1d,ph3d_)
282 : end if
283 :
284 : ! Allocate memory for projected scalars
285 11306802 : ABI_MALLOC(gx,(cplex,nlmn,nincat,nspinor*ndat))
286 13191269 : ABI_MALLOC(dgxdt,(cplex,ndgxdt,nlmn,nincat,nspinor*ndat))
287 13191269 : ABI_MALLOC(d2gxdt,(cplex,nd2gxdt,nlmn,nincat,nspinor*ndat))
288 5653401 : ABI_MALLOC(cplex_dgxdt,(ndgxdt))
289 5653401 : ABI_MALLOC(cplex_d2gxdt,(nd2gxdt))
290 :
291 : ! Retrieve eventually <p_i|c> coeffs
292 1884467 : if (cpopt==1) then
293 429096 : do ispinor=1,nspinor*ndat
294 754824 : do ia=1,nincat
295 8544588 : gx(1:cplex,1:nlmn,ia,ispinor)=cwaveprj(iatm+ia,ispinor)%cp(1:cplex,1:nlmn)
296 : end do
297 : end do
298 : end if
299 :
300 : ! Compute <p_i|c> scalars (and derivatives) for this block of atoms
301 1884467 : if (abs(choice_)>1.or.no_opernla_mv) then
302 3768934 : do idat=1,ndat
303 1884467 : call timab(1291,1,tsec)
304 : call opernla_ylm(choice_,cplex,cplex_dgxdt,cplex_d2gxdt,dimffnl,&
305 : & d2gxdt(:,:,:,:,1+nspinor*(idat-1):nspinor*idat),&
306 : & dgxdt(:,:,:,:,1+nspinor*(idat-1):nspinor*idat),ffnl_typ,&
307 : & gx(:,:,:,1+nspinor*(idat-1):nspinor*idat),&
308 : & ia3,idir,indlmn_typ,istwf_k,kpg_,matblk,mpi_enreg,nd2gxdt,ndgxdt,nincat,nkpg_,nlmn,&
309 1884467 : & nloalg,npw_k,nspinor,ph3d_,signs,ucvol,cwavef(:,1+npw_k*nspinor*(idat-1):npw_k*nspinor*idat))
310 3768934 : call timab(1291,2,tsec)
311 : end do
312 : else
313 0 : call timab(1292,1,tsec)
314 0 : do idat=1,ndat
315 : call opernla_ylm_mv(choice_,cplex,dimffnl,ffnl_typ,gx(:,:,:,1+nspinor*(idat-1):nspinor*idat),&
316 : & ia3,indlmn_typ,istwf_k,matblk,mpi_enreg,nincat,nlmn,&
317 0 : & nloalg,npw_k,nspinor,ph3d_,ucvol,cwavef(:,1+npw_k*nspinor*(idat-1):npw_k*nspinor*idat))
318 : end do
319 0 : call timab(1292,2,tsec)
320 : end if
321 :
322 : ! Transfer result to output variable cwaveprj
323 1884467 : if (cpopt==0) then
324 3356298 : do ispinor=1,nspinor*ndat
325 5288946 : do ia=1,nincat
326 1932648 : cwaveprj(iatm+ia,ispinor)%nlmn=nlmn
327 50680410 : cwaveprj(iatm+ia,ispinor)%cp(1:cplex,1:nlmn)=gx(1:cplex,1:nlmn,ia,ispinor)
328 3619027 : if(cplex==1) cwaveprj(iatm+ia,ispinor)%cp(2,1:nlmn)=zero
329 : end do
330 : end do
331 : end if
332 1884467 : if (cpopt>=0.and.choice>1) then
333 460344 : ishift=0
334 460344 : if ((idir>0).and.(cwaveprj(1,1)%ncpgr>ndgxdt)) ishift=idir-1
335 460344 : if(cplex==2)then
336 927228 : do ispinor=1,nspinor*ndat
337 1518908 : do ia=1,nincat
338 : ! cwaveprj(iatm+ia,ispinor)%ncpgr=ndgxdt+nd2gxdt
339 591680 : if (ndgxdt>0) cwaveprj(iatm+ia,ispinor)%dcp(1:2,1+ishift:ndgxdt+ishift,1:nlmn)=&
340 35976576 : & dgxdt(1:2,1:ndgxdt,1:nlmn,ia,ispinor)
341 591680 : if (nd2gxdt>0)cwaveprj(iatm+ia,ispinor)%dcp(1:2,ndgxdt+1+ishift:ndgxdt+nd2gxdt+ishift,1:nlmn)=&
342 466884 : & d2gxdt(1:2,1:nd2gxdt,1:nlmn,ia,ispinor)
343 : end do
344 : end do
345 : else
346 : ! cplex_dgxdt(i) = 1 if dgxdt(1,i,:,:) is real, 2 if it is pure imaginary
347 : ! cplex_d2gxdt(i) = 1 if d2gxdt(1,i,:,:) is real, 2 if it is pure imaginary
348 0 : do ispinor=1,nspinor*ndat
349 0 : do ia=1,nincat
350 : ! cwaveprj(iatm+ia,ispinor)%ncpgr=ndgxdt+nd2gxdt
351 0 : if (ndgxdt>0) then
352 0 : do ilmn =1,nlmn
353 0 : do ii = 1,ndgxdt
354 0 : ic = cplex_dgxdt(ii) ; jc = 3 - ic
355 0 : cwaveprj(iatm+ia,ispinor)%dcp(ic,ii+ishift,ilmn)=dgxdt(1,ii,ilmn,ia,ispinor)
356 0 : cwaveprj(iatm+ia,ispinor)%dcp(jc,ii+ishift,ilmn)=zero
357 : end do
358 : end do
359 : end if
360 0 : if (nd2gxdt>0) then
361 0 : do ilmn =1,nlmn
362 0 : do ii = 1,nd2gxdt
363 0 : ic = cplex_d2gxdt(ii) ; jc = 3 - ic
364 0 : cwaveprj(iatm+ia,ispinor)%dcp(ic,ndgxdt+ii+ishift,ilmn)=d2gxdt(1,ii,ilmn,ia,ispinor)
365 0 : cwaveprj(iatm+ia,ispinor)%dcp(jc,ndgxdt+ii+ishift,ilmn)=zero
366 : end do
367 : end do
368 : end if
369 : end do
370 : end do
371 : end if
372 : end if
373 :
374 : ! End loop inside block of atoms
375 1884467 : iatm=iatm+nincat
376 1884467 : ABI_FREE(gx)
377 1884467 : ABI_FREE(dgxdt)
378 1884467 : ABI_FREE(d2gxdt)
379 1884467 : ABI_FREE(cplex_dgxdt)
380 1884467 : ABI_FREE(cplex_d2gxdt)
381 : end do
382 :
383 : ! End loop over atom types
384 1884467 : ia1=ia2+1
385 1884467 : ABI_FREE(indlmn_typ)
386 2950572 : ABI_FREE(ffnl_typ)
387 : end do
388 :
389 : else
390 :
391 : ! Batched GEMM call : proceed "ndat" bands at once with GEMM opernla
392 :
393 0 : if(cplex==1) ABI_BUG("toto")
394 0 : if (nloalg(2)<=0) ABI_BUG("toto")
395 :
396 0 : nprojs = 0
397 0 : do itypat=1,ntypat
398 0 : nprojs = nprojs + count(indlmn(3,:,itypat)>0)*nattyp(itypat)
399 : end do
400 :
401 0 : ABI_MALLOC(cplex_dgxdt,(ndgxdt))
402 0 : ABI_MALLOC(cplex_d2gxdt,(nd2gxdt))
403 0 : ABI_MALLOC(vgx,(cplex,nprojs,nspinor*ndat))
404 0 : ABI_MALLOC(vdgxdt,(cplex,ndgxdt*nprojs,nspinor*ndat))
405 0 : vgx(:,:,:) = zero
406 0 : vdgxdt(:,:,:) = zero
407 : #ifdef HAVE_OPENMP_OFFLOAD
408 : !$OMP TARGET ENTER DATA MAP(to:vgx,vdgxdt,kpg_) IF (l_gpu_option==ABI_GPU_OPENMP)
409 : #endif
410 :
411 : ! Retrieve eventually <p_i|c> coeffs
412 0 : if (cpopt==1) then
413 0 : do idat=1, ndat*nspinor
414 : shift = 0
415 0 : do iatom = 1, natom
416 0 : nlmn = cwaveprj(iatom, idat)%nlmn
417 0 : vgx(1:cplex, shift+1:shift+nlmn, idat) = cwaveprj(iatom, idat)%cp(1:cplex, 1:nlmn)
418 0 : shift = shift + nlmn
419 : end do
420 : end do
421 : end if
422 :
423 0 : if (cplex /= 2) then
424 0 : ABI_MALLOC(temp_realvec_r,(npw_k*nspinor*ndat))
425 0 : ABI_MALLOC(temp_realvec_i,(npw_k*nspinor*ndat))
426 : #ifdef HAVE_OPENMP_OFFLOAD
427 : !$OMP TARGET ENTER DATA MAP(alloc:temp_realvec_r,temp_realvec_i) IF (l_gpu_option==ABI_GPU_OPENMP)
428 : #endif
429 : end if
430 :
431 : call opernla_gemm(choice,cplex,cplex_dgxdt,cplex_d2gxdt,dimffnl,&
432 : & d2gxdt_dum_in,vdgxdt,ffnl,vgx,&
433 : & idir,indlmn,istwf_k,kpg_,matblk,mpi_enreg,nd2gxdt,ndgxdt,nkpg_,&
434 : & npw_k,nspinor,ph3d,signs,ucvol,ndat,ntypat,lmnmax,nattyp,l_is_kprime,&
435 : & -1,0,cpopt,&
436 : & nprojs,&
437 : & cwavef,&
438 : & temp_realvec_r,temp_realvec_i,&
439 0 : & l_gpu_option,.false.)
440 :
441 0 : if (cplex /= 2) then
442 : #ifdef HAVE_OPENMP_OFFLOAD
443 : !$OMP TARGET EXIT DATA MAP(delete:temp_realvec_r,temp_realvec_i) IF (l_gpu_option==ABI_GPU_OPENMP)
444 : #endif
445 0 : ABI_FREE(temp_realvec_r)
446 0 : ABI_FREE(temp_realvec_i)
447 : end if
448 :
449 : #ifdef HAVE_OPENMP_OFFLOAD
450 : !$OMP TARGET UPDATE FROM(vgx,vdgxdt) IF (l_gpu_option==ABI_GPU_OPENMP)
451 : #endif
452 :
453 :
454 : ! Transfer result to output variable cwaveprj
455 0 : if (cpopt==0) then
456 0 : do idat=1, ndat*nspinor
457 : shift = 0
458 0 : do iatom = 1, natom
459 0 : nlmn = cwaveprj(iatom, idat)%nlmn
460 0 : cwaveprj(iatom, idat)%cp(1:cplex, 1:nlmn) = vgx(1:cplex, shift+1:shift+nlmn, idat)
461 0 : shift = shift + nlmn
462 : end do
463 : end do
464 : end if
465 :
466 0 : if (cpopt>=0.and.choice>1) then
467 0 : ishift=0
468 : if ((idir>0).and.(cwaveprj(1,1)%ncpgr>ndgxdt)) ishift=idir-1
469 0 : if(cplex==2)then
470 0 : do idat=1, ndat*nspinor
471 : shift = 0
472 0 : do iatom = 1, natom
473 0 : nlmn = cwaveprj(iatom, idat)%nlmn
474 0 : do igrad=1,ndgxdt
475 : cwaveprj(iatom, idat)%dcp(1:cplex,igrad,1:nlmn) = &
476 0 : & vdgxdt(1:cplex, shift+1:shift+nlmn, idat)
477 0 : shift = shift + nlmn
478 : end do
479 : end do
480 : end do
481 : else
482 0 : ABI_BUG("toto")
483 : end if
484 : end if
485 :
486 :
487 : #ifdef HAVE_OPENMP_OFFLOAD
488 : !$OMP TARGET EXIT DATA MAP(delete:vgx,vdgxdt,kpg_) IF (l_gpu_option==ABI_GPU_OPENMP)
489 : #endif
490 0 : ABI_FREE(vgx)
491 0 : ABI_FREE(vdgxdt)
492 0 : ABI_FREE(cplex_dgxdt)
493 0 : ABI_FREE(cplex_d2gxdt)
494 :
495 : end if
496 :
497 1066105 : if (nkpg==0) then
498 888721 : ABI_FREE(kpg_)
499 : end if
500 1066105 : if (nloalg(2)<=0) then
501 0 : ABI_FREE(ph3d_)
502 : end if
503 :
504 1066105 : call timab(1290,2,tsec)
505 :
506 : DBG_EXIT('COLL')
507 :
508 2132210 : end subroutine getcprj
509 : !!***
510 :
511 : !!****f* ABINIT/ctocprj
512 : !! NAME
513 : !! ctocprj
514 : !!
515 : !! FUNCTION
516 : !! Compute all <Proj_i|Cnk> for every wave function |Cnk> expressed in reciprocal space.
517 : !! |Proj_i> are non-local projectors (for each atom and each l,m,n)
518 : !! Can also compute derivatives of <Proj_i|Cnk> wrt to several parameters
519 : !!
520 : !! INPUTS
521 : !! atindx(natom)=index table for atoms
522 : !! cg(2,mcg)=planewave coefficients of wavefunctions
523 : !! choice: chooses derivatives to compute:
524 : !! =1 => no derivatives
525 : !! =2 => 1st derivatives with respect to atomic position(s)
526 : !! =3 => 1st derivatives with respect to strain(s)
527 : !! =23=> 1st derivatives with respect to strain(s) and atm pos
528 : !! =4 => 2nd derivatives with respect to atomic pos.
529 : !! =24=> 1st and 2nd derivatives with respect to atomic pos.
530 : !! =5 => derivatives with respect to k wavevector
531 : !! =6 => 2nd derivatives with respect to strain and atm. pos.
532 : !! gmet(3,3)=reciprocal space metric tensor in bohr**-2
533 : !! gprimd(3,3)=dimensional reciprocal space primitive translations
534 : !! iatom= if <=0, cprj=<p_i|Cnk> are computed for all atoms 1...natom
535 : !! if >0 cprj=<p_i|Cnk> are computed only for atom with index iatom
536 : !! idir=direction of the derivative, i.e. dir. of - atom to be moved in the case choice=2
537 : !! - strain component in the case choice=3
538 : !! - k point direction in the case choice=5
539 : !! Compatible only with choice=2,3,5; if idir=0, all derivatives are computed
540 : !! iorder_cprj=0 if output cprj=<p_i|Cnk> are sorted by atom type
541 : !! (first all elements of atom type 1, followed by those of atom type 2 and so on).
542 : !! 1 if output cprj=<p_i|Cnk> are sorted according to
543 : !! the variable typat in the main input file
544 : !! istwfk(nkpt)=option parameter that describes the storage of wfs
545 : !! kg(3,mpw*mkmem)=reduced planewave coordinates
546 : !! kpt(3,nkpt)=reduced coordinates of k points.
547 : !! mcg=size of wave-functions array (cg) =mpw*nspinor*mband_mem*mkmem*nsppol
548 : !! mcprj=size of projected wave-functions array (cprj) =nspinor*mband*mkmem*nsppol
549 : !! mgfft=maximum size of 1D FFTs
550 : !! mkmem=number of k points treated by this node.
551 : !! mpi_enreg=information about MPI parallelization
552 : !! mpsang=1+maximum angular momentum for nonlocal pseudopotentials
553 : !! mpw=maximum dimensioned size of npw
554 : !! my_nsppol=number of spin components in memory for current MPI process
555 : !! natom=number of atoms in cell
556 : !! nattyp(ntypat)= # atoms of each type
557 : !! nband(nkpt*nsppol)=number of bands at this k point for that spin polarization
558 : !! mband_mem=max number of bands for this processor (in case of band parallelism)
559 : !! ncprj=1st dim. of cprj array (natom if iatom<=0, 1 if iatom>0)
560 : !! ngfft(18)=contain all needed information about 3D FFT, see ~ABINIT/Infos/vargs.htm#ngfft
561 : !! nkpt=number of k points
562 : !! nloalg(3)=governs the choice of the algorithm for nonlocal operator
563 : !! npwarr(nkpt)=number of planewaves in basis at this k point
564 : !! nspinor=number of spinorial components of the wavefunctions
565 : !! nsppol=1 for unpolarized, 2 for spin-polarized
566 : !! ntypat=number of types of atoms in unit cell
567 : !! paral_kgb= 1 if kpt-band-FFT is activated
568 : !! ph1d(2,3*(2*mgfft+1)*natom)=1-dim structure factor phase information
569 : !! psps <type(pseudopotential_type)>=variables related to pseudopotentials
570 : !! rmet(3,3)=real space metric (bohr**2)
571 : !! tim_ctocprj=timing code of the calling routine
572 : !! typat(natom)= types of atoms
573 : !! uncp=unit number for <P_lmn|Cnk> data (if used)
574 : !! xred(3,natom)=reduced dimensionless atomic coordinates
575 : !! ylm(mpw*mkmem,mpsang*mpsang)=real spherical harmonics for each G and k point
576 : !! ylmgr(npw*mkmem,nylmgr,mpsang*mpsang*useylmgr)=gradients of real spherical harmonics wrt (k+G)
577 : !!!
578 : !! OUTPUT
579 : !! cprj(ncprj,mcprj) <type(pawcprj_type)>= projected input wave functions <Proj_i|Cnk> with NL projectors
580 : !! Usually ncprj=natom
581 : !!
582 : !! SOURCE
583 :
584 4458 : subroutine ctocprj(atindx,cg,choice,cprj,gmet,gprimd,iatom,idir,&
585 2229 : & iorder_cprj,istwfk,kg,kpt,mcg,mcprj,mgfft,mkmem,mpi_enreg,mpsang,&
586 2229 : & mpw,natom,nattyp,nband,ncprj,ngfft,nkpt,nloalg,npwarr,nspinor,&
587 2229 : & nsppol,my_nsppol,ntypat,paral_kgb,ph1d,psps,rmet,typat,ucvol,uncp,xred,ylm,ylmgr)
588 :
589 : !Arguments -------------------------------
590 : !scalars
591 : integer,intent(in) :: choice,iatom,idir,iorder_cprj,mcg,mcprj,mgfft,mkmem,mpsang,mpw
592 : integer,intent(in) :: my_nsppol,natom,ncprj,nkpt,nspinor,nsppol,ntypat,paral_kgb,uncp
593 : !TODO : distribute cprj over bands as well
594 : real(dp),intent(in) :: ucvol
595 : type(MPI_type),intent(in) :: mpi_enreg
596 : type(pseudopotential_type),target,intent(in) :: psps
597 : !arrays
598 : integer,intent(in) :: istwfk(nkpt),nband(nkpt*nsppol)
599 : integer,intent(in) :: ngfft(18),nloalg(3),npwarr(nkpt),kg(3,mpw*mkmem),typat(natom)
600 : integer,intent(in),target :: atindx(natom),nattyp(ntypat)
601 : real(dp),intent(in) :: cg(2,mcg)
602 : real(dp),intent(in) :: gmet(3,3),gprimd(3,3),kpt(3,nkpt),rmet(3,3)
603 : real(dp),intent(in) :: xred(3,natom),ylm(:,:),ylmgr(:,:,:)
604 : real(dp),intent(in),target :: ph1d(2,3*(2*mgfft+1)*natom)
605 : type(pawcprj_type),intent(inout) :: cprj(ncprj,mcprj)
606 :
607 : !Local variables-------------------------------
608 : !scalars
609 : integer :: blocksz,cg_bandpp,counter,cpopt,cprj_bandpp,dimffnl,ia,iatm,iatom1,iatom2
610 : integer :: iband_max,iband_min,iband_start,ibg,ibgb,iblockbd,ibp,icg,icgb,icp1,icp2
611 : integer :: ider,idir0,iend,ierr,ig,ii,ikg,ikpt,ilm,ipw,isize,isppol,istart,istwf_k,itypat,iwf1,iwf2,jdir
612 : integer :: matblk,me_distrb,my_nspinor,n1,n1_2p1,n2,n2_2p1,n3,n3_2p1,kk,nlmn
613 : integer :: mband,mband_cg,mband_cprj,npband_dfpt
614 : integer :: nband_k,nband_cprj_k,nblockbd,ncpgr,nkpg,npband_bandfft,npws,npw_k,npw_nk,ntypat0
615 : integer :: nband_cg_k
616 : integer :: shift1,shift1b,shift2,shift2b,shift3,shift3b
617 : integer :: spaceComm,spaceComm_band,spaceComm_fft,useylmgr
618 : logical :: cg_band_distributed,cprj_band_distributed,one_atom
619 : real(dp) :: arg
620 : character(len=500) :: msg
621 : !arrays
622 2229 : integer,allocatable :: bufsize(:),bufsize_wf(:),bufdisp(:),bufdisp_wf(:)
623 2229 : integer,allocatable :: dimlmn(:),kg_k(:,:),kg_k_loc(:,:)
624 2229 : integer,allocatable :: npw_block(:),npw_disp(:)
625 2229 : integer,pointer :: atindx_atm(:),indlmn_atm(:,:,:),nattyp_atm(:),pspso_atm(:)
626 : real(dp) :: kpoint(3),work(6),tsec(2)
627 2229 : real(dp),allocatable :: cwavef(:,:),cwavef_tmp(:,:)
628 2229 : real(dp),allocatable :: ffnl(:,:,:,:),ffnl_npw(:,:,:,:),ffnl_tmp(:,:,:,:),ffnl_tmp_npw(:,:,:,:)
629 2229 : real(dp),allocatable :: kpg_k(:,:)
630 2229 : real(dp),allocatable :: ph3d(:,:,:),ph3d_npw(:,:,:),ph3d_tmp(:,:,:),ph3d_tmp_npw(:,:,:)
631 2229 : real(dp),allocatable :: phkxred(:,:),ylm_k(:,:),ylmgr_k(:,:,:)
632 2229 : real(dp),ABI_CONTIGUOUS pointer :: ekb_atm(:,:),ffspl_atm(:,:,:,:),ph1d_atm(:,:)
633 2229 : type(pawcprj_type),allocatable :: cwaveprj(:,:)
634 :
635 : ! *********************************************************************
636 :
637 : DBG_ENTER('COLL')
638 :
639 : !Nothing to do if current MPI process does treat kpoints or plane-waves
640 2229 : if (mcg==0.or.mcprj==0) return
641 :
642 : !Preliminary tests
643 2229 : if (psps%useylm==0) then
644 0 : msg='Not available for useylm=0!'
645 0 : ABI_ERROR(msg)
646 : end if
647 2229 : if ((choice<1.or.choice>6).and.choice/=23.and.choice/=24) then
648 0 : msg='Bad choice!'
649 0 : ABI_BUG(msg)
650 : end if
651 2229 : if (idir>0.and.choice/=2.and.choice/=3.and.choice/=5) then
652 0 : msg='Does not support idir>0 for that choice!'
653 0 : ABI_BUG(msg)
654 : end if
655 6687 : if (size(ylm)/=mpw*mkmem*mpsang*mpsang) then
656 0 : msg='Wrong size for Ylm!'
657 0 : ABI_BUG(msg)
658 : end if
659 8916 : useylmgr=merge(0,1,(size(ylmgr)==0))
660 158 : if (useylmgr==0.and.(choice==3.or.choice==5.or.choice==23)) then
661 0 : msg=' Ylm gradients have to be in memory for choice=3, 5, or 23!'
662 0 : ABI_BUG(msg)
663 : end if
664 :
665 : !Init parallelism
666 2229 : npband_dfpt = 1
667 2229 : if (paral_kgb==1) then
668 832 : me_distrb=mpi_enreg%me_kpt
669 832 : spaceComm=mpi_enreg%comm_kpt
670 832 : spaceComm_fft=mpi_enreg%comm_fft
671 832 : npband_bandfft=mpi_enreg%nproc_band
672 832 : cg_bandpp=mpi_enreg%bandpp
673 832 : cprj_bandpp=mpi_enreg%bandpp
674 832 : spaceComm_band=mpi_enreg%comm_band
675 832 : cg_band_distributed=.true.
676 832 : cprj_band_distributed=(mpi_enreg%nproc_band>1)
677 832 : npband_dfpt=1
678 : else
679 1397 : me_distrb=mpi_enreg%me_kpt
680 1397 : spaceComm=mpi_enreg%comm_cell
681 1397 : spaceComm_fft=xmpi_comm_self
682 1397 : npband_bandfft=1
683 1397 : cg_bandpp=1
684 1397 : cprj_bandpp=1
685 1397 : if (mpi_enreg%paralbd==1) then
686 415 : spaceComm_band=mpi_enreg%comm_band ! not actually used as npband_bandfft=1
687 415 : cg_band_distributed=.true.
688 415 : cprj_band_distributed=.true.
689 : !npband_dfpt=1
690 415 : npband_dfpt=mpi_enreg%nproc_band
691 : else
692 982 : spaceComm_band=xmpi_comm_self
693 982 : cg_band_distributed=.false.
694 982 : cprj_band_distributed=.false.
695 : end if
696 : end if
697 2229 : if (cg_bandpp/=cprj_bandpp) then
698 0 : ABI_BUG('cg_bandpp must be equal to cprj_bandpp!')
699 : end if
700 :
701 : !Manage parallelization over bands
702 87150 : mband=maxval(nband(1:nkpt*my_nsppol))
703 2229 : mband_cg=mband/npband_bandfft/npband_dfpt
704 2229 : mband_cprj=mband_cg
705 :
706 : !Manage parallelization over spinors
707 2229 : my_nspinor=max(1,nspinor/mpi_enreg%nproc_spinor)
708 :
709 : !Check sizes
710 2229 : if (mpw*mband_cg*my_nspinor*mkmem*my_nsppol>mcg) then
711 0 : ABI_BUG(sjoin('Bad mcg value!', itoa(mcg)))
712 : end if
713 2229 : if (mband_cprj*my_nspinor*mkmem*my_nsppol>mcprj) then
714 0 : ABI_BUG(sjoin('Bad mcprj value!', itoa(mcprj)))
715 : end if
716 :
717 : !Check sizes for cprj (distribution is tricky)
718 2229 : one_atom=(iatom>0)
719 2229 : if (one_atom.and.ncprj/=1) then
720 0 : ABI_BUG('Bad value for ncprj dimension (should be 1) !')
721 : end if
722 2229 : if (.not.one_atom.and.ncprj/=natom) then
723 0 : ABI_BUG('Bad value for ncprj dimension (should be natom) !')
724 : end if
725 :
726 : !Initialize some variables
727 2229 : n1=ngfft(1);n2=ngfft(2);n3=ngfft(3)
728 2229 : n1_2p1=2*n1+1;n2_2p1=2*n2+1;n3_2p1=2*n3+1
729 2229 : ibg=0;icg=0;cpopt=0
730 2229 : ider=0;idir0=0;istart=idir;iend=idir
731 2229 : if (choice==3.or.choice==5.or.choice==23) ider=1
732 2229 : if (idir>0) then
733 156 : if (choice==3) idir0=-idir
734 156 : if (choice==5) idir0=idir
735 : else
736 : ! if (choice==23) idir0=-7
737 2073 : if (choice==3) idir0=-7
738 2073 : if (choice==5) idir0=4
739 : end if
740 2229 : if (idir0==0.or.idir0==4) then
741 2057 : dimffnl=1+3*ider
742 172 : else if (idir0/=-7) then
743 156 : dimffnl=1+ider
744 : else
745 16 : dimffnl=1+6*ider
746 16 : if(choice==3)then
747 2229 : istart=ider;iend=6*ider
748 : end if
749 : end if
750 2229 : nkpg=0
751 2229 : if (choice==3.or.choice==2.or.choice==23) nkpg=3*nloalg(3)
752 2229 : if (choice==4.or.choice==24) nkpg=9*nloalg(3)
753 :
754 : !Set number of gradients for <p_i|Cnk>
755 2229 : ncpgr=0
756 2229 : if (idir==0) then
757 2073 : if (choice==2) ncpgr=3
758 2073 : if (choice==3) ncpgr=6
759 2073 : if (choice==23)ncpgr=9
760 2073 : if (choice==4) ncpgr=6
761 2073 : if (choice==24)ncpgr=9
762 2073 : if (choice==5) ncpgr=3
763 2073 : if (choice==6) ncpgr=63
764 : else
765 156 : ncpgr=1
766 : end if
767 : !Test cprj gradients dimension (just to be sure)
768 2229 : if (cprj(1,1)%ncpgr/=ncpgr) then
769 0 : ABI_BUG('cprj are badly allocated !')
770 : end if
771 :
772 :
773 : !Extract data for treated atom(s)
774 2229 : if (one_atom) then
775 0 : iatom1=iatom;iatom2=iatom
776 0 : ntypat0=1;itypat=typat(iatom)
777 0 : ABI_MALLOC(nattyp_atm,(ntypat0))
778 0 : nattyp_atm(1)=1
779 0 : ABI_MALLOC(atindx_atm,(ntypat0))
780 0 : atindx_atm(1)=atindx(iatom)
781 0 : ABI_MALLOC(ph1d_atm,(2,(n1_2p1+n2_2p1+n3_2p1)*ntypat0))
782 0 : shift1=(atindx(iatom)-1)*n1_2p1
783 0 : shift2=(atindx(iatom)-1)*n2_2p1+natom*n1_2p1
784 0 : shift3=(atindx(iatom)-1)*n3_2p1+natom*(n1_2p1+n2_2p1)
785 0 : shift1b=0;shift2b=n1_2p1;shift3b=n1_2p1+n2_2p1
786 0 : ph1d_atm(:,shift1b+1:shift1b+n1_2p1)=ph1d(:,shift1+1:shift1+n1_2p1)
787 0 : ph1d_atm(:,shift2b+1:shift2b+n2_2p1)=ph1d(:,shift2+1:shift2+n2_2p1)
788 0 : ph1d_atm(:,shift3b+1:shift3b+n3_2p1)=ph1d(:,shift3+1:shift3+n3_2p1)
789 0 : ABI_MALLOC(ekb_atm,(psps%dimekb,ntypat0))
790 0 : ABI_MALLOC(indlmn_atm,(6,psps%lmnmax,ntypat0))
791 0 : ABI_MALLOC(ffspl_atm,(psps%mqgrid_ff,2,psps%lnmax,ntypat0))
792 0 : ABI_MALLOC(pspso_atm,(ntypat0))
793 0 : ekb_atm(:,1)=psps%ekb(:,itypat)
794 0 : indlmn_atm(:,:,1)=psps%indlmn(:,:,itypat)
795 0 : ffspl_atm(:,:,:,1)=psps%ffspl(:,:,:,itypat)
796 0 : pspso_atm(1)=psps%pspso(itypat)
797 : else
798 2229 : iatom1=1;iatom2=natom
799 2229 : ntypat0=ntypat
800 2229 : atindx_atm => atindx
801 2229 : nattyp_atm => nattyp
802 2229 : ph1d_atm => ph1d
803 2229 : ekb_atm => psps%ekb
804 2229 : indlmn_atm => psps%indlmn
805 2229 : ffspl_atm => psps%ffspl
806 2229 : pspso_atm => psps%pspso
807 : end if
808 :
809 : !Dimensioning and allocation of <p_i|Cnk>
810 6687 : ABI_MALLOC(dimlmn,(ncprj))
811 7139 : dimlmn=0 ! Type-sorted cprj
812 2229 : if (one_atom) then
813 0 : itypat=typat(iatom)
814 0 : dimlmn(1:nattyp(itypat))=count(indlmn_atm(3,:,itypat)>0)
815 : else
816 2229 : ia=0
817 4957 : do itypat=1,ntypat0
818 30779 : dimlmn(ia+1:ia+nattyp(itypat))=count(indlmn_atm(3,:,itypat)>0)
819 4957 : ia=ia+nattyp(itypat)
820 : end do
821 : end if
822 23647 : ABI_MALLOC(cwaveprj,(ncprj,my_nspinor*cprj_bandpp))
823 2229 : call pawcprj_alloc(cwaveprj,ncpgr,dimlmn)
824 :
825 : !Additional statements if band-fft parallelism
826 2229 : if (npband_bandfft>1) then
827 2496 : ABI_MALLOC(npw_block,(npband_bandfft))
828 1664 : ABI_MALLOC(npw_disp,(npband_bandfft))
829 : !FB ABI_MALLOC(bufsize,(npband_bandfft*cg_bandpp))
830 1664 : ABI_MALLOC(bufsize,(npband_bandfft))
831 : !FB ABI_MALLOC(bufdisp,(npband_bandfft*cg_bandpp))
832 1664 : ABI_MALLOC(bufdisp,(npband_bandfft))
833 : !FB ABI_MALLOC(bufsize_wf,(npband_bandfft*cg_bandpp))
834 1664 : ABI_MALLOC(bufsize_wf,(npband_bandfft))
835 : !FB ABI_MALLOC(bufdisp_wf,(npband_bandfft*cg_bandpp))
836 1664 : ABI_MALLOC(bufdisp_wf,(npband_bandfft))
837 : end if
838 :
839 : !Set output datastructure to zero
840 2229 : call pawcprj_set_zero(cprj)
841 :
842 : !LOOP OVER SPINS
843 4551 : do isppol=1,my_nsppol
844 2322 : ikg=0
845 :
846 : ! BIG FAT k POINT LOOP
847 89472 : do ikpt=1,nkpt
848 84921 : counter=100*ikpt+isppol
849 :
850 : ! Select k point to be treated by this proc
851 84921 : nband_k=nband(ikpt+(isppol-1)*nkpt)
852 84921 : if (proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,isppol,me_distrb)) cycle
853 :
854 : ! Retrieve k-point
855 212768 : kpoint(:)=kpt(:,ikpt)
856 53192 : istwf_k=istwfk(ikpt)
857 :
858 : ! Retrieve number of plane waves
859 53192 : npw_k=npwarr(ikpt)
860 53192 : if (npband_bandfft>1) then
861 : ! Special treatment for band-fft //
862 6656 : call xmpi_allgather(npw_k,npw_block,spaceComm_band,ierr)
863 19968 : npw_nk=sum(npw_block);npw_disp(1)=0
864 13312 : do ii=2,npband_bandfft
865 13312 : npw_disp(ii)=npw_disp(ii-1)+npw_block(ii-1)
866 : end do
867 : else
868 46536 : npw_nk=npw_k
869 : end if
870 :
871 : ! Retrieve (k+G) points and spherical harmonics
872 212768 : ABI_MALLOC(ylm_k,(npw_k,mpsang*mpsang))
873 212768 : ABI_MALLOC(ylmgr_k,(npw_k,3,mpsang*mpsang*useylmgr))
874 159576 : ABI_MALLOC(kg_k,(3,npw_nk))
875 53192 : if (npband_bandfft>1) then
876 : ! Special treatment for band-fft //
877 19968 : ABI_MALLOC(kg_k_loc,(3,npw_k))
878 1953536 : kg_k_loc(:,1:npw_k)=kg(:,1+ikg:npw_k+ikg)
879 33280 : bufsize(:)=3*npw_block(:);bufdisp(:)=3*npw_disp(:)
880 6656 : call xmpi_allgatherv(kg_k_loc,3*npw_k,kg_k,bufsize,bufdisp,spaceComm_band,ierr)
881 : else
882 30547868 : kg_k(:,1:npw_k)=kg(:,1+ikg:npw_k+ikg)
883 : end if
884 275330 : do ilm=1,mpsang*mpsang
885 34024285 : ylm_k(1:npw_k,ilm)=ylm(1+ikg:npw_k+ikg,ilm)
886 88615205 : if (useylmgr>0) ylmgr_k(1:npw_k,1:3,ilm)=ylmgr(1+ikg:npw_k+ikg,1:3,ilm)
887 : end do
888 :
889 : ! Compute (k+G) vectors
890 212768 : ABI_MALLOC(kpg_k,(npw_nk,nkpg))
891 53192 : if (nkpg>0) then
892 12174 : call mkkpg(kg_k,kpg_k,kpoint,nkpg,npw_nk)
893 : end if
894 : ! Allocate and compute the arrays phkxred and ph3d
895 159576 : ABI_MALLOC(phkxred,(2,ncprj))
896 158906 : do ia=iatom1,iatom2
897 105714 : iatm=min(atindx_atm(ia),ncprj)
898 105714 : arg=two_pi*(kpoint(1)*xred(1,ia)+kpoint(2)*xred(2,ia)+kpoint(3)*xred(3,ia))
899 158906 : phkxred(1,iatm)=cos(arg);phkxred(2,iatm)=sin(arg)
900 : end do
901 53192 : matblk=ncprj;if (nloalg(2)<=0) matblk=0
902 212768 : ABI_MALLOC(ph3d,(2,npw_nk,matblk))
903 53192 : if (matblk>0)then
904 : ! Here, precomputation of ph3d
905 53192 : if (npband_bandfft>1) then
906 : ! Special treatment for band-fft //
907 26624 : ABI_MALLOC(ph3d_tmp,(2,npw_k,matblk))
908 6656 : call ph1d3d(1,ncprj,kg_k_loc,matblk,ncprj,npw_k,n1,n2,n3,phkxred,ph1d_atm,ph3d_tmp)
909 26624 : ABI_MALLOC(ph3d_tmp_npw,(2,matblk,npw_k))
910 26624 : ABI_MALLOC(ph3d_npw,(2,matblk,npw_nk))
911 33280 : isize=2*matblk;bufsize(:)=isize*npw_block(:);bufdisp(:)=isize*npw_disp(:)
912 493376 : do ipw=1,npw_k
913 3413696 : ph3d_tmp_npw(:,:,ipw)=ph3d_tmp(:,ipw,:)
914 : end do
915 6656 : call xmpi_allgatherv(ph3d_tmp_npw,isize*npw_k,ph3d_npw,bufsize,bufdisp,spaceComm_band,ierr)
916 980096 : do ipw=1,npw_nk
917 6820736 : ph3d(:,ipw,:)=ph3d_npw(:,:,ipw)
918 : end do
919 6656 : ABI_FREE(ph3d_npw)
920 6656 : ABI_FREE(ph3d_tmp_npw)
921 6656 : ABI_FREE(ph3d_tmp)
922 : else
923 46536 : call ph1d3d(1,ncprj,kg_k,matblk,ncprj,npw_k,n1,n2,n3,phkxred,ph1d_atm,ph3d)
924 : end if
925 0 : else if (npband_bandfft>1) then
926 0 : ABI_ERROR('Band-fft parallelism +nloag(1)<0 forbidden !')
927 : end if
928 :
929 : ! Compute nonlocal form factors ffnl at all (k+G)
930 319152 : ABI_MALLOC(ffnl,(npw_nk,dimffnl,psps%lmnmax,ntypat0))
931 53192 : if (npband_bandfft>1) then
932 : ! Special treatment for band-fft //
933 39936 : ABI_MALLOC(ffnl_tmp,(npw_k,dimffnl,psps%lmnmax,ntypat0))
934 : call mkffnl(psps%dimekb,dimffnl,ekb_atm,ffnl_tmp,ffspl_atm,&
935 : & gmet,gprimd,ider,idir0,indlmn_atm,kg_k_loc,kpg_k,kpoint,psps%lmnmax,&
936 : & psps%lnmax,psps%mpsang,psps%mqgrid_ff,nkpg,npw_k,ntypat0,&
937 6656 : & pspso_atm,psps%qgrid_ff,rmet,psps%usepaw,psps%useylm,ylm_k,ylmgr_k)
938 39936 : ABI_MALLOC(ffnl_tmp_npw,(dimffnl,psps%lmnmax,ntypat0,npw_k))
939 39936 : ABI_MALLOC(ffnl_npw,(dimffnl,psps%lmnmax,ntypat0,npw_nk))
940 6656 : isize=dimffnl*psps%lmnmax*ntypat0
941 33280 : bufsize(:)=isize*npw_block(:);bufdisp(:)=isize*npw_disp(:)
942 493376 : do ipw=1,npw_k
943 8767616 : ffnl_tmp_npw(:,:,:,ipw)=ffnl_tmp(ipw,:,:,:)
944 : end do
945 6656 : call xmpi_allgatherv(ffnl_tmp_npw,isize*npw_k,ffnl_npw,bufsize,bufdisp,spaceComm_band,ierr)
946 980096 : do ipw=1,npw_nk
947 17528576 : ffnl(ipw,:,:,:)=ffnl_npw(:,:,:,ipw)
948 : end do
949 6656 : ABI_FREE(ffnl_npw)
950 6656 : ABI_FREE(ffnl_tmp_npw)
951 6656 : ABI_FREE(ffnl_tmp)
952 : else
953 : call mkffnl(psps%dimekb,dimffnl,ekb_atm,ffnl,ffspl_atm,&
954 : & gmet,gprimd,ider,idir0,indlmn_atm,kg_k,kpg_k,kpoint,psps%lmnmax,&
955 : & psps%lnmax,psps%mpsang,psps%mqgrid_ff,nkpg,npw_k,ntypat0,&
956 46536 : & pspso_atm,psps%qgrid_ff,rmet,psps%usepaw,psps%useylm,ylm_k,ylmgr_k)
957 : end if
958 :
959 : ! No more need of kg_g_tmp
960 53192 : if (npband_bandfft>1) then
961 6656 : ABI_FREE(kg_k_loc)
962 : end if
963 :
964 : ! Allocate arrays for a wave-function (or a block of WFs)
965 159576 : ABI_MALLOC(cwavef,(2,npw_nk*my_nspinor*cg_bandpp))
966 53192 : if (npband_bandfft>1) then
967 33280 : isize=2*my_nspinor*cg_bandpp;bufsize(:)=isize*npw_block(:);bufdisp(:)=isize*npw_disp(:)
968 19968 : isize=2*my_nspinor*npw_k*cg_bandpp;bufsize_wf(:)=isize
969 : !FB do ii=1,npband_bandfft*cg_bandpp
970 19968 : do ii=1,npband_bandfft
971 19968 : bufdisp_wf(ii)=(ii-1)*isize
972 : end do
973 : end if
974 :
975 : ! Loop over bands or blocks of bands
976 53192 : icgb=icg ; ibgb=ibg ; iband_start=1
977 53192 : blocksz=npband_bandfft*cg_bandpp
978 53192 : nblockbd=nband_k/blocksz
979 53192 : nband_cprj_k=merge(nband_k/npband_bandfft/npband_dfpt,nband_k,cprj_band_distributed)
980 : nband_cg_k=merge(nband_k/npband_bandfft/npband_dfpt,nband_k,cg_band_distributed)
981 377854 : do iblockbd=1,nblockbd
982 324662 : iband_min=1+(iblockbd-1)*blocksz
983 324662 : iband_max=iblockbd*blocksz
984 :
985 324662 : if (proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,iband_min,iband_max,isppol,me_distrb)) then
986 216 : if (.not.cg_band_distributed) icgb=icgb+npw_k*my_nspinor*blocksz
987 216 : if (.not.cprj_band_distributed) ibgb=ibgb+my_nspinor*blocksz
988 : cycle
989 : end if
990 :
991 : ! Extract wavefunction information
992 : ! Special treatment for band-fft parallelism
993 324446 : if (npband_bandfft>1) then
994 : !Transpose WF to get them in "FFT" representation
995 19968 : ABI_MALLOC(cwavef_tmp,(2,npw_k*my_nspinor*blocksz))
996 3900416 : do ig=1,npw_k*my_nspinor*blocksz
997 3893760 : cwavef_tmp(1,ig)=cg(1,ig+icgb)
998 3900416 : cwavef_tmp(2,ig)=cg(2,ig+icgb)
999 : end do
1000 6656 : call xmpi_alltoallv(cwavef_tmp,bufsize_wf,bufdisp_wf,cwavef,bufsize,bufdisp,spaceComm_band,ierr)
1001 6656 : ABI_FREE(cwavef_tmp)
1002 : !Reorder WF according to cg_bandpp and/or spinor
1003 13312 : if (cg_bandpp>1.or.my_nspinor>1) then
1004 : !FB ABI_MALLOC(cwavef_tmp,(2,npw_nk*my_nspinor*blocksz))
1005 13312 : ABI_MALLOC(cwavef_tmp,(2,npw_nk*my_nspinor*cg_bandpp))
1006 : !FB do ig=1,npw_nk*my_nspinor*blocksz
1007 3900416 : do ig=1,npw_nk*my_nspinor*cg_bandpp
1008 11687936 : cwavef_tmp(:,ig)=cwavef(:,ig)
1009 : end do
1010 : shift1=0
1011 33280 : do iwf2=1,cg_bandpp
1012 59904 : do ig=1,my_nspinor
1013 : shift2=0
1014 106496 : do iwf1=1,npband_bandfft
1015 53248 : npws=npw_block(iwf1)
1016 53248 : ipw=shift2+(iwf2-1)*my_nspinor*npws+(ig-1)*npws
1017 11734528 : cwavef(:,shift1+1:shift1+npws)=cwavef_tmp(:,ipw+1:ipw+npws)
1018 79872 : shift1=shift1+npws ; shift2=shift2+cg_bandpp*my_nspinor*npws
1019 : end do
1020 : end do
1021 : end do
1022 6656 : ABI_FREE(cwavef_tmp)
1023 : end if
1024 : else
1025 60986987 : do ig=1,npw_k*my_nspinor*cg_bandpp
1026 60669197 : cwavef(1,ig)=cg(1,ig+icgb)
1027 60986987 : cwavef(2,ig)=cg(2,ig+icgb)
1028 : end do
1029 : end if
1030 :
1031 : ! Compute scalar product of wavefunction with all NL projectors
1032 668860 : do ibp=1,cg_bandpp ! Note: we suppose cp_bandpp=cprj_bandpp
1033 344414 : iwf1=1+(ibp-1)*npw_nk*my_nspinor;iwf2=ibp*npw_nk*my_nspinor
1034 344414 : icp1=1+(ibp-1)*my_nspinor;icp2=ibp*my_nspinor
1035 1013704 : do jdir=istart,iend
1036 344844 : call timab(1294,1,tsec)
1037 : call getcprj(choice,cpopt,cwavef(:,iwf1:iwf2),cwaveprj(:,icp1:icp2),&
1038 : & ffnl,jdir,indlmn_atm,istwf_k,kg_k,kpg_k,kpoint,psps%lmnmax,&
1039 : & mgfft,mpi_enreg,1,ncprj,nattyp_atm,ngfft,nloalg,&
1040 344844 : & npw_nk,my_nspinor,ntypat0,phkxred,ph1d_atm,ph3d,ucvol,psps%useylm,.false.)
1041 689258 : call timab(1294,2,tsec)
1042 : end do
1043 : end do
1044 : ! Export cwaveprj to big array cprj
1045 : call pawcprj_put(atindx_atm,cwaveprj,cprj,ncprj,iband_start,ibgb,ikpt,iorder_cprj,isppol,&
1046 : & mband_cprj,mkmem,natom,cprj_bandpp,nband_cprj_k,dimlmn,my_nspinor,nsppol,uncp,&
1047 324446 : & mpi_comm_band=spaceComm_band,to_be_gathered=(cg_band_distributed.and.(.not.cprj_band_distributed)))
1048 :
1049 324446 : iband_start=iband_start+merge(cg_bandpp,blocksz,cprj_band_distributed)
1050 :
1051 : ! End loop over bands
1052 377638 : icgb=icgb+npw_k*my_nspinor*blocksz
1053 : end do
1054 :
1055 : ! Shift array memory (if mkmem/=0)
1056 53192 : if (mkmem/=0) then
1057 53192 : ibg=ibg+my_nspinor*nband_cprj_k
1058 : !FB icg=icg+my_nspinor*nband_cg_k*npw_k
1059 53192 : icg=icg+my_nspinor*nband_k*npw_k
1060 53192 : ikg=ikg+npw_k
1061 : end if
1062 :
1063 : ! End big k point loop
1064 53192 : ABI_FREE(ffnl)
1065 53192 : ABI_FREE(ph3d)
1066 53192 : ABI_FREE(phkxred)
1067 53192 : ABI_FREE(kg_k)
1068 53192 : ABI_FREE(kpg_k)
1069 53192 : ABI_FREE(ylm_k)
1070 53192 : ABI_FREE(ylmgr_k)
1071 87243 : ABI_FREE(cwavef)
1072 : end do
1073 : ! End loop over spins
1074 : end do
1075 :
1076 2229 : if ((iatom<=0).and.(choice==23)) then
1077 261 : do iatom1=1,ncprj
1078 41285 : do ii=1,mcprj
1079 41024 : nlmn=cprj(iatom1,ii)%nlmn
1080 369390 : do kk=1,nlmn
1081 2297344 : work(1:6)=cprj(iatom1,ii)%dcp(1,1:6,kk)
1082 328192 : call strconv(work,gprimd,work)
1083 2297344 : cprj(iatom1,ii)%dcp(1,1:6,kk)=work(1:6)
1084 2297344 : work(1:6)=cprj(iatom1,ii)%dcp(2,1:6,kk)
1085 328192 : call strconv(work,gprimd,work)
1086 2338368 : cprj(iatom1,ii)%dcp(2,1:6,kk)=work(1:6)
1087 : end do
1088 : end do
1089 : end do
1090 : end if
1091 :
1092 : !If needed, gather computed scalars
1093 2229 : if (.not.(cg_band_distributed .and. cprj_band_distributed)) then
1094 982 : call pawcprj_mpi_sum(cprj,spaceComm_band,ierr)
1095 : end if
1096 :
1097 : !Deallocate temporary storage
1098 2229 : if (one_atom) then
1099 0 : ABI_FREE(atindx_atm)
1100 0 : ABI_FREE(nattyp_atm)
1101 0 : ABI_FREE(ph1d_atm)
1102 0 : ABI_FREE(ekb_atm)
1103 0 : ABI_FREE(indlmn_atm)
1104 0 : ABI_FREE(ffspl_atm)
1105 0 : ABI_FREE(pspso_atm)
1106 : end if
1107 2229 : nullify(atindx_atm,nattyp_atm,ph1d_atm,ekb_atm,indlmn_atm,ffspl_atm,pspso_atm)
1108 2229 : call pawcprj_free(cwaveprj)
1109 12200 : ABI_FREE(cwaveprj)
1110 2229 : ABI_FREE(dimlmn)
1111 2229 : if (npband_bandfft>1) then
1112 832 : ABI_FREE(npw_block)
1113 832 : ABI_FREE(npw_disp)
1114 832 : ABI_FREE(bufsize)
1115 832 : ABI_FREE(bufdisp)
1116 832 : ABI_FREE(bufsize_wf)
1117 832 : ABI_FREE(bufdisp_wf)
1118 : end if
1119 :
1120 : DBG_EXIT('COLL')
1121 :
1122 4458 : end subroutine ctocprj
1123 : !!***
1124 :
1125 : !!****f* ABINIT/cprj_rotate
1126 : !! NAME
1127 : !! cprj_rotate
1128 : !!
1129 : !! FUNCTION
1130 : !! Compute cprj_nk = \sum_m z_m cprj_mk
1131 : !! where z_m is an array of complex values.
1132 : !! The input is overwritten.
1133 : !!
1134 : !! INPUTS
1135 : !!
1136 : !! SIDE EFFECTS
1137 : !!
1138 : !! SOURCE
1139 :
1140 0 : subroutine cprj_rotate(cprj_in,evec,dimcprj,natom,nband,nspinor)
1141 :
1142 : !Arguments -------------------------------
1143 : !scalars
1144 : integer,intent(in) :: natom,nband,nspinor
1145 : !arrays
1146 : integer,intent(in) :: dimcprj(:)
1147 : real(dp) :: evec(:,:)
1148 : type(pawcprj_type),intent(inout) :: cprj_in(natom,nspinor*nband)
1149 :
1150 : !Local variables-------------------------------
1151 : !scalars
1152 : integer :: iband,ncpgr
1153 : !arrays
1154 : ! real(dp) :: tsec(2)
1155 0 : real(dp) :: z_tmp(2,nband)
1156 0 : type(pawcprj_type),pointer :: cprj_iband(:,:)
1157 : type(pawcprj_type),allocatable,target :: cprj_tmp(:,:)
1158 :
1159 :
1160 : ! *********************************************************************
1161 :
1162 : DBG_ENTER('COLL')
1163 :
1164 : ! call timab(1211,1,tsec)
1165 :
1166 0 : ncpgr=cprj_in(1,1)%ncpgr
1167 0 : ABI_MALLOC(cprj_tmp,(natom,nspinor*nband))
1168 0 : call pawcprj_alloc(cprj_tmp,ncpgr,dimcprj)
1169 :
1170 0 : do iband=1,nband
1171 0 : z_tmp = reshape(evec(:,iband),(/2,nband/))
1172 0 : cprj_iband => cprj_tmp(:,nspinor*(iband-1)+1:nspinor*iband)
1173 0 : call pawcprj_lincom(z_tmp,cprj_in,cprj_iband,nband)
1174 : end do
1175 :
1176 0 : call pawcprj_copy(cprj_tmp,cprj_in)
1177 0 : call pawcprj_free(cprj_tmp)
1178 0 : ABI_FREE(cprj_tmp)
1179 :
1180 : ! call timab(1211,2,tsec)
1181 :
1182 : DBG_EXIT('COLL')
1183 :
1184 0 : end subroutine cprj_rotate
1185 : !!***
1186 :
1187 23072 : subroutine xg_cprj_copy(cprj,xg_cprj,xg_nonlop,option)
1188 :
1189 : integer, intent(in) :: option
1190 : type(pawcprj_type),intent(inout) :: cprj(:,:)
1191 : type(xgBlock_t), intent(inout) :: xg_cprj
1192 : type(xg_nonlop_t), intent(in) :: xg_nonlop
1193 :
1194 23072 : real(dp),pointer :: cprj_contiguous(:,:)
1195 : integer :: cplex,iatom,iband,iband_spin,ilmn,ispinor
1196 : integer :: natom,nband_cprj,nlmn,nspinor
1197 : integer :: start,end,space_cprj
1198 23072 : integer,allocatable :: nlmn_shift(:)
1199 :
1200 23072 : if (option/=XG_TO_CPRJ.and.option/=CPRJ_TO_XG) then
1201 0 : ABI_ERROR('Bad option')
1202 : end if
1203 :
1204 23072 : natom = xg_nonlop%natom
1205 23072 : nspinor = xg_nonlop%nspinor
1206 :
1207 23072 : nband_cprj = cols(xg_cprj)/nspinor
1208 :
1209 23072 : if (size(cprj,1)/=natom) then
1210 0 : ABI_ERROR('Bad size for cprj (for dim=1)')
1211 : end if
1212 23072 : if (size(cprj,2)/=nband_cprj*nspinor) then
1213 0 : ABI_ERROR('Bad size for cprj (for dim=2)')
1214 : end if
1215 :
1216 23072 : space_cprj = space(xg_cprj)
1217 23072 : cplex=2;if (space_cprj==SPACE_R) cplex=1
1218 :
1219 23072 : call xgBlock_reverseMap(xg_cprj,cprj_contiguous)
1220 :
1221 23072 : if (size(cprj_contiguous,1)/=cplex*xg_nonlop%cprjdim) then
1222 0 : ABI_ERROR('Bad size for cprj_contiguous (for dim=1)')
1223 : end if
1224 :
1225 69216 : ABI_MALLOC(nlmn_shift,(natom))
1226 23072 : nlmn_shift(1) = 0
1227 : ! Check nlm for every atom and fill nlmn_shift
1228 101224 : do iband=1,nband_cprj
1229 201656 : do ispinor=1,nspinor
1230 100432 : iband_spin = (iband-1)*nspinor+ispinor
1231 385208 : do iatom=1,natom
1232 206624 : nlmn=xg_nonlop%nlmn_natom(iatom)
1233 619872 : if (size(cprj(iatom,iband_spin)%cp)/=2*nlmn) then ! NOTE: cprj%cp size is always (2,nlmn) even in the real case
1234 0 : ABI_ERROR('Bad size for cprj_cwavef_bands%cp')
1235 : end if
1236 307056 : if (iatom<natom) then
1237 106192 : nlmn_shift(iatom+1) = nlmn_shift(iatom) + nlmn
1238 : end if
1239 : end do
1240 : end do
1241 : end do
1242 :
1243 : ! Do the copy
1244 23072 : if (option==CPRJ_TO_XG) then
1245 : !$omp parallel default (none) &
1246 : !$omp& shared(xg_nonlop,nlmn_shift,cprj_contiguous,cprj), &
1247 : !$omp& firstprivate(cplex,nband_cprj,nspinor,natom), &
1248 : !$omp& private(iband,ispinor,iatom,iband_spin,nlmn,start,end)
1249 : !$omp do collapse(3)
1250 0 : do iband=1,nband_cprj
1251 0 : do ispinor=1,nspinor
1252 0 : do iatom=1,natom
1253 0 : iband_spin = (iband-1)*nspinor+ispinor
1254 0 : nlmn=xg_nonlop%nlmn_natom(iatom)
1255 0 : do ilmn=1,nlmn
1256 0 : start = 1+cplex*(ilmn-1)+cplex*nlmn_shift(iatom)
1257 0 : end = cplex*ilmn +cplex*nlmn_shift(iatom)
1258 : cprj_contiguous(start:end,iband_spin) = &
1259 0 : cprj(iatom,iband_spin)%cp(1:cplex,ilmn)
1260 : end do
1261 : end do
1262 : end do
1263 : end do
1264 : !$omp end do
1265 : !$omp end parallel
1266 23072 : else if (option==XG_TO_CPRJ) then
1267 : !$omp parallel default (none) &
1268 : !$omp& shared(xg_nonlop,nlmn_shift,cprj_contiguous,cprj), &
1269 : !$omp& firstprivate(cplex,nband_cprj,nspinor,natom), &
1270 : !$omp& private(iband,ispinor,iatom,iband_spin,nlmn,start,end)
1271 : !$omp do collapse(3)
1272 101224 : do iband=1,nband_cprj
1273 201656 : do ispinor=1,nspinor
1274 385208 : do iatom=1,natom
1275 206624 : iband_spin = (iband-1)*nspinor+ispinor
1276 206624 : nlmn=xg_nonlop%nlmn_natom(iatom)
1277 2052848 : do ilmn=1,nlmn
1278 1745792 : start = 1+cplex*(ilmn-1)+cplex*nlmn_shift(iatom)
1279 1745792 : end = cplex*ilmn +cplex*nlmn_shift(iatom)
1280 : cprj(iatom,iband_spin)%cp(1:cplex,ilmn) = &
1281 5039120 : cprj_contiguous(start:end,iband_spin)
1282 : end do
1283 : end do
1284 : end do
1285 : end do
1286 : !$omp end do
1287 : !$omp end parallel
1288 : end if
1289 :
1290 23072 : ABI_FREE(nlmn_shift)
1291 :
1292 23072 : end subroutine xg_cprj_copy
1293 :
1294 : end module m_cgprj
1295 : !!***
|
