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