TABLE OF CONTENTS


ABINIT/pawuenergy [ Functions ]

[ Top ] [ Functions ]

NAME

 pawuenergy

FUNCTION

 Compute contributions to energy for PAW+U calculations

COPYRIGHT

 Copyright (C) 1998-2018 ABINIT group (BA,FJ,MT)
 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/Infos/contributors.

INPUTS

  iatom=index of current atom (note: this is the absolute index, not the index on current proc)
  pawprtvol=control print volume and debugging output for PAW
  pawtab <type(pawtab_type)>=paw tabulated starting data:
     %lpawu=l used for lda+u
     %vee(2*lpawu+1*4)=screened coulomb matrix
  paw_ij <type(paw_ij_type)>=paw arrays given on (i,j) channels
     %noccmmp(2*lpawu+1,2*lpawu+1,nspden)=density matrix in the PAW augm. region
     %nocctot(nspden)=number of electrons in the correlated subspace
  dmft_dc,e_ee,e_dc,e_dcdc,u_dmft,j_dmft= optional arguments for DMFT

OUTPUT

  eldaumdc= PAW+U contribution to total energy
  eldaumdcdc= PAW+U contribution to double-counting total energy

PARENTS

      m_energy,pawdenpot

CHILDREN

      wrtout

SOURCE

 40 #if defined HAVE_CONFIG_H
 41 #include "config.h"
 42 #endif
 43 
 44 #include "abi_common.h"
 45 
 46  subroutine pawuenergy(iatom,eldaumdc,eldaumdcdc,pawprtvol,pawtab,paw_ij,&
 47  &                     dmft_dc,e_ee,e_dc,e_dcdc,u_dmft,j_dmft) ! optional arguments (DMFT)
 48 
 49  use defs_basis
 50  use m_profiling_abi
 51  use m_errors
 52 
 53  use m_pawtab, only : pawtab_type
 54  use m_paw_ij, only : paw_ij_type
 55 
 56 !This section has been created automatically by the script Abilint (TD).
 57 !Do not modify the following lines by hand.
 58 #undef ABI_FUNC
 59 #define ABI_FUNC 'pawuenergy'
 60  use interfaces_14_hidewrite
 61 !End of the abilint section
 62 
 63  implicit none
 64 
 65 !Arguments ---------------------------------------------
 66 !scalars
 67  integer,intent(in) :: iatom,pawprtvol
 68  integer,optional,intent(in) :: dmft_dc
 69  real(dp),intent(inout) :: eldaumdc,eldaumdcdc
 70  real(dp),optional,intent(inout) :: e_ee,e_dc,e_dcdc
 71  real(dp),optional,intent(in) :: j_dmft,u_dmft
 72  type(paw_ij_type),intent(in) :: paw_ij
 73  type(pawtab_type),intent(in) :: pawtab
 74 
 75 !Local variables ---------------------------------------
 76 !scalars
 77 !Option for interaction energy in case of non-collinear magnetism:
 78 !           1: E_int=-U/4.N.(N-2)
 79 !           2: E_int=-U/2.(Nup.(Nup-1)+Ndn.(Ndn-1))
 80  integer,parameter :: option_interaction=1
 81 
 82  integer :: cplex_dij,dmftdc,ispden,jspden,lpawu,m1,m11,m2,m21,m3,m31,m4,m41
 83  real(dp) :: eks_opt3,edcdc_opt3,edcdctemp,edctemp,eldautemp,jpawu,jpawu_dc,mnorm,mx,my,mz
 84  real(dp) :: n_sig,n_sigs,n_msig,n_msigs,n_dndn,n_tot,n_upup
 85  real(dp) :: n12_ud_im,n12_du_im
 86  real(dp) :: n12_ud_re,n12_du_re
 87  real(dp) :: n34_ud_im,n34_du_im
 88  real(dp) :: n34_ud_re,n34_du_re
 89  real(dp) :: upawu
 90  real(dp),allocatable :: n12_sig(:),n34_msig(:),n34_sig(:)
 91  character(len=500) :: message
 92 
 93 ! *****************************************************
 94 
 95  if(present(dmft_dc))  then
 96    dmftdc=dmft_dc
 97    if(pawtab%usepawu<10) then
 98      write(message,'(4x,2a,i5)') "Error, usepawu should be equal to 10 if", &
 99 &     " dmft_dc is an argument of pawuenergy",pawtab%usepawu
