TABLE OF CONTENTS
ABINIT/dfptff_gradberry [ Functions ]
NAME
dfptff_gradberry
FUNCTION
Calculation of the gradient of Berry-phase term in finite electric field.
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.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 finite electric field calculation ikpt = the index of the current k point isppol=1 for unpolarized, 2 for spin-polarized 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_vtorho
CHILDREN
overlap_g
SOURCE
55 #if defined HAVE_CONFIG_H 56 #include "config.h" 57 #endif 58 59 #include "abi_common.h" 60 61 62 subroutine dfptff_gradberry(cg,cg1,dtefield,grad_berry,ikpt,isppol,mband,mpw,mpw1,mkmem,mk1mem,nkpt,& 63 & npwarr,npwar1,nspinor,nsppol,qmat,pwindall) 64 65 66 use defs_basis 67 use m_profiling_abi 68 use m_efield 69 70 use m_cgtools, only : overlap_g 71 72 !This section has been created automatically by the script Abilint (TD). 73 !Do not modify the following lines by hand. 74 #undef ABI_FUNC 75 #define ABI_FUNC 'dfptff_gradberry' 76 !End of the abilint section 77 78 implicit none 79 80 !Arguments ---------------------------------------- 81 !scalars 82 integer,intent(in) :: ikpt,isppol,mband,mk1mem,mkmem,mpw,mpw1,nkpt,nspinor 83 integer,intent(in) :: nsppol 84 type(efield_type),intent(in) :: dtefield 85 !arrays 86 integer,intent(in) :: npwar1(nkpt),npwarr(nkpt) 87 integer,intent(in) :: pwindall(max(mpw,mpw1)*mkmem,8,3) 88 real(dp),intent(in) :: cg(2,mpw*nspinor*mband*mkmem*nsppol) 89 real(dp),intent(in) :: cg1(2,mpw1*nspinor*mband*mk1mem*nsppol) 90 real(dp),intent(in) :: qmat(2,dtefield%mband_occ,dtefield%mband_occ,nkpt,2,3) 91 real(dp),intent(out) :: grad_berry(2,mpw1,dtefield%mband_occ) 92 93 !Local variables ------------------------- 94 !scalars 95 integer :: iband,icg,icg1,idir,ikpt1 96 integer :: ikpt1m,ikptn,ikptnm,ikptnp1,ipw,jband,jpw,kband 97 integer :: mpw_tmp,npw_k1,npw_k2,pwmax,pwmin 98 real(dp) :: doti,dotr,fac,wfi,wfr 99 !arrays 100 integer,allocatable :: pwind_tmp(:) 101 real(dp) :: z1(2),z2(2) 102 real(dp),allocatable :: Amat(:,:,:),Bmat(:,:,:),s1mat(:,:,:),vect1(:,:) 103 real(dp),allocatable :: vect2(:,:) 104 105 ! ************************************************************************* 106 107 mpw_tmp=max(mpw,mpw1) 108 ABI_ALLOCATE(vect1,(2,0:mpw_tmp)) 109 ABI_ALLOCATE(vect2,(2,0:mpw_tmp)) 110 ABI_ALLOCATE(s1mat,(2,dtefield%mband_occ,dtefield%mband_occ)) 111 ABI_ALLOCATE(pwind_tmp,(mpw_tmp)) 112 ABI_ALLOCATE(Amat,(2,dtefield%mband_occ,dtefield%mband_occ)) 113 ABI_ALLOCATE(Bmat,(2,dtefield%mband_occ,dtefield%mband_occ)) 114 vect1(:,0) = zero ; vect2(:,0) = zero 115 s1mat(:,:,:)=zero 116 grad_berry(:,:,:) = zero 117 118 do idir=1,3 119 fac = dtefield%efield_dot(idir)*dble(nkpt)/& 120 & (dble(dtefield%nstr(idir))*four_pi) 121 122 ! prepare 123 ikpt1 = dtefield%ikpt_dk(ikpt,1,idir) 124 icg1 = dtefield%cgindex(ikpt1,isppol+nsppol) 125 npw_k1 = npwar1(ikpt) 126 npw_k2 = npwar1(ikpt1) 127 pwind_tmp(1:npw_k1) = pwindall((ikpt-1)*mpw_tmp+1:(ikpt-1)*mpw_tmp+npw_k1,3,idir) 128 129 do ipw = 1, npw_k1 130 jpw = pwind_tmp(ipw) 131 132 if (jpw > 0) then 133 do iband = 1, dtefield%mband_occ 134 wfr = cg1(1,icg1 + (iband - 1)*npw_k2*nspinor + jpw) 135 wfi = cg1(2,icg1 + (iband - 1)*npw_k2*nspinor + jpw) 136 137 do jband = 1, dtefield%mband_occ 138 139 grad_berry(1,ipw,jband) = & 140 & grad_berry(1,ipw,jband) + & 141 & fac*qmat(1,iband,jband,ikpt,1,idir)*wfr - fac*qmat(2,iband,jband,ikpt,1,idir)*wfi 142 143 grad_berry(2,ipw,jband) = & 144 & grad_berry(2,ipw,jband) + & 145 & fac*qmat(1,iband,jband,ikpt,1,idir)*wfi + fac*qmat(2,iband,jband,ikpt,1,idir)*wfr 146 147 end do 148 end do 149 end if 150 151 end do 152 153 ! compute <u^(0)_{k_j+n}|u^(1)_{k_j+1,q}> matrix---------------------------------------------------- 154 155 ! prepare to calculate overlap matrix 156 ikptn = dtefield%ikpt_dk(ikpt,7,idir) 157 ikpt1 = dtefield%ikpt_dk(ikpt,1,idir) 158 icg = dtefield%cgindex(ikptn,isppol) 159 icg1 = dtefield%cgindex(ikpt1,isppol+nsppol) 160 npw_k1 = npwarr(ikptn) 161 npw_k2 = npwar1(ikpt1) 162 pwind_tmp(1:npw_k1) = pwindall((ikptn-1)*mpw_tmp+1:(ikptn-1)*mpw_tmp+npw_k1,7,idir) 163 164 165 vect1(:,0) = zero ; vect2(:,0) = zero 166 do jband = 1, dtefield%mband_occ 167 vect2(:,1:npw_k2) = & 168 & cg1(:,icg1 + 1 + (jband-1)*npw_k2*nspinor:icg1 + jband*npw_k2*nspinor) 169 if (npw_k2 < mpw_tmp) vect2(:,npw_k2+1:mpw_tmp) = zero 170 171 do iband = 1, dtefield%mband_occ 172 173 pwmin = (iband-1)*npw_k1*nspinor 174 pwmax = pwmin + npw_k1*nspinor 175 vect1(:,1:npw_k1) = & 176 & cg(:,icg + 1 + pwmin:icg + pwmax) 177 if (npw_k1 < mpw_tmp) vect1(:,npw_k1+1:mpw_tmp) = zero 178 call overlap_g(doti,dotr,mpw_tmp,npw_k1,npw_k2,nspinor,pwind_tmp,& 179 & vect1,vect2) 180 s1mat(1,iband,jband) = dotr 181 s1mat(2,iband,jband) = doti 182 183 end do ! iband 184 end do !jband 185 186 ! compute <u^(0)_{-k_j+1}|u^(1)_{-k_j+n},q> matrix-------------------- 187 188 ! prepare to calculate overlap matrix 189 ikptn = dtefield%ikpt_dk(ikpt,7,idir) 190 ikptnm= dtefield%ikpt_dk(ikptn,9,idir) 191 ikpt1 = dtefield%ikpt_dk(ikpt,1,idir) 192 ikpt1m= dtefield%ikpt_dk(ikpt1,9,idir) 193 194 icg = dtefield%cgindex(ikpt1m,isppol) 195 icg1 = dtefield%cgindex(ikptnm,isppol+nsppol) 196 npw_k1 = npwarr(ikpt1m) 197 npw_k2 = npwar1(ikptnm) 198 pwind_tmp(1:npw_k1) = pwindall((ikpt1m-1)*mpw_tmp+1:(ikpt1m-1)*mpw_tmp+npw_k1,7,idir) 199 200 vect1(:,0) = zero ; vect2(:,0) = zero 201 do jband = 1, dtefield%mband_occ 202 vect2(:,1:npw_k2) = & 203 & cg1(:,icg1 + 1 + (jband-1)*npw_k2*nspinor:icg1 + jband*npw_k2*nspinor) 204 if (npw_k2 < mpw_tmp) vect2(:,npw_k2+1:mpw_tmp) = zero 205 206 do iband = 1, dtefield%mband_occ 207 208 pwmin = (iband-1)*npw_k1*nspinor 209 pwmax = pwmin + npw_k1*nspinor 210 vect1(:,1:npw_k1) = & 211 & cg(:,icg + 1 + pwmin:icg + pwmax) 212 if (npw_k1 < mpw_tmp) vect1(:,npw_k1+1:mpw_tmp) = zero 213 call overlap_g(doti,dotr,mpw_tmp,npw_k1,npw_k2,nspinor,pwind_tmp,& 214 & vect1,vect2) 215 216 s1mat(1,jband,iband) = s1mat(1,jband,iband) + dotr 217 s1mat(2,jband,iband) = s1mat(2,jband,iband) + doti 218 219 end do ! iband 220 end do !jband 221 222 Amat(:,:,:)=zero 223 224 ! calculate Amat 225 do iband=1, dtefield%mband_occ 226 do jband=1, dtefield%mband_occ 227 do kband=1, dtefield%mband_occ 228 Amat(1,iband,jband) = Amat(1,iband,jband) + s1mat(1,iband,kband)*& 229 & qmat(1,kband,jband,ikpt,1,idir)& 230 & - s1mat(2,iband,kband)*qmat(2,kband,jband,ikpt,1,idir) 231 Amat(2,iband,jband) = Amat(2,iband,jband) + s1mat(1,iband,kband)*& 232 & qmat(2,kband,jband,ikpt,1,idir)& 233 & + s1mat(2,iband,kband)*qmat(1,kband,jband,ikpt,1,idir) 234 end do 235 end do 236 end do 237 238 Bmat(:,:,:)=zero 239 240 ! calculate Bmat 241 ikptn = dtefield%ikpt_dk(ikpt,7,idir) 242 do iband=1, dtefield%mband_occ 243 do jband=1, dtefield%mband_occ 244 do kband=1, dtefield%mband_occ 245 Bmat(1,jband,kband) = Bmat(1,jband,kband) + Amat(1,iband,kband)*& 246 & qmat(1,jband,iband,ikptn,1,idir)& 247 & - Amat(2,iband,kband)*qmat(2,jband,iband,ikptn,1,idir) 248 Bmat(2,jband,kband) = Bmat(2,jband,kband) + Amat(1,iband,kband)*& 249 & qmat(2,jband,iband,ikptn,1,idir)& 250 & + Amat(2,iband,kband)*qmat(1,jband,iband,ikptn,1,idir) 251 end do 252 end do 253 end do 254 255 ! calc. the second term of gradient------------------------------ 256 257 ! preparation 258 259 ikptnp1 = dtefield%ikpt_dk(ikpt,3,idir) 260 icg = dtefield%cgindex(ikptnp1,isppol) 261 npw_k1 = npwar1(ikpt) 262 npw_k2 = npwarr(ikptnp1) 263 pwind_tmp(1:npw_k1) = pwindall((ikpt-1)*mpw_tmp+1:(ikpt-1)*mpw_tmp+npw_k1,5,idir) 264 265 z1(:) = zero 266 z2(:) = zero 267 268 do ipw = 1, npw_k1 269 270 jpw = pwind_tmp(ipw) 271 272 if (jpw > 0) then 273 274 do iband = 1, dtefield%mband_occ 275 wfr = cg(1,icg + (iband - 1)*npw_k2*nspinor + jpw) 276 wfi = cg(2,icg + (iband - 1)*npw_k2*nspinor + jpw) 277 278 do jband=1, dtefield%mband_occ 279 280 grad_berry(1,ipw,jband) = grad_berry(1,ipw,jband) & 281 & - fac*(Bmat(1,iband,jband)*wfr - Bmat(2,iband,jband)*wfi) 282 grad_berry(2,ipw,jband) = grad_berry(2,ipw,jband) & 283 & - fac*(Bmat(1,iband,jband)*wfi + Bmat(2,iband,jband)*wfr) 284 285 end do 286 end do 287 end if 288 end do 289 290 ! Second part of gradient of Berry phase++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 291 292 vect1(:,0) = zero ; vect2(:,0) = zero 293 294 ! prepare 295 ikpt1 = dtefield%ikpt_dk(ikpt,2,idir) 296 icg1 = dtefield%cgindex(ikpt1,isppol+nsppol) 297 npw_k1 = npwar1(ikpt) 298 npw_k2 = npwar1(ikpt1) 299 pwind_tmp(1:npw_k1) = pwindall((ikpt-1)*mpw_tmp+1:(ikpt-1)*mpw_tmp+npw_k1,4,idir) 300 301 ! write(std_out,*)'dfpt_cgwf:pwind_tmp',pwind_tmp 302 ! stop 303 304 do ipw = 1, npw_k1 305 306 jpw = pwind_tmp(ipw) 307 308 if (jpw > 0) then 309 do iband = 1, dtefield%mband_occ 310 wfr = cg1(1,icg1 + (iband - 1)*npw_k2*nspinor + jpw) 311 wfi = cg1(2,icg1 + (iband - 1)*npw_k2*nspinor + jpw) 312 313 do jband = 1, dtefield%mband_occ 314 315 grad_berry(1,ipw,jband) = & 316 & grad_berry(1,ipw,jband) - & 317 & fac*qmat(1,iband,jband,ikpt,2,idir)*wfr + fac*qmat(2,iband,jband,ikpt,2,idir)*wfi 318 319 grad_berry(2,ipw,jband) = & 320 & grad_berry(2,ipw,jband) - & 321 & fac*qmat(1,iband,jband,ikpt,2,idir)*wfi - fac*qmat(2,iband,jband,ikpt,2,idir)*wfr 322 323 end do 324 end do 325 end if 326 end do 327 328 329 330 331 ! compute <u^(0)_{k_j+n}|u^(1)_{k_j-1,q}> matrix---------------------------------------------------- 332 333 ! prepare to calculate overlap matrix 334 ikptn = dtefield%ikpt_dk(ikpt,7,idir) 335 ikpt1 = dtefield%ikpt_dk(ikpt,2,idir) 336 icg = dtefield%cgindex(ikptn,isppol) 337 icg1 = dtefield%cgindex(ikpt1,isppol+nsppol) 338 npw_k1 = npwarr(ikptn) 339 npw_k2 = npwar1(ikpt1) 340 pwind_tmp(1:npw_k1) =pwindall((ikptn-1)*mpw_tmp+1:(ikptn-1)*mpw_tmp+npw_k1,8,idir) 341 342 vect1(:,0) = zero ; vect2(:,0) = zero 343 do jband = 1, dtefield%mband_occ 344 vect2(:,1:npw_k2) = & 345 & cg1(:,icg1 + 1 + (jband-1)*npw_k2*nspinor:icg1 + jband*npw_k2*nspinor) 346 if (npw_k2 < mpw_tmp) vect2(:,npw_k2+1:mpw_tmp) = zero 347 348 do iband = 1, dtefield%mband_occ 349 350 pwmin = (iband-1)*npw_k1*nspinor 351 pwmax = pwmin + npw_k1*nspinor 352 vect1(:,1:npw_k1) = & 353 & cg(:,icg + 1 + pwmin:icg + pwmax) 354 if (npw_k1 < mpw_tmp) vect1(:,npw_k1+1:mpw_tmp) = zero 355 call overlap_g(doti,dotr,mpw_tmp,npw_k1,npw_k2,nspinor,pwind_tmp,& 356 & vect1,vect2) 357 s1mat(1,iband,jband) = dotr 358 s1mat(2,iband,jband) = doti 359 360 end do ! iband 361 end do !jband 362 363 ! compute <u^(0)_{-k_j-1}|u^(1)_{-k_j+n,q}> matrix----------------------------------------------------- 364 365 ! prepare to calculate overlap matrix 366 ikptn = dtefield%ikpt_dk(ikpt,7,idir) 367 ikptnm= dtefield%ikpt_dk(ikptn,9,idir) 368 ikpt1 = dtefield%ikpt_dk(ikpt,2,idir) 369 ikpt1m= dtefield%ikpt_dk(ikpt1,9,idir) 370 icg = dtefield%cgindex(ikpt1m,isppol) 371 icg1 = dtefield%cgindex(ikptnm,isppol+nsppol) 372 npw_k1 = npwarr(ikpt1m) 373 npw_k2 = npwar1(ikptnm) 374 pwind_tmp(1:npw_k1) =pwindall((ikpt1m-1)*mpw_tmp+1:(ikpt1m-1)*mpw_tmp+npw_k1,8,idir) 375 376 377 378 vect1(:,0) = zero ; vect2(:,0) = zero 379 do jband = 1, dtefield%mband_occ 380 vect2(:,1:npw_k2) = & 381 & cg1(:,icg1 + 1 + (jband-1)*npw_k2*nspinor:icg1 + jband*npw_k2*nspinor) 382 if (npw_k2 < mpw_tmp) vect2(:,npw_k2+1:mpw_tmp) = zero 383 do iband = 1, dtefield%mband_occ 384 pwmin = (iband-1)*npw_k1*nspinor 385 pwmax = pwmin + npw_k1*nspinor 386 vect1(:,1:npw_k1) = & 387 & cg(:,icg + 1 + pwmin:icg + pwmax) 388 if (npw_k1 < mpw_tmp) vect1(:,npw_k1+1:mpw_tmp) = zero 389 call overlap_g(doti,dotr,mpw_tmp,npw_k1,npw_k2,nspinor,pwind_tmp,& 390 & vect1,vect2) 391 s1mat(1,jband,iband) = s1mat(1,jband,iband) + dotr 392 s1mat(2,jband,iband) = s1mat(2,jband,iband) + doti 393 394 end do ! iband 395 end do !jband 396 397 Amat(:,:,:)=zero 398 399 ! calculate Amat 400 do iband=1, dtefield%mband_occ 401 do jband=1, dtefield%mband_occ 402 do kband=1, dtefield%mband_occ 403 Amat(1,iband,jband) = Amat(1,iband,jband) + s1mat(1,iband,kband)*& 404 & qmat(1,kband,jband,ikpt,2,idir)& 405 & - s1mat(2,iband,kband)*qmat(2,kband,jband,ikpt,2,idir) 406 Amat(2,iband,jband) = Amat(2,iband,jband) + s1mat(1,iband,kband)*& 407 & qmat(2,kband,jband,ikpt,2,idir)& 408 & + s1mat(2,iband,kband)*qmat(1,kband,jband,ikpt,2,idir) 409 end do 410 end do 411 end do 412 413 Bmat(:,:,:)=zero 414 415 ! calculate Bmat 416 ikptn = dtefield%ikpt_dk(ikpt,7,idir) 417 do iband=1, dtefield%mband_occ 418 do jband=1, dtefield%mband_occ 419 do kband=1, dtefield%mband_occ 420 Bmat(1,jband,kband) = Bmat(1,jband,kband) + Amat(1,iband,kband)*& 421 & qmat(1,jband,iband,ikptn,2,idir)& 422 & - Amat(2,iband,kband)*qmat(2,jband,iband,ikptn,2,idir) 423 Bmat(2,jband,kband) = Bmat(2,jband,kband) + Amat(1,iband,kband)*& 424 & qmat(2,jband,iband,ikptn,2,idir)& 425 + Amat(2,iband,kband)*qmat(1,jband,iband,ikptn,2,idir) 426 end do 427 end do 428 end do 429 430 ! calc. the second term of gradient------------------------------ 431 432 ! preparation 433 434 ikptnp1 = dtefield%ikpt_dk(ikpt,4,idir) 435 icg = dtefield%cgindex(ikptnp1,isppol) 436 npw_k1 = npwar1(ikpt) 437 npw_k2 = npwarr(ikptnp1) 438 pwind_tmp(1:npw_k1) =pwindall((ikpt-1)*mpw_tmp+1:(ikpt-1)*mpw_tmp+npw_k1,6,idir) 439 440 z1(:) = zero 441 z2(:) = zero 442 do ipw = 1, npw_k1 443 444 jpw = pwind_tmp(ipw) 445 446 if (jpw > 0) then 447 448 do iband = 1, dtefield%mband_occ 449 wfr = cg(1,icg + (iband - 1)*npw_k2*nspinor + jpw) 450 wfi = cg(2,icg + (iband - 1)*npw_k2*nspinor + jpw) 451 452 do jband=1, dtefield%mband_occ 453 454 grad_berry(1,ipw,jband) = grad_berry(1,ipw,jband) + fac*(Bmat(1,iband,jband)*wfr - Bmat(2,iband,jband)*wfi) 455 grad_berry(2,ipw,jband) = grad_berry(2,ipw,jband) + fac*(Bmat(1,iband,jband)*wfi + Bmat(2,iband,jband)*wfr) 456 457 end do 458 end do 459 end if 460 end do 461 462 end do !idir 463 464 ABI_DEALLOCATE(vect1) 465 ABI_DEALLOCATE(vect2) 466 ABI_DEALLOCATE(s1mat) 467 ABI_DEALLOCATE(Amat) 468 ABI_DEALLOCATE(Bmat) 469 ABI_DEALLOCATE(pwind_tmp) 470 471 end subroutine dfptff_gradberry