TABLE OF CONTENTS
ABINIT/dfpt_prtene [ Functions ]
NAME
dfpt_prtene
FUNCTION
Print components of second derivative of total energy in nice format
COPYRIGHT
Copyright (C) 1999-2018 ABINIT group (XG, DRH, 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/doc/developers/contributors.txt .
INPUTS
eberry=energy associated with Berry phase edocc=correction to 2nd-order total energy coming from changes of occupation eeig0=0th-order eigenenergies part of 2nd-order total energy eew=Ewald part of 2nd-order total energy efrhar=hartree frozen-wavefunction part of 2nd-order tot. en. efrkin=kinetic frozen-wavefunction part of 2nd-order tot. en. efrloc=local psp. frozen-wavefunction part of 2nd-order tot. en. efrnl=nonlocal psp. frozen-wavefunction part of 2nd-order tot. en efrx1=xc core corr.(1) frozen-wavefunction part of 2nd-order tot. en efrx2=xc core corr.(2) frozen-wavefunction part of 2nd-order tot. en ehart01=inhomogeneous 1st-order Hartree part of 2nd-order total energy for strain perturbation only (zero otherwise, and not used) ehart1=1st-order Hartree part of 2nd-order total energy eii=pseudopotential core part of 2nd-order total energy ek0=0th-order kinetic energy part of 2nd-order total energy. ek1=1st-order kinetic energy part of 2nd-order total energy. eloc0=0th-order local (psp+vxc+Hart) part of 2nd-order total energy elpsp1=1st-order local pseudopot. part of 2nd-order total energy. enl0=0th-order nonlocal pseudopot. part of 2nd-order total energy. enl1=1st-order nonlocal pseudopot. part of 2nd-order total energy. eovl1=1st-order change of wave-functions overlap, part of 2nd-order energy PAW only - Eq(79) and Eq(80) of PRB 78, 035105 (2008) epaw1=1st-order PAW on-site part of 2nd-order total energy. evdw=DFT-D semi-empirical part of 2nd-order total energy exc1=1st-order exchange-correlation part of 2nd-order total energy iout=unit number to which output is written ipert=type of the perturbation natom=number of atoms in unit cell usepaw= 0 for non paw calculation; =1 for paw calculation usevdw= flag set to 1 if vdw DFT-D semi-empirical potential is in use
OUTPUT
(only writing)
NOTES
all energies in Hartree
PARENTS
dfpt_looppert
CHILDREN
wrtout
SOURCE
63 #if defined HAVE_CONFIG_H 64 #include "config.h" 65 #endif 66 67 #include "abi_common.h" 68 69 70 subroutine dfpt_prtene(berryopt,eberry,edocc,eeig0,eew,efrhar,efrkin,efrloc,efrnl,efrx1,efrx2,& 71 & ehart01,ehart1,eii,ek0,ek1,eloc0,elpsp1,enl0,enl1,eovl1,epaw1,evdw,exc1,iout,ipert,natom,& 72 & usepaw,usevdw) 73 74 use defs_basis 75 use m_profiling_abi 76 77 !This section has been created automatically by the script Abilint (TD). 78 !Do not modify the following lines by hand. 79 #undef ABI_FUNC 80 #define ABI_FUNC 'dfpt_prtene' 81 use interfaces_14_hidewrite 82 !End of the abilint section 83 84 implicit none 85 86 !Arguments ------------------------------- 87 !scalars 88 integer,intent(in) :: berryopt,iout,ipert,natom,usepaw,usevdw 89 real(dp),intent(in) :: eberry,edocc,eeig0,eew,efrhar,efrkin,efrloc,efrnl,efrx1 90 real(dp),intent(in) :: efrx2,ehart01,ehart1,eii,ek0,ek1,eloc0,elpsp1,enl0,enl1 91 real(dp),intent(in) :: eovl1,epaw1,evdw,exc1 92 93 !Local variables ------------------------- 94 !scalars 95 integer :: nn 96 logical :: berry_activated 97 real(dp) :: enl1_effective,erelax,etotal 98 character(len=10) :: numb 99 character(len=10),parameter :: numbstr(20) = & 100 & (/'One ','Two ','Three ','Four ','Five ', & 101 & 'Six ','Seven ','Eight ','Nine ','Ten ', & 102 & 'Eleven ','Twelve ','Thirteen ','Fourteen ','Fifteen ', & 103 & 'Sixteen ','Seventeen ','Eighteen ','Nineteen ','Twenty '/) 104 character(len=500) :: message 105 106 ! ********************************************************************* 107 108 !Count and print the number of components of 2nd-order energy 109 !MT feb 2015: this number is wrong! Should change it but 110 ! need to change a lot of ref. files 111 berry_activated=(berryopt== 4.or.berryopt== 6.or.berryopt== 7.or. & 112 & berryopt==14.or.berryopt==16.or.berryopt==17) 113 if (ipert==natom+1) nn=8 114 if (ipert==natom+5) nn=8 115 if (ipert==natom+2) nn=7 116 if (ipert>=1.and.ipert<=natom) nn=13 117 if (ipert==natom+3.or.ipert==natom+4) nn=17 118 if (ipert==natom+2.and.berry_activated) nn=nn+1 119 if (usepaw==1) nn=nn+1 120 if (usevdw==1) nn=nn+1 121 if (ipert==natom+10.or.ipert==natom+11) nn=1 ! means nothing, 122 ! because we do not compute derivatives of the energy in this case 123 write(message, '(4a)' ) ch10,& 124 & ' ',trim(numbstr(nn)),' components of 2nd-order total energy (hartree) are ' 125 call wrtout(iout,message,'COLL') 126 call wrtout(std_out,message,'COLL') 127 128 numb='1,2,3' 129 write(message, '(3a)' )& 130 & ' ',trim(numb),': 0th-order hamiltonian combined with 1st-order wavefunctions' 131 call wrtout(iout,message,'COLL') 132 call wrtout(std_out,message,'COLL') 133 write(message, '(a,es17.8,a,es17.8,a,es17.8)' )& 134 & ' kin0=',ek0, ' eigvalue=',eeig0,' local=',eloc0 135 call wrtout(iout,message,'COLL') 136 call wrtout(std_out,message,'COLL') 137 138 numb='4,5,6';if( ipert==natom+3.or.ipert==natom+4) numb='4,5,6,7' 139 write(message, '(3a)' )& 140 & ' ',trim(numb),': 1st-order hamiltonian combined with 1st and 0th-order wfs' 141 call wrtout(iout,message,'COLL') 142 call wrtout(std_out,message,'COLL') 143 if(ipert/=natom+1.and.ipert/=natom+2)then 144 write(message, '(a,es17.8,a,es17.8,a,es17.8,a,a)' ) & 145 & ' loc psp =',elpsp1,' Hartree=',ehart1,' xc=',exc1,ch10,& 146 & ' note that "loc psp" includes a xc core correction that could be resolved' 147 else if(ipert==natom+1) then 148 write(message, '(a,es17.8,a,es17.8,a,es17.8)' ) & 149 & ' kin1=',ek1, ' Hartree=',ehart1,' xc=',exc1 150 else if(ipert==natom+2) then 151 write(message, '(a,es17.8,a,es17.8,a,es17.8)' ) & 152 & ' dotwf=',enl1, ' Hartree=',ehart1,' xc=',exc1 153 end if 154 if(ipert==natom+3 .or. ipert==natom+4) then 155 write(message, '(a,es17.8,a,es17.8,a,es17.8,a,a,es17.8)' ) & 156 & ' loc psp =',elpsp1,' Hartree=',ehart1,' xc=',exc1,ch10,& 157 & ' kin1=',ek1 158 end if 159 call wrtout(iout,message,'COLL') 160 call wrtout(std_out,message,'COLL') 161 162 enl1_effective=enl1;if (ipert==natom+2) enl1_effective=zero 163 numb='7,8,9';if( ipert==natom+3.or.ipert==natom+4) numb='8,9,10' 164 write(message, '(5a,es17.8,a,es17.8,a,es17.8)' )& 165 & ' ',trim(numb),': eventually, occupation + non-local contributions',ch10,& 166 & ' edocc=',edocc,' enl0=',enl0,' enl1=',enl1_effective 167 call wrtout(iout,message,'COLL') 168 call wrtout(std_out,message,'COLL') 169 170 if (usepaw==1) then 171 numb='10';if( ipert==natom+3.or.ipert==natom+4) numb='11' 172 write(message, '(3a,es17.8)' )& 173 & ' ',trim(numb),': eventually, PAW "on-site" Hxc contribution: epaw1=',epaw1 174 call wrtout(iout,message,'COLL') 175 call wrtout(std_out,message,'COLL') 176 end if 177 178 if(ipert/=natom+10 .and.ipert/=natom+11) then 179 erelax=0.0_dp 180 if(ipert>=1.and.ipert<=natom)then 181 erelax=ek0+edocc+eeig0+eloc0+elpsp1+ehart1+exc1+enl0+enl1+epaw1 182 else if(ipert==natom+1.or.ipert==natom+2)then 183 erelax=ek0+edocc+eeig0+eloc0+ek1+ehart1+exc1+enl0+enl1+epaw1 184 else if(ipert==natom+3.or.ipert==natom+4)then 185 erelax=ek0+edocc+eeig0+eloc0+ek1+elpsp1+ehart1+exc1+enl0+enl1+epaw1 186 else if(ipert==natom+5)then 187 erelax=ek0+edocc+eeig0+eloc0+ek1+elpsp1+ehart1+exc1+enl0+enl1+epaw1 188 end if 189 enl1_effective=enl1 190 if (ipert==natom+1.or.ipert==natom+2) then 191 if (1.0_dp+enl1/10.0_dp==1.0_dp) enl1_effective=zero 192 end if 193 194 numb='1-9';if (usepaw==1) numb='1-10' 195 if( ipert==natom+3.or.ipert==natom+4) then 196 numb='1-10';if (usepaw==1) numb='1-11' 197 end if 198 write(message, '(5a,es17.8)' )& 199 & ' ',trim(numb),' gives the relaxation energy (to be shifted if some occ is /=2.0)',& 200 & ch10,' erelax=',erelax 201 call wrtout(iout,message,'COLL') 202 call wrtout(std_out,message,'COLL') 203 end if 204 205 if(ipert>=1.and.ipert<=natom)then 206 207 numb='10,11,12';if (usepaw==1) numb='11,12,13' 208 write(message, '(4a)' )& 209 & ' ',trim(numb),' Non-relaxation contributions : ',& 210 & 'frozen-wavefunctions and Ewald' 211 call wrtout(iout,message,'COLL') 212 call wrtout(std_out,message,'COLL') 213 write(message, '(a,es17.8,a,es17.8,a,es17.8)' ) & 214 & ' fr.local=',efrloc,' fr.nonlo=',efrnl,' Ewald=',eew 215 call wrtout(iout,message,'COLL') 216 call wrtout(std_out,message,'COLL') 217 218 write(message, '(a,es16.6)' )' dfpt_prtene : non-relax=',efrloc+efrnl+eew 219 call wrtout(std_out,message,'COLL') 220 221 numb='13,14';if (usepaw==1) numb='14,15' 222 write(message, '(3a)' )& 223 & ' ',trim(numb),' Frozen wf xc core corrections (1) and (2)' 224 call wrtout(iout,message,'COLL') 225 call wrtout(std_out,message,'COLL') 226 write(message, '(a,es17.8,a,es17.8)' ) & 227 & ' frxc 1 =',efrx1,' frxc 2 =',efrx2 228 call wrtout(iout,message,'COLL') 229 call wrtout(std_out,message,'COLL') 230 if (usepaw==1) then 231 numb='16' 232 write(message, '(5a,es17.8)' )& 233 & ' ',trim(numb),' Contribution from 1st-order change of wavefunctions overlap',& 234 & ch10,' eovl1 =',eovl1 235 call wrtout(iout,message,'COLL') 236 call wrtout(std_out,message,'COLL') 237 end if 238 if (usevdw==1) then 239 numb='15';if (usepaw==1) numb='17' 240 write(message, '(3a,es17.8)' )& 241 & ' ',trim(numb),' Contribution from van der Waals DFT-D: evdw =',evdw 242 call wrtout(iout,message,'COLL') 243 call wrtout(std_out,message,'COLL') 244 end if 245 246 write(message, '(a)' )' Resulting in : ' 247 call wrtout(iout,message,'COLL') 248 call wrtout(std_out,message,'COLL') 249 etotal=erelax+eew+efrloc+efrnl+efrx1+efrx2+evdw 250 write(message, '(a,e20.10,a,e22.12,a)' ) & 251 & ' 2DEtotal=',etotal,' Ha. Also 2DEtotal=',etotal*Ha_eV,' eV' 252 call wrtout(iout,message,'COLL') 253 call wrtout(std_out,message,'COLL') 254 write(message, '(a,es20.10,a,es20.10,a)' ) & 255 & ' (2DErelax=',erelax,' Ha. 2DEnonrelax=',etotal-erelax,' Ha)' 256 call wrtout(iout,message,'COLL') 257 call wrtout(std_out,message,'COLL') 258 write(message, '(a,es20.10,a,a)' ) & 259 & ' ( non-var. 2DEtotal :',& 260 & 0.5_dp*(elpsp1+enl1)+eovl1+eew+efrloc+efrnl+efrx1+efrx2+evdw,' Ha)',ch10 261 call wrtout(iout,message,'COLL') 262 call wrtout(std_out,message,'COLL') 263 264 else if(ipert==natom+1.or.ipert==natom+2)then 265 if (ipert==natom+1.and.usepaw==1) then 266 numb='11' 267 write(message, '(5a,es17.8)' )& 268 & ' ',trim(numb),' Contribution from 1st-order change of wavefunctions overlap',& 269 & ch10,' eovl1 =',eovl1 270 call wrtout(iout,message,'COLL') 271 call wrtout(std_out,message,'COLL') 272 end if 273 write(message,*)' No Ewald or frozen-wf contrib.:',& 274 & ' the relaxation energy is the total one' 275 if(berry_activated) then 276 numb='10'; 277 write(message,'(3a,es20.10)')' ',trim(numb),' Berry phase energy :',eberry 278 end if 279 call wrtout(iout,message,'COLL') 280 call wrtout(std_out,message,'COLL') 281 etotal=erelax 282 write(message, '(a,e20.10,a,e22.12,a)' ) & 283 & ' 2DEtotal=',etotal,' Ha. Also 2DEtotal=',etotal*Ha_eV,' eV' 284 call wrtout(iout,message,'COLL') 285 call wrtout(std_out,message,'COLL') 286 write(message, '(a,es20.10,a)' ) & 287 & ' ( non-var. 2DEtotal :',0.5_dp*(ek1+enl1_effective)+eovl1,' Ha)' 288 call wrtout(iout,message,'COLL') 289 call wrtout(std_out,message,'COLL') 290 291 else if(ipert==natom+3 .or. ipert==natom+4) then 292 numb='11,12,13';if (usepaw==1) numb='12,13,14' 293 write(message, '(4a)' )& 294 & ' ',trim(numb),' Non-relaxation contributions : ',& 295 & 'frozen-wavefunctions and Ewald' 296 call wrtout(iout,message,'COLL') 297 call wrtout(std_out,message,'COLL') 298 write(message, '(a,es17.8,a,es17.8,a,es17.8)' ) & 299 & ' fr.hart=',efrhar,' fr.kin=',efrkin,' fr.loc=',efrloc 300 call wrtout(iout,message,'COLL') 301 call wrtout(std_out,message,'COLL') 302 303 numb='14,15,16';if (usepaw==1) numb='15,16,17' 304 write(message, '(4a)' )& 305 & ' ',trim(numb),' Non-relaxation contributions : ',& 306 & 'frozen-wavefunctions and Ewald' 307 call wrtout(iout,message,'COLL') 308 call wrtout(std_out,message,'COLL') 309 write(message, '(a,es17.8,a,es17.8,a,es17.8)' ) & 310 & ' fr.nonl=',efrnl,' fr.xc=',efrx1,' Ewald=',eew 311 call wrtout(iout,message,'COLL') 312 call wrtout(std_out,message,'COLL') 313 314 numb='17';if (usepaw==1) numb='18' 315 write(message, '(4a)' )& 316 & ' ',trim(numb),' Non-relaxation contributions : ',& 317 & 'pseudopotential core energy' 318 call wrtout(iout,message,'COLL') 319 call wrtout(std_out,message,'COLL') 320 write(message, '(a,es17.8)' ) & 321 & ' pspcore=',eii 322 call wrtout(iout,message,'COLL') 323 call wrtout(std_out,message,'COLL') 324 if (usepaw==1) then 325 numb='19' 326 write(message, '(5a,es17.8)' )& 327 & ' ',trim(numb),' Contribution from 1st-order change of wavefunctions overlap',& 328 & ch10,' eovl1 =',eovl1 329 call wrtout(iout,message,'COLL') 330 call wrtout(std_out,message,'COLL') 331 end if 332 if (usevdw==1) then 333 numb='18';if (usepaw==1) numb='20' 334 write(message, '(3a,es17.8)' )& 335 & ' ',trim(numb),' Contribution from van der Waals DFT-D: evdw =',evdw 336 call wrtout(iout,message,'COLL') 337 call wrtout(std_out,message,'COLL') 338 end if 339 340 write(message, '(a,es16.6)' )' dfpt_prtene : non-relax=',& 341 & efrhar+efrkin+efrloc+efrnl+efrx1+eew+evdw 342 call wrtout(std_out,message,'COLL') 343 write(message, '(a)' )' Resulting in : ' 344 call wrtout(iout,message,'COLL') 345 call wrtout(std_out,message,'COLL') 346 etotal=erelax+efrhar+efrkin+efrloc+efrnl+efrx1+eew+eii+evdw 347 write(message, '(a,e20.10,a,e22.12,a)' ) & 348 & ' 2DEtotal=',etotal,' Ha. Also 2DEtotal=',etotal*Ha_eV,' eV' 349 call wrtout(iout,message,'COLL') 350 call wrtout(std_out,message,'COLL') 351 write(message, '(a,es20.10,a,es20.10,a)' ) & 352 & ' (2DErelax=',erelax,' Ha. 2DEnonrelax=',etotal-erelax,' Ha)' 353 call wrtout(iout,message,'COLL') 354 call wrtout(std_out,message,'COLL') 355 write(message, '(a,es20.10,a,a)' ) & 356 & ' ( non-var. 2DEtotal :',& 357 & 0.5_dp*(elpsp1+enl1+ek1+ehart01)+eovl1+& 358 & efrhar+efrkin+efrloc+efrnl+efrx1+eew+eii+evdw,' Ha)',ch10 359 call wrtout(iout,message,'COLL') 360 call wrtout(std_out,message,'COLL') 361 end if 362 363 end subroutine dfpt_prtene