TABLE OF CONTENTS


ABINIT/dfptff_die [ Functions ]

[ Top ] [ Functions ]

NAME

 dfptff_die

FUNCTION

 calculate electric susceptibility tensor in Eq.(28) in PRB 75, 115116(2007).

COPYRIGHT

 Copyright (C) 2004-2018 ABINIT group (XW).
 This file is distributed under the terms of the
 GNU General Public License, see ~abinit/COPYING
 or http://www.gnu.org/copyleft/gpl.txt .
 For the initials of contributors, see ~abinit/doc/developers/contributors .

INPUTS

 cg(2,mpw*nspinor*mband*mkmem*nsppol) = planewave coefficients of wavefunctions
 cg1(2,mpw1*nspinor*mband*mk1mem*nsppol) = pw coefficients of
 RF wavefunctions at k,q.
 dtefield = variables related to response Berry-phase calculation
 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

 55 #if defined HAVE_CONFIG_H
 56 #include "config.h"
 57 #endif
 58 
 59 #include "abi_common.h"
 60 
 61 subroutine dfptff_die(cg,cg1,dtefield,d2lo,idirpert,ipert,mband,mkmem,&
 62 &               mpw,mpw1,mpert,nkpt,npwarr,npwar1,nsppol,nspinor,pwindall,qmat,rprimd)
 63 
 64  use defs_basis
 65  use m_profiling_abi
 66  use m_efield
 67 
 68  use m_cgtools,   only : overlap_g
 69 
 70 !This section has been created automatically by the script Abilint (TD).
 71 !Do not modify the following lines by hand.
 72 #undef ABI_FUNC
 73 #define ABI_FUNC 'dfptff_die'
 74 !End of the abilint section
 75 
 76  implicit none
 77 
 78 !Arguments ----------------------------------------
 79 !scalars
 80  integer,intent(in) :: idirpert,ipert,mband,mkmem,mpert,mpw,mpw1,nkpt,nspinor
 81  integer,intent(in) :: nsppol
 82  type(efield_type),intent(in) :: dtefield
 83 !arrays
 84  integer,intent(in) :: npwar1(nkpt),npwarr(nkpt)
 85  integer,intent(in) :: pwindall(max(mpw,mpw1)*mkmem,8,3)
 86  real(dp),intent(in) :: cg(2,mpw*mband*mkmem*nspinor*nsppol)
 87  real(dp),intent(in) :: cg1(2,mpw1*mband*mkmem*nspinor*nsppol)
 88  real(dp),intent(in) :: qmat(2,dtefield%mband_occ,dtefield%mband_occ,nkpt,2,3)
 89  real(dp),intent(in) :: rprimd(3,3)
 90  real(dp),intent(inout) :: d2lo(2,3,mpert,3,mpert) !vz_i
 91 
 92 !Local variables ----------------------------------
 93 !scalars
 94  integer :: ialpha,iband,icg,icg1,idir,ikpt,ikpt1,jband,mpw_tmp,npw_k1
 95  integer :: npw_k2,pwmax,pwmin
 96  real(dp) :: doti,dotr,e0,fac
 97 !arrays
 98  integer,allocatable :: pwind_tmp(:)
 99  real(dp) :: edir(3)
