TABLE OF CONTENTS
ABINIT/pawuenergy [ 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