TABLE OF CONTENTS


ABINIT/dfptff_bec [ Functions ]

[ Top ] [ Functions ]

NAME

 dfptff_bec

FUNCTION

 calculate Born effective charge tensor in Eq.(33) in PRB 75, 115116(2007) [[cite:Wang2007]].

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
 idirpert = the current coloumn of the dielectric permittivity tensor
 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>
 rprimd(3,3)=dimensional primitive translations in real space (bohr)

OUTPUT

 d2lo(1,1:3,natom+5,1:3,1:natom) = Born effective charge tensor

PARENTS

      dfpt_scfcv

CHILDREN

      overlap_g

SOURCE

2469 subroutine dfptff_bec(cg,cg1,dtefield,natom,d2lo,idirpert,ipert,mband,mkmem,&
2470 &               mpw,mpw1,mpert,nkpt,npwarr,npwar1,nsppol,nspinor,pwindall,qmat,rprimd)
2471 
2472 
2473 !This section has been created automatically by the script Abilint (TD).
2474 !Do not modify the following lines by hand.
2475 #undef ABI_FUNC
2476 #define ABI_FUNC 'dfptff_bec'
2477 !End of the abilint section
2478 
2479  implicit none
2480 
2481 !Arguments ----------------------------------------
2482 !scalars
2483  integer,intent(in) :: idirpert,ipert,mband,mkmem,mpert,mpw,mpw1,natom,nkpt
2484  integer,intent(in) :: nspinor,nsppol
2485  type(efield_type),intent(in) :: dtefield
2486 !arrays
2487  integer,intent(in) :: npwar1(nkpt),npwarr(nkpt)
2488  integer,intent(in) :: pwindall(max(mpw,mpw1)*mkmem,8,3)
2489  real(dp),intent(in) :: cg(2,mpw*mband*mkmem*nspinor*nsppol)
2490  real(dp),intent(in) :: cg1(2,mpw1*mband*mkmem*nspinor*nsppol)
2491  real(dp),intent(in) :: qmat(2,dtefield%mband_occ,dtefield%mband_occ,nkpt,2,3)
2492  real(dp),intent(in) :: rprimd(3,3)
2493  real(dp),intent(inout) :: d2lo(2,3,mpert,3,mpert) !vz_i
2494 
2495 !Local variables ----------------------------------
2496 !scalars
2497  integer :: ialpha,iband,icg,icg1,idir,ikpt,ikpt1,jband,mpw_tmp,npw_k1
2498  integer :: npw_k2,pwmax,pwmin
2499  real(dp) :: doti,dotr,e0,fac
2500 !arrays
2501  integer,allocatable :: pwind_tmp(:)
2502  real(dp) :: edir(3)
2503  real(dp),allocatable :: s1mat(:,:,:),vect1(:,:),vect2(:,:)
2504 
2505 ! *************************************************************************
2506 
2507 !calculate s1 matrices -----------------------------
2508  mpw_tmp=max(mpw,mpw1)
2509  ABI_ALLOCATE(s1mat,(2,dtefield%mband_occ,dtefield%mband_occ))
2510  ABI_ALLOCATE(vect1,(2,0:mpw_tmp))
2511  ABI_ALLOCATE(vect2,(2,0:mpw_tmp))
2512  ABI_ALLOCATE(pwind_tmp,(mpw_tmp))
2513  vect1(:,0) = zero ; vect2(:,0) = zero
2514 
2515  edir(:)=zero
2516 
2517  do ikpt=1,nkpt
2518    do idir=1,3
2519 
2520 !    compute <u^(0)_{k_j}|u^(1)_{k_j+1,q}> matrix--- q=0 ----------------------------------------
2521 
2522 !    prepare to calculate overlap matrix
2523      ikpt1 = dtefield%ikpt_dk(ikpt,1,idir)
2524      icg = dtefield%cgindex(ikpt,1)
2525      icg1 = dtefield%cgindex(ikpt1,1+nsppol)
2526      npw_k1 = npwarr(ikpt)
2527      npw_k2 = npwar1(ikpt1)
2528      pwind_tmp(1:npw_k1) = pwindall((ikpt-1)*mpw_tmp+1:(ikpt-1)*mpw_tmp+npw_k1,7,idir)
2529      vect1(:,0) = zero ; vect2(:,0) = zero
2530      do jband = 1, dtefield%mband_occ
2531        vect2(:,1:npw_k2) = &
2532 &       cg1(:,icg1 + 1 + (jband-1)*npw_k2*nspinor:icg1 + jband*npw_k2*nspinor)
2533        if (npw_k2 < mpw_tmp) vect2(:,npw_k2+1:mpw_tmp) = zero
2534        do iband = 1, dtefield%mband_occ
2535          pwmin = (iband-1)*npw_k1*nspinor
2536          pwmax = pwmin + npw_k1*nspinor
2537          vect1(:,1:npw_k1) = &
2538 &         cg(:,icg + 1 + pwmin:icg + pwmax)
2539          if (npw_k1 < mpw_tmp) vect1(:,npw_k1+1:mpw_tmp) = zero
2540          call overlap_g(doti,dotr,mpw_tmp,npw_k1,npw_k2,nspinor,pwind_tmp,&
2541 &         vect1,vect2)
2542          s1mat(1,iband,jband) = dotr
2543          s1mat(2,iband,jband) = doti
2544        end do    ! iband
2545      end do    !jband
2546 
2547 !    compute <u^(1)_{k_j,q}|u^(0)_{k_j+1}> matrix-- q=0 -------------------------------------
2548 
2549 !    prepare to calculate overlap matrix
2550      ikpt1 = dtefield%ikpt_dk(ikpt,1,idir)
2551      icg = dtefield%cgindex(ikpt,1+nsppol)
2552      icg1 = dtefield%cgindex(ikpt1,1)
2553      npw_k1 = npwar1(ikpt)
2554      npw_k2 = npwarr(ikpt1)
2555      pwind_tmp(1:npw_k1) = pwindall((ikpt-1)*mpw_tmp+1:(ikpt-1)*mpw_tmp+npw_k1,5,idir)
2556      vect1(:,0) = zero ; vect2(:,0) = zero
2557      do jband = 1, dtefield%mband_occ
2558        vect2(:,1:npw_k2) = &
2559 &       cg(:,icg1 + 1 + (jband-1)*npw_k2*nspinor:icg1 + jband*npw_k2*nspinor)
2560        if (npw_k2 < mpw_tmp) vect2(:,npw_k2+1:mpw_tmp) = zero
2561        do iband = 1, dtefield%mband_occ
2562          pwmin = (iband-1)*npw_k1*nspinor
2563          pwmax = pwmin + npw_k1*nspinor
2564          vect1(:,1:npw_k1) = &
2565 &         cg1(:,icg + 1 + pwmin:icg + pwmax)
2566          if (npw_k1 < mpw_tmp) vect1(:,npw_k1+1:mpw_tmp) = zero
2567          call overlap_g(doti,dotr,mpw_tmp,npw_k1,npw_k2,nspinor,pwind_tmp,&
2568 &         vect1,vect2)
2569          s1mat(1,iband,jband) = s1mat(1,iband,jband) + dotr
2570          s1mat(2,iband,jband) = s1mat(2,iband,jband) + doti
2571        end do    ! iband
2572      end do    !jband
2573 
2574 !    sum over the whole------------------------------------------------------------
2575 
2576      e0=zero
2577 
2578      do iband=1,dtefield%mband_occ
2579        do jband=1,dtefield%mband_occ
2580          e0 = e0 + (s1mat(1,iband,jband)*qmat(2,jband,iband,ikpt,1,idir)&
2581 &         +    s1mat(2,iband,jband)*qmat(1,jband,iband,ikpt,1,idir))
2582        end do
2583      end do
2584 
2585      do ialpha=1,3
2586        fac = rprimd(ialpha,idir)/&
2587 &       (dble(dtefield%nstr(idir))*pi)
2588 
2589        edir(ialpha)=edir(ialpha)+ e0*fac
2590      end do
2591 
2592    end do
2593  end do
2594 
2595  d2lo(1,1:3,natom+2,idirpert,ipert)=edir(:)
2596 
2597  ABI_DEALLOCATE(s1mat)
2598  ABI_DEALLOCATE(vect1)
2599  ABI_DEALLOCATE(vect2)
2600  ABI_DEALLOCATE(pwind_tmp)
2601 
2602 end subroutine dfptff_bec

ABINIT/dfptff_die [ Functions ]

[ Top ] [ Functions ]

NAME

 dfptff_die

FUNCTION

 calculate electric susceptibility tensor in Eq.(28) in PRB 75, 115116(2007) [[cite:Wang2007]].

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
 idirpert = the current coloumn of the dielectric permittivity tensor
 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>
 rprimd(3,3)=dimensional primitive translations in real space (bohr)

OUTPUT

 diet = electric susceptibility tensor

PARENTS

      dfpt_scfcv

CHILDREN

      overlap_g

SOURCE

2287 subroutine dfptff_die(cg,cg1,dtefield,d2lo,idirpert,ipert,mband,mkmem,&
2288 &               mpw,mpw1,mpert,nkpt,npwarr,npwar1,nsppol,nspinor,pwindall,qmat,rprimd)
2289 
2290 
2291 !This section has been created automatically by the script Abilint (TD).
2292 !Do not modify the following lines by hand.
2293 #undef ABI_FUNC
2294 #define ABI_FUNC 'dfptff_die'
2295 !End of the abilint section
2296 
2297  implicit none
2298 
2299 !Arguments ----------------------------------------
2300 !scalars
2301  integer,intent(in) :: idirpert,ipert,mband,mkmem,mpert,mpw,mpw1,nkpt,nspinor
2302  integer,intent(in) :: nsppol
2303  type(efield_type),intent(in) :: dtefield
2304 !arrays
2305  integer,intent(in) :: npwar1(nkpt),npwarr(nkpt)
2306  integer,intent(in) :: pwindall(max(mpw,mpw1)*mkmem,8,3)
2307  real(dp),intent(in) :: cg(2,mpw*mband*mkmem*nspinor*nsppol)
2308  real(dp),intent(in) :: cg1(2,mpw1*mband*mkmem*nspinor*nsppol)
2309  real(dp),intent(in) :: qmat(2,dtefield%mband_occ,dtefield%mband_occ,nkpt,2,3)
2310  real(dp),intent(in) :: rprimd(3,3)
2311  real(dp),intent(inout) :: d2lo(2,3,mpert,3,mpert) !vz_i
2312 
2313 !Local variables ----------------------------------
2314 !scalars
2315  integer :: ialpha,iband,icg,icg1,idir,ikpt,ikpt1,jband,mpw_tmp,npw_k1
2316  integer :: npw_k2,pwmax,pwmin
2317  real(dp) :: doti,dotr,e0,fac
2318 !arrays
2319  integer,allocatable :: pwind_tmp(:)
2320  real(dp) :: edir(3)
2321  real(dp),allocatable :: s1mat(:,:,:),vect1(:,:),vect2(:,:)
2322 
2323 ! *************************************************************************
2324 
2325 !calculate s1 matrices -----------------------------
2326  mpw_tmp=max(mpw,mpw1)
2327  ABI_ALLOCATE(s1mat,(2,dtefield%mband_occ,dtefield%mband_occ))
2328  ABI_ALLOCATE(vect1,(2,0:mpw_tmp))
2329  ABI_ALLOCATE(vect2,(2,0:mpw_tmp))
2330  ABI_ALLOCATE(pwind_tmp,(mpw_tmp))
2331  vect1(:,0) = zero ; vect2(:,0) = zero
2332 
2333  edir(:)=zero
2334 
2335  do ikpt=1,nkpt
2336    do idir=1,3
2337 !    compute <u^(0)_{k_j}|u^(1)_{k_j+1,q}> matrix--- q=0 ----------------------------------------
2338 
2339 !    prepare to calculate overlap matrix
2340      ikpt1 = dtefield%ikpt_dk(ikpt,1,idir)
2341      icg = dtefield%cgindex(ikpt,1)
2342      icg1 = dtefield%cgindex(ikpt1,1+nsppol)
2343      npw_k1 = npwarr(ikpt)
2344      npw_k2 = npwar1(ikpt1)
2345      pwind_tmp(1:npw_k1) = pwindall((ikpt-1)*mpw_tmp+1:(ikpt-1)*mpw_tmp+npw_k1,7,idir)
2346 
2347      vect1(:,0) = zero ; vect2(:,0) = zero
2348      do jband = 1, dtefield%mband_occ
2349        vect2(:,1:npw_k2) = &
2350 &       cg1(:,icg1 + 1 + (jband-1)*npw_k2*nspinor:icg1 + jband*npw_k2*nspinor)
2351        if (npw_k2 < mpw_tmp) vect2(:,npw_k2+1:mpw_tmp) = zero
2352        do iband = 1, dtefield%mband_occ
2353          pwmin = (iband-1)*npw_k1*nspinor
2354          pwmax = pwmin + npw_k1*nspinor
2355          vect1(:,1:npw_k1) = &
2356 &         cg(:,icg + 1 + pwmin:icg + pwmax)
2357          if (npw_k1 < mpw_tmp) vect1(:,npw_k1+1:mpw_tmp) = zero
2358          call overlap_g(doti,dotr,mpw_tmp,npw_k1,npw_k2,nspinor,pwind_tmp,&
2359 &         vect1,vect2)
2360          s1mat(1,iband,jband) = dotr
2361          s1mat(2,iband,jband) = doti
2362        end do    ! iband
2363      end do    !jband
2364 
2365 !    compute <u^(1)_{k_j,q}|u^(0)_{k_j+1}> matrix-- q=0 -------------------------------------
2366 
2367 !    prepare to calculate overlap matrix
2368      ikpt1 = dtefield%ikpt_dk(ikpt,1,idir)
2369      icg = dtefield%cgindex(ikpt,1+nsppol)
2370      icg1 = dtefield%cgindex(ikpt1,1)
2371      npw_k1 = npwar1(ikpt)
2372      npw_k2 = npwarr(ikpt1)
2373      pwind_tmp(1:npw_k1) = pwindall((ikpt-1)*mpw_tmp+1:(ikpt-1)*mpw_tmp+npw_k1,5,idir)
2374      vect1(:,0) = zero ; vect2(:,0) = zero
2375      do jband = 1, dtefield%mband_occ
2376        vect2(:,1:npw_k2) = &
2377 &       cg(:,icg1 + 1 + (jband-1)*npw_k2*nspinor:icg1 + jband*npw_k2*nspinor)
2378        if (npw_k2 < mpw_tmp) vect2(:,npw_k2+1:mpw_tmp) = zero
2379        do iband = 1, dtefield%mband_occ
2380          pwmin = (iband-1)*npw_k1*nspinor
2381          pwmax = pwmin + npw_k1*nspinor
2382          vect1(:,1:npw_k1) = &
2383 &         cg1(:,icg + 1 + pwmin:icg + pwmax)
2384          if (npw_k1 < mpw_tmp) vect1(:,npw_k1+1:mpw_tmp) = zero
2385          call overlap_g(doti,dotr,mpw_tmp,npw_k1,npw_k2,nspinor,pwind_tmp,&
2386 &         vect1,vect2)
2387          s1mat(1,iband,jband) = s1mat(1,iband,jband) + dotr
2388          s1mat(2,iband,jband) = s1mat(2,iband,jband) + doti
2389        end do    ! iband
2390      end do    !jband
2391 
2392 !    sum over the whole------------------------------------------------------------
2393 
2394      e0=zero
2395 
2396      do iband=1,dtefield%mband_occ
2397        do jband=1,dtefield%mband_occ
2398          e0 = e0 + (s1mat(1,iband,jband)*qmat(2,jband,iband,ikpt,1,idir)&
2399 &         +    s1mat(2,iband,jband)*qmat(1,jband,iband,ikpt,1,idir))
2400 
2401        end do
2402      end do
2403 
2404      do ialpha=1,3
2405        fac = rprimd(ialpha,idir)/&
2406 &       (dble(dtefield%nstr(idir))*pi)
2407        edir(ialpha)=edir(ialpha)+ e0*fac
2408      end do
2409 
2410    end do !idir
2411  end do !ikpt
2412 
2413  d2lo(1,1:3,ipert,idirpert,ipert)=edir(:)
2414 
2415  ABI_DEALLOCATE(s1mat)
2416  ABI_DEALLOCATE(vect1)
2417  ABI_DEALLOCATE(vect2)
2418  ABI_DEALLOCATE(pwind_tmp)
2419 
2420 end subroutine dfptff_die