100  real(dp),allocatable :: s1mat(:,:,:),vect1(:,:),vect2(:,:)
101 
102 ! *************************************************************************
103 
104 !calculate s1 matrices -----------------------------
105  mpw_tmp=max(mpw,mpw1)
106  ABI_ALLOCATE(s1mat,(2,dtefield%mband_occ,dtefield%mband_occ))
107  ABI_ALLOCATE(vect1,(2,0:mpw_tmp))
108  ABI_ALLOCATE(vect2,(2,0:mpw_tmp))
109  ABI_ALLOCATE(pwind_tmp,(mpw_tmp))
110  vect1(:,0) = zero ; vect2(:,0) = zero
111 
112  edir(:)=zero
113 
114  do ikpt=1,nkpt
115    do idir=1,3
116 !    compute <u^(0)_{k_j}|u^(1)_{k_j+1,q}> matrix--- q=0 ----------------------------------------
117 
118 !    prepare to calculate overlap matrix
119      ikpt1 = dtefield%ikpt_dk(ikpt,1,idir)
120      icg = dtefield%cgindex(ikpt,1)
121      icg1 = dtefield%cgindex(ikpt1,1+nsppol)
122      npw_k1 = npwarr(ikpt)
123      npw_k2 = npwar1(ikpt1)
124      pwind_tmp(1:npw_k1) = pwindall((ikpt-1)*mpw_tmp+1:(ikpt-1)*mpw_tmp+npw_k1,7,idir)
125 
126      vect1(:,0) = zero ; vect2(:,0) = zero
127      do jband = 1, dtefield%mband_occ
128        vect2(:,1:npw_k2) = &
129 &       cg1(:,icg1 + 1 + (jband-1)*npw_k2*nspinor:icg1 + jband*npw_k2*nspinor)
130        if (npw_k2 < mpw_tmp) vect2(:,npw_k2+1:mpw_tmp) = zero
131        do iband = 1, dtefield%mband_occ
132          pwmin = (iband-1)*npw_k1*nspinor
133          pwmax = pwmin + npw_k1*nspinor
134          vect1(:,1:npw_k1) = &
135 &         cg(:,icg + 1 + pwmin:icg + pwmax)
136          if (npw_k1 < mpw_tmp) vect1(:,npw_k1+1:mpw_tmp) = zero
137          call overlap_g(doti,dotr,mpw_tmp,npw_k1,npw_k2,nspinor,pwind_tmp,&
138 &         vect1,vect2)
139          s1mat(1,iband,jband) = dotr
140          s1mat(2,iband,jband) = doti
141        end do    ! iband
142      end do    !jband
143 
144 !    compute <u^(1)_{k_j,q}|u^(0)_{k_j+1}> matrix-- q=0 -------------------------------------
145 
146 !    prepare to calculate overlap matrix
147      ikpt1 = dtefield%ikpt_dk(ikpt,1,idir)
148      icg = dtefield%cgindex(ikpt,1+nsppol)
149      icg1 = dtefield%cgindex(ikpt1,1)
150      npw_k1 = npwar1(ikpt)
151      npw_k2 = npwarr(ikpt1)
152      pwind_tmp(1:npw_k1) = pwindall((ikpt-1)*mpw_tmp+1:(ikpt-1)*mpw_tmp+npw_k1,5,idir)
153      vect1(:,0) = zero ; vect2(:,0) = zero
154      do jband = 1, dtefield%mband_occ
155        vect2(:,1:npw_k2) = &
156 &       cg(:,icg1 + 1 + (jband-1)*npw_k2*nspinor:icg1 + jband*npw_k2*nspinor)
157        if (npw_k2 < mpw_tmp) vect2(:,npw_k2+1:mpw_tmp) = zero
158        do iband = 1, dtefield%mband_occ
159          pwmin = (iband-1)*npw_k1*nspinor
160          pwmax = pwmin + npw_k1*nspinor
161          vect1(:,1:npw_k1) = &
162 &         cg1(:,icg + 1 + pwmin:icg + pwmax)
163          if (npw_k1 < mpw_tmp) vect1(:,npw_k1+1:mpw_tmp) = zero
164          call overlap_g(doti,dotr,mpw_tmp,npw_k1,npw_k2,nspinor,pwind_tmp,&
165 &         vect1,vect2)
166          s1mat(1,iband,jband) = s1mat(1,iband,jband) + dotr
167          s1mat(2,iband,jband) = s1mat(2,iband,jband) + doti
168        end do    ! iband
169      end do    !jband
170 
171 !    sum over the whole------------------------------------------------------------
172 
173      e0=zero
174 
175      do iband=1,dtefield%mband_occ
176        do jband=1,dtefield%mband_occ
177          e0 = e0 + (s1mat(1,iband,jband)*qmat(2,jband,iband,ikpt,1,idir)&
178 &         +    s1mat(2,iband,jband)*qmat(1,jband,iband,ikpt,1,idir))
179 
180        end do
181      end do
182 
183      do ialpha=1,3
184        fac = rprimd(ialpha,idir)/&
185 &       (dble(dtefield%nstr(idir))*pi)
186        edir(ialpha)=edir(ialpha)+ e0*fac
187      end do
188 
189    end do !idir
190  end do !ikpt
191 
192  d2lo(1,1:3,ipert,idirpert,ipert)=edir(:)
193 
194  ABI_DEALLOCATE(s1mat)
195  ABI_DEALLOCATE(vect1)
196  ABI_DEALLOCATE(vect2)
197  ABI_DEALLOCATE(pwind_tmp)
198 end subroutine dfptff_die