1c     Calculates the steric energy as defined in the reference below
2c     for a set of grid points and integration weights
3
4      subroutine steric_energy(tol_rho, fac, rho,
5     &       delrho, nq, qwght, ipol, energy)
6c
7      implicit none
8c
9#include "dft2drv.fh"
10c
11      double precision tol_rho, fac
12      integer nq, ipol
13      double precision qwght(nq) ! integration weights
14      double precision energy
15c
16c     Charge Density
17c
18      double precision rho(nq,ipol*(ipol+1)/2)
19c
20c     Charge Density Gradient
21c
22      double precision delrho(nq,3,ipol)
23c
24c     Sampling Matrices for the XC Potential
25c
26      double precision BETA
27      Parameter (BETA = 0.05D0)
28c
29c     References:
30c     S.B. Liu, J. Chem. Phys. 126, 244103 (2007).
31c     S.B. Liu, N. Govind, J. Phys. Chem. A. 112, 6690 (2008).
32c
33c***************************************************************************
34c
35      integer n
36      !double precision arcsinh
37      !double precision rho13, rho43, gamma, x, g, gdenom
38      double precision gamma, x
39c
40c     arcsinh(x)=log(x+dsqrt(1d0+x*x))
41c
42      if (ipol.eq.1) then
43c
44c        ======> SPIN-RESTRICTED <======
45c
46         energy = 0.d0
47         do 10 n = 1, nq
48            if (rho(n,1).lt.tol_rho) goto 10
49c
50c           Spin alpha:
51c
52c           rho13 = (0.5d0*rho(n,1))**(1.d0/3.d0)
53c           rho43 = rho13**4
54            gamma = delrho(n,1,1)*delrho(n,1,1) +
55     &              delrho(n,2,1)*delrho(n,2,1) +
56     &              delrho(n,3,1)*delrho(n,3,1)
57            if (gamma.gt.tol_rho)then
58               x = gamma/rho(n,1)
59            else
60               x = 0d0
61            endif
62            energy = energy + x*qwght(n)
63 10      continue
64c
65      else
66c
67c        ======> SPIN-UNRESTRICTED <======
68c
69         energy = 0.d0
70         do 20 n = 1, nq
71            if (rho(n,1).lt.tol_rho) goto 20
72            if (rho(n,2).lt.tol_rho) goto 25
73c
74c           Spin alpha:
75c
76c           rho13 = abs(rho(n,2))**(1.d0/3.d0)*sign(1d0,rho(n,2))
77c           rho43 = rho13**4
78            gamma = delrho(n,1,1)*delrho(n,1,1) +
79     &              delrho(n,2,1)*delrho(n,2,1) +
80     &              delrho(n,3,1)*delrho(n,3,1)
81            if (gamma.gt.tol_rho)then
82               x = gamma/rho(n,2)
83            else
84               x = 0d0
85            endif
86            energy = energy + x*qwght(n)
87 25         continue
88c
89c           Spin beta:
90c
91            if (rho(n,3).lt.tol_rho) goto 20
92c
93c           rho13 = abs(rho(n,3))**(1.d0/3.d0)*sign(1d0,rho(n,3))
94c           rho43 = rho13**4
95            gamma = delrho(n,1,2)*delrho(n,1,2) +
96     &              delrho(n,2,2)*delrho(n,2,2) +
97     &              delrho(n,3,2)*delrho(n,3,2)
98            if (gamma.gt.tol_rho)then
99               x = gamma/rho(n,3)
100            else
101               x = 0d0
102            endif
103            energy = energy + x*qwght(n)
104 20      continue
105c
106      endif
107c
108      energy = energy/8.d0
109
110      return
111      end
112c $Id$
113