1#ifndef SECOND_DERIV 2 Subroutine xc_xpw91(tol_rho, fac, lfac, nlfac, rho, delrho, 3 & Amat, Cmat, nq, ipol, Ex, qwght,ldew,func) 4#else 5 Subroutine xc_xpw91_d2(tol_rho, fac, lfac, nlfac, rho, delrho, 6 & Amat, Amat2, Cmat, Cmat2, nq, ipol, Ex, 7 & qwght,ldew,func) 8#endif 9c 10C$Id$ 11c 12 implicit none 13c 14#include "dft2drv.fh" 15c 16 double precision fac, Ex 17 integer nq, ipol 18 logical lfac, nlfac,ldew 19 double precision func(*) ! value of the functional [output] 20c 21c Charge Density & Its Cube Root 22c 23 double precision rho(nq,ipol*(ipol+1)/2) 24c 25c Charge Density Gradient 26c 27 double precision delrho(nq,3,ipol) 28c 29c Quadrature Weights 30c 31 double precision qwght(nq) 32c 33c Sampling Matrices for the XC Potential & Energy 34c 35 double precision Amat(nq,ipol), Cmat(nq,*) 36c 37c 38c Compute the partial derivatives of the exchange functional of Perdew91. 39c 40c Becke & Perdew Parameters 41c 42 double precision DPOW 43 double precision BETA, tol_rho, CPW91 44 45 Parameter (BETA = 0.0042D0, CPW91=1.6455D0,DPOW=4) ! pw91 paper 46C 47C 48 49#ifdef SECOND_DERIV 50c 51c Second Derivatives of the Exchange Energy Functional 52c 53 double precision Amat2(nq,NCOL_AMAT2), Cmat2(nq,NCOL_CMAT2) 54#endif 55c 56c References: 57c 58c 59c*************************************************************************** 60 61 integer n 62 double precision gamma 63c 64 if (ipol.eq.1 )then 65c 66c ======> SPIN-RESTRICTED <====== 67c 68 do 10 n = 1, nq 69 if (rho(n,1).lt.tol_rho) goto 10 70 gamma = delrho(n,1,1)*delrho(n,1,1) + 71 & delrho(n,2,1)*delrho(n,2,1) + 72 & delrho(n,3,1)*delrho(n,3,1) 73 gamma=0.25d0*gamma 74#ifndef SECOND_DERIV 75 call xc_xpw91core(DPOW,BETA,n,1, 76 & rho(n,1),gamma,qwght(n),func(n), 77 & tol_rho, fac, lfac, nlfac, 78 & Amat, Cmat, nq, ipol, Ex, ldew) 79#else 80 call xc_xpw91core_d2(DPOW,BETA,n,1, 81 & rho(n,1),gamma,qwght(n),func(n), 82 & tol_rho, fac, lfac, nlfac, 83 & Amat, Amat2, Cmat, Cmat2, nq, ipol, Ex, 84 & ldew) 85#endif 86 87 10 continue 88c 89 else 90c 91c ======> SPIN-UNRESTRICTED <====== 92c 93 do 20 n = 1, nq 94 if (rho(n,1).lt.tol_rho) goto 20 95c 96c Spin alpha: 97c 98 if (rho(n,2).gt.tol_rho) then 99 gamma = delrho(n,1,1)*delrho(n,1,1) + 100 & delrho(n,2,1)*delrho(n,2,1) + 101 & delrho(n,3,1)*delrho(n,3,1) 102#ifndef SECOND_DERIV 103 call xc_xpw91core(DPOW,BETA,n,1, 104 & rho(n,2),gamma,qwght(n),func(n), 105 & tol_rho, fac, lfac, nlfac, 106 & Amat, Cmat, nq, ipol, Ex, ldew) 107#else 108 call xc_xpw91core_d2(DPOW,BETA,n,1, 109 & rho(n,2),gamma,qwght(n),func(n), 110 & tol_rho, fac, lfac, nlfac, 111 & Amat, Amat2, Cmat, Cmat2, nq, ipol, Ex, 112 & ldew) 113#endif 114 endif 115c 116c Spin beta: 117c 118 if (rho(n,3).gt.tol_rho) then 119 120 gamma = delrho(n,1,2)*delrho(n,1,2) + 121 & delrho(n,2,2)*delrho(n,2,2) + 122 & delrho(n,3,2)*delrho(n,3,2) 123#ifndef SECOND_DERIV 124 call xc_xpw91core(DPOW,BETA,n,2, 125 & rho(n,3),gamma,qwght(n),func(n), 126 & tol_rho, fac, lfac, nlfac, 127 & Amat, Cmat, nq, ipol, Ex, ldew) 128#else 129 call xc_xpw91core_d2(DPOW,BETA,n,2, 130 & rho(n,3),gamma,qwght(n),func(n), 131 & tol_rho, fac, lfac, nlfac, 132 & Amat, Amat2, Cmat, Cmat2, nq, ipol, Ex, 133 & ldew) 134#endif 135 endif 136c 137 20 continue 138c 139 endif 140c 141 return 142 end 143#ifndef SECOND_DERIV 144#define SECOND_DERIV 145c 146c Compile source again for the 2nd derivative case 147c 148#include "xc_xpw91.F" 149#endif 150