1c 2c Modified to handle second derivatives while reusing code 3c 4c BGJ - 8/98 5c 6#if !defined(SECOND_DERIV) && !defined(THIRD_DERIV) 7 Subroutine xc_camxlsd(tol_rho, fac, lfac, nlfac, rho, Amat, nq, 8 & ipol, Ex, qwght, ldew, func) 9#elif defined(SECOND_DERIV) && !defined(THIRD_DERIV) 10c For locations of 2nd derivatives of functionals in array 11#include "dft2drv.fh" 12 Subroutine xc_camxlsd_d2(tol_rho, fac, lfac, nlfac, rho, Amat, 13 & Amat2, nq, ipol, Ex, qwght, ldew, func) 14#else 15c For locations of 3rd derivatives of functionals in array 16#include "dft3drv.fh" 17 Subroutine xc_camxlsd_d3(tol_rho, fac, lfac, nlfac, rho, Amat, 18 1 Amat2, Amat3, nq, ipol, Ex, qwght, ldew, 19 2 func) 20#endif 21c 22C$Id$ 23c 24 Implicit none 25#include "errquit.fh" 26c 27#include "stdio.fh" 28c 29 integer nq, ipol 30 double precision fac, Ex 31 logical ldew, lfac, nlfac 32 double precision func(*) ! value of the functional [output] 33c 34c Charge Density 35c 36 double precision rho(nq,(ipol*(ipol+1))/2) 37c 38c Quadrature Weights 39c 40 double precision qwght(nq) 41c 42c Partial First Derivatives of the Exchange Energy Functional 43c 44 double precision Amat(nq,ipol) 45 double precision Etmp,Atmp,Ctmp,A2tmp,C2tmp,C3tmp 46c 47#if defined(SECOND_DERIV) || defined(THIRD_DERIV) 48c 49c Partial Second Derivatives of the Exchange Energy Functional 50c 51c double precision Amat2(nq,*) 52 double precision Amat2(nq,NCOL_AMAT2) 53#endif 54#ifdef THIRD_DERIV 55c 56c Partial Third Derivatives of the Exchange Energy Functional 57c 58 double precision Amat3(nq,NCOL_AMAT3) 59 double precision A3tmp, C4tmp, C5tmp, C6tmp 60 double precision rhom23 61#endif 62c 63c Compute the partial derivatives of the exchange functional of Dirac. 64c 65 double precision P1, P2, P3, P4, tol_rho 66c 67c P1 = -(3/PI)**(1/3) 68c P2 = -(3/4)*(3/PI)**(1/3) 69c P3 = -(6/PI)**(1/3) 70c P4 = -(3/4)*(6/PI)**(1/3) 71c 72 Parameter (P1 = -0.9847450218426959D+00) 73 Parameter (P2 = -0.7385587663820219D+00) 74 Parameter (P3 = -0.1240700981798799D+01) 75 Parameter (P4 = -0.9305257363490993D+00) 76 double precision rho13, rho32, rho33, one_third 77 Parameter (one_third = 1.d0/3.d0) 78 double precision two_ninth 79 Parameter (two_ninth = 2.0d0/9.0d0) 80 integer n 81c 82 if (ipol.eq.1)then 83c 84c ======> SPIN-RESTRICTED <====== 85c 86 do 10 n = 1, nq 87 if (rho(n,1).gt.tol_rho)then 88 rho13=rho(n,1)**one_third 89 Etmp = rho(n,1)*rho13*P2*fac 90 if(ldew)func(n) = func(n) + rho(n,1)*rho13*fac*P2 91 Atmp = rho13*P1*fac 92 Ctmp = 0.0d0 93#if defined(SECOND_DERIV) || defined(THIRD_DERIV) 94 A2tmp = (rho13/rho(n,1))*2.0d0*one_third*P1*fac 95 C2tmp = 0.0d0 96 C3tmp = 0.0d0 97#endif 98#ifdef THIRD_DERIV 99 rhom23 = rho13/rho(n,1) 100 A3tmp = (rhom23/rho(n,1))*(-4.0d0)*two_ninth*P1*fac 101 C4tmp = 0.0d0 102 C5tmp = 0.0d0 103 C6tmp = 0.0d0 104#endif 105#ifdef THIRD_DERIV 106 call xc_att_xc_d3(rho(n,1),ipol,Etmp,Atmp,Ctmp,A2tmp, 107 & C2tmp,C3tmp,A3tmp,C4tmp,C5tmp,C6tmp) 108c 109 Amat2(n,D2_RA_RA) = Amat2(n,D2_RA_RA) + A2tmp 110c 111 Amat3(n,D3_RA_RA_RA) = Amat3(n,D3_RA_RA_RA) + A3tmp 112#elif defined(SECOND_DERIV) 113 call xc_att_xc_d2(rho(n,1),ipol,Etmp,Atmp,Ctmp,A2tmp, 114 & C2tmp,C3tmp) 115c 116 Amat2(n,D2_RA_RA) = Amat2(n,D2_RA_RA) + A2tmp 117#else 118 call xc_att_xc(rho(n,1),ipol,Etmp,Atmp,Ctmp) 119#endif 120 Ex = Ex + qwght(n)*Etmp 121 Amat(n,1) = Amat(n,1) + Atmp 122 endif 123 10 continue 124c 125 else 126c 127c ======> SPIN-UNRESTRICTED <====== 128c 129 do 20 n = 1,nq 130 if (rho(n,1).gt.tol_rho)then 131 rho32=0.0d0 132 rho33=0.0d0 133 if (rho(n,2).gt.tol_rho) rho32=rho(n,2)**one_third 134 if (rho(n,3).gt.tol_rho) rho33=rho(n,3)**one_third 135c ---- alpha ---- 136 Etmp = rho32*rho(n,2)*P4*fac 137c 138 Atmp = P3*rho32*fac 139 Ctmp = 0.0d0 140#if defined(SECOND_DERIV) || defined(THIRD_DERIV) 141 A2tmp = 0.0d0 142 C2tmp = 0.0d0 143 C3tmp = 0.0d0 144c 145 if (rho(n,2).gt.tol_rho) then 146 A2tmp = one_third*P3*rho32/rho(n,2)*fac 147 endif 148#endif 149#ifdef THIRD_DERIV 150 A3tmp = 0.0d0 151 C4tmp = 0.0d0 152 C5tmp = 0.0d0 153 C6tmp = 0.0d0 154c 155 if (rho(n,2).gt.tol_rho) then 156 A3tmp = -two_ninth*P3*rho32/(rho(n,2)**2)*fac 157 endif 158#endif 159#ifdef THIRD_DERIV 160 if (rho(n,2).gt.tol_rho) then 161 call xc_att_xc_d3(rho(n,2),ipol,Etmp,Atmp,Ctmp,A2tmp, 162 & C2tmp,C3tmp,A3tmp,C4tmp,C5tmp,C6tmp) 163 endif 164c 165 Amat2(n,D2_RA_RA) = Amat2(n,D2_RA_RA) + A2tmp 166c 167 Amat3(n,D3_RA_RA_RA) = Amat3(n,D3_RA_RA_RA) + A3tmp 168#elif defined(SECOND_DERIV) 169 if (rho(n,2).gt.tol_rho) then 170 call xc_att_xc_d2(rho(n,2),ipol,Etmp,Atmp,Ctmp,A2tmp, 171 & C2tmp,C3tmp) 172 endif 173c 174 Amat2(n,D2_RA_RA) = Amat2(n,D2_RA_RA) + A2tmp 175#else 176 if (rho(n,2).gt.tol_rho) then 177 call xc_att_xc(rho(n,2),ipol,Etmp,Atmp,Ctmp) 178 endif 179#endif 180 Amat(n,1) = Amat(n,1) + Atmp 181 Ex = Ex + qwght(n)*Etmp 182c ---- beta ---- 183 Etmp = rho33*rho(n,3)*P4*fac 184 Atmp = P3*rho33*fac 185 Ctmp = 0.0d0 186#if defined(SECOND_DERIV) || defined(THIRD_DERIV) 187 A2tmp = 0.0d0 188 C2tmp = 0.0d0 189 C3tmp = 0.0d0 190c 191 if (rho(n,3).gt.tol_rho) then 192 A2tmp = one_third*P3*rho33/rho(n,3)*fac 193 end if 194#endif 195#ifdef THIRD_DERIV 196 A3tmp = 0.0d0 197 C4tmp = 0.0d0 198 C5tmp = 0.0d0 199 C6tmp = 0.0d0 200c 201 if (rho(n,3).gt.tol_rho) then 202 A3tmp = -two_ninth*P3*rho33/(rho(n,3)**2)*fac 203 endif 204#endif 205#ifdef THIRD_DERIV 206 if (rho(n,3).gt.tol_rho) then 207 call xc_att_xc_d3(rho(n,3),ipol,Etmp,Atmp,Ctmp,A2tmp, 208 & C2tmp,C3tmp,A3tmp,C4tmp,C5tmp,C6tmp) 209 endif 210c 211 Amat2(n,D2_RB_RB) = Amat2(n,D2_RB_RB) + A2tmp 212c 213 Amat3(n,D3_RB_RB_RB) = Amat3(n,D3_RB_RB_RB) + A3tmp 214#elif defined(SECOND_DERIV) 215 if (rho(n,3).gt.tol_rho) then 216 call xc_att_xc_d2(rho(n,3),ipol,Etmp,Atmp,Ctmp,A2tmp, 217 & C2tmp,C3tmp) 218 end if 219c 220 Amat2(n,D2_RB_RB) = Amat2(n,D2_RB_RB) + A2tmp 221#else 222 if (rho(n,3).gt.tol_rho) then 223 call xc_att_xc(rho(n,3),ipol,Etmp,Atmp,Ctmp) 224 endif 225#endif 226 Amat(n,2) = Amat(n,2) + Atmp 227 Ex = Ex + qwght(n)*Etmp 228c 229 if (ldew)func(n) = func(n) + ( rho32*rho(n,2) + 230 & rho33*rho(n,3) )*P4*fac 231 endif 232 20 continue 233c 234 endif 235c 236 return 237 end 238c 239#ifndef SECOND_DERIV 240#define SECOND_DERIV 241c Compile source again for the 2nd derivative case 242#include "xc_camxlsd.F" 243#endif 244c 245#ifndef THIRD_DERIV 246#define THIRD_DERIV 247c Compile source again for the 3rd derivative case 248#include "xc_camxlsd.F" 249#endif 250