ABINIT/dfptff_ebp [ Functions ]

[ Top ] [ Functions ]

NAME

 dfptff_ebp

FUNCTION

 calculation of the energy from the term \Omega E \cdot P

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

1994 subroutine dfptff_ebp(cg,cg1,dtefield,eberry,mband,mkmem,&
1995 &               mpw,mpw1,nkpt,npwarr,npwar1,nsppol,nspinor,pwindall,qmat)
1996 
1997 
1998 !This section has been created automatically by the script Abilint (TD).
1999 !Do not modify the following lines by hand.
2000 #undef ABI_FUNC
2001 #define ABI_FUNC 'dfptff_ebp'
2002 !End of the abilint section
2003 
2004  implicit none
2005 
2006 !Arguments ----------------------------------------
2007 !scalars
2008  integer,intent(in) :: mband,mkmem,mpw,mpw1,nkpt,nspinor,nsppol
2009  real(dp),intent(out) :: eberry
2010  type(efield_type),intent(in) :: dtefield
2011 !arrays
2012  integer,intent(in) :: npwar1(nkpt),npwarr(nkpt)
2013  integer,intent(in) :: pwindall(max(mpw,mpw1)*mkmem,8,3)
2014  real(dp),intent(in) :: cg(2,mpw*mband*mkmem*nspinor*nsppol)
2015  real(dp),intent(in) :: cg1(2,mpw1*mband*mkmem*nspinor*nsppol)
2016  real(dp),intent(in) :: qmat(2,dtefield%mband_occ,dtefield%mband_occ,nkpt,2,3)
2017 
2018 !Local variables ----------------------------------
2019 !scalars
2020  integer :: iband,icg,icg1,idir
2021  integer :: ikpt,ikpt1,ikptn,ikptnm
2022  integer :: jband,kband,mpw_tmp,npw_k1,npw_k2,pwmax,pwmin
2023  real(dp) :: doti,dotr,e0,fac
2024 !arrays
2025  integer,allocatable :: pwind_tmp(:)
2026  real(dp) :: z1(2)
2027  real(dp),allocatable :: Amat(:,:,:),umat(:,:,:,:),vect1(:,:),vect2(:,:)
2028 
2029 ! *************************************************************************
2030 
2031 !calculate 4 matrices -----------------------------
2032  mpw_tmp=max(mpw,mpw1)
2033  ABI_ALLOCATE(umat,(2,dtefield%mband_occ,dtefield%mband_occ,4))
2034  ABI_ALLOCATE(vect1,(2,0:mpw_tmp))
2035  ABI_ALLOCATE(vect2,(2,0:mpw_tmp))
2036  ABI_ALLOCATE(pwind_tmp,(mpw_tmp))
2037  ABI_ALLOCATE(Amat,(2,dtefield%mband_occ,dtefield%mband_occ))
2038  vect1(:,0) = zero ; vect2(:,0) = zero
2039  eberry=zero
2040 
2041  do ikpt=1,nkpt
2042 
2043    do idir=1,3
2044 
2045      fac = dtefield%efield_dot(idir)/&
2046 &     (dble(dtefield%nstr(idir))*four_pi)
2047 
2048 !    compute <u^(1)_{k_j,q}|u^(1)_{k_j+1,q}> matrix---------
2049 
2050 !    prepare to calculate overlap matrix
2051      ikpt1 = dtefield%ikpt_dk(ikpt,1,idir)
2052      icg = dtefield%cgindex(ikpt,1+nsppol)
2053      icg1 = dtefield%cgindex(ikpt1,1+nsppol)
2054      npw_k1 = npwar1(ikpt)
2055      npw_k2 = npwar1(ikpt1)
2056      pwind_tmp(1:npw_k1) = pwindall((ikpt-1)*mpw_tmp+1:(ikpt-1)*mpw_tmp+npw_k1,3,idir)
2057      vect1(:,0) = zero ; vect2(:,0) = zero
2058      do jband = 1, dtefield%mband_occ
2059        vect2(:,1:npw_k2) = &
2060 &       cg1(:,icg1 + 1 + (jband-1)*npw_k2*nspinor:icg1 + jband*npw_k2*nspinor)
2061 
2062        if (npw_k2 < mpw_tmp) vect2(:,npw_k2+1:mpw_tmp) = zero
2063 
2064        do iband = 1, dtefield%mband_occ
2065 
2066          pwmin = (iband-1)*npw_k1*nspinor
2067          pwmax = pwmin + npw_k1*nspinor
2068          vect1(:,1:npw_k1) = &
2069 &         cg1(:,icg + 1 + pwmin:icg + pwmax)
2070          if (npw_k1 < mpw_tmp) vect1(:,npw_k1+1:mpw_tmp) = zero
2071          call overlap_g(doti,dotr,mpw_tmp,npw_k1,npw_k2,nspinor,pwind_tmp,&
2072 &         vect1,vect2)
2073          umat(1,iband,jband,1) = dotr
2074          umat(2,iband,jband,1) = doti
2075 
2076        end do    ! iband
2077      end do    !jband
2078 
2079 !    compute <u^(0)_{k_j}|u^(1)_{k_j-n+1,q}> matrix----------------------------------------------------
2080 
2081 !    prepare to calculate overlap matrix
2082      ikpt1 = dtefield%ikpt_dk(ikpt,5,idir)
2083      icg = dtefield%cgindex(ikpt,1)
2084      icg1 = dtefield%cgindex(ikpt1,1+nsppol)
2085      npw_k1 = npwarr(ikpt)
2086      npw_k2 = npwar1(ikpt1)
2087      pwind_tmp(1:npw_k1) = pwindall((ikpt-1)*mpw_tmp+1:(ikpt-1)*mpw_tmp+npw_k1,7,idir)
2088 
2089      vect1(:,0) = zero ; vect2(:,0) = zero
2090      do jband = 1, dtefield%mband_occ
2091        vect2(:,1:npw_k2) = &
2092 &       cg1(:,icg1 + 1 + (jband-1)*npw_k2*nspinor:icg1 + jband*npw_k2*nspinor)
2093        if (npw_k2 < mpw_tmp) vect2(:,npw_k2+1:mpw_tmp) = zero
2094 
2095        do iband = 1, dtefield%mband_occ
2096 
2097          pwmin = (iband-1)*npw_k1*nspinor
2098          pwmax = pwmin + npw_k1*nspinor
2099          vect1(:,1:npw_k1) = &
2100 &         cg(:,icg + 1 + pwmin:icg + pwmax)
2101          if (npw_k1 < mpw_tmp) vect1(:,npw_k1+1:mpw_tmp) = zero
2102          call overlap_g(doti,dotr,mpw_tmp,npw_k1,npw_k2,nspinor,pwind_tmp,&
2103 &         vect1,vect2)
2104          umat(1,iband,jband,2) = dotr
2105          umat(2,iband,jband,2) = doti
2106 
2107        end do    ! iband
2108      end do    !jband
2109 
2110 !    compute <u^(1)_{k_j-n,q}|u^(0)_{k_j+1}> matrix----------------------------------------------------
2111 
2112 !    prepare to calculate overlap matrix
2113      ikptn  = dtefield%ikpt_dk(ikpt,8,idir)
2114      ikpt1 = dtefield%ikpt_dk(ikpt,1,idir)
2115      icg = dtefield%cgindex(ikptn,1+nsppol)
2116      icg1 = dtefield%cgindex(ikpt1,1)
2117      npw_k1 = npwar1(ikptn)
2118      npw_k2 = npwarr(ikpt1)
2119      pwind_tmp(1:npw_k1) = pwindall((ikptn-1)*mpw_tmp+1:(ikptn-1)*mpw_tmp+npw_k1,5,idir)
2120 
2121      vect1(:,0) = zero ; vect2(:,0) = zero
2122      do jband = 1, dtefield%mband_occ
2123        vect2(:,1:npw_k2) = &
2124 &       cg(:,icg1 + 1 + (jband-1)*npw_k2*nspinor:icg1 + jband*npw_k2*nspinor)
2125        if (npw_k2 < mpw_tmp) vect2(:,npw_k2+1:mpw_tmp) = zero
2126 
2127        do iband = 1, dtefield%mband_occ
2128 
2129          pwmin = (iband-1)*npw_k1*nspinor
2130          pwmax = pwmin + npw_k1*nspinor
2131          vect1(:,1:npw_k1) = &
2132 &         cg1(:,icg + 1 + pwmin:icg + pwmax)
2133          if (npw_k1 < mpw_tmp) vect1(:,npw_k1+1:mpw_tmp) = zero
2134          call overlap_g(doti,dotr,mpw_tmp,npw_k1,npw_k2,nspinor,pwind_tmp,&
2135 &         vect1,vect2)
2136          umat(1,iband,jband,3) = dotr
2137          umat(2,iband,jband,3) = doti
2138 
2139        end do    ! iband
2140      end do    !jband
2141 
2142 !    compute <u^(0)_{-k_j-n+1}|u^(1)_{-k_j,q}> matrix----------------------------------------------------
2143 
2144 !    prepare to calculate overlap matrix
2145      ikptn = dtefield%ikpt_dk(ikpt,5,idir)
2146      ikptnm = dtefield%ikpt_dk(ikptn,9,idir)
2147      ikpt1 = dtefield%ikpt_dk(ikpt,9,idir)
2148      icg = dtefield%cgindex(ikptnm,1)
2149      icg1 = dtefield%cgindex(ikpt1,1+nsppol)
2150      npw_k1 = npwarr(ikptnm)
2151      npw_k2 = npwar1(ikpt1)
2152      pwind_tmp(1:npw_k1) = pwindall((ikptnm-1)*mpw_tmp+1:(ikptnm-1)*mpw_tmp+npw_k1,7,idir)
2153 
2154      vect1(:,0) = zero ; vect2(:,0) = zero
2155      do jband = 1, dtefield%mband_occ
2156        vect2(:,1:npw_k2) = &
2157 &       cg1(:,icg1 + 1 + (jband-1)*npw_k2*nspinor:icg1 + jband*npw_k2*nspinor)
2158        if (npw_k2 < mpw_tmp) vect2(:,npw_k2+1:mpw_tmp) = zero
2159 
2160        do iband = 1, dtefield%mband_occ
2161 
2162          pwmin = (iband-1)*npw_k1*nspinor
2163          pwmax = pwmin + npw_k1*nspinor
2164          vect1(:,1:npw_k1) = &
2165 &         cg(:,icg + 1 + pwmin:icg + pwmax)
2166          if (npw_k1 < mpw_tmp) vect1(:,npw_k1+1:mpw_tmp) = zero
2167          call overlap_g(doti,dotr,mpw_tmp,npw_k1,npw_k2,nspinor,pwind_tmp,&
2168 &         vect1,vect2)
2169          umat(1,iband,jband,4) = dotr
2170          umat(2,iband,jband,4) = doti
2171 
2172        end do    ! iband
2173      end do    !jband
2174 
2175 !    sum over the whole------------------------------------------------------------
2176 
2177      e0=zero
2178      do iband=1,dtefield%mband_occ
2179        do jband=1,dtefield%mband_occ
2180          e0 = e0 + 4_dp*(umat(1,iband,jband,1)*qmat(2,jband,iband,ikpt,1,idir)&
2181 &         +       umat(2,iband,jband,1)*qmat(1,jband,iband,ikpt,1,idir))
2182 
2183        end do
2184      end do
2185 
2186      eberry = eberry - e0*fac
2187 
2188      e0=zero
2189 
2190      ikptn=dtefield%ikpt_dk(ikpt,8,idir)
2191 
2192      Amat(:,:,:)=zero
2193 
2194 !    calculate Amat
2195      do iband=1, dtefield%mband_occ
2196        do jband=1, dtefield%mband_occ
2197          do kband=1, dtefield%mband_occ
2198            Amat(1,iband,jband) = Amat(1,iband,jband) + (umat(1,iband,kband,3))*&
2199 &           qmat(1,kband,jband,ikpt,1,idir)&
2200 &           - (umat(2,iband,kband,3))*qmat(2,kband,jband,ikpt,1,idir)
2201            Amat(2,iband,jband) = Amat(2,iband,jband) + (umat(1,iband,kband,3))*&
2202 &           qmat(2,kband,jband,ikpt,1,idir)&
2203 &           + (umat(2,iband,kband,3))*qmat(1,kband,jband,ikpt,1,idir)
2204          end do
2205        end do
2206      end do
2207 
2208      do iband=1, dtefield%mband_occ
2209        do jband=1, dtefield%mband_occ
2210          do kband=1, dtefield%mband_occ
2211 
2212            z1(1) = (umat(1,jband,iband,4)+umat(1,iband,jband,2))*&
2213 &           qmat(1,jband,kband,ikptn,1,idir)&
2214 &           -    (umat(2,jband,iband,4)+umat(2,iband,jband,2))*&
2215 &           qmat(2,jband,kband,ikptn,1,idir)
2216            z1(2) = (umat(1,jband,iband,4)+umat(1,iband,jband,2))*&
2217 &           qmat(2,jband,kband,ikptn,1,idir)&
2218 &           +    (umat(2,jband,iband,4)+umat(2,iband,jband,2))*&
2219 &           qmat(1,jband,kband,ikptn,1,idir)
2220 
2221            e0 = e0 - 4_dp*(z1(1)*Amat(2,kband,iband)+z1(2)*Amat(1,kband,iband))
2222 
2223          end do
2224        end do
2225      end do
2226 
2227      eberry = eberry - e0*fac
2228 
2229    end do !end idir
2230  end do !end ikpt
2231 
2232  ABI_DEALLOCATE(umat)
2233  ABI_DEALLOCATE(vect1)
2234  ABI_DEALLOCATE(vect2)
2235  ABI_DEALLOCATE(pwind_tmp)
2236  ABI_DEALLOCATE(Amat)
2237 
2238 end subroutine dfptff_ebp

ABINIT/dfptff_edie [ Functions ]

[ Top ] [ 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) [[cite:Wang2007]].

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

