TABLE OF CONTENTS
ABINIT/ctocprj [ Functions ]
NAME
ctocprj
FUNCTION
Compute all <Proj_i|Cnk> for every wave function |Cnk> expressed in reciprocal space. |Proj_i> are non-local projectors (for each atom and each l,m,n) Can also compute derivatives of <Proj_i|Cnk> wrt to several parameters
COPYRIGHT
Copyright (C) 1998-2018 ABINIT group (MT) This file is distributed under the terms of the GNU General Public License, see ~ABINIT/COPYING or http://www.gnu.org/copyleft/gpl.txt . For the initials of contributors, see ~ABINIT/Infos/contributors .
INPUTS
atindx(natom)=index table for atoms
cg(2,mcg)=planewave coefficients of wavefunctions
choice: chooses derivatives to compute:
=1 => no derivatives
=2 => 1st derivatives with respect to atomic position(s)
=3 => 1st derivatives with respect to strain(s)
=23=> 1st derivatives with respect to strain(s) and atm pos
=4 => 2nd derivatives with respect to atomic pos.
=24=> 1st and 2nd derivatives with respect to atomic pos.
=5 => derivatives with respect to k wavevector
=6 => 2nd derivatives with respect to strain and atm. pos.
gmet(3,3)=reciprocal space metric tensor in bohr**-2
gprimd(3,3)=dimensional reciprocal space primitive translations
iatom= if <=0, cprj=<p_i|Cnk> are computed for all atoms 1...natom
if >0 cprj=<p_i|Cnk> are computed only for atom with index iatom
idir=direction of the derivative, i.e. dir. of - atom to be moved in the case choice=2
- strain component in the case choice=3
- k point direction in the case choice=5
Compatible only with choice=2,3,5; if idir=0, all derivatives are computed
iorder_cprj=0 if output cprj=<p_i|Cnk> are sorted by atom type
(first all elements of atom type 1, followed by those of atom type 2 and so on).
1 if output cprj=<p_i|Cnk> are sorted according to
the variable typat in the main input file
istwfk(nkpt)=option parameter that describes the storage of wfs
kg(3,mpw*mkmem)=reduced planewave coordinates
kpt(3,nkpt)=reduced coordinates of k points.
mcg=size of wave-functions array (cg) =mpw*nspinor*mband*mkmem*nsppol
mcprj=size of projected wave-functions array (cprj) =nspinor*mband*mkmem*nsppol
mgfft=maximum size of 1D FFTs
mkmem=number of k points treated by this node.
mpi_enreg=information about MPI parallelization
mpsang=1+maximum angular momentum for nonlocal pseudopotentials
mpw=maximum dimensioned size of npw
natom=number of atoms in cell
nattyp(ntypat)= # atoms of each type
nband(nkpt*nsppol)=number of bands at this k point for that spin polarization
ncprj=1st dim. of cprj array (natom if iatom<=0, 1 if iatom>0)
ngfft(18)=contain all needed information about 3D FFT, see ~ABINIT/Infos/vargs.htm#ngfft
nkpt=number of k points
nloalg(3)=governs the choice of the algorithm for nonlocal operator
npwarr(nkpt)=number of planewaves in basis at this k point
nspinor=number of spinorial components of the wavefunctions
nsppol=1 for unpolarized, 2 for spin-polarized
ntypat=number of types of atoms in unit cell
paral_kgb= 1 if kpt-band-FFT is activated
ph1d(2,3*(2*mgfft+1)*natom)=1-dim structure factor phase information
psps <type(pseudopotential_type)>=variables related to pseudopotentials
rmet(3,3)=real space metric (bohr**2)
tim_ctocprj=timing code of the calling routine
typat(natom)= types of atoms
uncp=unit number for <P_lmn|Cnk> data (if used)
xred(3,natom)=reduced dimensionless atomic coordinates
ylm(mpw*mkmem,mpsang*mpsang)=real spherical harmonics for each G and k point
ylmgr(npw*mkmem,nylmgr,mpsang*mpsang*useylmgr)=gradients of real spherical harmonics wrt (k+G)
OUTPUT
cprj(ncprj,mcprj) <type(pawcprj_type)>= projected input wave functions <Proj_i|Cnk> with NL projectors Usually ncprj=natom
PARENTS
dfpt_looppert,extrapwf,forstr,scfcv,update_e_field_vars,vtorho
CHILDREN
getcprj,mkffnl,mkkpg,pawcprj_alloc,pawcprj_free,pawcprj_mpi_sum pawcprj_put,pawcprj_set_zero,ph1d3d,strconv,xmpi_allgather xmpi_allgatherv,xmpi_alltoallv
SOURCE
88 #if defined HAVE_CONFIG_H 89 #include "config.h" 90 #endif 91 92 #include "abi_common.h" 93 94 subroutine ctocprj(atindx,cg,choice,cprj,gmet,gprimd,iatom,idir,& 95 & iorder_cprj,istwfk,kg,kpt,mcg,mcprj,mgfft,mkmem,mpi_enreg,mpsang,& 96 & mpw,natom,nattyp,nband,ncprj,ngfft,nkpt,nloalg,npwarr,nspinor,& 97 & nsppol,ntypat,paral_kgb,ph1d,psps,rmet,typat,ucvol,uncp,xred,ylm,ylmgr) 98 99 use defs_basis 100 use defs_datatypes 101 use defs_abitypes 102 use m_profiling_abi 103 use m_xmpi 104 use m_errors 105 use m_hdr 106 107 use m_kg, only : ph1d3d 108 use m_pawcprj, only : pawcprj_type, pawcprj_alloc, pawcprj_put, pawcprj_free, & 109 & pawcprj_set_zero, pawcprj_mpi_sum 110 111 !This section has been created automatically by the script Abilint (TD). 112 !Do not modify the following lines by hand. 113 #undef ABI_FUNC 114 #define ABI_FUNC 'ctocprj' 115 use interfaces_32_util 116 use interfaces_41_geometry 117 use interfaces_66_nonlocal, except_this_one => ctocprj 118 !End of the abilint section 119 120 implicit none 121 122 !Arguments ------------------------------- 123 !scalars 124 integer,intent(in) :: choice,iatom,idir,iorder_cprj,mcg,mcprj,mgfft,mkmem,mpsang,mpw 125 integer,intent(in) :: natom,ncprj,nkpt,nspinor,nsppol,ntypat,paral_kgb,uncp 126 real(dp),intent(in) :: ucvol 127 type(MPI_type),intent(in) :: mpi_enreg 128 type(pseudopotential_type),target,intent(in) :: psps 129 !arrays 130 integer,intent(in) :: istwfk(nkpt),nband(nkpt*nsppol) 131 integer,intent(in) :: ngfft(18),nloalg(3),npwarr(nkpt),kg(3,mpw*mkmem),typat(natom) 132 integer,intent(in),target :: atindx(natom),nattyp(ntypat) 133 real(dp),intent(in) :: cg(2,mcg) 134 real(dp),intent(in) :: gmet(3,3),gprimd(3,3),kpt(3,nkpt),rmet(3,3) 135 real(dp),intent(in) :: xred(3,natom),ylm(:,:),ylmgr(:,:,:) 136 real(dp),intent(in),target :: ph1d(2,3*(2*mgfft+1)*natom) 137 type(pawcprj_type),intent(inout) :: cprj(ncprj,mcprj) 138 139 !Local variables------------------------------- 140 !scalars 141 integer :: blocksz,cg_bandpp,counter,cpopt,cprj_bandpp,dimffnl,ia,iatm,iatom1,iatom2 142 integer :: iband_max,iband_min,iband_start,ibg,ibgb,iblockbd,ibp,icg,icgb,icp1,icp2 143 integer :: ider,idir0,iend,ierr,ig,ii,ikg,ikpt,ilm,ipw,isize,isppol,istart,istwf_k,itypat,iwf1,iwf2,jdir 144 integer :: matblk,mband_cprj,me_distrb,my_nspinor,n1,n1_2p1,n2,n2_2p1,n3,n3_2p1,kk,nlmn 145 integer :: nband_k,nband_cprj_k,nblockbd,ncpgr,nkpg,npband_bandfft,npws,npw_k,npw_nk,ntypat0 146 integer :: shift1,shift1b,shift2,shift2b,shift3,shift3b 147 integer :: spaceComm,spaceComm_band,spaceComm_fft,useylmgr 148 logical :: cg_band_distributed,cprj_band_distributed,one_atom 149 real(dp) :: arg 150 character(len=500) :: msg 151 !arrays 152 integer,allocatable :: bufsize(:),bufsize_wf(:),bufdisp(:),bufdisp_wf(:) 153 integer,allocatable :: dimlmn(:),kg_k(:,:),kg_k_loc(:,:) 154 integer,allocatable :: npw_block(:),npw_disp(:) 155 integer,pointer :: atindx_atm(:),indlmn_atm(:,:,:),nattyp_atm(:),pspso_atm(:) 156 real(dp) :: kpoint(3),work(6) 157 real(dp),allocatable :: cwavef(:,:),cwavef_tmp(:,:) 158 real(dp),allocatable :: ffnl(:,:,:,:),ffnl_npw(:,:,:,:),ffnl_tmp(:,:,:,:),ffnl_tmp_npw(:,:,:,:) 159 real(dp),allocatable :: kpg_k(:,:) 160 real(dp),allocatable :: ph3d(:,:,:),ph3d_npw(:,:,:),ph3d_tmp(:,:,:),ph3d_tmp_npw(:,:,:) 161 real(dp),allocatable :: phkxred(:,:),ylm_k(:,:),ylmgr_k(:,:,:) 162 real(dp),ABI_CONTIGUOUS pointer :: ekb_atm(:,:),ffspl_atm(:,:,:,:),ph1d_atm(:,:) 163 type(pawcprj_type),allocatable :: cwaveprj(:,:) 164 165 ! ********************************************************************* 166 167 DBG_ENTER('COLL') 168 169 !Nothing to do if current MPI process does treat kpoints or plane-waves 170 if (mcg==0.or.mcprj==0) return 171 172 !Preliminary tests 173 if (psps%useylm==0) then 174 msg='Not available for useylm=0!' 175 MSG_ERROR(msg) 176 end if 177 if ((choice<1.or.choice>6).and.choice/=23.and.choice/=24) then 178 msg='Bad choice!' 179 MSG_BUG(msg) 180 end if 181 if (idir>0.and.choice/=2.and.choice/=3.and.choice/=5) then 182 msg='Does not support idir>0 for that choice!' 183 MSG_BUG(msg) 184 end if 185 if (size(ylm)/=mpw*mkmem*mpsang*mpsang) then 186 msg='Wrong size for Ylm!' 187 MSG_BUG(msg) 188 end if 189 useylmgr=merge(0,1,(size(ylmgr)==0)) 190 if (useylmgr==0.and.(choice==3.or.choice==5.or.choice==23)) then 191 msg=' Ylm gradients have to be in memory for choice=3, 5, or 23!' 192 MSG_BUG(msg) 193 end if 194 195 !Init parallelism 196 if (paral_kgb==1) then 197 me_distrb=mpi_enreg%me_kpt 198 spaceComm=mpi_enreg%comm_kpt 199 spaceComm_fft=mpi_enreg%comm_fft 200 npband_bandfft=mpi_enreg%nproc_band 201 cg_bandpp=mpi_enreg%bandpp 202 cprj_bandpp=mpi_enreg%bandpp 203 spaceComm_band=mpi_enreg%comm_band 204 cg_band_distributed=.true. 205 cprj_band_distributed=(mcprj/=mcg/mpw) 206 else 207 me_distrb=mpi_enreg%me_kpt 208 spaceComm=mpi_enreg%comm_cell 209 spaceComm_fft=xmpi_comm_self 210 npband_bandfft=1 211 cg_bandpp=1;cprj_bandpp=1 212 cg_band_distributed=.false. 213 cprj_band_distributed=.false. 214 if (mpi_enreg%paralbd==1) then 215 spaceComm_band=mpi_enreg%comm_band 216 else 217 spaceComm_band=xmpi_comm_self 218 end if 219 end if 220 if (cg_bandpp/=cprj_bandpp) then 221 MSG_BUG('cg_bandpp must be equal to cprj_bandpp !') 222 end if 223 224 !Manage parallelization over spinors 225 my_nspinor=max(1,nspinor/mpi_enreg%nproc_spinor) 226 !Check sizes for cprj (distribution is tricky) 227 one_atom=(iatom>0) 228 if (one_atom.and.ncprj/=1) then 229 MSG_BUG('Bad value for ncprj dimension (should be 1) !') 230 end if 231 if (.not.one_atom.and.ncprj/=natom) then 232 MSG_BUG('Bad value for ncprj dimension (should be natom) !') 233 end if 234 235 !Initialize some variables 236 mband_cprj=mcprj/(my_nspinor*mkmem*nsppol) 237 n1=ngfft(1);n2=ngfft(2);n3=ngfft(3) 238 n1_2p1=2*n1+1;n2_2p1=2*n2+1;n3_2p1=2*n3+1 239 ibg=0;icg=0;cpopt=0 240 ider=0;idir0=0;istart=idir;iend=idir 241 if (choice==3.or.choice==5.or.choice==23) ider=1 242 if (idir>0) then 243 if (choice==3) idir0=-idir 244 if (choice==5) idir0=idir 245 else 246 ! if (choice==23) idir0=-7 247 if (choice==3) idir0=-7 248 if (choice==5) idir0=4 249 end if 250 if (idir0==0.or.idir0==4) then 251 dimffnl=1+3*ider 252 else if (idir0/=-7) then 253 dimffnl=1+ider 254 else 255 dimffnl=1+6*ider 256 if(choice==3)then 257 istart=ider;iend=6*ider 258 end if 259 end if 260 nkpg=0 261 if (choice==3.or.choice==2.or.choice==23) nkpg=3*nloalg(3) 262 if (choice==4.or.choice==24) nkpg=9*nloalg(3) 263 264 !Set number of gradients for <p_i|Cnk> 265 ncpgr=0 266 if (idir==0) then 267 if (choice==2) ncpgr=3 268 if (choice==3) ncpgr=6 269 if (choice==23)ncpgr=9 270 if (choice==4) ncpgr=6 271 if (choice==24)ncpgr=9 272 if (choice==5) ncpgr=3 273 if (choice==6) ncpgr=63 274 else 275 ncpgr=1 276 end if 277 !Test cprj gradients dimension (just to be sure) 278 if (cprj(1,1)%ncpgr/=ncpgr) then 279 MSG_BUG('cprj are badly allocated !') 280 end if 281 282 283 !Extract data for treated atom(s) 284 if (one_atom) then 285 iatom1=iatom;iatom2=iatom 286 ntypat0=1;itypat=typat(iatom) 287 ABI_ALLOCATE(nattyp_atm,(ntypat0)) 288 nattyp_atm(1)=1 289 ABI_ALLOCATE(atindx_atm,(ntypat0)) 290 atindx_atm(1)=atindx(iatom) 291 ABI_ALLOCATE(ph1d_atm,(2,(n1_2p1+n2_2p1+n3_2p1)*ntypat0)) 292 shift1=(atindx(iatom)-1)*n1_2p1 293 shift2=(atindx(iatom)-1)*n2_2p1+natom*n1_2p1 294 shift3=(atindx(iatom)-1)*n3_2p1+natom*(n1_2p1+n2_2p1) 295 shift1b=0;shift2b=n1_2p1;shift3b=n1_2p1+n2_2p1 296 ph1d_atm(:,shift1b+1:shift1b+n1_2p1)=ph1d(:,shift1+1:shift1+n1_2p1) 297 ph1d_atm(:,shift2b+1:shift2b+n2_2p1)=ph1d(:,shift2+1:shift2+n2_2p1) 298 ph1d_atm(:,shift3b+1:shift3b+n3_2p1)=ph1d(:,shift3+1:shift3+n3_2p1) 299 ABI_ALLOCATE(ekb_atm,(psps%dimekb,ntypat0)) 300 ABI_ALLOCATE(indlmn_atm,(6,psps%lmnmax,ntypat0)) 301 ABI_ALLOCATE(ffspl_atm,(psps%mqgrid_ff,2,psps%lnmax,ntypat0)) 302 ABI_ALLOCATE(pspso_atm,(ntypat0)) 303 ekb_atm(:,1)=psps%ekb(:,itypat) 304 indlmn_atm(:,:,1)=psps%indlmn(:,:,itypat) 305 ffspl_atm(:,:,:,1)=psps%ffspl(:,:,:,itypat) 306 pspso_atm(1)=psps%pspso(itypat) 307 else 308 iatom1=1;iatom2=natom 309 ntypat0=ntypat 310 atindx_atm => atindx 311 nattyp_atm => nattyp 312 ph1d_atm => ph1d 313 ekb_atm => psps%ekb 314 indlmn_atm => psps%indlmn 315 ffspl_atm => psps%ffspl 316 pspso_atm => psps%pspso 317 end if 318 319 !Dimensioning and allocation of <p_i|Cnk> 320 ABI_ALLOCATE(dimlmn,(ncprj)) 321 dimlmn=0 ! Type-sorted cprj 322 if (one_atom) then 323 itypat=typat(iatom) 324 dimlmn(ia+1:ia+nattyp(itypat))=count(indlmn_atm(3,:,itypat)>0) 325 else 326 ia=0 327 do itypat=1,ntypat0 328 dimlmn(ia+1:ia+nattyp(itypat))=count(indlmn_atm(3,:,itypat)>0) 329 ia=ia+nattyp(itypat) 330 end do 331 end if 332 ABI_DATATYPE_ALLOCATE(cwaveprj,(ncprj,my_nspinor*cprj_bandpp)) 333 call pawcprj_alloc(cwaveprj,ncpgr,dimlmn) 334 335 !Additional statements if band-fft parallelism 336 if (npband_bandfft>1) then 337 ABI_ALLOCATE(npw_block,(npband_bandfft)) 338 ABI_ALLOCATE(npw_disp,(npband_bandfft)) 339 ABI_ALLOCATE(bufsize,(npband_bandfft*cg_bandpp)) 340 ABI_ALLOCATE(bufdisp,(npband_bandfft*cg_bandpp)) 341 ABI_ALLOCATE(bufsize_wf,(npband_bandfft*cg_bandpp)) 342 ABI_ALLOCATE(bufdisp_wf,(npband_bandfft*cg_bandpp)) 343 end if 344 345 !Set output datastructure to zero 346 call pawcprj_set_zero(cprj) 347 348 !LOOP OVER SPINS 349 do isppol=1,nsppol 350 ikg=0 351 352 ! BIG FAT k POINT LOOP 353 do ikpt=1,nkpt 354 counter=100*ikpt+isppol 355 356 ! Select k point to be treated by this proc 357 nband_k=nband(ikpt+(isppol-1)*nkpt) 358 if (proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,isppol,me_distrb)) cycle 359 360 ! Retrieve k-point 361 kpoint(:)=kpt(:,ikpt) 362 istwf_k=istwfk(ikpt) 363 364 ! Retrieve number of plane waves 365 npw_k=npwarr(ikpt) 366 if (npband_bandfft>1) then 367 ! Special treatment for band-fft // 368 call xmpi_allgather(npw_k,npw_block,spaceComm_band,ierr) 369 npw_nk=sum(npw_block);npw_disp(1)=0 370 do ii=2,npband_bandfft 371 npw_disp(ii)=npw_disp(ii-1)+npw_block(ii-1) 372 end do 373 else 374 npw_nk=npw_k 375 end if 376 377 ! Retrieve (k+G) points and spherical harmonics 378 ABI_ALLOCATE(ylm_k,(npw_k,mpsang*mpsang)) 379 ABI_ALLOCATE(ylmgr_k,(npw_k,3,mpsang*mpsang*useylmgr)) 380 ABI_ALLOCATE(kg_k,(3,npw_nk)) 381 if (npband_bandfft>1) then 382 ! Special treatment for band-fft // 383 ABI_ALLOCATE(kg_k_loc,(3,npw_k)) 384 kg_k_loc(:,1:npw_k)=kg(:,1+ikg:npw_k+ikg) 385 bufsize(:)=3*npw_block(:);bufdisp(:)=3*npw_disp(:) 386 call xmpi_allgatherv(kg_k_loc,3*npw_k,kg_k,bufsize,bufdisp,spaceComm_band,ierr) 387 else 388 kg_k(:,1:npw_k)=kg(:,1+ikg:npw_k+ikg) 389 end if 390 do ilm=1,mpsang*mpsang 391 ylm_k(1:npw_k,ilm)=ylm(1+ikg:npw_k+ikg,ilm) 392 if (useylmgr>0) ylmgr_k(1:npw_k,1:3,ilm)=ylmgr(1+ikg:npw_k+ikg,1:3,ilm) 393 end do 394 395 ! Compute (k+G) vectors 396 ABI_ALLOCATE(kpg_k,(npw_nk,nkpg)) 397 if (nkpg>0) then 398 call mkkpg(kg_k,kpg_k,kpoint,nkpg,npw_nk) 399 end if 400 ! Allocate and compute the arrays phkxred and ph3d 401 ABI_ALLOCATE(phkxred,(2,ncprj)) 402 do ia=iatom1,iatom2 403 iatm=min(atindx_atm(ia),ncprj) 404 arg=two_pi*(kpoint(1)*xred(1,ia)+kpoint(2)*xred(2,ia)+kpoint(3)*xred(3,ia)) 405 phkxred(1,iatm)=cos(arg);phkxred(2,iatm)=sin(arg) 406 end do 407 matblk=ncprj;if (nloalg(2)<=0) matblk=0 408 ABI_ALLOCATE(ph3d,(2,npw_nk,matblk)) 409 if (matblk>0)then 410 ! Here, precomputation of ph3d 411 if (npband_bandfft>1) then 412 ! Special treatment for band-fft // 413 ABI_ALLOCATE(ph3d_tmp,(2,npw_k,matblk)) 414 call ph1d3d(1,ncprj,kg_k_loc,matblk,ncprj,npw_k,n1,n2,n3,phkxred,ph1d_atm,ph3d_tmp) 415 ABI_ALLOCATE(ph3d_tmp_npw,(2,matblk,npw_k)) 416 ABI_ALLOCATE(ph3d_npw,(2,matblk,npw_nk)) 417 isize=2*matblk;bufsize(:)=isize*npw_block(:);bufdisp(:)=isize*npw_disp(:) 418 do ipw=1,npw_k 419 ph3d_tmp_npw(:,:,ipw)=ph3d_tmp(:,ipw,:) 420 end do 421 call xmpi_allgatherv(ph3d_tmp_npw,isize*npw_k,ph3d_npw,bufsize,bufdisp,spaceComm_band,ierr) 422 do ipw=1,npw_nk 423 ph3d(:,ipw,:)=ph3d_npw(:,:,ipw) 424 end do 425 ABI_DEALLOCATE(ph3d_npw) 426 ABI_DEALLOCATE(ph3d_tmp_npw) 427 ABI_DEALLOCATE(ph3d_tmp) 428 else 429 call ph1d3d(1,ncprj,kg_k,matblk,ncprj,npw_k,n1,n2,n3,phkxred,ph1d_atm,ph3d) 430 end if 431 else if (npband_bandfft>1) then 432 MSG_ERROR('Band-fft parallelism +nloag(1)<0 forbidden !') 433 end if 434 435 ! Compute nonlocal form factors ffnl at all (k+G) 436 ABI_ALLOCATE(ffnl,(npw_nk,dimffnl,psps%lmnmax,ntypat0)) 437 if (npband_bandfft>1) then 438 ! Special treatment for band-fft // 439 ABI_ALLOCATE(ffnl_tmp,(npw_k,dimffnl,psps%lmnmax,ntypat0)) 440 call mkffnl(psps%dimekb,dimffnl,ekb_atm,ffnl_tmp,ffspl_atm,& 441 & gmet,gprimd,ider,idir0,indlmn_atm,kg_k_loc,kpg_k,kpoint,psps%lmnmax,& 442 & psps%lnmax,psps%mpsang,psps%mqgrid_ff,nkpg,npw_k,ntypat0,& 443 & pspso_atm,psps%qgrid_ff,rmet,psps%usepaw,psps%useylm,ylm_k,ylmgr_k) 444 ABI_ALLOCATE(ffnl_tmp_npw,(dimffnl,psps%lmnmax,ntypat0,npw_k)) 445 ABI_ALLOCATE(ffnl_npw,(dimffnl,psps%lmnmax,ntypat0,npw_nk)) 446 isize=dimffnl*psps%lmnmax*ntypat0 447 bufsize(:)=isize*npw_block(:);bufdisp(:)=isize*npw_disp(:) 448 do ipw=1,npw_k 449 ffnl_tmp_npw(:,:,:,ipw)=ffnl_tmp(ipw,:,:,:) 450 end do 451 call xmpi_allgatherv(ffnl_tmp_npw,isize*npw_k,ffnl_npw,bufsize,bufdisp,spaceComm_band,ierr) 452 do ipw=1,npw_nk 453 ffnl(ipw,:,:,:)=ffnl_npw(:,:,:,ipw) 454 end do 455 ABI_DEALLOCATE(ffnl_npw) 456 ABI_DEALLOCATE(ffnl_tmp_npw) 457 ABI_DEALLOCATE(ffnl_tmp) 458 else 459 call mkffnl(psps%dimekb,dimffnl,ekb_atm,ffnl,ffspl_atm,& 460 & gmet,gprimd,ider,idir0,indlmn_atm,kg_k,kpg_k,kpoint,psps%lmnmax,& 461 & psps%lnmax,psps%mpsang,psps%mqgrid_ff,nkpg,npw_k,ntypat0,& 462 & pspso_atm,psps%qgrid_ff,rmet,psps%usepaw,psps%useylm,ylm_k,ylmgr_k) 463 end if 464 465 ! No more need of kg_g_tmp 466 if (npband_bandfft>1) then 467 ABI_DEALLOCATE(kg_k_loc) 468 end if 469 470 ! Allocate arrays for a wave-function (or a block of WFs) 471 ABI_ALLOCATE(cwavef,(2,npw_nk*my_nspinor*cg_bandpp)) 472 if (npband_bandfft>1) then 473 isize=2*my_nspinor*cg_bandpp;bufsize(:)=isize*npw_block(:);bufdisp(:)=isize*npw_disp(:) 474 isize=2*my_nspinor*npw_k*cg_bandpp;bufsize_wf(:)=isize 475 do ii=1,npband_bandfft*cg_bandpp 476 bufdisp_wf(ii)=(ii-1)*isize 477 end do 478 end if 479 480 ! Loop over bands or blocks of bands 481 icgb=icg ; ibgb=ibg ; iband_start=1 482 blocksz=npband_bandfft*cg_bandpp 483 nblockbd=nband_k/blocksz 484 nband_cprj_k=merge(nband_k/npband_bandfft,nband_k,cprj_band_distributed) 485 do iblockbd=1,nblockbd 486 iband_min=1+(iblockbd-1)*blocksz 487 iband_max=iblockbd*blocksz 488 489 if (proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,iband_min,iband_max,isppol,me_distrb)) then 490 if (.not.cg_band_distributed) icgb=icgb+npw_k*my_nspinor*blocksz 491 if (.not.cprj_band_distributed) ibgb=ibgb+my_nspinor*blocksz 492 cycle 493 end if 494 495 ! Extract wavefunction information 496 ! Special treatment for band-fft parallelism 497 if (npband_bandfft>1) then 498 !Transpose WF to get them in "FFT" representation 499 ABI_ALLOCATE(cwavef_tmp,(2,npw_k*my_nspinor*blocksz)) 500 do ig=1,npw_k*my_nspinor*blocksz 501 cwavef_tmp(1,ig)=cg(1,ig+icgb) 502 cwavef_tmp(2,ig)=cg(2,ig+icgb) 503 end do 504 call xmpi_alltoallv(cwavef_tmp,bufsize_wf,bufdisp_wf,cwavef,bufsize,bufdisp,spaceComm_band,ierr) 505 ABI_DEALLOCATE(cwavef_tmp) 506 !Reorder WF according to cg_bandpp and/or spinor 507 if (cg_bandpp>1.or.my_nspinor>1) then 508 ABI_ALLOCATE(cwavef_tmp,(2,npw_nk*my_nspinor*blocksz)) 509 do ig=1,npw_nk*my_nspinor*blocksz 510 cwavef_tmp(:,ig)=cwavef(:,ig) 511 end do 512 shift1=0 513 do iwf2=1,cg_bandpp 514 do ig=1,my_nspinor 515 shift2=0 516 do iwf1=1,npband_bandfft 517 npws=npw_block(iwf1) 518 ipw=shift2+(iwf2-1)*my_nspinor*npws+(ig-1)*npws 519 cwavef(:,shift1+1:shift1+npws)=cwavef_tmp(:,ipw+1:ipw+npws) 520 shift1=shift1+npws ; shift2=shift2+cg_bandpp*my_nspinor*npws 521 end do 522 end do 523 end do 524 ABI_DEALLOCATE(cwavef_tmp) 525 end if 526 else 527 do ig=1,npw_k*my_nspinor*cg_bandpp 528 cwavef(1,ig)=cg(1,ig+icgb) 529 cwavef(2,ig)=cg(2,ig+icgb) 530 end do 531 end if 532 533 ! Compute scalar product of wavefunction with all NL projectors 534 do ibp=1,cg_bandpp ! Note: we suppose cp_bandpp=cprj_bandpp 535 iwf1=1+(ibp-1)*npw_nk*my_nspinor;iwf2=ibp*npw_nk*my_nspinor 536 icp1=1+(ibp-1)*my_nspinor;icp2=ibp*my_nspinor 537 do jdir=istart,iend 538 call getcprj(choice,cpopt,cwavef(:,iwf1:iwf2),cwaveprj(:,icp1:icp2),& 539 & ffnl,jdir,indlmn_atm,istwf_k,kg_k,kpg_k,kpoint,psps%lmnmax,& 540 & mgfft,mpi_enreg,ncprj,nattyp_atm,ngfft,nloalg,& 541 & npw_nk,my_nspinor,ntypat0,phkxred,ph1d_atm,ph3d,ucvol,psps%useylm) 542 end do 543 end do 544 ! Export cwaveprj to big array cprj 545 call pawcprj_put(atindx_atm,cwaveprj,cprj,ncprj,iband_start,ibgb,ikpt,iorder_cprj,isppol,& 546 & mband_cprj,mkmem,natom,cprj_bandpp,nband_cprj_k,dimlmn,my_nspinor,nsppol,uncp,& 547 & mpi_comm_band=spaceComm_band,to_be_gathered=(cg_band_distributed.and.(.not.cprj_band_distributed))) 548 549 iband_start=iband_start+merge(cg_bandpp,blocksz,cprj_band_distributed) 550 551 ! End loop over bands 552 icgb=icgb+npw_k*my_nspinor*blocksz 553 end do 554 555 ! Shift array memory (if mkmem/=0) 556 if (mkmem/=0) then 557 ibg=ibg+my_nspinor*nband_cprj_k 558 icg=icg+my_nspinor*nband_k*npw_k 559 ikg=ikg+npw_k 560 end if 561 562 ! End big k point loop 563 ABI_DEALLOCATE(ffnl) 564 ABI_DEALLOCATE(ph3d) 565 ABI_DEALLOCATE(phkxred) 566 ABI_DEALLOCATE(kg_k) 567 ABI_DEALLOCATE(kpg_k) 568 ABI_DEALLOCATE(ylm_k) 569 ABI_DEALLOCATE(ylmgr_k) 570 ABI_DEALLOCATE(cwavef) 571 end do 572 ! End loop over spins 573 end do 574 575 if ((iatom<=0).and.(choice==23)) then 576 do iatom1=1,ncprj 577 do ii=1,mcprj 578 nlmn=cprj(iatom1,ii)%nlmn 579 do kk=1,nlmn 580 work(1:6)=cprj(iatom1,ii)%dcp(1,1:6,kk) 581 call strconv(work,gprimd,work) 582 cprj(iatom1,ii)%dcp(1,1:6,kk)=work(1:6) 583 work(1:6)=cprj(iatom1,ii)%dcp(2,1:6,kk) 584 call strconv(work,gprimd,work) 585 cprj(iatom1,ii)%dcp(2,1:6,kk)=work(1:6) 586 end do 587 end do 588 end do 589 end if 590 591 !If needed, gather computed scalars 592 if (.not.cg_band_distributed) then 593 call pawcprj_mpi_sum(cprj,spaceComm_band,ierr) 594 end if 595 596 !Deallocate temporary storage 597 if (one_atom) then 598 ABI_DEALLOCATE(atindx_atm) 599 ABI_DEALLOCATE(nattyp_atm) 600 ABI_DEALLOCATE(ph1d_atm) 601 ABI_DEALLOCATE(ekb_atm) 602 ABI_DEALLOCATE(indlmn_atm) 603 ABI_DEALLOCATE(ffspl_atm) 604 ABI_DEALLOCATE(pspso_atm) 605 end if 606 nullify(atindx_atm,nattyp_atm,ph1d_atm,ekb_atm,indlmn_atm,ffspl_atm,pspso_atm) 607 call pawcprj_free(cwaveprj) 608 ABI_DATATYPE_DEALLOCATE(cwaveprj) 609 ABI_DEALLOCATE(dimlmn) 610 if (npband_bandfft>1) then 611 ABI_DEALLOCATE(npw_block) 612 ABI_DEALLOCATE(npw_disp) 613 ABI_DEALLOCATE(bufsize) 614 ABI_DEALLOCATE(bufdisp) 615 ABI_DEALLOCATE(bufsize_wf) 616 ABI_DEALLOCATE(bufdisp_wf) 617 end if 618 619 DBG_EXIT('COLL') 620 621 end subroutine ctocprj