TABLE OF CONTENTS
ABINIT/dfptff_edie [ Functions ]
NAME
dfptff_edie
FUNCTION
calculate the second order energy from the contribution of \Omega E \cdot P term, Eq.(6) in PRB 75, 115116(2007).
COPYRIGHT
Copyright (C) 2004-2018 ABINIT group (XW). 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 .
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. dtefield = variables related to response Berry-phase calculation ikpt = the index of the current k point isppol = the index of the spin component mband = maximum number of bands mkmem = maximum number of k-points in core memory mpw = maximum number of plane waves mpw1 = maximum number of plane waves for response wavefunctions nkpt = number of k points npwarr(nkpt) = number of planewaves in basis and boundary at this k point npwar1(nkpt) = number of planewaves in basis and boundary for response wfs nspinor = 1 for scalar wfs, 2 for spinor wfs nsppol = 1 for unpolarized, 2 for spin-polarized qmat(2,dtefield%mband_occ,dtefield%mband_occ,nkpt,2,3) = inverse of the overlap matrix pwindall(max(mpw,mpw1)*mkmem,8,3) = array used to compute the overlap matrices pwindall(:,1,:) <- <u^(0)_i|u^(0)_i+1> pwindall(:,2,:) <- <u^(0)_i|u^(0)_i-1> pwindall(:,3,:) <- <u^(1)_i|u^(1)_i+1> pwindall(:,4,:) <- <u^(1)_i|u^(1)_i-1> pwindall(:,5,:) <- <u^(1)_i|u^(0)_i+n+1> pwindall(:,6,:) <- <u^(1)_i|u^(0)_i+n-1> pwindall(:,7,:) <- <u^(0)_i|u^(1)_i-n+1> pwindall(:,8,:) <- <u^(0)_i|u^(1)_i-n-1>
OUTPUT
eberry = the energy of the Berry phase term
PARENTS
dfpt_scfcv
CHILDREN
overlap_g
SOURCE
56 #if defined HAVE_CONFIG_H 57 #include "config.h" 58 #endif 59 60 #include "abi_common.h" 61 62 63 subroutine dfptff_edie(cg,cg1,dtefield,eberry,idir_efield,mband,mkmem,& 64 & mpw,mpw1,nkpt,npwarr,npwar1,nsppol,nspinor,pwindall,qmat,rprimd) 65 66 67 use defs_basis 68 use m_profiling_abi 69 use m_efield 70 71 use m_cgtools, only : overlap_g 72 73 !This section has been created automatically by the script Abilint (TD). 74 !Do not modify the following lines by hand. 75 #undef ABI_FUNC 76 #define ABI_FUNC 'dfptff_edie' 77 !End of the abilint section 78 79 implicit none 80 81 !Arguments ---------------------------------------- 82 !scalars 83 integer,intent(in) :: idir_efield,mband,mkmem,mpw,mpw1,nkpt,nspinor,nsppol 84 real(dp),intent(out) :: eberry 85 type(efield_type),intent(in) :: dtefield 86 !arrays 87 integer,intent(in) :: npwar1(nkpt),npwarr(nkpt) 88 integer,intent(in) :: pwindall(max(mpw,mpw1)*mkmem,8,3) 89 real(dp),intent(in) :: cg(2,mpw*mband*mkmem*nspinor*nsppol) 90 real(dp),intent(in) :: cg1(2,mpw1*mband*mkmem*nspinor*nsppol) 91 real(dp),intent(in) :: qmat(2,dtefield%mband_occ,dtefield%mband_occ,nkpt,2,3) 92 real(dp),intent(in) :: rprimd(3,3) 93 94 !Local variables ---------------------------------- 95 !scalars 96 integer :: iband,icg,icg1,idir 97 integer :: ikpt,ikpt1,ikptn,ikptnm 98 integer :: jband,kband,mpw_tmp,npw_k1,npw_k2,pwmax,pwmin 99 real(dp) :: doti,dotr,e0,fac 100 !arrays 101 integer,allocatable :: pwind_tmp(:) 102 real(dp) :: z1(2) 103 real(dp),allocatable :: Amat(:,:,:),umat(:,:,:,:),vect1(:,:),vect2(:,:) 104 105 ! ************************************************************************* 106 107 !calculate 4 matrices ----------------------------- 108 mpw_tmp=max(mpw,mpw1) 109 ABI_ALLOCATE(umat,(2,dtefield%mband_occ,dtefield%mband_occ,4)) 110 ABI_ALLOCATE(vect1,(2,0:mpw_tmp)) 111 ABI_ALLOCATE(vect2,(2,0:mpw_tmp)) 112 ABI_ALLOCATE(pwind_tmp,(mpw_tmp)) 113 ABI_ALLOCATE(Amat,(2,dtefield%mband_occ,dtefield%mband_occ)) 114 vect1(:,0) = zero ; vect2(:,0) = zero 115 eberry=zero 116 117 do ikpt=1,nkpt 118 do idir=1,3 119 fac = dtefield%efield_dot(idir)/& 120 & (dble(dtefield%nstr(idir))*four_pi) 121 122 ! compute <u^(1)_{k_j,q}|u^(1)_{k_j+1,q}> matrix---------------------------------------------------- 123 124 ! prepare to calculate overlap matrix 125 ikpt1 = dtefield%ikpt_dk(ikpt,1,idir) 126 icg = dtefield%cgindex(ikpt,1+nsppol) 127 icg1 = dtefield%cgindex(ikpt1,1+nsppol) 128 npw_k1 = npwar1(ikpt) 129 npw_k2 = npwar1(ikpt1) 130 pwind_tmp(1:npw_k1) = pwindall((ikpt-1)*mpw_tmp+1:(ikpt-1)*mpw_tmp+npw_k1,3,idir) 131 vect1(:,0) = zero ; vect2(:,0) = zero 132 do jband = 1, dtefield%mband_occ 133 vect2(:,1:npw_k2) = & 134 & cg1(:,icg1 + 1 + (jband-1)*npw_k2*nspinor:icg1 + jband*npw_k2*nspinor) 135 if (npw_k2 < mpw_tmp) vect2(:,npw_k2+1:mpw_tmp) = zero 136 do iband = 1, dtefield%mband_occ 137 pwmin = (iband-1)*npw_k1*nspinor 138 pwmax = pwmin + npw_k1*nspinor 139 vect1(:,1:npw_k1) = & 140 & cg1(:,icg + 1 + pwmin:icg + pwmax) 141 if (npw_k1 < mpw_tmp) vect1(:,npw_k1+1:mpw_tmp) = zero 142 call overlap_g(doti,dotr,mpw_tmp,npw_k1,npw_k2,nspinor,pwind_tmp,& 143 & vect1,vect2) 144 umat(1,iband,jband,1) = dotr 145 umat(2,iband,jband,1) = doti 146 end do ! iband 147 end do !jband 148 149 ! compute <u^(0)_{k_j}|u^(1)_{k_j-n+1,q}> matrix---------------------------------------------------- 150 151 ! prepare to calculate overlap matrix 152 ikpt1 = dtefield%ikpt_dk(ikpt,5,idir) 153 icg = dtefield%cgindex(ikpt,1) 154 icg1 = dtefield%cgindex(ikpt1,1+nsppol) 155 npw_k1 = npwarr(ikpt) 156 npw_k2 = npwar1(ikpt1) 157 pwind_tmp(1:npw_k1) = pwindall((ikpt-1)*mpw_tmp+1:(ikpt-1)*mpw_tmp+npw_k1,7,idir) 158 159 vect1(:,0) = zero ; vect2(:,0) = zero 160 do jband = 1, dtefield%mband_occ 161 vect2(:,1:npw_k2) = & 162 & cg1(:,icg1 + 1 + (jband-1)*npw_k2*nspinor:icg1 + jband*npw_k2*nspinor) 163 if (npw_k2 < mpw_tmp) vect2(:,npw_k2+1:mpw_tmp) = zero 164 do iband = 1, dtefield%mband_occ 165 pwmin = (iband-1)*npw_k1*nspinor 166 pwmax = pwmin + npw_k1*nspinor 167 vect1(:,1:npw_k1) = & 168 & cg(:,icg + 1 + pwmin:icg + pwmax) 169 if (npw_k1 < mpw_tmp) vect1(:,npw_k1+1:mpw_tmp) = zero 170 call overlap_g(doti,dotr,mpw_tmp,npw_k1,npw_k2,nspinor,pwind_tmp,& 171 & vect1,vect2) 172 umat(1,iband,jband,2) = dotr 173 umat(2,iband,jband,2) = doti 174 end do ! iband 175 end do !jband 176 177 ! compute <u^(1)_{k_j-n,q}|u^(0)_{k_j+1}> matrix---------------------------------------------------- 178 179 ! prepare to calculate overlap matrix 180 ikptn = dtefield%ikpt_dk(ikpt,8,idir) 181 ikpt1 = dtefield%ikpt_dk(ikpt,1,idir) 182 icg = dtefield%cgindex(ikptn,1+nsppol) 183 icg1 = dtefield%cgindex(ikpt1,1) 184 npw_k1 = npwar1(ikptn) 185 npw_k2 = npwarr(ikpt1) 186 pwind_tmp(1:npw_k1) = pwindall((ikptn-1)*mpw_tmp+1:(ikptn-1)*mpw_tmp+npw_k1,5,idir) 187 188 vect1(:,0) = zero ; vect2(:,0) = zero 189 do jband = 1, dtefield%mband_occ 190 vect2(:,1:npw_k2) = & 191 & cg(:,icg1 + 1 + (jband-1)*npw_k2*nspinor:icg1 + jband*npw_k2*nspinor) 192 if (npw_k2 < mpw_tmp) vect2(:,npw_k2+1:mpw_tmp) = zero 193 do iband = 1, dtefield%mband_occ 194 pwmin = (iband-1)*npw_k1*nspinor 195 pwmax = pwmin + npw_k1*nspinor 196 vect1(:,1:npw_k1) = & 197 & cg1(:,icg + 1 + pwmin:icg + pwmax) 198 if (npw_k1 < mpw_tmp) vect1(:,npw_k1+1:mpw_tmp) = zero 199 call overlap_g(doti,dotr,mpw_tmp,npw_k1,npw_k2,nspinor,pwind_tmp,& 200 & vect1,vect2) 201 umat(1,iband,jband,3) = dotr 202 umat(2,iband,jband,3) = doti 203 end do ! iband 204 end do !jband 205 206 ! compute <u^(0)_{-k_j-n+1}|u^(1)_{-k_j,q}> matrix---------------------------------------------------- 207 208 ! prepare to calculate overlap matrix 209 ikptn = dtefield%ikpt_dk(ikpt,5,idir) 210 ikptnm = dtefield%ikpt_dk(ikptn,9,idir) 211 ikpt1 = dtefield%ikpt_dk(ikpt,9,idir) 212 icg = dtefield%cgindex(ikptnm,1) 213 icg1 = dtefield%cgindex(ikpt1,1+nsppol) 214 npw_k1 = npwarr(ikptnm) 215 npw_k2 = npwar1(ikpt1) 216 pwind_tmp(1:npw_k1) = pwindall((ikptnm-1)*mpw_tmp+1:(ikptnm-1)*mpw_tmp+npw_k1,7,idir) 217 vect1(:,0) = zero ; vect2(:,0) = zero 218 do jband = 1, dtefield%mband_occ 219 vect2(:,1:npw_k2) = & 220 & cg1(:,icg1 + 1 + (jband-1)*npw_k2*nspinor:icg1 + jband*npw_k2*nspinor) 221 if (npw_k2 < mpw_tmp) vect2(:,npw_k2+1:mpw_tmp) = zero 222 do iband = 1, dtefield%mband_occ 223 pwmin = (iband-1)*npw_k1*nspinor 224 pwmax = pwmin + npw_k1*nspinor 225 vect1(:,1:npw_k1) = & 226 & cg(:,icg + 1 + pwmin:icg + pwmax) 227 if (npw_k1 < mpw_tmp) vect1(:,npw_k1+1:mpw_tmp) = zero 228 call overlap_g(doti,dotr,mpw_tmp,npw_k1,npw_k2,nspinor,pwind_tmp,& 229 & vect1,vect2) 230 umat(1,iband,jband,4) = dotr 231 umat(2,iband,jband,4) = doti 232 end do ! iband 233 end do !jband 234 235 ! sum over the whole------------------------------------------------------------ 236 237 e0=zero 238 do iband=1,dtefield%mband_occ 239 do jband=1,dtefield%mband_occ 240 e0 = e0 + 2_dp*(umat(1,iband,jband,1)*qmat(2,jband,iband,ikpt,1,idir)& 241 & + umat(2,iband,jband,1)*qmat(1,jband,iband,ikpt,1,idir)) 242 end do 243 end do 244 eberry = eberry - e0*fac 245 e0=zero 246 ikptn=dtefield%ikpt_dk(ikpt,8,idir) 247 Amat(:,:,:)=zero 248 ! calculate Amat 249 do iband=1, dtefield%mband_occ 250 do jband=1, dtefield%mband_occ 251 do kband=1, dtefield%mband_occ 252 Amat(1,iband,jband) = Amat(1,iband,jband) + (umat(1,iband,kband,3))*qmat(1,kband,jband,ikpt,1,idir)& 253 & - (umat(2,iband,kband,3))*qmat(2,kband,jband,ikpt,1,idir) 254 Amat(2,iband,jband) = Amat(2,iband,jband) + (umat(1,iband,kband,3))*qmat(2,kband,jband,ikpt,1,idir)& 255 & + (umat(2,iband,kband,3))*qmat(1,kband,jband,ikpt,1,idir) 256 end do 257 end do 258 end do 259 260 do iband=1, dtefield%mband_occ 261 do jband=1, dtefield%mband_occ 262 do kband=1, dtefield%mband_occ 263 z1(1) = (umat(1,jband,iband,4)+umat(1,iband,jband,2))*qmat(1,jband,kband,ikptn,1,idir)& 264 & - (umat(2,jband,iband,4)+umat(2,iband,jband,2))*qmat(2,jband,kband,ikptn,1,idir) 265 z1(2) = (umat(1,jband,iband,4)+umat(1,iband,jband,2))*qmat(2,jband,kband,ikptn,1,idir)& 266 & + (umat(2,jband,iband,4)+umat(2,iband,jband,2))*qmat(1,jband,kband,ikptn,1,idir) 267 268 e0 = e0 - 4_dp*(z1(1)*Amat(2,kband,iband)+z1(2)*Amat(1,kband,iband)) 269 end do 270 end do 271 end do 272 273 eberry = eberry - e0*fac 274 275 ! !---------------------------------last part--------------------------------------------- 276 277 fac = rprimd(idir_efield,idir)/& 278 & (dble(dtefield%nstr(idir))*two_pi) 279 280 ! compute <u^(1)_{k_j-n,q}|u^(0)_{k_j+1}> matrix---------------------------------------------------- 281 282 ! prepare to calculate overlap matrix 283 ikptn = dtefield%ikpt_dk(ikpt,8,idir) 284 ikpt1 = dtefield%ikpt_dk(ikpt,1,idir) 285 icg = dtefield%cgindex(ikptn,1+nsppol) 286 icg1 = dtefield%cgindex(ikpt1,1) 287 npw_k1 = npwar1(ikptn) 288 npw_k2 = npwarr(ikpt1) 289 pwind_tmp(1:npw_k1) = pwindall((ikptn-1)*mpw_tmp+1:(ikptn-1)*mpw_tmp+npw_k1,5,idir) 290 291 vect1(:,0) = zero ; vect2(:,0) = zero 292 do jband = 1, dtefield%mband_occ 293 vect2(:,1:npw_k2) = & 294 & cg(:,icg1 + 1 + (jband-1)*npw_k2*nspinor:icg1 + jband*npw_k2*nspinor) 295 if (npw_k2 < mpw_tmp) vect2(:,npw_k2+1:mpw_tmp) = zero 296 do iband = 1, dtefield%mband_occ 297 pwmin = (iband-1)*npw_k1*nspinor 298 pwmax = pwmin + npw_k1*nspinor 299 vect1(:,1:npw_k1) = & 300 & cg1(:,icg + 1 + pwmin:icg + pwmax) 301 if (npw_k1 < mpw_tmp) vect1(:,npw_k1+1:mpw_tmp) = zero 302 call overlap_g(doti,dotr,mpw_tmp,npw_k1,npw_k2,nspinor,pwind_tmp,& 303 & vect1,vect2) 304 umat(1,iband,jband,1) = dotr 305 umat(2,iband,jband,1) = doti 306 end do ! iband 307 end do !jband 308 309 ! compute <u^(0)_{k_j}|u^(1)_{k_j-n+1,q}> matrix---------------------------------------------------- 310 311 ! prepare to calculate overlap matrix 312 ikpt1 = dtefield%ikpt_dk(ikpt,5,idir) 313 icg = dtefield%cgindex(ikpt,1) 314 icg1 = dtefield%cgindex(ikpt1,1+nsppol) 315 npw_k1 = npwarr(ikpt) 316 npw_k2 = npwar1(ikpt1) 317 pwind_tmp(1:npw_k1) = pwindall((ikpt-1)*mpw_tmp+1:(ikpt-1)*mpw_tmp+npw_k1,7,idir) 318 vect1(:,0) = zero ; vect2(:,0) = zero 319 do jband = 1, dtefield%mband_occ 320 vect2(:,1:npw_k2) = & 321 & cg1(:,icg1 + 1 + (jband-1)*npw_k2*nspinor:icg1 + jband*npw_k2*nspinor) 322 if (npw_k2 < mpw_tmp) vect2(:,npw_k2+1:mpw_tmp) = zero 323 do iband = 1, dtefield%mband_occ 324 pwmin = (iband-1)*npw_k1*nspinor 325 pwmax = pwmin + npw_k1*nspinor 326 vect1(:,1:npw_k1) = & 327 & cg(:,icg + 1 + pwmin:icg + pwmax) 328 if (npw_k1 < mpw_tmp) vect1(:,npw_k1+1:mpw_tmp) = zero 329 call overlap_g(doti,dotr,mpw_tmp,npw_k1,npw_k2,nspinor,pwind_tmp,& 330 & vect1,vect2) 331 umat(1,iband,jband,1) = umat(1,iband,jband,1) + dotr 332 umat(2,iband,jband,1) = umat(2,iband,jband,1) + doti 333 end do ! iband 334 end do !jband 335 336 e0=zero 337 338 do iband=1,dtefield%mband_occ 339 do jband=1,dtefield%mband_occ 340 e0 = e0 + (umat(1,iband,jband,1)*qmat(2,jband,iband,ikpt,1,idir)& 341 & + umat(2,iband,jband,1)*qmat(1,jband,iband,ikpt,1,idir)) 342 end do 343 end do 344 345 eberry = eberry - e0*fac 346 347 end do !end idir 348 end do !end ikpt 349 350 ABI_DEALLOCATE(umat) 351 ABI_DEALLOCATE(vect1) 352 ABI_DEALLOCATE(vect2) 353 ABI_DEALLOCATE(pwind_tmp) 354 ABI_DEALLOCATE(Amat) 355 end subroutine dfptff_edie