1658 subroutine dfptff_edie(cg,cg1,dtefield,eberry,idir_efield,mband,mkmem,&
1659 &                mpw,mpw1,nkpt,npwarr,npwar1,nsppol,nspinor,pwindall,qmat,rprimd)
1660 
1661 
1662 !This section has been created automatically by the script Abilint (TD).
1663 !Do not modify the following lines by hand.
1664 #undef ABI_FUNC
1665 #define ABI_FUNC 'dfptff_edie'
1666 !End of the abilint section
1667 
1668  implicit none
1669 
1670 !Arguments ----------------------------------------
1671 !scalars
1672  integer,intent(in) :: idir_efield,mband,mkmem,mpw,mpw1,nkpt,nspinor,nsppol
1673  real(dp),intent(out) :: eberry
1674  type(efield_type),intent(in) :: dtefield
1675 !arrays
1676  integer,intent(in) :: npwar1(nkpt),npwarr(nkpt)
1677  integer,intent(in) :: pwindall(max(mpw,mpw1)*mkmem,8,3)
1678  real(dp),intent(in) :: cg(2,mpw*mband*mkmem*nspinor*nsppol)
1679  real(dp),intent(in) :: cg1(2,mpw1*mband*mkmem*nspinor*nsppol)
1680  real(dp),intent(in) :: qmat(2,dtefield%mband_occ,dtefield%mband_occ,nkpt,2,3)
1681  real(dp),intent(in) :: rprimd(3,3)
1682 
1683 !Local variables ----------------------------------
1684 !scalars
1685  integer :: iband,icg,icg1,idir
1686  integer :: ikpt,ikpt1,ikptn,ikptnm
1687  integer :: jband,kband,mpw_tmp,npw_k1,npw_k2,pwmax,pwmin
1688  real(dp) :: doti,dotr,e0,fac
1689 !arrays
1690  integer,allocatable :: pwind_tmp(:)
1691  real(dp) :: z1(2)
1692  real(dp),allocatable :: Amat(:,:,:),umat(:,:,:,:),vect1(:,:),vect2(:,:)
1693 
1694 ! *************************************************************************
1695 
1696 !calculate 4 matrices -----------------------------
1697  mpw_tmp=max(mpw,mpw1)
1698  ABI_ALLOCATE(umat,(2,dtefield%mband_occ,dtefield%mband_occ,4))
1699  ABI_ALLOCATE(vect1,(2,0:mpw_tmp))
1700  ABI_ALLOCATE(vect2,(2,0:mpw_tmp))
1701  ABI_ALLOCATE(pwind_tmp,(mpw_tmp))
1702  ABI_ALLOCATE(Amat,(2,dtefield%mband_occ,dtefield%mband_occ))
1703  vect1(:,0) = zero ; vect2(:,0) = zero
1704  eberry=zero
1705 
1706  do ikpt=1,nkpt
1707    do idir=1,3
1708      fac = dtefield%efield_dot(idir)/&
1709 &     (dble(dtefield%nstr(idir))*four_pi)
1710 
1711 !    compute <u^(1)_{k_j,q}|u^(1)_{k_j+1,q}> matrix----------------------------------------------------
1712 
1713 !    prepare to calculate overlap matrix
1714      ikpt1 = dtefield%ikpt_dk(ikpt,1,idir)
1715      icg = dtefield%cgindex(ikpt,1+nsppol)
1716      icg1 = dtefield%cgindex(ikpt1,1+nsppol)
1717      npw_k1 = npwar1(ikpt)
1718      npw_k2 = npwar1(ikpt1)
1719      pwind_tmp(1:npw_k1) = pwindall((ikpt-1)*mpw_tmp+1:(ikpt-1)*mpw_tmp+npw_k1,3,idir)
1720      vect1(:,0) = zero ; vect2(:,0) = zero
1721      do jband = 1, dtefield%mband_occ
1722        vect2(:,1:npw_k2) = &
1723 &       cg1(:,icg1 + 1 + (jband-1)*npw_k2*nspinor:icg1 + jband*npw_k2*nspinor)
1724        if (npw_k2 < mpw_tmp) vect2(:,npw_k2+1:mpw_tmp) = zero
1725        do iband = 1, dtefield%mband_occ
1726          pwmin = (iband-1)*npw_k1*nspinor
1727          pwmax = pwmin + npw_k1*nspinor
1728          vect1(:,1:npw_k1) = &
1729 &         cg1(:,icg + 1 + pwmin:icg + pwmax)
1730          if (npw_k1 < mpw_tmp) vect1(:,npw_k1+1:mpw_tmp) = zero
1731          call overlap_g(doti,dotr,mpw_tmp,npw_k1,npw_k2,nspinor,pwind_tmp,&
1732 &         vect1,vect2)
1733          umat(1,iband,jband,1) = dotr
1734          umat(2,iband,jband,1) = doti
1735        end do    ! iband
1736      end do    !jband
1737 
1738 !    compute <u^(0)_{k_j}|u^(1)_{k_j-n+1,q}> matrix----------------------------------------------------
1739 
1740 !    prepare to calculate overlap matrix
1741      ikpt1 = dtefield%ikpt_dk(ikpt,5,idir)
1742      icg = dtefield%cgindex(ikpt,1)
1743      icg1 = dtefield%cgindex(ikpt1,1+nsppol)
1744      npw_k1 = npwarr(ikpt)
1745      npw_k2 = npwar1(ikpt1)
1746      pwind_tmp(1:npw_k1) = pwindall((ikpt-1)*mpw_tmp+1:(ikpt-1)*mpw_tmp+npw_k1,7,idir)
1747 
1748      vect1(:,0) = zero ; vect2(:,0) = zero
1749      do jband = 1, dtefield%mband_occ
1750        vect2(:,1:npw_k2) = &
1751 &       cg1(:,icg1 + 1 + (jband-1)*npw_k2*nspinor:icg1 + jband*npw_k2*nspinor)
1752        if (npw_k2 < mpw_tmp) vect2(:,npw_k2+1:mpw_tmp) = zero
1753        do iband = 1, dtefield%mband_occ
1754          pwmin = (iband-1)*npw_k1*nspinor
1755          pwmax = pwmin + npw_k1*nspinor
1756          vect1(:,1:npw_k1) = &
1757 &         cg(:,icg + 1 + pwmin:icg + pwmax)
1758          if (npw_k1 < mpw_tmp) vect1(:,npw_k1+1:mpw_tmp) = zero
1759          call overlap_g(doti,dotr,mpw_tmp,npw_k1,npw_k2,nspinor,pwind_tmp,&
1760 &         vect1,vect2)
1761          umat(1,iband,jband,2) = dotr
1762          umat(2,iband,jband,2) = doti
1763        end do    ! iband
1764      end do    !jband
1765 
1766 !    compute <u^(1)_{k_j-n,q}|u^(0)_{k_j+1}> matrix----------------------------------------------------
1767 
1768 !    prepare to calculate overlap matrix
1769      ikptn  = dtefield%ikpt_dk(ikpt,8,idir)
1770      ikpt1 = dtefield%ikpt_dk(ikpt,1,idir)
1771      icg = dtefield%cgindex(ikptn,1+nsppol)
1772      icg1 = dtefield%cgindex(ikpt1,1)
1773      npw_k1 = npwar1(ikptn)
1774      npw_k2 = npwarr(ikpt1)
1775      pwind_tmp(1:npw_k1) = pwindall((ikptn-1)*mpw_tmp+1:(ikptn-1)*mpw_tmp+npw_k1,5,idir)
1776 
1777      vect1(:,0) = zero ; vect2(:,0) = zero
1778      do jband = 1, dtefield%mband_occ
1779        vect2(:,1:npw_k2) = &
1780 &       cg(:,icg1 + 1 + (jband-1)*npw_k2*nspinor:icg1 + jband*npw_k2*nspinor)
1781        if (npw_k2 < mpw_tmp) vect2(:,npw_k2+1:mpw_tmp) = zero
1782        do iband = 1, dtefield%mband_occ
1783          pwmin = (iband-1)*npw_k1*nspinor
1784          pwmax = pwmin + npw_k1*nspinor
1785          vect1(:,1:npw_k1) = &
1786 &         cg1(:,icg + 1 + pwmin:icg + pwmax)
1787          if (npw_k1 < mpw_tmp) vect1(:,npw_k1+1:mpw_tmp) = zero
1788          call overlap_g(doti,dotr,mpw_tmp,npw_k1,npw_k2,nspinor,pwind_tmp,&
1789 &         vect1,vect2)
1790          umat(1,iband,jband,3) = dotr
1791          umat(2,iband,jband,3) = doti
1792        end do    ! iband
1793      end do    !jband
1794 
1795 !    compute <u^(0)_{-k_j-n+1}|u^(1)_{-k_j,q}> matrix----------------------------------------------------
1796 
1797 !    prepare to calculate overlap matrix
1798      ikptn = dtefield%ikpt_dk(ikpt,5,idir)
1799      ikptnm = dtefield%ikpt_dk(ikptn,9,idir)
1800      ikpt1 = dtefield%ikpt_dk(ikpt,9,idir)
1801      icg = dtefield%cgindex(ikptnm,1)
1802      icg1 = dtefield%cgindex(ikpt1,1+nsppol)
1803      npw_k1 = npwarr(ikptnm)
1804      npw_k2 = npwar1(ikpt1)
1805      pwind_tmp(1:npw_k1) = pwindall((ikptnm-1)*mpw_tmp+1:(ikptnm-1)*mpw_tmp+npw_k1,7,idir)
1806      vect1(:,0) = zero ; vect2(:,0) = zero
1807      do jband = 1, dtefield%mband_occ
1808        vect2(:,1:npw_k2) = &
1809 &       cg1(:,icg1 + 1 + (jband-1)*npw_k2*nspinor:icg1 + jband*npw_k2*nspinor)
1810        if (npw_k2 < mpw_tmp) vect2(:,npw_k2+1:mpw_tmp) = zero
1811        do iband = 1, dtefield%mband_occ
1812          pwmin = (iband-1)*npw_k1*nspinor
1813          pwmax = pwmin + npw_k1*nspinor
1814          vect1(:,1:npw_k1) = &
1815 &         cg(:,icg + 1 + pwmin:icg + pwmax)
1816          if (npw_k1 < mpw_tmp) vect1(:,npw_k1+1:mpw_tmp) = zero
1817          call overlap_g(doti,dotr,mpw_tmp,npw_k1,npw_k2,nspinor,pwind_tmp,&
1818 &         vect1,vect2)
1819          umat(1,iband,jband,4) = dotr
1820          umat(2,iband,jband,4) = doti
1821        end do    ! iband
1822      end do    !jband
1823 
1824 !    sum over the whole------------------------------------------------------------
1825 
1826      e0=zero
1827      do iband=1,dtefield%mband_occ
1828        do jband=1,dtefield%mband_occ
1829          e0 = e0 + 2_dp*(umat(1,iband,jband,1)*qmat(2,jband,iband,ikpt,1,idir)&
1830 &         +       umat(2,iband,jband,1)*qmat(1,jband,iband,ikpt,1,idir))
1831        end do
1832      end do
1833      eberry = eberry - e0*fac
1834      e0=zero
1835      ikptn=dtefield%ikpt_dk(ikpt,8,idir)
1836      Amat(:,:,:)=zero
1837 !    calculate Amat
1838      do iband=1, dtefield%mband_occ
1839        do jband=1, dtefield%mband_occ
1840          do kband=1, dtefield%mband_occ
1841            Amat(1,iband,jband) = Amat(1,iband,jband) + (umat(1,iband,kband,3))*qmat(1,kband,jband,ikpt,1,idir)&
1842 &           - (umat(2,iband,kband,3))*qmat(2,kband,jband,ikpt,1,idir)
1843            Amat(2,iband,jband) = Amat(2,iband,jband) + (umat(1,iband,kband,3))*qmat(2,kband,jband,ikpt,1,idir)&
1844 &           + (umat(2,iband,kband,3))*qmat(1,kband,jband,ikpt,1,idir)
1845          end do
1846        end do
1847      end do
1848 
1849      do iband=1, dtefield%mband_occ
1850        do jband=1, dtefield%mband_occ
1851          do kband=1, dtefield%mband_occ
1852            z1(1) = (umat(1,jband,iband,4)+umat(1,iband,jband,2))*qmat(1,jband,kband,ikptn,1,idir)&
1853 &           - (umat(2,jband,iband,4)+umat(2,iband,jband,2))*qmat(2,jband,kband,ikptn,1,idir)
1854            z1(2) = (umat(1,jband,iband,4)+umat(1,iband,jband,2))*qmat(2,jband,kband,ikptn,1,idir)&
1855 &           + (umat(2,jband,iband,4)+umat(2,iband,jband,2))*qmat(1,jband,kband,ikptn,1,idir)
1856 
1857            e0 = e0 - 4_dp*(z1(1)*Amat(2,kband,iband)+z1(2)*Amat(1,kband,iband))
1858          end do
1859        end do
1860      end do
1861 
1862      eberry = eberry - e0*fac
1863 
1864 !    !---------------------------------last part---------------------------------------------
1865 
1866      fac = rprimd(idir_efield,idir)/&
1867 &     (dble(dtefield%nstr(idir))*two_pi)
1868 
1869 !    compute <u^(1)_{k_j-n,q}|u^(0)_{k_j+1}> matrix----------------------------------------------------
1870 
1871 !    prepare to calculate overlap matrix
1872      ikptn  = dtefield%ikpt_dk(ikpt,8,idir)
1873      ikpt1 = dtefield%ikpt_dk(ikpt,1,idir)
1874      icg = dtefield%cgindex(ikptn,1+nsppol)
1875      icg1 = dtefield%cgindex(ikpt1,1)
1876      npw_k1 = npwar1(ikptn)
1877      npw_k2 = npwarr(ikpt1)
1878      pwind_tmp(1:npw_k1) = pwindall((ikptn-1)*mpw_tmp+1:(ikptn-1)*mpw_tmp+npw_k1,5,idir)
1879 
1880      vect1(:,0) = zero ; vect2(:,0) = zero
1881      do jband = 1, dtefield%mband_occ
1882        vect2(:,1:npw_k2) = &
1883 &       cg(:,icg1 + 1 + (jband-1)*npw_k2*nspinor:icg1 + jband*npw_k2*nspinor)
1884        if (npw_k2 < mpw_tmp) vect2(:,npw_k2+1:mpw_tmp) = zero
1885        do iband = 1, dtefield%mband_occ
1886          pwmin = (iband-1)*npw_k1*nspinor
1887          pwmax = pwmin + npw_k1*nspinor
1888          vect1(:,1:npw_k1) = &
1889 &         cg1(:,icg + 1 + pwmin:icg + pwmax)
1890          if (npw_k1 < mpw_tmp) vect1(:,npw_k1+1:mpw_tmp) = zero
1891          call overlap_g(doti,dotr,mpw_tmp,npw_k1,npw_k2,nspinor,pwind_tmp,&
1892 &         vect1,vect2)
1893          umat(1,iband,jband,1) = dotr
1894          umat(2,iband,jband,1) = doti
1895        end do    ! iband
1896      end do    !jband
1897 
1898 !    compute <u^(0)_{k_j}|u^(1)_{k_j-n+1,q}> matrix----------------------------------------------------
1899 
1900 !    prepare to calculate overlap matrix
1901      ikpt1 = dtefield%ikpt_dk(ikpt,5,idir)
1902      icg = dtefield%cgindex(ikpt,1)
1903      icg1 = dtefield%cgindex(ikpt1,1+nsppol)
1904      npw_k1 = npwarr(ikpt)
1905      npw_k2 = npwar1(ikpt1)
1906      pwind_tmp(1:npw_k1) = pwindall((ikpt-1)*mpw_tmp+1:(ikpt-1)*mpw_tmp+npw_k1,7,idir)
1907      vect1(:,0) = zero ; vect2(:,0) = zero
1908      do jband = 1, dtefield%mband_occ
1909        vect2(:,1:npw_k2) = &
1910 &       cg1(:,icg1 + 1 + (jband-1)*npw_k2*nspinor:icg1 + jband*npw_k2*nspinor)
1911        if (npw_k2 < mpw_tmp) vect2(:,npw_k2+1:mpw_tmp) = zero
1912        do iband = 1, dtefield%mband_occ
1913          pwmin = (iband-1)*npw_k1*nspinor
1914          pwmax = pwmin + npw_k1*nspinor
1915          vect1(:,1:npw_k1) = &
1916 &         cg(:,icg + 1 + pwmin:icg + pwmax)
1917          if (npw_k1 < mpw_tmp) vect1(:,npw_k1+1:mpw_tmp) = zero
1918          call overlap_g(doti,dotr,mpw_tmp,npw_k1,npw_k2,nspinor,pwind_tmp,&
1919 &         vect1,vect2)
1920          umat(1,iband,jband,1) = umat(1,iband,jband,1) + dotr
1921          umat(2,iband,jband,1) = umat(2,iband,jband,1) + doti
1922        end do    ! iband
1923      end do    !jband
1924 
1925      e0=zero
1926 
1927      do iband=1,dtefield%mband_occ
1928        do jband=1,dtefield%mband_occ
1929          e0 = e0 +    (umat(1,iband,jband,1)*qmat(2,jband,iband,ikpt,1,idir)&
1930 &         +       umat(2,iband,jband,1)*qmat(1,jband,iband,ikpt,1,idir))
1931        end do
1932      end do
1933 
1934      eberry = eberry - e0*fac
1935 
1936    end do !end idir
1937  end do !end ikpt
1938 
1939  ABI_DEALLOCATE(umat)
1940  ABI_DEALLOCATE(vect1)
1941  ABI_DEALLOCATE(vect2)
1942  ABI_DEALLOCATE(pwind_tmp)
1943  ABI_DEALLOCATE(Amat)
1944 
1945 end subroutine dfptff_edie

ABINIT/dfptff_gbefd [ Functions ]