100      call wrtout(std_out,message,'COLL')
101    end if
102  else 
103    dmftdc=0
104  end if
105 
106  DBG_ENTER("COLL")
107 
108  lpawu=pawtab%lpawu
109  upawu=pawtab%upawu;if(present(u_dmft)) upawu=u_dmft
110  jpawu=pawtab%jpawu;if(present(j_dmft)) jpawu=j_dmft
111 
112 !======================================================
113 !Compute LDA+U Energy
114 !-----------------------------------------------------
115 
116  eldautemp=zero
117  edcdc_opt3=zero
118  eks_opt3=zero
119  cplex_dij=paw_ij%cplex_dij
120  ABI_ALLOCATE(n12_sig,(cplex_dij))
121  ABI_ALLOCATE(n34_msig,(cplex_dij))
122  ABI_ALLOCATE(n34_sig,(cplex_dij))
123  do ispden=1,min(paw_ij%ndij,2)
124    jspden=min(paw_ij%ndij,2)-ispden+1
125 !  compute n_sigs and n_msigs for pawtab%usepawu=3
126    if (paw_ij%nspden<=2) then
127      n_sig =paw_ij%nocctot(ispden)
128      n_msig=paw_ij%nocctot(jspden)
129      n_tot=n_sig+n_msig
130    else
131      n_tot=paw_ij%nocctot(1)
132      mx=paw_ij%nocctot(2)
133      my=paw_ij%nocctot(3)
134      mz=paw_ij%nocctot(4)
135      mnorm=sqrt(mx*mx+my*my+mz*mz)
136      if (ispden==1) then
137 !      n_sig =half*(n_tot+mnorm)
138 !      n_msig=half*(n_tot-mnorm)
139        n_sig =half*(n_tot+sign(mnorm,mz))
140        n_msig=half*(n_tot-sign(mnorm,mz))
141      else
142 !      n_sig =half*(n_tot-mnorm)
143 !      n_msig=half*(n_tot+mnorm)
144        n_sig =half*(n_tot-sign(mnorm,mz))
145        n_msig=half*(n_tot+sign(mnorm,mz))
146      end if
147    end if
148    n_sigs =n_sig/(float(2*lpawu+1))
149    n_msigs =n_msig/(float(2*lpawu+1))
150 !  if(pawtab%usepawu==3) then
151 !  write(message,fmt=12) "noccmmp11 ",ispden,paw_ij%noccmmp(1,1,1,ispden)
152 !  call wrtout(std_out,message,'COLL')
153 !  write(message,fmt=12) "noccmmp11 ",jspden,paw_ij%noccmmp(1,1,1,jspden)
154 !  call wrtout(std_out,message,'COLL')
155 !  write(message,fmt=12) "n_sig      ",ispden,n_sig
156 !  call wrtout(std_out,message,'COLL')
157 !  write(message,fmt=12) "n_msig     ",jspden,n_msig
158 !  call wrtout(std_out,message,'COLL')
159 !  write(message,fmt=12) "n_sigs     ",ispden,n_sigs
160 !  call wrtout(std_out,message,'COLL')
161 !  write(message,fmt=12) "n_msigs    ",jspden,n_msigs
162 !  call wrtout(std_out,message,'COLL')
163 !  endif
164 !  12 format(a,i4,e20.10)
165 !  compute interaction energy E_{ee}
166    do m1=-lpawu,lpawu
167      m11=m1+lpawu+1
168      do m2=-lpawu,lpawu
169        m21=m2+lpawu+1
170        n12_sig(:)=paw_ij%noccmmp(:,m11,m21,ispden)
171        if(m21==m11.and.(pawtab%usepawu==3.or.dmftdc==3)) n12_sig(1)=n12_sig(1)-n_sigs
172        do m3=-lpawu,lpawu
173          m31=m3+lpawu+1
174          do m4=-lpawu,lpawu
175            m41=m4+lpawu+1
176            n34_sig(:) =paw_ij%noccmmp(:,m31,m41,ispden)
177            n34_msig(:)=paw_ij%noccmmp(:,m31,m41,jspden)
178            if(m31==m41.and.(pawtab%usepawu==3.or.dmftdc==3)) then
179              n34_sig(1)= n34_sig(1) - n_sigs
180              n34_msig(1)= n34_msig(1) - n_msigs
181            end if
182            eldautemp=eldautemp &
183 &           + n12_sig(1)*n34_msig(1)*pawtab%vee(m11,m31,m21,m41) &
184 &           + n12_sig(1)*n34_sig(1) *(pawtab%vee(m11,m31,m21,m41)-pawtab%vee(m11,m31,m41,m21))
185            if(cplex_dij==2) then
186              eldautemp=eldautemp &
187 &             - n12_sig(2)*n34_msig(2)*pawtab%vee(m11,m31,m21,m41) &
188 &             - n12_sig(2)*n34_sig(2) *(pawtab%vee(m11,m31,m21,m41)-pawtab%vee(m11,m31,m41,m21))
189            end if
190            if (pawtab%usepawu==3.or.dmftdc==3) then
191              edcdc_opt3=edcdc_opt3 &
192 &             + n_sigs*n34_msig(1)*pawtab%vee(m11,m31,m21,m41) &
193 &             + n_sigs*n34_sig(1) *(pawtab%vee(m11,m31,m21,m41)-pawtab%vee(m11,m31,m41,m21))
194              eks_opt3=eks_opt3 &
195 &             + paw_ij%noccmmp(1,m11,m21,ispden)*n34_msig(1)*pawtab%vee(m11,m31,m21,m41) &
196 &             + paw_ij%noccmmp(1,m11,m21,ispden)*n34_sig(1) *(pawtab%vee(m11,m31,m21,m41)-pawtab%vee(m11,m31,m41,m21))
197              if(cplex_dij==2) then
198                eks_opt3=eks_opt3 &
199 &               - paw_ij%noccmmp(2,m11,m21,ispden)*n34_msig(2)*pawtab%vee(m11,m31,m21,m41) &
200 &               - paw_ij%noccmmp(2,m11,m21,ispden)*n34_sig(2) *(pawtab%vee(m11,m31,m21,m41)-pawtab%vee(m11,m31,m41,m21))
201              end if
202            end if
203          end do ! m4
204        end do ! m3
205      end do ! m2
206    end do ! m1
207  end do ! ispden
208  if (paw_ij%ndij==1) eldautemp=two*eldautemp ! Non-magn. system: sum up and dn energies
209  ABI_DEALLOCATE(n12_sig)
210  ABI_DEALLOCATE(n34_msig)
211  ABI_DEALLOCATE(n34_sig)
212 
213 !Non-collinear magnetism: add non-diagonal term; see (Eq 3) in PRB 72, 024458 (2005)
214  if (paw_ij%ndij==4) then
215    do m1=-lpawu,lpawu
216      m11=m1+lpawu+1
217      do m2=-lpawu,lpawu
218        m21=m2+lpawu+1
219        n12_ud_re=paw_ij%noccmmp(1,m11,m21,3) ! updn
220        n12_ud_im=paw_ij%noccmmp(2,m11,m21,3) ! updn
221        n12_du_re=paw_ij%noccmmp(1,m11,m21,4) ! dnup
222        n12_du_im=paw_ij%noccmmp(2,m11,m21,4) ! dnup
223        do m3=-lpawu,lpawu
224          m31=m3+lpawu+1
225          do m4=-lpawu,lpawu
226            m41=m4+lpawu+1
227            n34_ud_re=paw_ij%noccmmp(1,m31,m41,3)  ! updn
228            n34_ud_im=paw_ij%noccmmp(2,m31,m41,3)  ! updn
229            n34_du_re=paw_ij%noccmmp(1,m31,m41,4)  ! dnup
230            n34_du_im=paw_ij%noccmmp(2,m31,m41,4)  ! dnup
231            eldautemp=eldautemp-pawtab%vee(m11,m31,m41,m21) &
232 &           *(n12_ud_re*n34_du_re-n12_ud_im*n34_du_im &
233 &           +n12_du_re*n34_ud_re-n12_du_im*n34_ud_im)
234            if (pawtab%usepawu==3.or.dmftdc==3) then
235              eks_opt3=eks_opt3-pawtab%vee(m11,m31,m41,m21) &
236 &             *(n12_ud_re*n34_du_re-n12_ud_im*n34_du_im &
237 &             +n12_du_re*n34_ud_re-n12_du_im*n34_ud_im)
238            end if
239          end do ! m4
240        end do ! m3
241      end do ! m2
242    end do ! m1
243  end if
244 
245 !Divide eldautemp by 2; see (Eq 1) in PRB 77, 155104 (2008)
246  eldautemp=half*eldautemp
247 
248 !if (paw_ij%ndij==1) then
249 !n_tot=two*paw_ij%nocctot(1)
250 !n_upup=paw_ij%nocctot(1)
251 !n_dndn=paw_ij%nocctot(1)
252 !else if (paw_ij%ndij==2) then
253 !n_tot=paw_ij%nocctot(1)+paw_ij%nocctot(2)
254 !n_upup=paw_ij%nocctot(1)
255 !n_dndn=paw_ij%nocctot(2)
256 !else if (paw_ij%ndij==4) then
257 !n_tot=paw_ij%nocctot(1)
258 !mx=paw_ij%nocctot(2)
259 !my=paw_ij%nocctot(3)
260 !mz=paw_ij%nocctot(4)
261 !mnorm=sqrt(mx*mx+my*my+mz*mz)
262 !n_upup=half*(n_tot+mnorm)
263 !n_dndn=half*(n_tot-mnorm)
264 !end if
265  n_upup=n_sig
266  n_dndn=n_msig
267 
268  edcdctemp=zero;edctemp=zero
269 
270 !Full localized limit
271  if((pawtab%usepawu==1.or.pawtab%usepawu==4).or.(dmftdc==1.or.dmftdc==4.or.dmftdc==5)) then
272    jpawu_dc=jpawu
273    if(dmftdc==4)  then
274      jpawu_dc=zero
275    end if
276    edcdctemp=edcdctemp-half*upawu*n_tot**2
277    edctemp  =edctemp  +half*upawu*(n_tot*(n_tot-one))
278    if (paw_ij%ndij/=4.or.option_interaction==2) then
279      if(dmftdc/=5) then
280        edcdctemp=edcdctemp+half*jpawu_dc*(n_upup**2+n_dndn**2)
281        edctemp  =edctemp  -half*jpawu_dc*(n_upup*(n_upup-one)+n_dndn*(n_dndn-one))
282      else if(dmftdc==5)  then
283        edcdctemp=edcdctemp+quarter*jpawu_dc*n_tot**2
284        edctemp  =edctemp  -quarter*jpawu_dc*(n_tot*(n_tot-two))
285      end if
286    else if (paw_ij%ndij==4.and.option_interaction==1) then
287 !    write(message,'(a)') "  warning: option_interaction==1 for test         "
288 !    call wrtout(std_out,message,'COLL')
289      edcdctemp=edcdctemp+quarter*jpawu_dc*n_tot**2
290      edctemp  =edctemp  -quarter*jpawu_dc*(n_tot*(n_tot-two))
291    else if (paw_ij%ndij==4.and.option_interaction==3) then
292 !    edcdctemp= \frac{J}/{4}[ N(N) + \vect{m}.\vect{m}]
293      edcdctemp=edcdctemp+quarter*jpawu_dc*(n_tot**2 + &
294 &     mx**2+my**2+mz**2)  ! +\frac{J}/{4}\vect{m}.\vect{m}
295 !    edctemp= -\frac{J}/{4}[ N(N-2) + \vect{m}.\vect{m}]
296      edctemp  =edctemp  -quarter*jpawu_dc*(  &
297 &     (n_tot*(n_tot-two)) +   &
298 &     mx**2+my**2+mz**2)  ! -\frac{J}/{4}\vect{m}.\vect{m}
299    end if
300 
301 !  Around mean field
302  else if(pawtab%usepawu==2.or.dmftdc==2) then
303    edctemp=edctemp+upawu*(n_upup*n_dndn)&
304 &   +half*(upawu-jpawu)*(n_upup**2+n_dndn**2) &
305 &   *(dble(2*lpawu)/dble(2*lpawu+1))
306    edcdctemp=-edctemp
307  else if(pawtab%usepawu==3.or.dmftdc==3) then
308    edcdctemp=edcdc_opt3
309    if(abs(pawprtvol)>=3) then
310      write(message,fmt=11) "edcdc_opt3          ",edcdc_opt3
311      call wrtout(std_out,message,'COLL')
312      write(message,fmt=11) "eks_opt3            ",eks_opt3
313      call wrtout(std_out,message,'COLL')
314      write(message,fmt=11) "eks+edcdc_opt3      ",eks_opt3+edcdc_opt3
315      call wrtout(std_out,message,'COLL')
316      write(message,fmt=11) "(eks+edcdc_opt3)/2  ",(eks_opt3+edcdc_opt3)/2.d0
317      call wrtout(std_out,message,'COLL')
318    end if
319  end if
320 
321 
322  eldaumdc  =eldaumdc  +eldautemp-edctemp
323  eldaumdcdc=eldaumdcdc-eldautemp-edcdctemp
324 
325 !if(pawtab%usepawu/=10.or.pawprtvol>=3) then
326  if(abs(pawprtvol)>=3) then
327    if(pawtab%usepawu<10) then
328      write(message, '(5a,i4)')ch10,'======= LDA+U Energy terms (in Hartree) ====',ch10,&
329 &     ch10,' For Atom ',iatom
330    else if (pawtab%usepawu >= 10) then
331      write(message, '(5a,i4)')ch10,'  ===   LDA+U Energy terms for the DMFT occupation matrix ==',ch10,&
332 &     ch10,' For Atom ',iatom
333    end if
334 
335    call wrtout(std_out,message,'COLL')
336    write(message, '(a)' )"   Contributions to the direct expression of energy:"
337    call wrtout(std_out,  message,'COLL')
338    write(message,fmt=11) "     Double counting  correction   =",edctemp
339    call wrtout(std_out,  message,'COLL')
340    write(message,fmt=11) "     Interaction energy            =",eldautemp
341    call wrtout(std_out,  message,'COLL')
342    write(message,fmt=11) "     Total LDA+U Contribution      =",eldautemp-edctemp
343    call wrtout(std_out,  message,'COLL')
344    write(message, '(a)' )' '
345    call wrtout(std_out,  message,'COLL')
346    write(message, '(a)' )"   For the ""Double-counting"" decomposition:"
347    call wrtout(std_out,  message,'COLL')
348    write(message,fmt=11) "     LDA+U Contribution            =",-eldautemp-edcdctemp
349    call wrtout(std_out,  message,'COLL')
350    11 format(a,e20.10)
351    if(abs(pawprtvol)>=2) then
352      write(message,fmt=11)"     edcdctemp                     =",edcdctemp
353      call wrtout(std_out,  message,'COLL')
354      write(message,fmt=11)"     eldaumdcdc for current atom   =",-eldautemp-edcdctemp
355      call wrtout(std_out,  message,'COLL')
356      write(message, '(a)' )' '
357      call wrtout(std_out,  message,'COLL')
358      write(message,fmt=11)"   pawuenergy: -VUKS pred          =",eldaumdcdc-eldaumdc
359      call wrtout(std_out,  message,'COLL')
360    end if
361    write(message, '(a)' )' '
362    call wrtout(std_out,  message,'COLL')
363  end if
364 
365 !For DMFT calculation
366  if(present(e_ee))   e_ee=e_ee+eldautemp
367  if(present(e_dc))   e_dc=e_dc+edctemp
368  if(present(e_dcdc)) e_dcdc=e_dcdc+edcdctemp
369 
370  DBG_EXIT("COLL")
371 
372  end subroutine pawuenergy