TABLE OF CONTENTS


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).

COPYRIGHT

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

INPUTS

 cg(2,mpw*nspinor*mband*mkmem*nsppol) = planewave coefficients of wavefunctions
 cg1(2,mpw1*nspinor*mband*mk1mem*nsppol) = pw coefficients of
 RF wavefunctions at k,q.
 dtefield = variables related to response Berry-phase calculation
 ikpt = the index of the current k point
 isppol = the index of the spin component
 mband =  maximum number of bands
 mkmem = maximum number of k-points in core memory
 mpw = maximum number of plane waves
 mpw1 = maximum number of plane waves for response wavefunctions
 nkpt = number of k points
 npwarr(nkpt) = number of planewaves in basis and boundary at this k point
 npwar1(nkpt) = number of planewaves in basis and boundary for response wfs
 nspinor = 1 for scalar wfs, 2 for spinor wfs
 nsppol = 1 for unpolarized, 2 for spin-polarized
 qmat(2,dtefield%mband_occ,dtefield%mband_occ,nkpt,2,3) =
 inverse of the overlap matrix
 pwindall(max(mpw,mpw1)*mkmem,8,3) = array used to compute the overlap matrices
 pwindall(:,1,:) <- <u^(0)_i|u^(0)_i+1>
 pwindall(:,2,:) <- <u^(0)_i|u^(0)_i-1>
 pwindall(:,3,:) <- <u^(1)_i|u^(1)_i+1>
 pwindall(:,4,:) <- <u^(1)_i|u^(1)_i-1>
 pwindall(:,5,:) <- <u^(1)_i|u^(0)_i+n+1>
 pwindall(:,6,:) <- <u^(1)_i|u^(0)_i+n-1>
 pwindall(:,7,:) <- <u^(0)_i|u^(1)_i-n+1>
 pwindall(:,8,:) <- <u^(0)_i|u^(1)_i-n-1>

OUTPUT

 eberry = the energy of the Berry phase term

PARENTS

      dfpt_scfcv

CHILDREN

      overlap_g

SOURCE

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