[ Top ] [ Functions ]

NAME

 dfptff_gbefd

FUNCTION

 calculate the gradient of the second order \Omega E \cdot P
 term, Eq.(23) in PRB 75, 115116(2007) [[cite:Wang2007]].

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 = the gradient of the Berry phase term

PARENTS

      dfpt_vtorho

CHILDREN

      overlap_g

SOURCE

1206 subroutine dfptff_gbefd(cg,cg1,dtefield,grad_berry,idir_efield,ikpt,isppol,mband,mpw,mpw1,mkmem,mk1mem,nkpt,&
1207 &                 npwarr,npwar1,nspinor,nsppol,qmat,pwindall,rprimd)
1208 
1209 
1210 !This section has been created automatically by the script Abilint (TD).
1211 !Do not modify the following lines by hand.
1212 #undef ABI_FUNC
1213 #define ABI_FUNC 'dfptff_gbefd'
1214 !End of the abilint section
1215 
1216  implicit none
1217 
1218 !Arguments ----------------------------------------
1219 !scalars
1220  integer,intent(in) :: idir_efield,ikpt,isppol,mband,mk1mem,mkmem,mpw,mpw1,nkpt
1221  integer,intent(in) :: nspinor,nsppol
1222  type(efield_type),intent(in) :: dtefield
1223 !arrays
1224  integer,intent(in) :: npwar1(nkpt),npwarr(nkpt)
1225  integer,intent(in) :: pwindall(max(mpw,mpw1)*mkmem,8,3)
1226  real(dp),intent(in) :: cg(2,mpw*nspinor*mband*mkmem*nsppol)
1227  real(dp),intent(in) :: cg1(2,mpw1*nspinor*mband*mk1mem*nsppol)
1228  real(dp),intent(in) :: qmat(2,dtefield%mband_occ,dtefield%mband_occ,nkpt,2,3)
1229  real(dp),intent(in) :: rprimd(3,3)
1230  real(dp),intent(out) :: grad_berry(2,mpw1,dtefield%mband_occ)
1231 
1232 !Local variables -------------------------
1233 !scalars
1234  integer :: iband,icg,icg1,idir,ikpt1
1235  integer :: ikptn,ikptnp1,ipw,jband,jpw,kband
1236  integer :: mpw_tmp,npw_k1,npw_k2,pwmax,pwmin
1237  real(dp) :: doti,dotr,fac,wfi,wfr
1238 !arrays
1239  integer,allocatable :: pwind_tmp(:)
1240  real(dp) :: z1(2),z2(2)
1241  real(dp),allocatable :: Amat(:,:,:),Bmat(:,:,:),s1mat(:,:,:),vect1(:,:)
1242  real(dp),allocatable :: vect2(:,:)
1243 
1244 ! *************************************************************************
1245 
1246  mpw_tmp=max(mpw,mpw1)
1247  ABI_ALLOCATE(vect1,(2,0:mpw_tmp))
1248  ABI_ALLOCATE(vect2,(2,0:mpw_tmp))
1249  ABI_ALLOCATE(s1mat,(2,dtefield%mband_occ,dtefield%mband_occ))
1250  ABI_ALLOCATE(pwind_tmp,(mpw_tmp))
1251  ABI_ALLOCATE(Amat,(2,dtefield%mband_occ,dtefield%mband_occ))
1252  ABI_ALLOCATE(Bmat,(2,dtefield%mband_occ,dtefield%mband_occ))
1253  vect1(:,0) = zero ; vect2(:,0) = zero
1254  s1mat(:,:,:)=zero
1255  grad_berry(:,:,:) = zero
1256 
1257  do idir=1,3
1258    fac = dtefield%efield_dot(idir)*dble(nkpt)/&
1259 &   (dble(dtefield%nstr(idir))*four_pi)
1260 
1261 !  prepare
1262    ikpt1 = dtefield%ikpt_dk(ikpt,1,idir)
1263    icg1 = dtefield%cgindex(ikpt1,isppol+nsppol)
1264    npw_k1 = npwar1(ikpt)
1265    npw_k2 = npwar1(ikpt1)
1266    pwind_tmp(1:npw_k1) = pwindall((ikpt-1)*mpw_tmp+1:(ikpt-1)*mpw_tmp+npw_k1,3,idir)
1267 
1268    do ipw = 1, npw_k1
1269      jpw = pwind_tmp(ipw)
1270      if (jpw > 0) then
1271        do iband = 1, dtefield%mband_occ
1272          wfr = cg1(1,icg1 + (iband - 1)*npw_k2*nspinor + jpw)
1273          wfi = cg1(2,icg1 + (iband - 1)*npw_k2*nspinor + jpw)
1274          do  jband = 1, dtefield%mband_occ
1275            grad_berry(1,ipw,jband) = &
1276 &           grad_berry(1,ipw,jband) + &
1277 &           fac*qmat(1,iband,jband,ikpt,1,idir)*wfr - fac*qmat(2,iband,jband,ikpt,1,idir)*wfi
1278 
1279            grad_berry(2,ipw,jband) = &
1280 &           grad_berry(2,ipw,jband) + &
1281 &           fac*qmat(1,iband,jband,ikpt,1,idir)*wfi + fac*qmat(2,iband,jband,ikpt,1,idir)*wfr
1282          end do
1283        end do
1284      end if
1285    end do
1286 
1287 !  compute <u^(0)_{k_j}|u^(1)_{k_j+1}> matrix----------------------------------------------------
1288 
1289 !  prepare to calculate overlap matrix
1290    ikptn = ikpt
1291    ikpt1 = dtefield%ikpt_dk(ikpt,1,idir)
1292    icg = dtefield%cgindex(ikptn,isppol)
1293    icg1 = dtefield%cgindex(ikpt1,isppol+nsppol)
1294    npw_k1 = npwarr(ikptn)
1295    npw_k2 = npwar1(ikpt1)
1296    pwind_tmp(1:npw_k1) = pwindall((ikptn-1)*mpw_tmp+1:(ikptn-1)*mpw_tmp+npw_k1,7,idir)
1297 
1298    vect1(:,0) = zero ; vect2(:,0) = zero
1299    do jband = 1, dtefield%mband_occ
1300      vect2(:,1:npw_k2) = &
1301 &     cg1(:,icg1 + 1 + (jband-1)*npw_k2*nspinor:icg1 + jband*npw_k2*nspinor)
1302      if (npw_k2 < mpw_tmp) vect2(:,npw_k2+1:mpw_tmp) = zero
1303      do iband = 1, dtefield%mband_occ
1304        pwmin = (iband-1)*npw_k1*nspinor
1305        pwmax = pwmin + npw_k1*nspinor
1306        vect1(:,1:npw_k1) = &
1307 &       cg(:,icg + 1 + pwmin:icg + pwmax)
1308        if (npw_k1 < mpw_tmp) vect1(:,npw_k1+1:mpw_tmp) = zero
1309        call overlap_g(doti,dotr,mpw_tmp,npw_k1,npw_k2,nspinor,pwind_tmp,&
1310 &       vect1,vect2)
1311        s1mat(1,iband,jband) = dotr
1312        s1mat(2,iband,jband) = doti
1313      end do    ! iband
1314    end do    !jband
1315 
1316 !  compute <u^(1)_{k_j}|u^(0)_{k_j+1}> matrix-----------------------------------------------------
1317 
1318 !  prepare to calculate overlap matrix
1319    ikpt1 = dtefield%ikpt_dk(ikpt,1,idir)
1320    icg = dtefield%cgindex(ikpt,isppol+nsppol)
1321    icg1 = dtefield%cgindex(ikpt1,isppol)
1322    npw_k1 = npwar1(ikpt)
1323    npw_k2 = npwarr(ikpt1)
1324    pwind_tmp(1:npw_k1) = pwindall((ikpt-1)*mpw_tmp+1:(ikpt-1)*mpw_tmp+npw_k1,5,idir)
1325 
1326    vect1(:,0) = zero ; vect2(:,0) = zero
1327    do jband = 1, dtefield%mband_occ
1328      vect2(:,1:npw_k2) = &
1329 &     cg(:,icg1 + 1 + (jband-1)*npw_k2*nspinor:icg1 + jband*npw_k2*nspinor)
1330      if (npw_k2 < mpw_tmp) vect2(:,npw_k2+1:mpw_tmp) = zero
1331      do iband = 1, dtefield%mband_occ
1332        pwmin = (iband-1)*npw_k1*nspinor
1333        pwmax = pwmin + npw_k1*nspinor
1334        vect1(:,1:npw_k1) = &
1335 &       cg1(:,icg + 1 + pwmin:icg + pwmax)
1336        if (npw_k1 < mpw_tmp) vect1(:,npw_k1+1:mpw_tmp) = zero
1337        call overlap_g(doti,dotr,mpw_tmp,npw_k1,npw_k2,nspinor,pwind_tmp,&
1338 &       vect1,vect2)
1339        s1mat(1,jband,iband) = s1mat(1,jband,iband) + dotr
1340        s1mat(2,jband,iband) = s1mat(2,jband,iband) + doti
1341      end do    ! iband
1342    end do    !jband
1343 
1344    Amat(:,:,:)=zero
1345 
1346 !  calculate Amat
1347    do iband=1, dtefield%mband_occ
1348      do jband=1, dtefield%mband_occ
1349        do kband=1, dtefield%mband_occ
1350          Amat(1,iband,jband) = Amat(1,iband,jband) + s1mat(1,iband,kband)*qmat(1,kband,jband,ikpt,1,idir)&
1351 &         - s1mat(2,iband,kband)*qmat(2,kband,jband,ikpt,1,idir)
1352          Amat(2,iband,jband) = Amat(2,iband,jband) + s1mat(1,iband,kband)*qmat(2,kband,jband,ikpt,1,idir)&
1353 &         + s1mat(2,iband,kband)*qmat(1,kband,jband,ikpt,1,idir)
1354        end do
1355      end do
1356    end do
1357 
1358    Bmat(:,:,:)=zero
1359 
1360 !  calculate Bmat
1361    ikptn = dtefield%ikpt_dk(ikpt,7,idir)
1362    do iband=1, dtefield%mband_occ
1363      do jband=1, dtefield%mband_occ
1364        do kband=1, dtefield%mband_occ
1365          Bmat(1,jband,kband) = Bmat(1,jband,kband) + Amat(1,iband,kband)*qmat(1,jband,iband,ikptn,1,idir)&
1366 &         - Amat(2,iband,kband)*qmat(2,jband,iband,ikptn,1,idir)
1367          Bmat(2,jband,kband) = Bmat(2,jband,kband) + Amat(1,iband,kband)*qmat(2,jband,iband,ikptn,1,idir)&
1368 &         + Amat(2,iband,kband)*qmat(1,jband,iband,ikptn,1,idir)
1369        end do
1370      end do
1371    end do
1372 
1373 !  calc. the second term of gradient------------------------------
1374 
1375 !  preparation
1376 
1377    ikptnp1 = dtefield%ikpt_dk(ikpt,3,idir)
1378    icg = dtefield%cgindex(ikptnp1,isppol)
1379    npw_k1 = npwar1(ikpt)
1380    npw_k2 = npwarr(ikptnp1)
1381    pwind_tmp(1:npw_k1) = pwindall((ikpt-1)*mpw_tmp+1:(ikpt-1)*mpw_tmp+npw_k1,5,idir)
1382 
1383    z1(:) = zero
1384    z2(:) = zero
1385 
1386    do ipw = 1, npw_k1
1387      jpw = pwind_tmp(ipw)
1388      if (jpw > 0) then
1389        do iband = 1, dtefield%mband_occ
1390          wfr = cg(1,icg + (iband - 1)*npw_k2*nspinor + jpw)
1391          wfi = cg(2,icg + (iband - 1)*npw_k2*nspinor + jpw)
1392          do jband=1, dtefield%mband_occ
1393            grad_berry(1,ipw,jband) = grad_berry(1,ipw,jband) - fac*(Bmat(1,iband,jband)*wfr - Bmat(2,iband,jband)*wfi)
1394            grad_berry(2,ipw,jband) = grad_berry(2,ipw,jband) - fac*(Bmat(1,iband,jband)*wfi + Bmat(2,iband,jband)*wfr)
1395          end do
1396        end do
1397      end if
1398    end do
1399 
1400 !  Second part of gradient of Berry phase++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1401 
1402    vect1(:,0) = zero ; vect2(:,0) = zero
1403 
1404 !  prepare
1405    ikpt1 = dtefield%ikpt_dk(ikpt,2,idir)
1406    icg1 = dtefield%cgindex(ikpt1,isppol+nsppol)
1407    npw_k1 = npwar1(ikpt)
1408    npw_k2 = npwar1(ikpt1)
1409    pwind_tmp(1:npw_k1) = pwindall((ikpt-1)*mpw_tmp+1:(ikpt-1)*mpw_tmp+npw_k1,4,idir)
1410 
1411    do ipw = 1, npw_k1
1412      jpw = pwind_tmp(ipw)
1413      if (jpw > 0) then
1414        do iband = 1, dtefield%mband_occ
1415          wfr = cg1(1,icg1 + (iband - 1)*npw_k2*nspinor + jpw)
1416          wfi = cg1(2,icg1 + (iband - 1)*npw_k2*nspinor + jpw)
1417          do  jband = 1, dtefield%mband_occ
1418            grad_berry(1,ipw,jband) = &
1419 &           grad_berry(1,ipw,jband) - &
1420 &           fac*qmat(1,iband,jband,ikpt,2,idir)*wfr + fac*qmat(2,iband,jband,ikpt,2,idir)*wfi
1421 
1422            grad_berry(2,ipw,jband) = &
1423 &           grad_berry(2,ipw,jband) - &
1424 &           fac*qmat(1,iband,jband,ikpt,2,idir)*wfi - fac*qmat(2,iband,jband,ikpt,2,idir)*wfr
1425          end do
1426        end do
1427      end if
1428    end do
1429 
1430 !  compute <u^(0)_{k_j}|u^(1)_{k_j-1}> matrix----------------------------------------------------
1431 
1432 !  prepare to calculate overlap matrix
1433    ikptn = ikpt
1434    ikpt1 = dtefield%ikpt_dk(ikpt,2,idir)
1435    icg = dtefield%cgindex(ikptn,isppol)
1436    icg1 = dtefield%cgindex(ikpt1,isppol+nsppol)
1437    npw_k1 = npwarr(ikptn)
1438    npw_k2 = npwar1(ikpt1)
1439    pwind_tmp(1:npw_k1) =pwindall((ikptn-1)*mpw_tmp+1:(ikptn-1)*mpw_tmp+npw_k1,8,idir)
1440 
1441    vect1(:,0) = zero ; vect2(:,0) = zero
1442    do jband = 1, dtefield%mband_occ
1443      vect2(:,1:npw_k2) = &
1444 &     cg1(:,icg1 + 1 + (jband-1)*npw_k2*nspinor:icg1 + jband*npw_k2*nspinor)
1445      if (npw_k2 < mpw_tmp) vect2(:,npw_k2+1:mpw_tmp) = zero
1446      do iband = 1, dtefield%mband_occ
1447        pwmin = (iband-1)*npw_k1*nspinor
1448        pwmax = pwmin + npw_k1*nspinor
1449        vect1(:,1:npw_k1) = &
1450 &       cg(:,icg + 1 + pwmin:icg + pwmax)
1451        if (npw_k1 < mpw_tmp) vect1(:,npw_k1+1:mpw_tmp) = zero
1452        call overlap_g(doti,dotr,mpw_tmp,npw_k1,npw_k2,nspinor,pwind_tmp,&
1453 &       vect1,vect2)
1454        s1mat(1,iband,jband) = dotr
1455        s1mat(2,iband,jband) = doti
1456      end do    ! iband
1457    end do    !jband
1458 
1459 !  compute <u^(1)_{k_j}|u^(0)_{k_j-1}> matrix-----------------------------------------------------
1460 
1461 !  prepare to calculate overlap matrix
1462    ikpt1 = dtefield%ikpt_dk(ikpt,2,idir)
1463    icg = dtefield%cgindex(ikpt,isppol+nsppol)
1464    icg1 = dtefield%cgindex(ikpt1,isppol)
1465    npw_k1 = npwarr(ikpt)
1466    npw_k2 = npwar1(ikpt1)
1467    pwind_tmp(1:npw_k1) =pwindall((ikpt-1)*mpw_tmp+1:(ikpt-1)*mpw_tmp+npw_k1,6,idir)
1468 
1469    vect1(:,0) = zero ; vect2(:,0) = zero
1470    do jband = 1, dtefield%mband_occ
1471      vect2(:,1:npw_k2) = &
1472 &     cg(:,icg1 + 1 + (jband-1)*npw_k2*nspinor:icg1 + jband*npw_k2*nspinor)
1473      if (npw_k2 < mpw_tmp) vect2(:,npw_k2+1:mpw_tmp) = zero
1474      do iband = 1, dtefield%mband_occ
1475        pwmin = (iband-1)*npw_k1*nspinor
1476        pwmax = pwmin + npw_k1*nspinor
1477        vect1(:,1:npw_k1) = &
1478 &       cg1(:,icg + 1 + pwmin:icg + pwmax)
1479        if (npw_k1 < mpw_tmp) vect1(:,npw_k1+1:mpw_tmp) = zero
1480        call overlap_g(doti,dotr,mpw_tmp,npw_k1,npw_k2,nspinor,pwind_tmp,&
1481 &       vect1,vect2)
1482        s1mat(1,jband,iband) = s1mat(1,jband,iband) + dotr
1483        s1mat(2,jband,iband) = s1mat(2,jband,iband) + doti
1484      end do    ! iband
1485    end do    !jband
1486 
1487    Amat(:,:,:)=zero
1488 
1489 !  calculate Amat
1490    do iband=1, dtefield%mband_occ
1491      do jband=1, dtefield%mband_occ
1492        do kband=1, dtefield%mband_occ
1493          Amat(1,iband,jband) = Amat(1,iband,jband) + s1mat(1,iband,kband)*qmat(1,kband,jband,ikpt,2,idir)&
1494 &         - s1mat(2,iband,kband)*qmat(2,kband,jband,ikpt,2,idir)
1495          Amat(2,iband,jband) = Amat(2,iband,jband) + s1mat(1,iband,kband)*qmat(2,kband,jband,ikpt,2,idir)&
1496 &         + s1mat(2,iband,kband)*qmat(1,kband,jband,ikpt,2,idir)
1497        end do
1498      end do
1499    end do
1500 
1501    Bmat(:,:,:)=zero
1502 
1503 !  calculate Bmat
1504    ikptn = dtefield%ikpt_dk(ikpt,7,idir)
1505    do iband=1, dtefield%mband_occ
1506      do jband=1, dtefield%mband_occ
1507        do kband=1, dtefield%mband_occ
1508          Bmat(1,jband,kband) = Bmat(1,jband,kband) + Amat(1,iband,kband)*qmat(1,jband,iband,ikptn,2,idir)&
1509 &         - Amat(2,iband,kband)*qmat(2,jband,iband,ikptn,2,idir)
1510          Bmat(2,jband,kband) = Bmat(2,jband,kband) + Amat(1,iband,kband)*qmat(2,jband,iband,ikptn,2,idir)&
1511          + Amat(2,iband,kband)*qmat(1,jband,iband,ikptn,2,idir)
1512        end do
1513      end do
1514    end do
1515 
1516 !  calc. the second term of gradient------------------------------
1517 
1518 !  preparation
1519 
1520    ikptnp1 = dtefield%ikpt_dk(ikpt,4,idir)
1521    icg = dtefield%cgindex(ikptnp1,isppol)
1522    npw_k1 = npwar1(ikpt)
1523    npw_k2 = npwarr(ikptnp1)
1524    pwind_tmp(1:npw_k1) =pwindall((ikpt-1)*mpw_tmp+1:(ikpt-1)*mpw_tmp+npw_k1,6,idir)
1525    z1(:) = zero
1526    z2(:) = zero
1527    do ipw = 1, npw_k1
1528      jpw = pwind_tmp(ipw)
1529      if (jpw > 0) then
1530        do iband = 1, dtefield%mband_occ
1531          wfr = cg(1,icg + (iband - 1)*npw_k2*nspinor + jpw)
1532          wfi = cg(2,icg + (iband - 1)*npw_k2*nspinor + jpw)
1533          do jband=1, dtefield%mband_occ
1534            grad_berry(1,ipw,jband) = grad_berry(1,ipw,jband) + fac*(Bmat(1,iband,jband)*wfr - Bmat(2,iband,jband)*wfi)
1535            grad_berry(2,ipw,jband) = grad_berry(2,ipw,jband) + fac*(Bmat(1,iband,jband)*wfi + Bmat(2,iband,jband)*wfr)
1536          end do
1537        end do
1538      end if
1539    end do
1540 
1541  end do !idir
1542 
1543 !!----------------------------------------third part of gradient------------------------------------------------------
1544  do idir=1,3
1545    fac = rprimd(idir_efield,idir)*dble(nkpt)/&
1546 &   (dble(dtefield%nstr(idir))*four_pi)
1547 
1548 !  prepare
1549    ikpt1 = dtefield%ikpt_dk(ikpt,1,idir)
1550    icg1 = dtefield%cgindex(ikpt1,isppol)
1551    npw_k1 = npwar1(ikpt)
1552    npw_k2 = npwarr(ikpt1)
1553    pwind_tmp(1:npw_k1) = pwindall((ikpt-1)*mpw_tmp+1:(ikpt-1)*mpw_tmp+npw_k1,5,idir)
1554    do ipw = 1, npw_k1
1555      jpw = pwind_tmp(ipw)
1556      if (jpw > 0) then
1557        do iband = 1, dtefield%mband_occ
1558          wfr = cg(1,icg1 + (iband - 1)*npw_k2*nspinor + jpw)
1559          wfi = cg(2,icg1 + (iband - 1)*npw_k2*nspinor + jpw)
1560          do  jband = 1, dtefield%mband_occ
1561            grad_berry(1,ipw,jband) = &
1562 &           grad_berry(1,ipw,jband) + &
1563 &           fac*qmat(1,iband,jband,ikpt,1,idir)*wfr - fac*qmat(2,iband,jband,ikpt,1,idir)*wfi
1564            grad_berry(2,ipw,jband) = &
1565 &           grad_berry(2,ipw,jband) + &
1566 &           fac*qmat(1,iband,jband,ikpt,1,idir)*wfi + fac*qmat(2,iband,jband,ikpt,1,idir)*wfr
1567          end do
1568        end do
1569      end if
1570    end do
1571 
1572 !  prepare
1573    ikpt1 = dtefield%ikpt_dk(ikpt,2,idir)
1574    icg1 = dtefield%cgindex(ikpt1,isppol)
1575    npw_k1 = npwar1(ikpt)
1576    npw_k2 = npwarr(ikpt1)
1577    pwind_tmp(1:npw_k1) = pwindall((ikpt-1)*mpw_tmp+1:(ikpt-1)*mpw_tmp+npw_k1,6,idir)
1578 
1579    do ipw = 1, npw_k1
1580      jpw = pwind_tmp(ipw)
1581      if (jpw > 0) then
1582        do iband = 1, dtefield%mband_occ
1583          wfr = cg(1,icg1 + (iband - 1)*npw_k2*nspinor + jpw)
1584          wfi = cg(2,icg1 + (iband - 1)*npw_k2*nspinor + jpw)
1585          do  jband = 1, dtefield%mband_occ
1586            grad_berry(1,ipw,jband) = &
1587 &           grad_berry(1,ipw,jband) - &
1588 &           fac*qmat(1,iband,jband,ikpt,2,idir)*wfr + fac*qmat(2,iband,jband,ikpt,2,idir)*wfi
1589 
1590            grad_berry(2,ipw,jband) = &
1591 &           grad_berry(2,ipw,jband) - &
1592 &           fac*qmat(1,iband,jband,ikpt,2,idir)*wfi - fac*qmat(2,iband,jband,ikpt,2,idir)*wfr
1593 
1594          end do
1595        end do
1596      end if
1597    end do
1598 
1599  end do !idir
1600 
1601  ABI_DEALLOCATE(vect1)
1602  ABI_DEALLOCATE(vect2)
1603  ABI_DEALLOCATE(s1mat)
1604  ABI_DEALLOCATE(Amat)
1605  ABI_DEALLOCATE(Bmat)
1606  ABI_DEALLOCATE(pwind_tmp)
1607 
1608 end subroutine dfptff_gbefd

