1c 2c Routines for calculating XC functional second derivatives 3c by finite difference 4c 5c BGJ - 1/02 6c 7c $Id$ 8c 9c Finite difference step size 10c 11#define STEP_SIZE 0.0001 12c 13 subroutine xc_setup_fd(tol_rho, rho, delrho, qwght, nq, ipol, GC, 14 & l_storage, i_prho, i_pdelrho, i_pAmat, i_pCmat, i_pfunc, 15 & i_qwght_copy) 16c 17 implicit none 18#include "errquit.fh" 19c 20 integer nq, ipol, l_storage, i_prho, i_pdelrho, i_pAmat, i_pCmat, 21 & i_pfunc, i_qwght_copy 22 double precision tol_rho, rho(nq,*), delrho(nq,3,*), qwght(nq) 23 logical GC 24c 25c Local variables 26c 27 integer npert, len_prho, len_pAmat, len_prho_tot, len_pAmat_tot 28 integer i 29c 30#include "mafdecls.fh" 31c 32 npert = 2 33 len_prho = 2*3*nq 34 len_pAmat = npert*npert*nq 35 if (GC) then 36 npert = 5 37 len_prho_tot = len_prho + 2*6*nq 38 len_pAmat_tot = npert*npert*nq 39 else 40 len_prho_tot = len_prho 41 len_pAmat_tot = len_pAmat 42 endif 43 if (.not. 44 & MA_Alloc_Get(MT_DBL,len_prho_tot+len_pAmat_tot+4*nq*npert, 45 & 'xc_fd',l_storage,i_prho)) 46 & call errquit('xc_setup_fd: cannot allocate storage',0, 47 & MA_ERR) 48 call dfill(len_prho_tot+len_pAmat_tot+4*nq*npert, 0d0, 49 & dbl_mb(i_prho), 1) 50 i_pAmat = i_prho + len_prho_tot 51 i_pdelrho = i_prho + len_prho 52 i_pCmat = i_pAmat + len_pAmat 53 i_pfunc = i_pAmat + len_pAmat_tot 54 i_qwght_copy = i_pfunc + 2*nq*npert 55c 56c Copy weights 57c 58 do i = 1, 2*npert 59 call dcopy(nq, qwght, 1, dbl_mb(i_qwght_copy+(i-1)*nq), 1) 60 enddo 61c 62c Set up perturbed density values 63c 64 call xc_pert_rho(tol_rho, rho, delrho, nq, ipol, GC, 65 & dbl_mb(i_prho), dbl_mb(i_pdelrho)) 66c 67 return 68 end 69c 70 subroutine xc_pert_rho(tol_rho, rho, delrho, nq, ipol, npert, 71 & prho, pdelrho) 72c 73 implicit none 74c 75 integer nq, ipol, npert 76 double precision tol_rho, rho(nq,*), delrho(nq,3,*), 77 & prho(nq,3,2,npert), pdelrho(nq,3,2,2,npert) 78c 79c Local variables 80c 81 double precision h 82 parameter (h = STEP_SIZE) 83 logical GC, DoA, DoB 84 integer i, ipert 85 double precision ra, drax, dray, draz, Gaa, Gab, Gbb, 86 & mod_drhoa, mod_drhob, cos_theta, x, mod_drhob_p, 87 & cos_theta_p, sin_theta_p 88c 89 GC = npert.eq.5 90c 91c Initial copy of density values - copy as unrestricted 92c 93 call dfill(nq*3, 0d0, prho, 1) 94 call dfill(nq*6, 0d0, pdelrho, 1) 95 if (ipol.eq.1) then 96 do i = 1, nq 97 if (rho(i,1).gt.tol_rho) then 98 ra = rho(i,1) * 0.5d0 99 prho(i,1,1,1) = rho(i,1) 100 prho(i,2,1,1) = ra 101 prho(i,3,1,1) = ra 102 endif 103 enddo 104 if (GC) then 105 do i = 1, nq 106 drax = delrho(i,1,1) * 0.5d0 107 dray = delrho(i,2,1) * 0.5d0 108 draz = delrho(i,3,1) * 0.5d0 109 pdelrho(i,1,1,1,1) = drax 110 pdelrho(i,2,1,1,1) = dray 111 pdelrho(i,3,1,1,1) = draz 112 pdelrho(i,1,2,1,1) = drax 113 pdelrho(i,2,2,1,1) = dray 114 pdelrho(i,3,2,1,1) = draz 115 enddo 116 endif 117 else 118 do i = 1, nq 119 if (rho(i,1).gt.tol_rho) then 120 prho(i,1,1,1) = rho(i,1) 121 prho(i,2,1,1) = rho(i,2) 122 prho(i,3,1,1) = rho(i,3) 123 endif 124 enddo 125 if (GC) then 126 call dcopy(nq*6, delrho, 1, pdelrho, 1) 127 endif 128 endif 129c 130c Copy to all other perturbation locations 131c 132 call dcopy(nq*3, prho, 1, prho(1,1,2,1), 1) 133 if (GC) then 134 call dcopy(nq*6, pdelrho, 1, pdelrho(1,1,1,2,1), 1) 135 endif 136 do ipert = 2, npert 137 do i = 1, 2 138 call dcopy(nq*3, prho, 1, prho(1,1,i,ipert), 1) 139 if (GC) then 140 call dcopy(nq*6, pdelrho, 1, pdelrho(1,1,1,i,ipert), 1) 141 endif 142 enddo 143 enddo 144c 145c Perturb the density parameter values - we aren't concerned if 146c some density values go negative as a result of this since the 147c functional implementation routines map negative densities to 148c a zero value for the functional and its derivatives (e.g. like 149c analytic continuation) 150c 151 do ipert = 1, 2 152 do i = 1, nq 153 if (prho(i,1,1,ipert).gt.tol_rho) then 154 prho(i,1,1,ipert) = prho(i,1,1,ipert) + h 155 prho(i,ipert+1,1,ipert) = prho(i,ipert+1,1,ipert) + h 156 prho(i,1,2,ipert) = prho(i,1,2,ipert) - h 157 prho(i,ipert+1,2,ipert) = prho(i,ipert+1,2,ipert) - h 158 endif 159 enddo 160 enddo 161 if (GC) then 162c 163c Perturb gamma values 164c 165 do i = 1, nq 166 Gaa = pdelrho(i,1,1,1,3)*pdelrho(i,1,1,1,3) 167 & + pdelrho(i,2,1,1,3)*pdelrho(i,2,1,1,3) 168 & + pdelrho(i,3,1,1,3)*pdelrho(i,3,1,1,3) 169 Gab = pdelrho(i,1,1,1,3)*pdelrho(i,1,2,1,3) 170 & + pdelrho(i,2,1,1,3)*pdelrho(i,2,2,1,3) 171 & + pdelrho(i,3,1,1,3)*pdelrho(i,3,2,1,3) 172 Gbb = pdelrho(i,1,2,1,3)*pdelrho(i,1,2,1,3) 173 & + pdelrho(i,2,2,1,3)*pdelrho(i,2,2,1,3) 174 & + pdelrho(i,3,2,1,3)*pdelrho(i,3,2,1,3) 175 mod_drhoa = sqrt(Gaa) 176 mod_drhob = sqrt(Gbb) 177c 178c Ensure perturbed gammas are always > 0, otherwise take derivative as 0 179c 180 DoA = Gaa.gt.h 181 DoB = DoA.and.Gbb.gt.h 182 if (DoB) then 183 cos_theta = Gab/(mod_drhoa*mod_drhob) 184 else 185 cos_theta = 1d0 186 endif 187c 188c Now that we have the vital parameters clear out locations 189c in the perturbed density gradient array - we only need three 190c locations to construct the proper perturbed gamma values 191c 192 pdelrho(i,1,1,1,3) = 0d0 193 pdelrho(i,2,1,1,3) = 0d0 194 pdelrho(i,3,1,1,3) = 0d0 195 pdelrho(i,1,2,1,3) = 0d0 196 pdelrho(i,2,2,1,3) = 0d0 197 pdelrho(i,3,2,1,3) = 0d0 198 pdelrho(i,1,1,1,4) = 0d0 199 pdelrho(i,2,1,1,4) = 0d0 200 pdelrho(i,3,1,1,4) = 0d0 201 pdelrho(i,1,2,1,4) = 0d0 202 pdelrho(i,2,2,1,4) = 0d0 203 pdelrho(i,3,2,1,4) = 0d0 204 pdelrho(i,1,1,1,5) = 0d0 205 pdelrho(i,2,1,1,5) = 0d0 206 pdelrho(i,3,1,1,5) = 0d0 207 pdelrho(i,1,2,1,5) = 0d0 208 pdelrho(i,2,2,1,5) = 0d0 209 pdelrho(i,3,2,1,5) = 0d0 210 if (DoA) then 211c 212c Perturb Gaa up 213c 214 x = mod_drhoa*cos_theta 215 pdelrho(i,1,1,1,3) = x 216 pdelrho(i,2,1,1,3) = sqrt(Gaa-x*x+h) 217 pdelrho(i,1,2,1,3) = mod_drhob 218 endif 219 if (DoB) then 220c 221c Perturb Gab up 222c 223 x = (Gab+h)/mod_drhob 224 pdelrho(i,1,1,1,4) = x 225 pdelrho(i,2,1,1,4) = sqrt(Gaa-x*x) 226 pdelrho(i,1,2,1,4) = mod_drhob 227c 228c Perturb Gbb up 229c 230 mod_drhob_p = sqrt(Gbb+h) 231 cos_theta_p = (mod_drhob/mod_drhob_p)*cos_theta 232 sin_theta_p = sqrt(1d0-cos_theta_p*cos_theta_p) 233 pdelrho(i,1,1,1,5) = mod_drhoa*cos_theta_p 234 pdelrho(i,2,1,1,5) = mod_drhoa*sin_theta_p 235 pdelrho(i,1,2,1,5) = mod_drhob_p 236 endif 237 if (DoA) then 238c 239c Perturb Gaa back 240c 241 x = mod_drhoa*cos_theta 242 pdelrho(i,1,1,2,3) = x 243 pdelrho(i,2,1,2,3) = sqrt(Gaa-x*x-h) 244 pdelrho(i,1,2,2,3) = mod_drhob 245 endif 246 if (DoB) then 247c 248c Perturb Gab back 249c 250 x = (Gab-h)/mod_drhob 251 pdelrho(i,1,1,2,4) = x 252 pdelrho(i,2,1,2,4) = sqrt(Gaa-x*x) 253 pdelrho(i,1,2,2,4) = mod_drhob 254c 255c Perturb Gbb back 256c 257 mod_drhob_p = sqrt(Gbb-h) 258 cos_theta_p = (mod_drhob/mod_drhob_p)*cos_theta 259 sin_theta_p = sqrt(1d0-cos_theta_p*cos_theta_p) 260 pdelrho(i,1,1,2,5) = mod_drhoa*cos_theta_p 261 pdelrho(i,2,1,2,5) = mod_drhoa*sin_theta_p 262 pdelrho(i,1,2,2,5) = mod_drhob_p 263 endif 264 enddo 265 endif 266c 267 return 268 end 269c 270 subroutine xc_make_fd(Amat2, Cmat2, nq, GC, pAmat, pCmat) 271c 272 implicit none 273c 274 integer nq 275 double precision Amat2(nq,*), Cmat2(nq,*), 276 & pAmat(nq,2,2,*), pCmat(nq,3,2,*) 277 logical GC 278c 279c Local variables 280c 281 double precision h, r2h 282 parameter (h = STEP_SIZE) 283 integer i, npert, ipert 284c 285#include "dft2drv.fh" 286c 287 if (GC) then 288 npert = 5 289 else 290 npert = 2 291 endif 292 r2h = 0.5d0*h 293c 294c Construct finite differences in the temporary arrays 295c 296 do ipert = 1, npert 297 do i = 1, nq 298 pAmat(i,1,1,ipert) = (pAmat(i,1,1,ipert) 299 & -pAmat(i,1,2,ipert))*r2h 300 pAmat(i,2,1,ipert) = (pAmat(i,2,1,ipert) 301 & -pAmat(i,2,2,ipert))*r2h 302 enddo 303 if (GC) then 304 do i = 1, nq 305 pCmat(i,1,1,ipert) = (pCmat(i,1,1,ipert) 306 & -pCmat(i,1,2,ipert))*r2h 307 pCmat(i,2,1,ipert) = (pCmat(i,2,1,ipert) 308 & -pCmat(i,2,2,ipert))*r2h 309 pCmat(i,3,1,ipert) = (pCmat(i,3,1,ipert) 310 & -pCmat(i,3,2,ipert))*r2h 311 enddo 312 endif 313 enddo 314c 315c Now scatter to Amat2 and Cmat2 316c 317 call dcopy(nq, pAmat(1,1,1,1), 1, Amat2(1,D2_RA_RA), 1) 318 call dcopy(nq, pAmat(1,2,1,1), 1, Amat2(1,D2_RA_RB), 1) 319c Or could be as below, etc., etc. 320c call dcopy(nq, pAmat(1,1,1,2), 1, Amat2(1,D2_RA_RB), 1) 321 call dcopy(nq, pAmat(1,2,1,2), 1, Amat2(1,D2_RB_RB), 1) 322 if (GC) then 323 call dcopy(nq, pCmat(1,D1_GAA,1,1), 1, Cmat2(1,D2_RA_GAA), 1) 324c Or could be as below, etc., etc. 325c call dcopy(nq, pAmat(1,1,1,3), 1, Cmat2(1,D2_RA_GAA), 1) 326 call dcopy(nq, pCmat(1,D1_GAA,1,2), 1, Cmat2(1,D2_RB_GAA), 1) 327 call dcopy(nq, pCmat(1,D1_GAA,1,3), 1, Cmat2(1,D2_GAA_GAA), 1) 328c 329 call dcopy(nq, pCmat(1,D1_GAB,1,1), 1, Cmat2(1,D2_RA_GAB), 1) 330 call dcopy(nq, pCmat(1,D1_GAB,1,2), 1, Cmat2(1,D2_RB_GAB), 1) 331 call dcopy(nq, pCmat(1,D1_GAB,1,3), 1, Cmat2(1,D2_GAA_GAB), 1) 332 call dcopy(nq, pCmat(1,D1_GAB,1,4), 1, Cmat2(1,D2_GAB_GAB), 1) 333c 334 call dcopy(nq, pCmat(1,D1_GBB,1,1), 1, Cmat2(1,D2_RA_GBB), 1) 335 call dcopy(nq, pCmat(1,D1_GBB,1,2), 1, Cmat2(1,D2_RB_GBB), 1) 336 call dcopy(nq, pCmat(1,D1_GBB,1,3), 1, Cmat2(1,D2_GAA_GBB), 1) 337 call dcopy(nq, pCmat(1,D1_GBB,1,4), 1, Cmat2(1,D2_GAB_GBB), 1) 338 call dcopy(nq, pCmat(1,D1_GBB,1,5), 1, Cmat2(1,D2_GBB_GBB), 1) 339 endif 340c 341 return 342 end 343