TABLE OF CONTENTS
ABINIT/getcprj [ Functions ]
NAME
getcprj
FUNCTION
Compute <Proj_i|Cnk> for one wave function |Cnk> expressed in reciprocal space. Compute also derivatives of <Proj_i|Cnk>. |Proj_i> are non-local projectors (for each atom and each l,m,n)
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
choice=chooses possible output:
In addition to projected wave function:
choice=1 => nothing else
=2 => 1st gradients with respect to atomic position(s)
=3 => 1st gradients with respect to strain(s)
=23=> 1st gradients 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 => 1st gradients with respect to k wavevector
=6 => 2nd derivatives with respect to strain and atm. pos.
cpopt=1 if <Proj_i|Cnk> are already in memory; see below (side effects).
cwavef(2,nspinor*npw_k)=input cmplx wavefunction coefficients <G|Cnk>
ffnl(npw_k,dimffnl,lmnmax,ntypat)=nonlocal form factors to be used for the application of the nl operator
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
indlmn(6,i,ntypat)= array giving l,m,n,lm,ln,s for i=lmn
istwf_k=option parameter that describes the storage of wfs
kg_k(3,npw_k)=reduced planewave coordinates
kpg(npw_k,npk)=(k+G) components and related data
kpoint(3)=k point in terms of recip. translations
lmnmax=max. number of (l,m,n) components over all types of atoms
mgfft=maximum size of 1D FFTs
mpi_enreg=informations about MPI parallelization
natom=number of atoms in cell
nattyp(ntypat)=number of atoms of each type
ngfft(18)=contain all needed information about 3D FFT, see ~ABINIT/Infos/vargs.htm#ngfft
nloalg(3)=governs the choice of the algorithm for nonlocal operator
npw_k=number of planewaves for given k point
nspinor=number of spinorial components of the wavefunctions (on current proc)
ntypat=number of types of atoms in unit cell
phkxred(2,natom)=phase factors exp(2 pi kpoint.xred)
ph1d(2,3*(2*mgfft+1)*natom)=1-dim structure factor phase information
ph3d(2,npw_k,natom)=3D structure factors, for each atom and plane wave
only used if nloalg(2)>0
ucvol= unit cell volume
useylm=governs the way the nonlocal operator is to be applied
SIDE EFFECTS
cwaveprj(natom,nspinor) <type(pawcprj_type)>=projected input wave function <Proj_i|Cnk> with all NL projectors (and derivatives) if cpopt=1 the projected scalars have been already been computed and only derivatives are computed here if cpopt=0 the projected scalars and derivatives are computed here
TODO
Spin-orbit
PARENTS
cgwf,ctocprj,debug_tools,dfpt_accrho,dfpt_nstpaw,ks_ddiago,m_wfd rf2_init
CHILDREN
mkkpg,opernla_ylm,ph1d3d
SOURCE
77 #if defined HAVE_CONFIG_H 78 #include "config.h" 79 #endif 80 81 #include "abi_common.h" 82 83 subroutine getcprj(choice,cpopt,cwavef,cwaveprj,ffnl,& 84 & idir,indlmn,istwf_k,kg_k,kpg,kpoint,lmnmax,mgfft,mpi_enreg,& 85 & natom,nattyp,ngfft,nloalg,npw_k,nspinor,ntypat,& 86 & phkxred,ph1d,ph3d,ucvol,useylm) 87 88 use defs_basis 89 use defs_abitypes 90 use m_profiling_abi 91 use m_errors 92 93 use m_kg, only : ph1d3d 94 use m_pawcprj, only : pawcprj_type 95 96 !This section has been created automatically by the script Abilint (TD). 97 !Do not modify the following lines by hand. 98 #undef ABI_FUNC 99 #define ABI_FUNC 'getcprj' 100 use interfaces_66_nonlocal, except_this_one => getcprj 101 !End of the abilint section 102 103 implicit none 104 105 !Arguments ------------------------------- 106 !scalars 107 integer,intent(in) :: choice,cpopt,idir,istwf_k,lmnmax 108 integer,intent(in) :: mgfft,natom,npw_k,nspinor,ntypat,useylm 109 real(dp),intent(in) :: ucvol 110 type(MPI_type),intent(in) :: mpi_enreg 111 !arrays 112 integer,intent(in) :: indlmn(6,lmnmax,ntypat),kg_k(3,npw_k),nattyp(ntypat) 113 integer,intent(in) :: ngfft(18),nloalg(3) 114 real(dp),intent(in) :: cwavef(2,npw_k*nspinor) 115 real(dp),intent(in),target :: ffnl(:,:,:,:),kpg(:,:),ph3d(:,:,:) 116 real(dp),intent(in) :: kpoint(3),ph1d(2,3*(2*mgfft+1)*natom),phkxred(2,natom) 117 type(pawcprj_type),intent(inout) :: cwaveprj(natom,nspinor) 118 119 !Local variables------------------------------- 120 !scalars 121 integer :: choice_,cplex,dimffnl,ia,ia1,ia2,ia3,ia4,iatm,ic,ii,ilmn,ishift,ispinor,itypat 122 integer :: jc,matblk,mincat,nd2gxdt,ndgxdt,nincat,nkpg,nkpg_,nlmn,signs 123 !arrays 124 integer,allocatable :: cplex_dgxdt(:),cplex_d2gxdt(:),indlmn_typ(:,:) 125 real(dp),allocatable :: d2gxdt(:,:,:,:,:),dgxdt(:,:,:,:,:),ffnl_typ(:,:,:) 126 real(dp),allocatable :: gx(:,:,:,:) 127 real(dp), pointer :: kpg_(:,:),ph3d_(:,:,:) 128 129 ! ********************************************************************* 130 131 DBG_ENTER('COLL') 132 133 !Nothing to do in that case 134 if (cpopt==1.and.choice==1) return 135 136 !Not available for useylm=0 137 if (useylm==0) then 138 MSG_ERROR('Not available for useylm=0 !') 139 end if 140 141 !Error on bad choice 142 if ((choice<1.or.choice>6).and.choice/=23.and.choice/=24) then 143 MSG_BUG('Does not presently support this choice !') 144 end if 145 146 !Error on bad idir 147 if (idir>0.and.choice/=2.and.choice/=3.and.choice/=5) then 148 MSG_BUG('Does not support idir>0 for this choice') 149 end if 150 151 !Error on sizes 152 nkpg=size(kpg,2) 153 if (nkpg>0) then 154 if( (choice==2.and.nkpg<3) .or. & 155 & ((choice==4.or.choice==24).and.nkpg<9) .or. & 156 & ((choice==6.or.choice==3.or.choice==23).and.nkpg<3) ) then 157 MSG_BUG('Incorrect size for kpg array !') 158 end if 159 end if 160 if (size(ffnl,1)/=npw_k.or.size(ffnl,3)/=lmnmax) then 161 MSG_BUG('Incorrect size for ffnl!') 162 end if 163 if (size(ph3d)>0) then 164 if (size(ph3d,2)/=npw_k) then 165 MSG_BUG('Incorrect size for ph3d!') 166 end if 167 end if 168 169 !Define dimensions of projected scalars 170 dimffnl=size(ffnl,2) 171 ndgxdt=0;nd2gxdt=0 172 if (idir==0) then 173 if (choice==2) ndgxdt=3 174 if (choice==3) ndgxdt=6 175 if (choice==23) ndgxdt=9 176 if (choice==4) nd2gxdt=6 177 if (choice==24) then 178 ndgxdt=3;nd2gxdt=6 179 end if 180 if (choice==5) ndgxdt=3 181 if (choice==6) then 182 ndgxdt=9;nd2gxdt=54 183 end if 184 else 185 ndgxdt=1 186 end if 187 if(cwaveprj(1,1)%ncpgr<ndgxdt+nd2gxdt) then 188 MSG_BUG('Incorrect size for ncpgr') 189 end if 190 !Eventually re-compute (k+G) vectors (and related data) 191 if (nkpg==0) then 192 nkpg_=0 193 if (choice==4.or.choice==24) nkpg_=9 194 if (choice==2.or.choice==3.or.choice==23) nkpg_=3 195 ABI_ALLOCATE(kpg_,(npw_k,nkpg_)) 196 if (nkpg_>0) then 197 call mkkpg(kg_k,kpg_,kpoint,nkpg_,npw_k) 198 end if 199 else 200 nkpg_=nkpg 201 kpg_ => kpg 202 end if 203 204 !Some other dims 205 mincat=min(NLO_MINCAT,maxval(nattyp)) 206 cplex=2;if (istwf_k>1) cplex=1 207 choice_=choice;if (cpopt==1) choice_=-choice 208 signs=1;if (idir>0) signs=2 209 !Eventually allocate temporary array for ph3d 210 if (nloalg(2)<=0) then 211 matblk=mincat 212 ABI_ALLOCATE(ph3d_,(2,npw_k,matblk)) 213 else 214 matblk=size(ph3d,3) 215 ph3d_ => ph3d 216 end if 217 218 !Loop over atom types 219 ia1=1;iatm=0 220 do itypat=1,ntypat 221 ia2=ia1+nattyp(itypat)-1;if (ia2<ia1) cycle 222 nlmn=count(indlmn(3,:,itypat)>0) 223 224 ! Retrieve some data for this type of atom 225 ABI_ALLOCATE(indlmn_typ,(6,nlmn)) 226 ABI_ALLOCATE(ffnl_typ,(npw_k,dimffnl,nlmn)) 227 indlmn_typ(:,1:nlmn)=indlmn(:,1:nlmn,itypat) 228 ffnl_typ(:,:,1:nlmn)=ffnl(:,:,1:nlmn,itypat) 229 230 ! Loop on blocks of atoms inside type 231 do ia3=ia1,ia2,mincat 232 ia4=min(ia2,ia3+mincat-1);nincat=ia4-ia3+1 233 ! Prepare the phase factors if they were not already computed 234 if (nloalg(2)<=0) then 235 call ph1d3d(ia3,ia4,kg_k,matblk,natom,npw_k,ngfft(1),ngfft(2),ngfft(3),& 236 & phkxred,ph1d,ph3d_) 237 end if 238 239 ! Allocate memory for projected scalars 240 ABI_ALLOCATE(gx,(cplex,nlmn,nincat,nspinor)) 241 ABI_ALLOCATE(dgxdt,(cplex,ndgxdt,nlmn,nincat,nspinor)) 242 ABI_ALLOCATE(d2gxdt,(cplex,nd2gxdt,nlmn,nincat,nspinor)) 243 ABI_ALLOCATE(cplex_dgxdt,(ndgxdt)) 244 ABI_ALLOCATE(cplex_d2gxdt,(nd2gxdt)) 245 246 ! Retrieve eventually <p_i|c> coeffs 247 if (cpopt==1) then 248 do ispinor=1,nspinor 249 do ia=1,nincat 250 gx(1:cplex,1:nlmn,ia,ispinor)=cwaveprj(iatm+ia,ispinor)%cp(1:cplex,1:nlmn) 251 end do 252 end do 253 end if 254 255 ! Compute <p_i|c> scalars (and derivatives) for this block of atoms 256 call opernla_ylm(choice_,cplex,cplex_dgxdt,cplex_d2gxdt,dimffnl,d2gxdt,dgxdt,ffnl_typ,gx,& 257 & ia3,idir,indlmn_typ,istwf_k,kpg_,matblk,mpi_enreg,nd2gxdt,ndgxdt,nincat,nkpg_,nlmn,& 258 & nloalg,npw_k,nspinor,ph3d_,signs,ucvol,cwavef) 259 260 ! Transfer result to output variable cwaveprj 261 if (cpopt==0) then 262 do ispinor=1,nspinor 263 do ia=1,nincat 264 cwaveprj(iatm+ia,ispinor)%nlmn=nlmn 265 cwaveprj(iatm+ia,ispinor)%cp(1:cplex,1:nlmn)=gx(1:cplex,1:nlmn,ia,ispinor) 266 if(cplex==1) cwaveprj(iatm+ia,ispinor)%cp(2,1:nlmn)=zero 267 end do 268 end do 269 end if 270 if (cpopt>=0.and.choice>1) then 271 ishift=0 272 if ((idir>0).and.(cwaveprj(1,1)%ncpgr>ndgxdt)) ishift=idir-1 273 if(cplex==2)then 274 do ispinor=1,nspinor 275 do ia=1,nincat 276 ! cwaveprj(iatm+ia,ispinor)%ncpgr=ndgxdt+nd2gxdt 277 if (ndgxdt>0) cwaveprj(iatm+ia,ispinor)%dcp(1:2,1+ishift:ndgxdt+ishift,1:nlmn)=& 278 & dgxdt(1:2,1:ndgxdt,1:nlmn,ia,ispinor) 279 if (nd2gxdt>0)cwaveprj(iatm+ia,ispinor)%dcp(1:2,ndgxdt+1+ishift:ndgxdt+nd2gxdt+ishift,1:nlmn)=& 280 & d2gxdt(1:2,1:nd2gxdt,1:nlmn,ia,ispinor) 281 end do 282 end do 283 else 284 ! cplex_dgxdt(i) = 1 if dgxdt(1,i,:,:) is real, 2 if it is pure imaginary 285 ! cplex_d2gxdt(i) = 1 if d2gxdt(1,i,:,:) is real, 2 if it is pure imaginary 286 do ispinor=1,nspinor 287 do ia=1,nincat 288 ! cwaveprj(iatm+ia,ispinor)%ncpgr=ndgxdt+nd2gxdt 289 if (ndgxdt>0) then 290 do ilmn =1,nlmn 291 do ii = 1,ndgxdt 292 ic = cplex_dgxdt(ii) ; jc = 3 - ic 293 cwaveprj(iatm+ia,ispinor)%dcp(ic,ii+ishift,ilmn)=dgxdt(1,ii,ilmn,ia,ispinor) 294 cwaveprj(iatm+ia,ispinor)%dcp(jc,ii+ishift,ilmn)=zero 295 end do 296 end do 297 end if 298 if (nd2gxdt>0) then 299 do ilmn =1,nlmn 300 do ii = 1,nd2gxdt 301 ic = cplex_d2gxdt(ii) ; jc = 3 - ic 302 cwaveprj(iatm+ia,ispinor)%dcp(ic,ndgxdt+ii+ishift,ilmn)=d2gxdt(1,ii,ilmn,ia,ispinor) 303 cwaveprj(iatm+ia,ispinor)%dcp(jc,ndgxdt+ii+ishift,ilmn)=zero 304 end do 305 end do 306 end if 307 end do 308 end do 309 end if 310 end if 311 312 ! End loop inside block of atoms 313 iatm=iatm+nincat 314 ABI_DEALLOCATE(gx) 315 ABI_DEALLOCATE(dgxdt) 316 ABI_DEALLOCATE(d2gxdt) 317 ABI_DEALLOCATE(cplex_dgxdt) 318 ABI_DEALLOCATE(cplex_d2gxdt) 319 end do 320 321 ! End loop over atom types 322 ia1=ia2+1 323 ABI_DEALLOCATE(indlmn_typ) 324 ABI_DEALLOCATE(ffnl_typ) 325 end do 326 327 if (nkpg==0) then 328 ABI_DEALLOCATE(kpg_) 329 end if 330 if (nloalg(2)<=0) then 331 ABI_DEALLOCATE(ph3d_) 332 end if 333 334 DBG_EXIT('COLL') 335 336 end subroutine getcprj