1! --- 2! Copyright (C) 1996-2016 The SIESTA group 3! This file is distributed under the terms of the 4! GNU General Public License: see COPYING in the top directory 5! or http://www.gnu.org/copyleft/gpl.txt . 6! See Docs/Contributors.txt for a list of contributors. 7! --- 8module basis_enthalpy 9 implicit none 10 private 11 public :: write_basis_enthalpy 12 13CONTAINS 14 15 subroutine write_basis_enthalpy(FE,FEharris) 16 use precision, only: dp 17 use fdf, only: fdf_physical 18 use units, only: GPa, eV 19 20 real(dp), intent(in) :: FE ! Electronic (Free)Energy 21 real(dp), intent(in) :: FEharris ! Electronic (Free)HarrisEnergy 22 23 real(dp) orb_vol, basis_enthalpy, basis_harris_enthalpy, basis_pressure 24 integer iu 25 26 call orb_volumes(orb_vol) 27 basis_pressure = GPa * fdf_physical("BasisPressure",0.2_dp,"GPa") 28 basis_enthalpy = FE + basis_pressure * orb_vol 29 basis_harris_enthalpy = FEHarris + basis_pressure * orb_vol 30 31 write(6,"(a,f18.6)") "(Free)E+ p_basis*V_orbitals = ", basis_enthalpy/eV 32 write(6,"(a,f18.6)") "(Free)Eharris+ p_basis*V_orbitals = ", basis_harris_enthalpy/eV 33 34 call io_assign(iu) 35 open(iu,file="BASIS_ENTHALPY",form="formatted",status="replace") 36 write(iu,*) basis_enthalpy/eV 37 write(iu,*) "The above number is the electronic (free)energy:", FE/eV 38 write(iu,*) "Plus the pressure : ", basis_pressure, & 39 " ( ", basis_pressure/GPa, " GPa)" 40 write(iu,*) " times the orbital volume (in Bohr**3): ", orb_vol 41 call io_close(iu) 42 43 call io_assign(iu) 44 open(iu,file="BASIS_HARRIS_ENTHALPY",form="formatted",status="replace") 45 write(iu,*) basis_harris_enthalpy/eV 46 write(iu,*) "The above number is the electronic (free) harris energy:", FEHarris/eV 47 write(iu,*) "Plus the pressure : ", basis_pressure, & 48 " ( ", basis_pressure/GPa, " GPa)" 49 write(iu,*) " times the orbital volume (in Bohr**3): ", orb_vol 50 call io_close(iu) 51 52 end subroutine write_basis_enthalpy 53!-------------------------------------------------------------------- 54 subroutine orb_volumes(orb_vol) 55 ! 56 ! Computes the total volume of the basis orbitals 57 ! in the unit cell, for the purposes of calculating 58 ! the "basis enthalpy" used in the optimization 59 ! procedure. See: 60 ! E. Anglada et al, PRB ... 61 ! 62 ! Notes: 63 ! 1. Only a single representative of each "nlz" shell is 64 ! used (e.g., just one and not 5 for a 'd' shell) 65 ! 2. All the species are treated equally. Presumably 66 ! one might want to give different weights to different 67 ! species (for example, to avoid the case in which 68 ! an impurity atom's basis details are drowned by the 69 ! much larger influence of the host orbitals). 70 ! This will be implemented via an fdf block: 71 ! 72 ! %block PAO.BasisEnthalpyWeights 73 ! Species_number_i Weigth_i 74 ! ... 75 ! %endblock PAO.BasisEnthalpyWeights. 76 77 use precision, only: dp 78 use atmfuncs, only: lofio, rcut 79 use atomlist, only: lasto, iphorb 80 use siesta_geom, only: isa, na_u 81 use units, only: pi 82 83 real(dp), intent(out) :: orb_vol 84 85 integer :: ia, ioa, is, l, io 86 real(dp) :: r 87 88 orb_vol = 0.0_dp 89 90 do ia = 1,na_u 91 is = isa(ia) 92 do io = lasto(ia-1)+1,lasto(ia) 93 ioa = iphorb(io) 94 l = lofio(is,ioa) 95 r = rcut(is,ioa) 96 orb_vol = orb_vol + ((4.0_dp/3.0_dp)*pi*r**3)/(2*l+1) 97 !Just one per nl shell 98 enddo 99 enddo 100 end subroutine orb_volumes 101 102end module basis_enthalpy 103