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