TABLE OF CONTENTS
ABINIT/dfpt_nsteltwf [ Functions ]
NAME
dfpt_nsteltwf
FUNCTION
This routine computes the non-local and kinetic contribution to the 2DTE matrix elements, in the non-stationary formulation
COPYRIGHT
Copyright (C) 1999-2018 ABINIT group (DRH,XG,AR,MB,MVer) 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/doc/developers/contributors.txt .
INPUTS
cg(2,mpw*nspinor*mband*mkmem*nsppol)=planewave coefficients of wavefunctions
cg1(2,mpw1*nspinor*mband*mk1mem*nsppol)=pw coefficients of RF
wavefunctions at k,q.
ecut=cut-off energy for plane wave basis sphere (Ha)
ecutsm=smearing energy for plane wave kinetic energy (Ha) (NOT NEEDED !)
effmass_free=effective mass for electrons (1. in common case)
gs_hamk <type(gs_hamiltonian_type)>=all data for the Hamiltonian at k
icg=shift to be applied on the location of data in the array cg
icg1=shift to be applied on the location of data in the array cg1
ikpt=number of the k-point
isppol=1 for unpolarized, 2 for spin-polarized
istwf_k=flag controlling the storage of WFs
kg_k(3,npw_k)=reduced planewave coordinates.
kg1_k(3,npw1_k)=reduced planewave coordinates at k+q, with RF k points
kpoint(3)=k-point in reduced coordinates
mband=maximum number of bands
mkmem =number of k points treated by this node.
mk1mem =number of k points treated by this node (RF data).
mpert =maximum number of ipert
mpi_enreg=information about MPI parallelization
mpw=maximum dimensioned size of npw or wfs at k
mpw1=maximum dimensioned size of npw for wfs at k+q (also for 1-order wfs).
natom=number of atoms in cell.
nband_k=number of bands at this k point for that spin polarization
npw_k=number of plane waves at this k point
npw1_k=number of plane waves at this k+q 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.
occ_k(nband_k)=occupation number for each band (usually 2) for each k.
psps <type(pseudopotential_type)>=variables related to pseudopotentials
rmet(3,3)=real space metric (bohr**2)
wtk_k=weight assigned to the k point.
ylm(npw_k,mpsang*mpsang)= real spherical harmonics for each G and k point
ylmgr(npw_k,3,mpsang*mpsang*useylm)= gradients of real spherical for each G and k point
OUTPUT
d2nl_k(2,3,mpert)=non-local contributions to non-stationary 2DTE, for the present k point, and perturbation idir, ipert
PARENTS
dfpt_nselt
CHILDREN
dotprod_g,kpgstr,load_k_hamiltonian,mkffnl,mkkin,nonlop
SOURCE
66 #if defined HAVE_CONFIG_H 67 #include "config.h" 68 #endif 69 70 #include "abi_common.h" 71 72 73 subroutine dfpt_nsteltwf(cg,cg1,d2nl_k,ecut,ecutsm,effmass_free,gs_hamk,icg,icg1,ikpt,isppol,& 74 & istwf_k,kg_k,kg1_k,kpoint,mband,mkmem,mk1mem,mpert,mpi_enreg,mpw,mpw1,natom,nband_k,& 75 & npw_k,npw1_k,nspinor,nsppol,ntypat,occ_k,psps,rmet,wtk_k,ylm,ylmgr) 76 77 use defs_basis 78 use defs_datatypes 79 use defs_abitypes 80 use m_wffile 81 use m_xmpi 82 use m_profiling_abi 83 use m_cgtools 84 85 use m_pawcprj, only : pawcprj_type 86 use m_hamiltonian, only : gs_hamiltonian_type,load_k_hamiltonian 87 use m_kg, only : mkkin, kpgstr 88 89 !This section has been created automatically by the script Abilint (TD). 90 !Do not modify the following lines by hand. 91 #undef ABI_FUNC 92 #define ABI_FUNC 'dfpt_nsteltwf' 93 use interfaces_66_nonlocal 94 !End of the abilint section 95 96 implicit none 97 98 !Arguments ------------------------------------ 99 !scalars 100 integer,intent(in) :: icg,icg1,ikpt,isppol,istwf_k,mband,mk1mem,mkmem,mpert,mpw,mpw1,natom 101 integer,intent(in) :: nspinor,nsppol,ntypat 102 integer,intent(inout) :: nband_k,npw1_k,npw_k 103 real(dp),intent(in) :: ecut,ecutsm,effmass_free,wtk_k 104 type(MPI_type),intent(in) :: mpi_enreg 105 type(pseudopotential_type),intent(in) :: psps 106 !arrays 107 integer,intent(in) :: kg1_k(3,npw1_k),kg_k(3,npw_k) 108 real(dp),intent(in) :: kpoint(3) 109 real(dp),intent(in) :: cg(2,mpw*nspinor*mband*mkmem*nsppol) 110 real(dp),intent(in) :: cg1(2,mpw1*nspinor*mband*mk1mem*nsppol) 111 real(dp),intent(in) :: occ_k(nband_k),rmet(3,3) 112 real(dp),intent(in) :: ylm(npw_k,psps%mpsang*psps%mpsang) 113 real(dp),intent(in) :: ylmgr(npw_k,3,psps%mpsang*psps%mpsang) 114 real(dp),intent(inout) :: d2nl_k(2,3,mpert) 115 116 !Local variables------------------------------- 117 !scalars 118 integer :: choice,cpopt,dimffnl,dimffnl2,iband 119 integer :: ider,idir0,idir1,ilmn,ipert1,ipw,ipws,ispinor,istr1,itypat 120 integer :: nkpg,nnlout,paw_opt,signs,tim_nonlop 121 real(dp) :: doti,dotr 122 type(gs_hamiltonian_type) :: gs_hamk 123 !arrays 124 real(dp) :: enlout(6),dum_svectout(1,1),dum(1),kpg_dum(0,0) 125 real(dp),allocatable :: cwave0(:,:),cwavef(:,:),dkinpw(:),eig2_k(:) 126 real(dp),allocatable :: ffnl(:,:,:,:),ffnl_ylm(:,:,:,:),ghc(:,:) 127 real(dp),allocatable :: gvnl1(:,:),gvnlc(:,:),kinpw1(:),ph3d(:,:,:) 128 type(pawcprj_type) :: cprj_dum(0,0) 129 130 ! ********************************************************************* 131 132 133 !DEBUG 134 !write(std_out,*)' dfpt_nstwf : enter ' 135 !ENDDEBUG 136 137 !Init me 138 ABI_ALLOCATE(ghc,(2,npw1_k*nspinor)) 139 ABI_ALLOCATE(gvnlc,(2,npw1_k*nspinor)) 140 ABI_ALLOCATE(gvnl1,(2,npw1_k*nspinor)) 141 ABI_ALLOCATE(eig2_k,(2*nsppol*mband**2)) 142 ABI_ALLOCATE(kinpw1,(npw1_k)) 143 ABI_ALLOCATE(dkinpw,(npw_k)) 144 nkpg=0 145 146 !Compute nonlocal form factors ffnl at (k+G), for all atoms 147 dimffnl=2 148 ABI_ALLOCATE(ffnl,(npw_k,dimffnl,psps%lmnmax,ntypat)) 149 if (psps%useylm==0) then 150 ider=1;idir0=0 151 call mkffnl(psps%dimekb,dimffnl,psps%ekb,ffnl,psps%ffspl,gs_hamk%gmet,gs_hamk%gprimd,ider,idir0,& 152 & psps%indlmn,kg_k,kpg_dum,kpoint,psps%lmnmax,psps%lnmax,psps%mpsang,psps%mqgrid_ff,nkpg,& 153 & npw_k,ntypat,psps%pspso,psps%qgrid_ff,rmet,psps%usepaw,psps%useylm,ylm,ylmgr) 154 else 155 ider=1;idir0=-7;dimffnl2=7 156 ABI_ALLOCATE(ffnl_ylm,(npw_k,dimffnl2,psps%lmnmax,ntypat)) 157 call mkffnl(psps%dimekb,dimffnl2,psps%ekb,ffnl_ylm,psps%ffspl,gs_hamk%gmet,gs_hamk%gprimd,& 158 & ider,idir0,psps%indlmn,kg_k,kpg_dum,kpoint,psps%lmnmax,psps%lnmax,psps%mpsang,psps%mqgrid_ff,& 159 & nkpg,npw_k,ntypat,psps%pspso,psps%qgrid_ff,rmet,psps%usepaw,psps%useylm,ylm,ylmgr) 160 do itypat=1,ntypat 161 do ilmn=1,psps%lmnmax 162 ffnl(:,1,ilmn,itypat)=ffnl_ylm(:,1,ilmn,itypat) 163 end do 164 end do 165 end if 166 167 !Compute kinetic contributions (1/2) (2 Pi)**2 (k+G)**2: 168 ! call mkkin(ecut,ecutsm,effmass_free,gs_hamk%gmet,kg1_k,kinpw1,kpoint,npw1_k) 169 call mkkin(ecut,ecutsm,effmass_free,gs_hamk%gmet,kg1_k,kinpw1,kpoint,npw1_k,0,0) 170 171 !Load k/k+q-dependent part in the Hamiltonian datastructure 172 ABI_ALLOCATE(ph3d,(2,npw_k,gs_hamk%matblk)) 173 call load_k_hamiltonian(gs_hamk,kpt_k=kpoint,npw_k=npw_k,istwf_k=istwf_k,kg_k=kg_k,ffnl_k=ffnl,& 174 & ph3d_k=ph3d,compute_ph3d=.true.) 175 176 ABI_ALLOCATE(cwave0,(2,npw_k*nspinor)) 177 ABI_ALLOCATE(cwavef,(2,npw1_k*nspinor)) 178 179 !Loop over bands 180 do iband=1,nband_k 181 182 if(mpi_enreg%proc_distrb(ikpt, iband, isppol) /= mpi_enreg%me_kpt) then 183 ! Skip the eigenvalue and the gvnl records of this band 184 cycle 185 end if 186 187 ! Get ground-state and first-order wavefunctions 188 cwave0(:,:)=cg(:,1+(iband-1)*npw_k*nspinor+icg:iband*npw_k*nspinor+icg) 189 cwavef(:,:)=cg1(:,1+(iband-1)*npw1_k*nspinor+icg1:iband*npw1_k*nspinor+icg1) 190 191 ! Double loop over strain perturbations 192 do ipert1=natom+3,natom+4 193 do idir1=1,3 194 if (ipert1==natom+3) istr1=idir1 195 if (ipert1==natom+4) istr1=idir1+3 196 197 ! Compute the derivative of the kinetic operator vs strain in dkinpw 198 call kpgstr(dkinpw,ecut,ecutsm,effmass_free,gs_hamk%gmet,gs_hamk%gprimd,istr1,& 199 & kg1_k,kpoint,npw1_k) 200 201 ! Get |vnon-locj1|u0> : 202 ! first-order non-local, applied to zero-order wavefunction 203 ! (??) this routine gives MINUS the non-local contribution 204 205 ! When using Ylms, load the correct ffnl derivative 206 if (psps%useylm==1) then 207 do itypat=1,ntypat 208 do ilmn=1,psps%lmnmax 209 ffnl(:,2,ilmn,itypat)=ffnl_ylm(:,1+istr1,ilmn,itypat) 210 end do 211 end do 212 end if 213 214 signs=2 ; choice=3 ; nnlout=6 ; paw_opt=0 ; cpopt=-1 ; tim_nonlop=5 215 call nonlop(choice,cpopt,cprj_dum,enlout,gs_hamk,istr1,dum,mpi_enreg,1,nnlout,paw_opt,& 216 & signs,dum_svectout,tim_nonlop,cwave0,gvnl1) 217 ! <G|Vnl1|Cnk> is contained in gvnl1 218 219 ! Kinetic contribution 220 do ispinor=1,nspinor 221 do ipw=1,npw1_k 222 ipws=ipw+npw1_k*(ispinor-1) 223 if(kinpw1(ipw)<huge(0.0_dp)*1.d-11)then 224 gvnl1(1,ipws)=gvnl1(1,ipws)+dkinpw(ipw)*cwave0(1,ipws) 225 gvnl1(2,ipws)=gvnl1(2,ipws)+dkinpw(ipw)*cwave0(2,ipws) 226 else 227 gvnl1(1,ipws)=0.0_dp 228 gvnl1(2,ipws)=0.0_dp 229 end if 230 end do 231 end do 232 233 ! construct the matrix element (<uj2|vj1|u0>)complex conjug. 234 ! and add it to the 2nd-order matrix 235 ! imaginary term should be zero for strain-strain 2nd derivatives, 236 ! but keep it as a test for now 237 call dotprod_g(dotr,doti,gs_hamk%istwf_k,npw1_k*nspinor,2,cwavef,gvnl1,& 238 & mpi_enreg%me_g0,mpi_enreg%comm_spinorfft) 239 240 d2nl_k(1,idir1,ipert1)= d2nl_k(1,idir1,ipert1)+wtk_k*occ_k(iband)*2.0_dp*dotr 241 d2nl_k(2,idir1,ipert1)= d2nl_k(2,idir1,ipert1)-wtk_k*occ_k(iband)*2.0_dp*doti 242 243 end do !idir1 244 end do !ipert1 245 246 ! UNTIL NOW, DO NOT TAKE INTO ACCOUNT istwf_k 247 248 ! End loop over bands 249 end do 250 251 ABI_DEALLOCATE(cwave0) 252 ABI_DEALLOCATE(cwavef) 253 254 !################################################################### 255 256 ABI_DEALLOCATE(eig2_k) 257 ABI_DEALLOCATE(ghc) 258 ABI_DEALLOCATE(gvnlc) 259 ABI_DEALLOCATE(gvnl1) 260 ABI_DEALLOCATE(kinpw1) 261 ABI_DEALLOCATE(dkinpw) 262 ABI_DEALLOCATE(ffnl) 263 ABI_DEALLOCATE(ph3d) 264 if (psps%useylm==1) then 265 ABI_DEALLOCATE(ffnl_ylm) 266 end if 267 268 end subroutine dfpt_nsteltwf