ABINIT/dfptff_gradberry [ Functions ]

[ Top ] [ 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

 755 subroutine dfptff_gradberry(cg,cg1,dtefield,grad_berry,ikpt,isppol,mband,mpw,mpw1,mkmem,mk1mem,nkpt,&
 756 &                     npwarr,npwar1,nspinor,nsppol,qmat,pwindall)
 757 
 758 
 759 !This section has been created automatically by the script Abilint (TD).
 760 !Do not modify the following lines by hand.
 761 #undef ABI_FUNC
 762 #define ABI_FUNC 'dfptff_gradberry'
 763 !End of the abilint section
 764 
 765  implicit none
 766 
 767 !Arguments ----------------------------------------
 768 !scalars
 769  integer,intent(in) :: ikpt,isppol,mband,mk1mem,mkmem,mpw,mpw1,nkpt,nspinor
 770  integer,intent(in) :: nsppol
 771  type(efield_type),intent(in) :: dtefield
 772 !arrays
 773  integer,intent(in) :: npwar1(nkpt),npwarr(nkpt)
 774  integer,intent(in) :: pwindall(max(mpw,mpw1)*mkmem,8,3)
 775  real(dp),intent(in) :: cg(2,mpw*nspinor*mband*mkmem*nsppol)
 776  real(dp),intent(in) :: cg1(2,mpw1*nspinor*mband*mk1mem*nsppol)
 777  real(dp),intent(in) :: qmat(2,dtefield%mband_occ,dtefield%mband_occ,nkpt,2,3)
 778  real(dp),intent(out) :: grad_berry(2,mpw1,dtefield%mband_occ)
 779 
 780 !Local variables -------------------------
 781 !scalars
 782  integer :: iband,icg,icg1,idir,ikpt1
 783  integer :: ikpt1m,ikptn,ikptnm,ikptnp1,ipw,jband,jpw,kband
 784  integer :: mpw_tmp,npw_k1,npw_k2,pwmax,pwmin
 785  real(dp) :: doti,dotr,fac,wfi,wfr
 786 !arrays
 787  integer,allocatable :: pwind_tmp(:)
 788  real(dp) :: z1(2),z2(2)
 789  real(dp),allocatable :: Amat(:,:,:),Bmat(:,:,:),s1mat(:,:,:),vect1(:,:)
 790  real(dp),allocatable :: vect2(:,:)
 791 
 792 ! *************************************************************************
 793 
 794  mpw_tmp=max(mpw,mpw1)
 795  ABI_ALLOCATE(vect1,(2,0:mpw_tmp))
 796  ABI_ALLOCATE(vect2,(2,0:mpw_tmp))
 797  ABI_ALLOCATE(s1mat,(2,dtefield%mband_occ,dtefield%mband_occ))
 798  ABI_ALLOCATE(pwind_tmp,(mpw_tmp))
 799  ABI_ALLOCATE(Amat,(2,dtefield%mband_occ,dtefield%mband_occ))
 800  ABI_ALLOCATE(Bmat,(2,dtefield%mband_occ,dtefield%mband_occ))
 801  vect1(:,0) = zero ; vect2(:,0) = zero
 802  s1mat(:,:,:)=zero
 803  grad_berry(:,:,:) = zero
 804 
 805  do idir=1,3
 806    fac = dtefield%efield_dot(idir)*dble(nkpt)/&
 807 &   (dble(dtefield%nstr(idir))*four_pi)
 808 
 809 !  prepare
 810    ikpt1 = dtefield%ikpt_dk(ikpt,1,idir)
 811    icg1 = dtefield%cgindex(ikpt1,isppol+nsppol)
 812    npw_k1 = npwar1(ikpt)
 813    npw_k2 = npwar1(ikpt1)
 814    pwind_tmp(1:npw_k1) = pwindall((ikpt-1)*mpw_tmp+1:(ikpt-1)*mpw_tmp+npw_k1,3,idir)
 815 
 816    do ipw = 1, npw_k1
 817      jpw = pwind_tmp(ipw)
 818 
 819      if (jpw > 0) then
 820        do iband = 1, dtefield%mband_occ
 821          wfr = cg1(1,icg1 + (iband - 1)*npw_k2*nspinor + jpw)
 822          wfi = cg1(2,icg1 + (iband - 1)*npw_k2*nspinor + jpw)
 823 
 824          do  jband = 1, dtefield%mband_occ
 825 
 826            grad_berry(1,ipw,jband) = &
 827 &           grad_berry(1,ipw,jband) + &
 828 &           fac*qmat(1,iband,jband,ikpt,1,idir)*wfr - fac*qmat(2,iband,jband,ikpt,1,idir)*wfi
 829 
 830            grad_berry(2,ipw,jband) = &
 831 &           grad_berry(2,ipw,jband) + &
 832 &           fac*qmat(1,iband,jband,ikpt,1,idir)*wfi + fac*qmat(2,iband,jband,ikpt,1,idir)*wfr
 833 
 834          end do
 835        end do
 836      end if
 837 
 838    end do
 839 
 840 !  compute <u^(0)_{k_j+n}|u^(1)_{k_j+1,q}> matrix----------------------------------------------------
 841 
 842 !  prepare to calculate overlap matrix
 843    ikptn = dtefield%ikpt_dk(ikpt,7,idir)
 844    ikpt1 = dtefield%ikpt_dk(ikpt,1,idir)
 845    icg = dtefield%cgindex(ikptn,isppol)
 846    icg1 = dtefield%cgindex(ikpt1,isppol+nsppol)
 847    npw_k1 = npwarr(ikptn)
 848    npw_k2 = npwar1(ikpt1)
 849    pwind_tmp(1:npw_k1) = pwindall((ikptn-1)*mpw_tmp+1:(ikptn-1)*mpw_tmp+npw_k1,7,idir)
 850 
 851 
 852    vect1(:,0) = zero ; vect2(:,0) = zero
 853    do jband = 1, dtefield%mband_occ
 854      vect2(:,1:npw_k2) = &
 855 &     cg1(:,icg1 + 1 + (jband-1)*npw_k2*nspinor:icg1 + jband*npw_k2*nspinor)
 856      if (npw_k2 < mpw_tmp) vect2(:,npw_k2+1:mpw_tmp) = zero
 857 
 858      do iband = 1, dtefield%mband_occ
 859 
 860        pwmin = (iband-1)*npw_k1*nspinor
 861        pwmax = pwmin + npw_k1*nspinor
 862        vect1(:,1:npw_k1) = &
 863 &       cg(:,icg + 1 + pwmin:icg + pwmax)
 864        if (npw_k1 < mpw_tmp) vect1(:,npw_k1+1:mpw_tmp) = zero
 865        call overlap_g(doti,dotr,mpw_tmp,npw_k1,npw_k2,nspinor,pwind_tmp,&
 866 &       vect1,vect2)
 867        s1mat(1,iband,jband) = dotr
 868        s1mat(2,iband,jband) = doti
 869 
 870      end do    ! iband
 871    end do    !jband
 872 
 873 !  compute <u^(0)_{-k_j+1}|u^(1)_{-k_j+n},q> matrix--------------------
 874 
 875 !  prepare to calculate overlap matrix
 876    ikptn = dtefield%ikpt_dk(ikpt,7,idir)
 877    ikptnm= dtefield%ikpt_dk(ikptn,9,idir)
 878    ikpt1 = dtefield%ikpt_dk(ikpt,1,idir)
 879    ikpt1m= dtefield%ikpt_dk(ikpt1,9,idir)
 880 
 881    icg = dtefield%cgindex(ikpt1m,isppol)
 882    icg1 = dtefield%cgindex(ikptnm,isppol+nsppol)
 883    npw_k1 = npwarr(ikpt1m)
 884    npw_k2 = npwar1(ikptnm)
 885    pwind_tmp(1:npw_k1) = pwindall((ikpt1m-1)*mpw_tmp+1:(ikpt1m-1)*mpw_tmp+npw_k1,7,idir)
 886 
 887    vect1(:,0) = zero ; vect2(:,0) = zero
 888    do jband = 1, dtefield%mband_occ
 889      vect2(:,1:npw_k2) = &
 890 &     cg1(:,icg1 + 1 + (jband-1)*npw_k2*nspinor:icg1 + jband*npw_k2*nspinor)
 891      if (npw_k2 < mpw_tmp) vect2(:,npw_k2+1:mpw_tmp) = zero
 892 
 893      do iband = 1, dtefield%mband_occ
 894 
 895        pwmin = (iband-1)*npw_k1*nspinor
 896        pwmax = pwmin + npw_k1*nspinor
 897        vect1(:,1:npw_k1) = &
 898 &       cg(:,icg + 1 + pwmin:icg + pwmax)
 899        if (npw_k1 < mpw_tmp) vect1(:,npw_k1+1:mpw_tmp) = zero
 900        call overlap_g(doti,dotr,mpw_tmp,npw_k1,npw_k2,nspinor,pwind_tmp,&
 901 &       vect1,vect2)
 902 
 903        s1mat(1,jband,iband) = s1mat(1,jband,iband) + dotr
 904        s1mat(2,jband,iband) = s1mat(2,jband,iband) + doti
 905 
 906      end do    ! iband
 907    end do    !jband
 908 
 909    Amat(:,:,:)=zero
 910 
 911 !  calculate Amat
 912    do iband=1, dtefield%mband_occ
 913      do jband=1, dtefield%mband_occ
 914        do kband=1, dtefield%mband_occ
 915          Amat(1,iband,jband) = Amat(1,iband,jband) + s1mat(1,iband,kband)*&
 916 &         qmat(1,kband,jband,ikpt,1,idir)&
 917 &         - s1mat(2,iband,kband)*qmat(2,kband,jband,ikpt,1,idir)
 918          Amat(2,iband,jband) = Amat(2,iband,jband) + s1mat(1,iband,kband)*&
 919 &         qmat(2,kband,jband,ikpt,1,idir)&
 920 &         + s1mat(2,iband,kband)*qmat(1,kband,jband,ikpt,1,idir)
 921        end do
 922      end do
 923    end do
 924 
 925    Bmat(:,:,:)=zero
 926 
 927 !  calculate Bmat
 928    ikptn = dtefield%ikpt_dk(ikpt,7,idir)
 929    do iband=1, dtefield%mband_occ
 930      do jband=1, dtefield%mband_occ
 931        do kband=1, dtefield%mband_occ
 932          Bmat(1,jband,kband) = Bmat(1,jband,kband) + Amat(1,iband,kband)*&
 933 &         qmat(1,jband,iband,ikptn,1,idir)&
 934 &         - Amat(2,iband,kband)*qmat(2,jband,iband,ikptn,1,idir)
 935          Bmat(2,jband,kband) = Bmat(2,jband,kband) + Amat(1,iband,kband)*&
 936 &         qmat(2,jband,iband,ikptn,1,idir)&
 937 &         + Amat(2,iband,kband)*qmat(1,jband,iband,ikptn,1,idir)
 938        end do
 939      end do
 940    end do
 941 
 942 !  calc. the second term of gradient------------------------------
 943 
 944 !  preparation
 945 
 946    ikptnp1 = dtefield%ikpt_dk(ikpt,3,idir)
 947    icg = dtefield%cgindex(ikptnp1,isppol)
 948    npw_k1 = npwar1(ikpt)
 949    npw_k2 = npwarr(ikptnp1)
 950    pwind_tmp(1:npw_k1) = pwindall((ikpt-1)*mpw_tmp+1:(ikpt-1)*mpw_tmp+npw_k1,5,idir)
 951 
 952    z1(:) = zero
 953    z2(:) = zero
 954 
 955    do ipw = 1, npw_k1
 956 
 957      jpw = pwind_tmp(ipw)
 958 
 959      if (jpw > 0) then
 960 
 961        do iband = 1, dtefield%mband_occ
 962          wfr = cg(1,icg + (iband - 1)*npw_k2*nspinor + jpw)
 963          wfi = cg(2,icg + (iband - 1)*npw_k2*nspinor + jpw)
 964 
 965          do jband=1, dtefield%mband_occ
 966 
 967            grad_berry(1,ipw,jband) = grad_berry(1,ipw,jband) &
 968 &           - fac*(Bmat(1,iband,jband)*wfr - Bmat(2,iband,jband)*wfi)
 969            grad_berry(2,ipw,jband) = grad_berry(2,ipw,jband) &
 970 &           - fac*(Bmat(1,iband,jband)*wfi + Bmat(2,iband,jband)*wfr)
 971 
 972          end do
 973        end do
 974      end if
 975    end do
 976 
 977 !  Second part of gradient of Berry phase++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 978 
 979    vect1(:,0) = zero ; vect2(:,0) = zero
 980 
 981 !  prepare
 982    ikpt1 = dtefield%ikpt_dk(ikpt,2,idir)
 983    icg1 = dtefield%cgindex(ikpt1,isppol+nsppol)
 984    npw_k1 = npwar1(ikpt)
 985    npw_k2 = npwar1(ikpt1)
 986    pwind_tmp(1:npw_k1) = pwindall((ikpt-1)*mpw_tmp+1:(ikpt-1)*mpw_tmp+npw_k1,4,idir)
 987 
 988 !  write(std_out,*)'dfpt_cgwf:pwind_tmp',pwind_tmp
 989 !  stop
 990 
 991    do ipw = 1, npw_k1
 992 
 993      jpw = pwind_tmp(ipw)
 994 
 995      if (jpw > 0) then
 996        do iband = 1, dtefield%mband_occ
 997          wfr = cg1(1,icg1 + (iband - 1)*npw_k2*nspinor + jpw)
 998          wfi = cg1(2,icg1 + (iband - 1)*npw_k2*nspinor + jpw)
 999 
1000          do  jband = 1, dtefield%mband_occ
1001 
1002            grad_berry(1,ipw,jband) = &
1003 &           grad_berry(1,ipw,jband) - &
1004 &           fac*qmat(1,iband,jband,ikpt,2,idir)*wfr + fac*qmat(2,iband,jband,ikpt,2,idir)*wfi
1005 
1006            grad_berry(2,ipw,jband) = &
1007 &           grad_berry(2,ipw,jband) - &
1008 &           fac*qmat(1,iband,jband,ikpt,2,idir)*wfi - fac*qmat(2,iband,jband,ikpt,2,idir)*wfr
1009 
1010          end do
1011        end do
1012      end if
1013    end do
1014 
1015 
1016 !  compute <u^(0)_{k_j+n}|u^(1)_{k_j-1,q}> matrix----------------------------------------------------
1017 
1018 !  prepare to calculate overlap matrix
1019    ikptn = dtefield%ikpt_dk(ikpt,7,idir)
1020    ikpt1 = dtefield%ikpt_dk(ikpt,2,idir)
1021    icg = dtefield%cgindex(ikptn,isppol)
1022    icg1 = dtefield%cgindex(ikpt1,isppol+nsppol)
1023    npw_k1 = npwarr(ikptn)
1024    npw_k2 = npwar1(ikpt1)
1025    pwind_tmp(1:npw_k1) =pwindall((ikptn-1)*mpw_tmp+1:(ikptn-1)*mpw_tmp+npw_k1,8,idir)
1026 
1027    vect1(:,0) = zero ; vect2(:,0) = zero
1028    do jband = 1, dtefield%mband_occ
1029      vect2(:,1:npw_k2) = &
1030 &     cg1(:,icg1 + 1 + (jband-1)*npw_k2*nspinor:icg1 + jband*npw_k2*nspinor)
1031      if (npw_k2 < mpw_tmp) vect2(:,npw_k2+1:mpw_tmp) = zero
1032 
1033      do iband = 1, dtefield%mband_occ
1034 
1035        pwmin = (iband-1)*npw_k1*nspinor
1036        pwmax = pwmin + npw_k1*nspinor
1037        vect1(:,1:npw_k1) = &
1038 &       cg(:,icg + 1 + pwmin:icg + pwmax)
1039        if (npw_k1 < mpw_tmp) vect1(:,npw_k1+1:mpw_tmp) = zero
1040        call overlap_g(doti,dotr,mpw_tmp,npw_k1,npw_k2,nspinor,pwind_tmp,&
1041 &       vect1,vect2)
1042        s1mat(1,iband,jband) = dotr
1043        s1mat(2,iband,jband) = doti
1044 
1045      end do    ! iband
1046    end do    !jband
1047 
1048 !  compute <u^(0)_{-k_j-1}|u^(1)_{-k_j+n,q}> matrix-----------------------------------------------------
1049 
1050 !  prepare to calculate overlap matrix
1051    ikptn = dtefield%ikpt_dk(ikpt,7,idir)
1052    ikptnm= dtefield%ikpt_dk(ikptn,9,idir)
1053    ikpt1 = dtefield%ikpt_dk(ikpt,2,idir)
1054    ikpt1m= dtefield%ikpt_dk(ikpt1,9,idir)
1055    icg = dtefield%cgindex(ikpt1m,isppol)
1056    icg1 = dtefield%cgindex(ikptnm,isppol+nsppol)
1057    npw_k1 = npwarr(ikpt1m)
1058    npw_k2 = npwar1(ikptnm)
1059    pwind_tmp(1:npw_k1) =pwindall((ikpt1m-1)*mpw_tmp+1:(ikpt1m-1)*mpw_tmp+npw_k1,8,idir)
1060 
1061 
1062 
1063    vect1(:,0) = zero ; vect2(:,0) = zero
1064    do jband = 1, dtefield%mband_occ
1065      vect2(:,1:npw_k2) = &
1066 &     cg1(:,icg1 + 1 + (jband-1)*npw_k2*nspinor:icg1 + jband*npw_k2*nspinor)
1067      if (npw_k2 < mpw_tmp) vect2(:,npw_k2+1:mpw_tmp) = zero
1068      do iband = 1, dtefield%mband_occ
1069        pwmin = (iband-1)*npw_k1*nspinor
1070        pwmax = pwmin + npw_k1*nspinor
1071        vect1(:,1:npw_k1) = &
1072 &       cg(:,icg + 1 + pwmin:icg + pwmax)
1073        if (npw_k1 < mpw_tmp) vect1(:,npw_k1+1:mpw_tmp) = zero
1074        call overlap_g(doti,dotr,mpw_tmp,npw_k1,npw_k2,nspinor,pwind_tmp,&
1075 &       vect1,vect2)
1076        s1mat(1,jband,iband) = s1mat(1,jband,iband) + dotr
1077        s1mat(2,jband,iband) = s1mat(2,jband,iband) + doti
1078 
1079      end do    ! iband
1080    end do    !jband
1081 
1082    Amat(:,:,:)=zero
1083 
1084 !  calculate Amat
1085    do iband=1, dtefield%mband_occ
1086      do jband=1, dtefield%mband_occ
1087        do kband=1, dtefield%mband_occ
1088          Amat(1,iband,jband) = Amat(1,iband,jband) + s1mat(1,iband,kband)*&
1089 &         qmat(1,kband,jband,ikpt,2,idir)&
1090 &         - s1mat(2,iband,kband)*qmat(2,kband,jband,ikpt,2,idir)
1091          Amat(2,iband,jband) = Amat(2,iband,jband) + s1mat(1,iband,kband)*&
1092 &         qmat(2,kband,jband,ikpt,2,idir)&
1093 &         + s1mat(2,iband,kband)*qmat(1,kband,jband,ikpt,2,idir)
1094        end do
1095      end do
1096    end do
1097 
1098    Bmat(:,:,:)=zero
1099 
1100 !  calculate Bmat
1101    ikptn = dtefield%ikpt_dk(ikpt,7,idir)
1102    do iband=1, dtefield%mband_occ
1103      do jband=1, dtefield%mband_occ
1104        do kband=1, dtefield%mband_occ
1105          Bmat(1,jband,kband) = Bmat(1,jband,kband) + Amat(1,iband,kband)*&
1106 &         qmat(1,jband,iband,ikptn,2,idir)&
1107 &         - Amat(2,iband,kband)*qmat(2,jband,iband,ikptn,2,idir)
1108          Bmat(2,jband,kband) = Bmat(2,jband,kband) + Amat(1,iband,kband)*&
1109 &         qmat(2,jband,iband,ikptn,2,idir)&
1110          + Amat(2,iband,kband)*qmat(1,jband,iband,ikptn,2,idir)
1111        end do
1112      end do
1113    end do
1114 
1115 !  calc. the second term of gradient------------------------------
1116 
1117 !  preparation
1118 
1119    ikptnp1 = dtefield%ikpt_dk(ikpt,4,idir)
1120    icg = dtefield%cgindex(ikptnp1,isppol)
1121    npw_k1 = npwar1(ikpt)
1122    npw_k2 = npwarr(ikptnp1)
1123    pwind_tmp(1:npw_k1) =pwindall((ikpt-1)*mpw_tmp+1:(ikpt-1)*mpw_tmp+npw_k1,6,idir)
1124 
1125    z1(:) = zero
1126    z2(:) = zero
1127    do ipw = 1, npw_k1
1128 
1129      jpw = pwind_tmp(ipw)
1130 
1131      if (jpw > 0) then
1132 
1133        do iband = 1, dtefield%mband_occ
1134          wfr = cg(1,icg + (iband - 1)*npw_k2*nspinor + jpw)
1135          wfi = cg(2,icg + (iband - 1)*npw_k2*nspinor + jpw)
1136 
1137          do jband=1, dtefield%mband_occ
1138 
1139            grad_berry(1,ipw,jband) = grad_berry(1,ipw,jband) + fac*(Bmat(1,iband,jband)*wfr - Bmat(2,iband,jband)*wfi)
1140            grad_berry(2,ipw,jband) = grad_berry(2,ipw,jband) + fac*(Bmat(1,iband,jband)*wfi + Bmat(2,iband,jband)*wfr)
1141 
1142          end do
1143        end do
1144      end if
1145    end do
1146 
1147  end do !idir
1148 
1149  ABI_DEALLOCATE(vect1)
1150  ABI_DEALLOCATE(vect2)
1151  ABI_DEALLOCATE(s1mat)
1152  ABI_DEALLOCATE(Amat)
1153  ABI_DEALLOCATE(Bmat)
1154  ABI_DEALLOCATE(pwind_tmp)
1155 
1156 end subroutine dfptff_gradberry

ABINIT/dfptff_initberry [ Functions ]

[ Top ] [ Functions ]

NAME

 dfptff_initberry

FUNCTION

 Initialization of response calculations in finite electric field.

INPUTS

 dtset <type(dataset_type)> = all input variables in this dataset
 gmet(3,3) = reciprocal space metric tensor in bohr**-2
 kg(3,mpw*mkmem) = reduced (integer) coordinates of G vecs in basis sphere
 kg1(3,mpw1*mkmem) = reduced (integer) coordinates of G vecs for response wfs
 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
 nsppol = 1 for unpolarized, 2 for spin-polarized
 occ(mband*nkpt*nsppol) = occup number for each band at each k point
 rprimd(3,3) = dimensional primitive vectors

OUTPUT

 dtefield=variables related to response Berry-phase calculation
 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>

NOTES

      this duplicates in part initberry for the init of the dtefield - should be made
      into a common constructor in m_dtefield or somethin

PARENTS

      dfpt_scfcv

CHILDREN

      kpgio,wrtout

SOURCE

104 subroutine dfptff_initberry(dtefield,dtset,gmet,kg,kg1,mband,mkmem,mpi_enreg,&
105 &                mpw,mpw1,nkpt,npwarr,npwar1,nsppol,occ,pwindall,rprimd)
106 
107 
108 !This section has been created automatically by the script Abilint (TD).
109 !Do not modify the following lines by hand.
110 #undef ABI_FUNC
111 #define ABI_FUNC 'dfptff_initberry'
112 !End of the abilint section
113 
114  implicit none
115 
116 !Arguments ----------------------------------------
117 !scalars
118  integer,intent(in) :: mband,mkmem,mpw,mpw1,nkpt,nsppol
119  type(MPI_type),intent(inout) :: mpi_enreg
120  type(dataset_type),intent(in) :: dtset
121  type(efield_type),intent(inout) :: dtefield !vz_i needs efield2
122 !arrays
123  integer,intent(in) :: kg(3,mpw*mkmem),kg1(3,mpw1*mkmem),npwar1(nkpt)
124  integer,intent(in) :: npwarr(nkpt)
125  integer,intent(out) :: pwindall(max(mpw,mpw1)*mkmem,8,3)
126  real(dp),intent(in) :: gmet(3,3),occ(mband*nkpt*nsppol),rprimd(3,3)
127 
128 !Local variables ----------------------------------
129 !scalars
130  integer :: flag,iband,icg,idir,ifor,ikg,ikg1,ikpt,ikpt1,ikpt2,ikstr
131  integer :: index,ipw,isppol,istr,iunmark,jpw,nband_k,mband_occ_k,nkstr,npw_k
132  integer :: npw_k1,orig
133  real(dp) :: ecut_eff,occ_val
134  character(len=500) :: message
135 !arrays
136  integer :: dg(3)
137  integer,allocatable :: kg_tmp(:,:),kpt_mark(:),npwarr_tmp(:),npwtot(:)
138  real(dp) :: diffk(3),dk(3)
139  real(dp),allocatable :: kpt1(:,:)
140 
141 ! *************************************************************************
142 
143 !-----------------------------------------------------------
144 !---------------- initilize dtefield //---------------------
145 !-----------------------------------------------------------
146 
147 !Initialization of efield_type variables
148  dtefield%efield_dot(:) = zero
149  dtefield%dkvecs(:,:) = zero
150  dtefield%maxnstr = 0    ; dtefield%maxnkstr  = 0
151  dtefield%nstr(:) = 0    ; dtefield%nkstr(:) = 0
152  ABI_ALLOCATE(dtefield%ikpt_dk,(nkpt,9,3))
153  ABI_ALLOCATE(dtefield%cgindex,(nkpt,nsppol*2))
154  ABI_ALLOCATE(dtefield%kgindex,(nkpt))
155  dtefield%ikpt_dk(:,:,:) = 0
156  dtefield%cgindex(:,:) = 0
157  dtefield%mband_occ = 0
158  ABI_ALLOCATE(dtefield%nband_occ,(nsppol))
159  dtefield%nband_occ = 0
160  pwindall(:,:,:) = 0
161 
162 !Compute the number of occupied bands and check that--------
163 !it is the same for each k-point----------------------------
164 
165  occ_val = two/(dtset%nsppol*one)
166 
167  index = 0
168  do isppol = 1, nsppol
169    do ikpt = 1, nkpt
170 
171      mband_occ_k = 0
172      nband_k = dtset%nband(ikpt)
173 
174      do iband = 1, nband_k
175        index = index + 1
176        if (abs(occ(index) - occ_val) < tol8) mband_occ_k = mband_occ_k + 1
177      end do
178 
179      if (ikpt > 1) then
180        if (dtefield%nband_occ(isppol) /= mband_occ_k) then
181          message = ' The number of valence bands is not the same for every k-point for present spin'
182          MSG_ERROR(message)
183        end if
184      else
185        dtefield%mband_occ = max(dtefield%mband_occ,mband_occ_k)
186        dtefield%nband_occ(isppol) = mband_occ_k
187      end if
188 
189    end do
190  end do
191 
192 !DEBUG
193 !write(std_out,*)'dfptff_initberry:nkpt',nkpt
194 !ENDDEBUG
195 
196 !Compute the location of zero-order wavefunction --------------
197  icg = 0
198  do isppol = 1, nsppol
199    do ikpt = 1, nkpt
200 
201      dtefield%cgindex(ikpt,isppol) = icg
202      nband_k = dtset%nband(ikpt)
203      npw_k = npwarr(ikpt)
204      icg = icg + dtset%nspinor*npw_k*nband_k
205 
206    end do
207  end do
208 
209 !Compute the location of kg --------------
210  ikg = 0
211  do ikpt = 1, nkpt
212 
213    dtefield%kgindex(ikpt) = ikg
214    npw_k = npwarr(ikpt)
215    ikg = ikg + npw_k
216 
217  end do
218 
219 !Compute the location of first-order wavefunction -------------------
220  icg = 0
221  do isppol = 1, nsppol
222    do ikpt = 1, nkpt
223 
224      dtefield%cgindex(ikpt,isppol+nsppol) = icg
225      nband_k = dtset%nband(ikpt)
226      npw_k = npwar1(ikpt)
227      icg = icg + dtset%nspinor*npw_k*nband_k
228 
229    end do
230  end do
231 
232 !Compute the reciprocal lattice coordinates of the electric field-----
233 
234  dtefield%efield_dot(1) = dot_product(dtset%efield(:),rprimd(:,1))
235  dtefield%efield_dot(2) = dot_product(dtset%efield(:),rprimd(:,2))
236  dtefield%efield_dot(3) = dot_product(dtset%efield(:),rprimd(:,3))
237 
238  write(message,'(a,a,a,a,3(2x,f16.9),a)')ch10,&
239 & ' initberry: Reciprocal lattice coordinates of the electric field',ch10,&
240 & '  efield_dot(1:3) = ',dtefield%efield_dot(1:3),ch10
241  call wrtout(std_out,message,'COLL')
242 
243 !find the related k points to every k point in full BZ-----------------
244 
245 !loop over three reciprocal directions
246  do idir = 1, 3
247 
248    if (dtset%rfdir(idir) == 1) then
249 
250 !    Compute dk(:), the vector between a k-point and its nearest
251 !    neighbour along the direction idir
252 
253      dk(:) = zero
254      dk(idir) = 1000_dp   ! initialize with a big number
255      do ikpt = 2, nkpt
256        diffk(:) = abs(dtset%kptns(:,ikpt) - dtset%kptns(:,1))
257        if ((diffk(1) < dk(1)+tol8).and.(diffk(2) < dk(2)+tol8).and.&
258 &       (diffk(3) < dk(3)+tol8)) dk(:) = diffk(:)
259      end do
260      dtefield%dkvecs(:,idir) = dk(:)
261 
262 !    For each k point, find k_prim such that k_prim= k + dk mod(G)
263 !    where G is a vector of the reciprocal lattice
264 
265      do ikpt = 1, nkpt
266 
267 !      First: k + dk
268        do ikpt1 = 1, nkpt
269          diffk(:) = abs(dtset%kptns(:,ikpt1) - dtset%kptns(:,ikpt) - dk(:))
270          if(sum(abs(diffk(:) - nint(diffk(:)))) < 3*tol8) then
271            dtefield%ikpt_dk(ikpt,1,idir) = ikpt1
272            exit
273          end if
274        end do
275 
276 !      Second: k - dk
277        do ikpt1 = 1, nkpt
278          diffk(:) = abs(dtset%kptns(:,ikpt1) - dtset%kptns(:,ikpt) + dk(:))
279          if(sum(abs(diffk(:) - nint(diffk(:)))) < 3*tol8) then
280            dtefield%ikpt_dk(ikpt,2,idir) = ikpt1
281            exit
282          end if
283        end do
284 
285 !      new
286 !      3rd: k + (n+1)dk
287 
288        do ikpt1 = 1, nkpt
289          diffk(:) = abs(dtset%kptns(:,ikpt1) - dtset%kptns(:,ikpt) - dk(:) - dtset%qptn(:))
290          if(sum(abs(diffk(:) - nint(diffk(:)))) < 3*tol8) then
291            dtefield%ikpt_dk(ikpt,3,idir) = ikpt1
292            exit
293          end if
294        end do
295 
296 !      6th: k - (n+1)dk
297 
298        do ikpt1 = 1, nkpt
299          diffk(:) = abs(dtset%kptns(:,ikpt1) - dtset%kptns(:,ikpt) + dk(:) + dtset%qptn(:))
300          if(sum(abs(diffk(:) - nint(diffk(:)))) < 3*tol8) then
301            dtefield%ikpt_dk(ikpt,6,idir) = ikpt1
302            exit
303          end if
304        end do
305 
306 !      4th: k + (n-1)dk
307 
308        do ikpt1 = 1, nkpt
309          diffk(:) = abs(dtset%kptns(:,ikpt1) - dtset%kptns(:,ikpt) + dk(:) - dtset%qptn(:))
310          if(sum(abs(diffk(:) - nint(diffk(:)))) < 3*tol8) then
311            dtefield%ikpt_dk(ikpt,4,idir) = ikpt1
312            exit
313          end if
314        end do
315 
316 !      5th: k - (n-1)dk
317 
318        do ikpt1 = 1, nkpt
319          diffk(:) = abs(dtset%kptns(:,ikpt1) - dtset%kptns(:,ikpt) - dk(:) + dtset%qptn(:))
320          if(sum(abs(diffk(:) - nint(diffk(:)))) < 3*tol8) then
321            dtefield%ikpt_dk(ikpt,5,idir) = ikpt1
322            exit
323          end if
324        end do
325 
326 !      7th: k+n dk
327 
328        do ikpt1 = 1, nkpt
329          diffk(:) = abs(dtset%kptns(:,ikpt1) - dtset%kptns(:,ikpt) - dtset%qptn(:))
330          if(sum(abs(diffk(:) - nint(diffk(:)))) < 3*tol8) then
331            dtefield%ikpt_dk(ikpt,7,idir) = ikpt1
332            exit
333          end if
334        end do
335 
336 !      8th: k-n dk
337 
338        do ikpt1 = 1, nkpt
339          diffk(:) = abs(dtset%kptns(:,ikpt1) - dtset%kptns(:,ikpt) + dtset%qptn(:))
340          if(sum(abs(diffk(:) - nint(diffk(:)))) < 3*tol8) then
341            dtefield%ikpt_dk(ikpt,8,idir) = ikpt1
342            exit
343          end if
344        end do
345 
346 !      9th: -k
347 
348        do ikpt1 = 1, nkpt
349          diffk(:) = abs(dtset%kptns(:,ikpt1) + dtset%kptns(:,ikpt))
350          if(sum(diffk(:))  < 3*tol8) then
351            dtefield%ikpt_dk(ikpt,9,idir) = ikpt1
352            exit
353          end if
354        end do
355 
356      end do     ! ikpt
357 
358 
359 
360 !    Find the string length, starting from k point 1
361 !    (all strings must have the same number of points)
362 
363      nkstr = 1
364      ikpt1 = 1
365      do ikpt = 1, nkpt
366        ikpt1 = dtefield%ikpt_dk(ikpt1,1,idir)
367        if (ikpt1 == 1) exit
368        nkstr = nkstr + 1
369      end do
370 
371 !    Check that the string length is a divisor of nkpt
372      if(mod(nkpt,nkstr) /= 0) then
373        write(message,'(a,i0,a,i0)')' The string length = ',nkstr,', is not a divisor of nkpt =',nkpt
374        MSG_BUG(message)
375      end if
376      dtefield%nkstr(idir) = nkstr
377      dtefield%nstr(idir)  = nkpt/nkstr
378 
379    end if      ! dtset%rfdir(idir) == 1
380 
381    write(message,'(a,i1,a,i3,a,i3)')&
382 &   '  dfptff_initberry: for direction ',idir,', nkstr = ',dtefield%nkstr(idir),&
383 &   ', nstr = ',dtefield%nstr(idir)
384    call wrtout(std_out,message,'COLL')
385 
386  end do     ! close loop over idir
387 
388  dtefield%maxnstr  = maxval(dtefield%nstr(:))
389  dtefield%maxnkstr = maxval(dtefield%nkstr(:))
390  ABI_ALLOCATE(dtefield%idxkstr,(dtefield%maxnkstr,dtefield%maxnstr,3))
391  dtefield%idxkstr(:,:,:) = 0
392 
393 
394 !Build the different strings------------------------------------------
395 
396  ABI_ALLOCATE(kpt_mark,(nkpt))
397  do idir = 1, 3
398 
399    if (dtset%rfdir(idir) == 1) then
400 
401      iunmark = 1
402      kpt_mark(:) = 0
403      do istr = 1, dtefield%nstr(idir)
404 
405        do while(kpt_mark(iunmark) /= 0)
406          iunmark = iunmark + 1
407        end do
408        dtefield%idxkstr(1,istr,idir) = iunmark
409        kpt_mark(iunmark) = 1
410        do ikstr = 2, dtefield%nkstr(idir)
411          ikpt1 = dtefield%idxkstr(ikstr-1,istr,idir)
412          ikpt2 = dtefield%ikpt_dk(ikpt1,1,idir)
413          dtefield%idxkstr(ikstr,istr,idir) = ikpt2
414          kpt_mark(ikpt2) = 1
415        end do
416 
417      end do    ! istr
418 
419    end if         ! rfdir(idir) == 1
420 
421  end do           ! close loop over idir
422 
423  ABI_DEALLOCATE(kpt_mark)
424 
425 
426 !Build the array pwindall that is needed to compute the different overlap matrices
427 !at k +- dk
428 
429  ABI_ALLOCATE(kg_tmp,(3,max(mpw,mpw1)*mkmem))
430  ABI_ALLOCATE(kpt1,(3,nkpt))
431  ABI_ALLOCATE(npwarr_tmp,(nkpt))
432  ABI_ALLOCATE(npwtot,(nkpt))
433  ecut_eff = dtset%ecut*(dtset%dilatmx)**2
434 
435  do idir = 1, 3
436 
437    if (dtset%rfdir(idir) == 1) then
438 
439      dk(:) = dtefield%dkvecs(:,idir)
440 
441      do ifor = 1, 2
442 
443        if (ifor == 2) dk(:) = -1_dp*dk(:)
444 
445 !      Build the array kpt1 = kptns + dk
446 !      all k-poins of kpt1 must be in the same BZ as those of kptns
447        kpt1(:,:) = zero
448        do ikpt = 1, nkpt
449          ikpt1 = dtefield%ikpt_dk(ikpt,ifor,idir)
450          kpt1(:,ikpt) = dtset%kptns(:,ikpt1)
451        end do  ! close loop over ikpt
452 
453 !      Set up the basis sphere of plane waves at kpt1
454        kg_tmp(:,:) = 0
455        call kpgio(ecut_eff,dtset%exchn2n3d,gmet,dtset%istwfk,kg_tmp,&
456 &       kpt1,dtset%mkmem,dtset%nband,nkpt,'PERS',mpi_enreg,mpw,&
457 &       npwarr_tmp,npwtot,dtset%nsppol)
458 
459        ikg = 0 ; ikg1 = 0
460        do ikpt = 1, nkpt
461 
462          nband_k = dtset%nband(ikpt)
463          ikpt1 = dtefield%ikpt_dk(ikpt,ifor,idir)
464 
465 
466          dg(:) = nint(dtset%kptns(:,ikpt) + dk(:) + tol10 - &
467 &         dtset%kptns(:,ikpt1))
468 
469          flag = 0; orig = 1
470          if (dg(1)*dg(1) + dg(2)*dg(2) + dg(3)*dg(3) > 0) flag = 1
471 
472          ikpt1 = dtefield%ikpt_dk(ikpt,ifor,idir)
473          npw_k  = npwarr(ikpt)
474          npw_k1 = npwarr_tmp(ikpt)
475 
476          do ipw = 1, npw_k
477            do jpw = orig, npw_k1
478              if ((kg(1,ikg + ipw) == kg_tmp(1,ikg1 + jpw) - dg(1)).and. &
479              (kg(2,ikg + ipw) == kg_tmp(2,ikg1 + jpw) - dg(2)).and. &
480              (kg(3,ikg + ipw) == kg_tmp(3,ikg1 + jpw) - dg(3))) then
481                pwindall((ikpt-1)*max(mpw,mpw1) + ipw,ifor,idir) = jpw
482                if (flag == 0) orig = jpw
483                exit
484              end if
485            end do
486          end do
487 
488          ikg  = ikg + npw_k
489          ikg1 = ikg1 + npw_k1
490 
491        end do    ! close loop over ikpt
492 
493      end do    ! close loop over ifor
494 
495    end if      ! rfdir(idir) == 1
496 
497  end do        ! close loop over idir
498 
499 !------------------------------------------------------------------------------
500 !<u_q|u_q>
501 !at k +- dk
502 
503  do idir = 1, 3
504 
505    if (dtset%rfdir(idir) == 1) then
506 
507      dk(:) = dtefield%dkvecs(:,idir)
508 
509      do ifor = 1, 2
510 
511        if (ifor == 2) dk(:) = -1_dp*dk(:)
512 
513 !      Build the array kpt1 = kptns + qptn + dk
514 !      all k-poins of kpt1 must be in the same BZ as those of kptns
515        kpt1(:,:) = zero
516        do ikpt = 1, nkpt
517          ikpt1 = dtefield%ikpt_dk(ikpt,ifor,idir)
518          kpt1(:,ikpt) = dtset%kptns(:,ikpt1)+dtset%qptn(:)
519        end do  ! close loop over ikpt
520 
521 !      Set UP THE BASIS SPHERE OF PLANE waves at kpt1
522        kg_tmp(:,:) = 0
523        call kpgio(ecut_eff,dtset%exchn2n3d,gmet,dtset%istwfk,kg_tmp,&
524 &       kpt1,dtset%mkmem,dtset%nband,nkpt,'PERS',mpi_enreg,mpw1,&
525 &       npwarr_tmp,npwtot,dtset%nsppol)
526 
527 
528        ikg = 0 ; ikg1 = 0
529        do ikpt = 1, nkpt
530 
531          nband_k = dtset%nband(ikpt)
532          ikpt1 = dtefield%ikpt_dk(ikpt,ifor,idir)
533 
534          if (proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,1,mpi_enreg%me)) cycle
535 
536          dg(:) = nint(dtset%kptns(:,ikpt) + dk(:) + tol10 - &
537 &         dtset%kptns(:,ikpt1))
538 
539          flag = 0; orig = 1
540          if (dg(1)*dg(1) + dg(2)*dg(2) + dg(3)*dg(3) > 0) flag = 1
541 
542          ikpt1 = dtefield%ikpt_dk(ikpt,ifor,idir)
543          npw_k  = npwar1(ikpt)
544          npw_k1 = npwarr_tmp(ikpt)
545 
546          do ipw = 1, npw_k
547            do jpw = orig, npw_k1
548              if ((kg1(1,ikg + ipw) == kg_tmp(1,ikg1 + jpw) - dg(1)).and. &
549              (kg1(2,ikg + ipw) == kg_tmp(2,ikg1 + jpw) - dg(2)).and. &
550              (kg1(3,ikg + ipw) == kg_tmp(3,ikg1 + jpw) - dg(3))) then
551                pwindall((ikpt-1)*max(mpw,mpw1) +ipw,ifor+2,idir) = jpw
552                if (flag == 0) orig = jpw
553                exit
554              end if
555            end do
556          end do
557 
558          ikg  = ikg + npw_k
559          ikg1 = ikg1 + npw_k1
560 
561        end do    ! close loop over ikpt
562      end do    ! close loop over ifor
563    end if      ! rfdir(idir) == 1
564  end do        ! close loop over idir
565 
566 
567 !---------------------------------------------------------------------------
568 
569  do idir = 1, 3
570 
571    if (dtset%rfdir(idir) == 1) then
572 
573      dk(:) = dtset%qptn(:) + dtefield%dkvecs(:,idir)
574 
575      do ifor = 1, 2
576 
577        if (ifor == 2) dk(:) = dtset%qptn(:) - dtefield%dkvecs(:,idir)
578 
579 !      Build the array kpt1 = kptns + qptn + dk
580 !      all k-poins of kpt1 must be in the same BZ as those of kptns
581        kpt1(:,:) = zero
582        do ikpt = 1, nkpt
583          ikpt1 = dtefield%ikpt_dk(ikpt,ifor+2,idir)
584          kpt1(:,ikpt) = dtset%kptns(:,ikpt1)
585        end do  ! close loop over ikpt
586 
587 !      Set UP THE BASIS SPHERE OF PLANE waves at kpt1
588        kg_tmp(:,:) = 0
589        call kpgio(ecut_eff,dtset%exchn2n3d,gmet,dtset%istwfk,kg_tmp,&
590 &       kpt1,dtset%mkmem,dtset%nband,nkpt,'PERS',mpi_enreg,mpw,&
591 &       npwarr_tmp,npwtot,dtset%nsppol)
592 
593        ikg = 0 ; ikg1 = 0
594        do ikpt = 1, nkpt
595 
596          nband_k = dtset%nband(ikpt)
597          ikpt1 = dtefield%ikpt_dk(ikpt,ifor+2,idir)
598 
599          if (proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,1,mpi_enreg%me)) cycle
600 
601          dg(:) = nint(dtset%kptns(:,ikpt) + dk(:) + tol10 - &
602 &         dtset%kptns(:,ikpt1))
603 
604          flag = 0; orig = 1
605          if (dg(1)*dg(1) + dg(2)*dg(2) + dg(3)*dg(3) > 0) flag = 1
606 
607          ikpt1 = dtefield%ikpt_dk(ikpt,ifor+2,idir)
608          npw_k  = npwar1(ikpt)
609          npw_k1 = npwarr_tmp(ikpt)
610 
611          do ipw = 1, npw_k
612            do jpw = orig, npw_k1
613              if ((kg1(1,ikg + ipw) == kg_tmp(1,ikg1 + jpw) - dg(1)).and. &
614              (kg1(2,ikg + ipw) == kg_tmp(2,ikg1 + jpw) - dg(2)).and. &
615              (kg1(3,ikg + ipw) == kg_tmp(3,ikg1 + jpw) - dg(3))) then
616                pwindall((ikpt-1)*max(mpw,mpw1)+ipw,ifor+4,idir) = jpw
617                if (flag == 0) orig = jpw
618                exit
619              end if
620            end do
621          end do
622          ikg  = ikg + npw_k
623          ikg1 = ikg1 + npw_k1
624 
625        end do    ! close loop over ikpt
626      end do    ! close loop over ifor
627    end if      ! rfdir(idir) == 1
628  end do        ! close loop over idir===============================================================
629 
630 
631 !Build the array pwind3 that is needed to compute the overlap matrices===============================
632 
633  do idir = 1, 3
634 
635    if (dtset%rfdir(idir) == 1) then
636 
637      dk(:) = - dtset%qptn(:) + dtefield%dkvecs(:,idir)
638 
639      do ifor = 1, 2
640 
641        if (ifor == 2) dk(:) = - dtset%qptn(:) -  dtefield%dkvecs(:,idir)
642 
643 !      Build the array kpt1 = kptns + qptn + dk
644 !      all k-poins of kpt1 must be in the same BZ as those of kptns
645        kpt1(:,:) = zero
646        do ikpt = 1, nkpt
647          ikpt1 = dtefield%ikpt_dk(ikpt,ifor+4,idir)
648          kpt1(:,ikpt) = dtset%kptns(:,ikpt1)+ dtset%qptn(:)
649        end do  ! close loop over ikpt
650 
651 !      Set UP THE BASIS SPHERE OF PLANE waves at kpt1
652        kg_tmp(:,:) = 0
653        call kpgio(ecut_eff,dtset%exchn2n3d,gmet,dtset%istwfk,kg_tmp,&
654 &       kpt1,dtset%mkmem,dtset%nband,nkpt,'PERS',mpi_enreg,mpw1,&
655 &       npwarr_tmp,npwtot,dtset%nsppol)
656 
657        ikg = 0 ; ikg1 = 0
658        do ikpt = 1, nkpt
659 
660          nband_k = dtset%nband(ikpt)
661          ikpt1 = dtefield%ikpt_dk(ikpt,ifor+4,idir)
662 
663          if (proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,1,mpi_enreg%me)) cycle
664 
665          dg(:) = nint(dtset%kptns(:,ikpt) + dk(:) + tol10 - &
666 &         dtset%kptns(:,ikpt1))
667 
668          flag = 0; orig = 1
669          if (dg(1)*dg(1) + dg(2)*dg(2) + dg(3)*dg(3) > 0) flag = 1
670 
671          ikpt1 = dtefield%ikpt_dk(ikpt,ifor+4,idir)
672          npw_k  = npwarr(ikpt)
673          npw_k1 = npwarr_tmp(ikpt)
674 
675          do ipw = 1, npw_k
676            do jpw = orig, npw_k1
677              if ((kg(1,ikg + ipw) == kg_tmp(1,ikg1 + jpw) - dg(1)).and. &
678              (kg(2,ikg + ipw) == kg_tmp(2,ikg1 + jpw) - dg(2)).and. &
679              (kg(3,ikg + ipw) == kg_tmp(3,ikg1 + jpw) - dg(3))) then
680                pwindall((ikpt-1)*max(mpw,mpw1) + ipw,ifor+6,idir) = jpw
681                if (flag == 0) orig = jpw
682                exit
683              end if
684            end do
685          end do
686          ikg  = ikg + npw_k
687          ikg1 = ikg1 + npw_k1
688 
689        end do    ! close loop over ikpt
690      end do    ! close loop over ifor
691    end if      ! rfdir(idir) == 1
692  end do        ! close loop over idir====================================================================
693 
694  ABI_DEALLOCATE(kg_tmp)
695  ABI_DEALLOCATE(kpt1)
696  ABI_DEALLOCATE(npwarr_tmp)
697  ABI_DEALLOCATE(npwtot)
698 
699 end subroutine dfptff_initberry

