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