ABINIT/m_dfpt_fef [ Modules ]

[ Top ] [ Modules ]

NAME

  m_dfpt_fef

FUNCTION

  Response calculations 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 .

PARENTS

CHILDREN

SOURCE

21 #if defined HAVE_CONFIG_H
22 #include "config.h"
23 #endif
24 
25 #include "abi_common.h"
26 
27 module m_dfpt_fef
28 
29  use defs_basis
30  use defs_abitypes
31  use m_abicore
32  use m_errors
33  use m_efield
34 
35  use m_kg,        only : kpgio
36  use m_cgtools,   only : overlap_g
37  use m_mpinfo,    only : proc_distrb_cycle
38 
39  implicit none
40 
41  private

ABINIT/qmatrix [ Functions ]

[ Top ] [ Functions ]

NAME

 qmatrix

FUNCTION

 calculation of the inverse of the overlap matrix

INPUTS

 cg(2,mpw*nspinor*mband*mkmem*nsppol) = planewave coefficients of wavefunctions
 RF wavefunctions at k,q.
 dtefield = variables related to response Berry-phase calculation
 ikpt = the index of the current k point
 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
 pwindall(max(mpw,mpw1)*mkmem,8,3) = array used to compute the overlap matrices

OUTPUT

 qmat(2,dtefield%nband_occ,dtefield%nband_occ,nkpt,2,3) = inverse of the overlap matrix

PARENTS

      dfpt_scfcv

CHILDREN

      dzgedi,dzgefa,overlap_g

SOURCE

2639 subroutine qmatrix(cg,dtefield,qmat,mpw,mpw1,mkmem,mband,npwarr,nkpt,nspinor,nsppol,pwindall)
2640 
2641  use m_hide_lapack, only : dzgedi, dzgefa
2642 
2643 !This section has been created automatically by the script Abilint (TD).
2644 !Do not modify the following lines by hand.
2645 #undef ABI_FUNC
2646 #define ABI_FUNC 'qmatrix'
2647 !End of the abilint section
2648 
2649  implicit none
2650 
2651 !Arguments ----------------------------------------
2652 !scalars
2653  integer,intent(in) :: mband,mkmem,mpw,mpw1,nkpt,nspinor,nsppol
2654  type(efield_type),intent(in) :: dtefield
2655 !arrays
2656  integer,intent(in) :: npwarr(nkpt),pwindall(max(mpw,mpw1)*mkmem,8,3)
2657  real(dp),intent(in) :: cg(2,mpw*nspinor*mband*mkmem*nsppol)
2658  real(dp),intent(out) :: qmat(2,dtefield%mband_occ,dtefield%mband_occ,nkpt,2,3)
2659 
2660 !Local variables -------------------------
2661 !scalars
2662  integer :: iband,icg,icg1,idir,ifor,ikpt,ikpt2,info,jband,job
2663  integer :: npw_k1,npw_k2,pwmax,pwmin
2664  integer :: isppol
2665  real(dp) :: doti,dotr
2666 !arrays
2667  integer,allocatable :: ipvt(:),pwind_k(:)
2668  real(dp) :: det(2,2)
2669  real(dp),allocatable :: sinv(:,:,:),smat_k(:,:,:),vect1(:,:),vect2(:,:)
2670  real(dp),allocatable :: zgwork(:,:)
2671 
2672 ! *************************************************************************
2673 
2674  ABI_ALLOCATE(ipvt,(dtefield%mband_occ))
2675  ABI_ALLOCATE(sinv,(2,dtefield%mband_occ,dtefield%mband_occ))
2676  ABI_ALLOCATE(zgwork,(2,dtefield%mband_occ))
2677  ABI_ALLOCATE(vect1,(2,0:mpw))
2678  ABI_ALLOCATE(vect2,(2,0:mpw))
2679  ABI_ALLOCATE(smat_k,(2,dtefield%mband_occ,dtefield%mband_occ))
2680  ABI_ALLOCATE(pwind_k,(max(mpw,mpw1)))
2681  vect1(:,0) = zero ; vect2(:,0) = zero
2682 
2683  job = 11
2684 
2685 !loop over k points
2686  do isppol = 1, nsppol
2687    do ikpt = 1, nkpt
2688      npw_k1 = npwarr(ikpt)
2689      icg  = dtefield%cgindex(ikpt,1)
2690      do idir = 1, 3
2691        do ifor = 1, 2
2692 
2693          ikpt2 = dtefield%ikpt_dk(ikpt,ifor,idir)
2694          npw_k2 = npwarr(ikpt2)
2695          icg1 = dtefield%cgindex(ikpt2,1)
2696          pwind_k(1:npw_k1) = pwindall((ikpt-1)*max(mpw,mpw1)+1:(ikpt-1)*max(mpw,mpw1)+npw_k1,ifor,idir)
2697 
2698          do jband = 1, dtefield%nband_occ(isppol)
2699            vect2(:,1:npw_k2) = cg(:,icg1 + 1 + (jband-1)*npw_k2*nspinor:icg1 + jband*npw_k2*nspinor)
2700            if (npw_k2 < mpw) vect2(:,npw_k2+1:mpw) = zero
2701 
2702            do iband = 1, dtefield%nband_occ(isppol)
2703 
2704              pwmin = (iband-1)*npw_k1*nspinor
2705              pwmax = pwmin + npw_k1*nspinor
2706              vect1(:,1:npw_k1) = cg(:,icg + 1 + pwmin:icg + pwmax)
2707              if (npw_k1 < mpw) vect1(:,npw_k1+1:mpw) = zero
2708              call overlap_g(doti,dotr,mpw,npw_k1,npw_k2,nspinor,pwind_k,vect1,vect2)
2709              smat_k(1,iband,jband) = dotr
2710              smat_k(2,iband,jband) = doti
2711            end do    ! iband
2712          end do    !jband
2713 
2714          sinv(:,:,:) = smat_k(:,:,:)
2715 
2716          call dzgefa(sinv,dtefield%mband_occ,dtefield%nband_occ(isppol),ipvt,info)
2717          call dzgedi(sinv,dtefield%mband_occ,dtefield%nband_occ(isppol),ipvt,det,zgwork,job)
2718 
2719          qmat(:,:,:,ikpt,ifor,idir) = sinv(:,:,:)
2720        end do
2721      end do
2722    end do  !end loop over k
2723  end do
2724 
2725  ABI_DEALLOCATE(ipvt)
2726  ABI_DEALLOCATE(sinv)
2727  ABI_DEALLOCATE(zgwork)
2728  ABI_DEALLOCATE(vect1)
2729  ABI_DEALLOCATE(vect2)
2730  ABI_DEALLOCATE(smat_k)
2731  ABI_DEALLOCATE(pwind_k)
2732 
2733 end subroutine qmatrix