1c $Id$ 2 subroutine xc_xdm_init(rtdb,iixdm_v,iixdm_ml) 3 implicit none 4 5#include "rtdb.fh" 6#include "global.fh" 7#include "mafdecls.fh" 8#include "cdft.fh" 9 10 integer rtdb ! input 11 integer natoms ! input 12 integer iixdm_v, iixdm_ml ! output 13 14 logical savedvml 15 16c xdm common 17 integer nxdm 18 integer ixdm_v, ixdm_ml 19 integer lxdm_v, lxdm_ml 20 common /xdmd/ nxdm, ixdm_v, lxdm_v, ixdm_ml, lxdm_ml 21 data ixdm_v / 0 / 22 data ixdm_ml / 0 / 23 24 if (.not. rtdb_get(rtdb, 'dft:xdm', mt_int, 1, 25 & lxdm)) lxdm = 0 26 if(lxdm.eq.0) return 27 28c are the volumes and moments available from a previous iteration? 29 if (.not. rtdb_get(rtdb,'dft:xdmsave', mt_log, 1, savedvml)) 30 & savedvml = .false. 31 32c get memory for volume, alpha, ml 33 nxdm = ntypes 34 35 if (savedvml) then 36c check if v and ml have been allocated 37 if(ixdm_v.eq.0.or.ixdm_ml.eq.0) then 38 if(ga_nodeid().eq.0) then 39 write(6,*) ' xdm: v and ml not allocated ' 40 write(6,*) ' xdm: resetting savedvml ' 41 endif 42 savedvml=.false. 43 if (.not. rtdb_put(rtdb,'dft:xdmsave', mt_log, 1, savedvml)) 44 c call errquit('xc_xdm_init: cant rtdb_put',0,0) 45 else 46 endif 47 endif 48 if (.not.savedvml) then 49 if (.not.MA_alloc_get(mt_dbl, ntypes, 50 & 'xdm_v', lxdm_v, ixdm_v)) then 51 call errquit('xc_xdm_init: cant alloc xdm_v',0,0) 52 endif 53 if (.not.MA_alloc_get(mt_dbl, 3*ntypes, 54 & 'xdm_ml', lxdm_ml, ixdm_ml)) then 55 call errquit('xc_xdm_init: cant alloc xdm_ml',0,0) 56 endif 57 endif 58 59 iixdm_v = ixdm_v 60 iixdm_ml = ixdm_ml 61 62 return 63 end 64 65 subroutine xc_xdm(rtdb,g_dens,g_vxc,n,nexc,exdm,fxdm,v,ml,what) 66 implicit none 67 68#include "errquit.fh" 69#include "rtdb.fh" 70#include "global.fh" 71#include "mafdecls.fh" 72#include "cdft.fh" 73#include "stdio.fh" 74#include "msgids.fh" 75#include "util_params.fh" 76 77 integer rtdb, g_dens(2), g_vxc(4) ! input 78 integer n, nexc ! input 79 double precision exdm, fxdm(3,n) ! output 80 double precision v(ntypes) ! inout, atomic volumes 81 double precision ml(3,ntypes) ! inout, moments 82 character*(*) what ! energy/forces flag 83 84 integer i, j, me, nmu 85 double precision rho_n, fxx 86 double precision dum1, ftemp(3,n) 87 logical setparam, savedvml 88 integer izz, ixx, jxx 89 90c coefficents and vdw radii 91 double precision c6(ntypes,ntypes), c8(ntypes,ntypes) 92 double precision c10(ntypes,ntypes), rc(ntypes,ntypes) 93 double precision rvdw(ntypes,ntypes) 94 double precision a(ntypes) 95 96c geometry get 97 logical geom_cent_get 98 external geom_cent_get 99 character*16 tag 100 double precision x1(3), x2(3), q, r 101 102c input parameters 103 double precision a1, a2 104 logical onlyc, varc 105 106c xdm common 107 integer nxdm 108 integer ixdm_v, ixdm_ml, lxdm_v, lxdm_ml 109 common /xdmd/ nxdm, ixdm_v, lxdm_v, ixdm_ml, lxdm_ml 110 111c free volumes and polarizabilities 112C DKH LSDA/UGBS Free Atomic Volumes 113 double precision frevol(103) 114 DATA (FREVOL(I),I=1,103)/ 115 + 9.194D0, 4.481D0, 91.957D0, 61.357D0, 49.813D0, 36.728D0, 116 + 27.633D0, 23.517D0, 19.322D0, 15.950D0, 109.359D0, 103.064D0, 117 + 120.419D0, 104.229D0, 86.782D0, 77.133D0, 66.372D0, 57.336D0, 118 + 203.093D0, 212.202D0, 183.101D0, 162.278D0, 143.250D0, 119 + 108.209D0, 123.098D0, 105.735D0, 92.944D0, 83.794D0, 75.750D0, 120 + 81.177D0, 118.371D0, 116.334D0, 107.474D0, 103.221D0, 95.111D0, 121 + 87.605D0, 248.772D0, 273.748D0, 249.211D0, 223.801D0, 122 + 175.809D0, 156.831D0, 160.042D0, 136.654D0, 127.754D0, 123 + 97.024D0, 112.778D0, 121.627D0, 167.906D0, 172.030D0, 124 + 165.500D0, 163.038D0, 153.972D0, 146.069D0, 341.992D0, 125 + 385.767D0, 343.377D0, 350.338D0, 334.905D0, 322.164D0, 126 + 310.337D0, 299.537D0, 289.567D0, 216.147D0, 268.910D0, 127 + 259.838D0, 251.293D0, 243.174D0, 235.453D0, 228.284D0, 128 + 229.617D0, 209.971D0, 197.541D0, 183.236D0, 174.685D0, 129 + 164.139D0, 150.441D0, 135.765D0, 125.297D0, 131.258D0, 130 + 185.769D0, 195.671D0, 193.036D0, 189.142D0, 185.919D0, 131 + 181.089D0, 357.787D0, 407.283D0, 383.053D0, 362.099D0, 132 + 346.565D0, 332.462D0, 319.591D0, 308.095D0, 297.358D0, 133 + 300.572D0, 275.792D0, 266.317D0, 257.429D0, 209.687D0, 134 + 203.250D0, 230.248D0, 236.878D0/ 135 double precision frepol(103) 136C free atom polariz.: CRC Handbook of Chemistry and Physics, 88th Ed. 137 DATA (FREPOL(I),I=1,102)/ 138 + 0.6668D0, 0.2051D0, 24.3300D0, 5.6000D0, 3.0300D0, 139 + 1.7600D0, 1.1000D0, 0.8020D0, 0.5570D0, 0.3956D0, 140 + 24.1100D0, 10.6000D0, 6.8000D0, 5.3800D0, 3.6300D0, 141 + 2.9000D0, 2.1800D0, 1.6411D0, 43.4000D0, 22.8000D0, 142 + 17.8000D0, 14.6000D0, 12.4000D0, 11.6000D0, 9.4000D0, 143 + 8.4000D0, 7.5000D0, 6.8000D0, 6.2000D0, 5.7500D0, 144 + 8.1200D0, 6.0700D0, 4.3100D0, 3.7700D0, 3.0500D0, 145 + 2.4844D0, 47.3000D0, 27.6000D0, 22.7000D0, 17.9000D0, 146 + 15.7000D0, 12.8000D0, 11.4000D0, 9.6000D0, 8.6000D0, 147 + 4.8000D0, 7.2000D0, 7.3600D0, 10.2000D0, 7.7000D0, 148 + 6.6000D0, 5.5000D0, 5.3500D0, 4.0440D0, 59.4200D0, 149 + 39.7000D0, 31.1000D0, 29.6000D0, 28.2000D0, 31.4000D0, 150 + 30.1000D0, 28.8000D0, 27.7000D0, 23.5000D0, 25.5000D0, 151 + 24.5000D0, 23.6000D0, 22.7000D0, 21.8000D0, 21.0000D0, 152 + 21.9000D0, 16.2000D0, 13.1000D0, 11.1000D0, 9.7000D0, 153 + 8.5000D0, 7.6000D0, 6.5000D0, 5.8000D0, 5.0200D0, 154 + 7.6000D0, 6.8000D0, 7.4000D0, 6.8000D0, 6.0000D0, 155 + 5.3000D0, 48.6000D0, 38.3000D0, 32.1000D0, 32.1000D0, 156 + 25.4000D0, 24.9000D0, 24.8000D0, 24.5000D0, 23.3000D0, 157 + 23.0000D0, 22.7000D0, 20.5000D0, 19.7000D0, 23.8000D0, 158 + 18.2000D0, 17.5000D0/ 159c 160c atomic densities 161 integer ngau, k, ii 162 double precision pi, vfree, rmid, qq, rr, h, rwei, rhofree 163 parameter (ngau = 251) 164c 165 pi=acos(-1d0) 166c 167 168c read input parameters 169 setparam = .not.rtdb_get(rtdb,'dft:xdm_a1',mt_dbl,1,a1) 170 setparam = setparam .or. 171 & .not.rtdb_get(rtdb,'dft:xdm_a2',mt_dbl,1,a2) 172 if (setparam) call xc_xdm_setdefaults(a1,a2) 173 174 if (.not.rtdb_get(rtdb,'dft:xdm_onlyc',mt_log,1,onlyc)) 175 & onlyc = .false. 176 if (.not.rtdb_get(rtdb,'dft:xdm_varc',mt_log,1,varc)) 177 & varc = .false. 178 179 me = ga_nodeid() 180 181c redo the volume and moment calculation? 182 if (.not. rtdb_get(rtdb,'dft:xdmsave', mt_log, 1, savedvml)) 183 & savedvml = .false. 184 if (varc .and. what.eq.'energy') then 185 savedvml = .false. 186 endif 187 188c header and convert a2 to au 189 if (me.eq.0 .and. what.eq.'energy' .and..not.savedvml) then 190 write (luout,*) 191 write (luout,'("+ XDM: a1 = ",F12.5)') a1 192 write (luout,'(" a2 (ang) = ",F12.5)') a2 193 write (luout,'(" onlyc? = ",L12)') onlyc 194 endif 195 a2 = a2 /cau2ang 196 197c do the v and ml calculation 198 if (.not.savedvml) then 199c initialize common values 200 call dfill(ntypes, 0d0, v, 1) 201 call dfill(3*ntypes, 0d0, ml, 1) 202 203c calculate volumes and moments using the dft grid 204 rho_n = 0 205 exdm = 0 206 207 call grid_quadv0_gen(rtdb,g_dens,g_vxc,nexc,rho_n,exdm,1,6, 208 & dum1,.false.,.false.) 209 call ga_dgop(msg_xdm1,v,ntypes,'+') 210 call ga_dgop(msg_xdm3,ml,3*ntypes,'+') 211 call ga_sync() 212 213 do i = 1, ntypes 214c skip ghost atoms 215 izz = znuc_atom_type(i) 216 if (izz.eq.0) cycle 217 218c Only half of it if spin-polarized calculation 219 if (ipol .eq. 2) then 220 v(i) = v(i) / 2d0 221 ml(1,i) = ml(1,i) / 2d0 222 ml(2,i) = ml(2,i) / 2d0 223 ml(3,i) = ml(3,i) / 2d0 224 endif 225 226c Divide by the multiplicity 227 nmu = 0 228 do j = 1, n 229 if (iatype(j).eq.i) nmu = nmu + 1 230 enddo 231 v(i) = v(i) / nmu 232 ml(1,i) = ml(1,i) / nmu 233 ml(2,i) = ml(2,i) / nmu 234 ml(3,i) = ml(3,i) / nmu 235 enddo 236 237c Calculate volumes and moments only once, then save for future use. 238 if (.not. rtdb_put(rtdb,'dft:xdmsave', mt_log, 1, .true.)) 239 & call errquit('xc_xdm: rtdb_put dft:xdmsave failed', 0, 240 & RTDB_ERR) 241 endif 242c fill polarizabilities vector 243 do i = 1, ntypes 244c skip ghost atoms 245 izz = znuc_atom_type(i) 246 if (izz.eq.0) cycle 247 a(i) = frepol(izz)/0.148184D0 * min(v(i)/frevol(izz),1d0) 248 enddo 249 250 if (me .eq. 0 .and..not.savedvml) then 251 write (luout,*) 252 write (luout,'("+ XDM: volume and moments")') 253 write (luout,'(" All results in atomic units.")') 254 write (luout,'(" i V alpha M1 M2 M3")') 255 do i = 1, n 256 ixx = iatype(i) 257 izz = znuc_atom_type(ixx) 258 if (izz.eq.0) cycle 259 write (luout,'(I3,5(X,F18.10))') i, v(ixx), a(ixx), 260 & ml(1,ixx), ml(2,ixx), ml(3,ixx) 261 end do 262 write (luout,*) 263 endif 264 265c calculate interaction coefficients 266 if (me .eq. 0 .and. what.eq.'energy' .and..not.savedvml) then 267 write (luout,'("+ XDM: dispersion coefficients")') 268 write (luout,'(" All results in atomic units.")') 269 write (luout,'(" i j C6 C8 C10 Rc Rvdw")') 270 endif 271 do i = 1, ntypes 272 izz = znuc_atom_type(i) 273 if (izz.eq.0) cycle 274 275 do j = 1, i 276 izz = znuc_atom_type(j) 277 if (izz.eq.0) cycle 278 c6(i,j) = a(i)*a(j)*ml(1,i)*ml(1,j) / 279 & (ml(1,i)*a(j) + ml(1,j)*a(i)) 280 c6(j,i) = c6(i,j) 281 c8(i,j) = 3d0/2d0 * (a(i)*a(j)*(ml(1,i)* 282 & ml(2,j)+ml(2,i)*ml(1,j))) / 283 & (ml(1,i)*a(j)+ml(1,j)*a(i)) 284 c8(j,i) = c8(i,j) 285 c10(i,j) = 2 * a(i)*a(j) * (ml(1,i)* 286 & ml(3,j) + ml(3,i)*ml(1,j)) / 287 & (ml(1,i)*a(j) + ml(1,j)*a(i)) + 288 & 21d0/5d0 * a(i)*a(j)* 289 & ml(2,i)*ml(2,j) / (a(j)*ml(1,i)+ 290 & a(i)*ml(1,j)) 291 c10(j,i) = c10(i,j) 292 rc(i,j) = (dsqrt(c8(i,j)/c6(i,j)) + 293 & dsqrt(c10(i,j)/c8(i,j)) + 294 & (c10(i,j)/c6(i,j))**(0.25d0)) / 3 295 rc(j,i) = rc(i,j) 296 rvdw(i,j) = a1 * rc(i,j) + a2 297 rvdw(j,i) = rvdw(i,j) 298 enddo 299 end do 300 301 if (what.eq.'energy'.and..not.savedvml) then 302 do i = 1, n 303 do j = 1, i 304 if (me .eq. 0) then 305 ixx = iatype(i) 306 izz = znuc_atom_type(ixx) 307 if (izz.eq.0) cycle 308 jxx = iatype(j) 309 izz = znuc_atom_type(jxx) 310 if (izz.eq.0) cycle 311 write (luout,'(I3,X,I3,1p,5(X,E18.10))') i, j, c6(ixx 312 & ,jxx),c8(ixx,jxx), c10(ixx,jxx), rc(ixx,jxx) 313 & ,rvdw(ixx,jxx) 314 endif 315 enddo 316 enddo 317 endif 318 319 if (onlyc) return 320 321c Sum the dispersion energy 322 exdm = 0d0 323 do i = 1, n 324 do j = 1,3 325 ftemp(j,i) = 0d0 326 enddo 327 enddo 328 329 do i = 1, n 330 ixx = iatype(i) 331 izz = znuc_atom_type(ixx) 332 if (izz.eq.0) cycle 333 if (.not.geom_cent_get(geom,i,tag,x1,q)) 334 & call errquit('xc_xdm: geom_cent_get failed',geom,GEOM_ERR) 335 336 do j = 1, i-1 337 jxx = iatype(j) 338 izz = znuc_atom_type(jxx) 339 if (izz.eq.0) cycle 340 if (.not.geom_cent_get(geom,j,tag,x2,q)) 341 & call errquit('xc_xdm: geom_cent_get failed',geom 342 & ,GEOM_ERR) 343 344 r = dsqrt((x1(1)-x2(1))**2+(x1(2)-x2(2))**2+ 345 & (x1(3)-x2(3))**2) 346 if (what.eq.'energy') then 347 exdm = exdm - c6(ixx,jxx) / (rvdw(ixx,jxx)**6 + r**6) 348 - - c8(ixx,jxx) / (rvdw(ixx,jxx)**8 + r**8) 349 - - c10(ixx,jxx) / (rvdw(ixx,jxx)**10 + r**10) 350 else 351 fxx = -( 352 + 6 * c6(ixx,jxx) * r**4 / (rvdw(ixx,jxx)**6 + r**6)**2+ 353 + 8 * c8(ixx,jxx) * r**6 / (rvdw(ixx,jxx)**8 + r**8)**2+ 354 + 10 * c10(ixx,jxx) * r**8 / (rvdw(ixx,jxx)**10 + r**10)**2) 355 356 do k = 1, 3 357 ftemp(k,i) = ftemp(k,i) + (x2(k) - x1(k)) * fxx 358 ftemp(k,j) = ftemp(k,j) - (x2(k) - x1(k)) * fxx 359 enddo 360 endif 361 enddo 362 enddo 363 if(what.eq.'forces') then 364c bit needed for forces 365 do i = 1, n 366 do k = 1, 3 367 fxdm(k,i) = fxdm(k,i) + ftemp(k,i) 368 end do 369 end do 370 endif 371 372 if (me .eq. 0) then 373 if (what.eq.'energy') then 374 write (luout, 375 & '("+ XDM dispersion energy (Hy) = ",1p,E22.12)') 376 & exdm 377 else 378 do i = 1, n 379 write (luout, 380 & '(" FXDM(",I2.2,")=",1p,3(E20.12,X),"a.u.")') 381 & i, ftemp(1:3,i) 382 end do 383 endif 384 write (luout,*) 385 endif 386 call ga_sync() 387 return 388 end 389 390 subroutine xc_eval_xdm(rho,delrho,lap,nq, 391 & qxyz,qwght,ttau, 392 & natoms,xyz,v,ml) 393 implicit none 394 395#include "stdio.fh" 396#include "errquit.fh" 397#include "global.fh" 398#include "cdft.fh" 399#include "mafdecls.fh" 400 401 double precision rho(nq,ipol*(ipol+1)/2) ! input, electron density 402 double precision delrho(nq,3,ipol) ! input, gradient of the electron density 403 double precision lap(nq,ipol*(ipol+1)/2) ! input, laplacian 404 integer nq ! input, no. of quadrature weights 405 double precision qxyz(3,nq) ! input, quadrature node positions 406 double precision qwght(nq) ! input, quadrature weights 407 double precision ttau(nq,ipol) ! input, ked 408 integer natoms ! input, no. of atoms 409 double precision xyz(3,natoms) ! input, atomic positions 410 double precision v(ntypes) ! inout, atomic volumes 411 double precision a(ntypes) ! inout, atomic polarizabilities 412 double precision ml(3,ntypes) ! inout, moments 413 414 integer iz, it 415 416c parameters 417 double precision pi 418c 419 integer i, j, k, iat, isp 420 double precision wei, rhoat, rhoi, r, x, y, z, vsum 421 double precision rhot, rhos, grho, taus, laps 422 double precision ds, qs, rhs, xroot, xshift, xold 423 double precision expx, gx, fx, ffx, db, db2, db3 424 double precision r2, r3 425 double precision eps,thresh,mtwo3rds 426 parameter(eps=1d-40,thresh=1d-12,mtwo3rds=-2d0/3d0) 427 428c atomic densities 429 double precision rhop(nq), b(nq) 430 double precision xdm_rhop 431 external xdm_rhop 432c 433 pi=acos(-1d0) 434 435 do i = 1, nq 436c calculate promolecular density on the grid 437 rhop(i) = 0d0 438 do j = 1, natoms 439 x = qxyz(1,i) - xyz(1,j) 440 y = qxyz(2,i) - xyz(2,j) 441 z = qxyz(3,i) - xyz(3,j) 442 r = dsqrt(x*x + y*y + z*z) 443 444 iz = znuc_atom_type(iatype(j)) 445c iz(j) == 0 is a ghost atom (bsse) 446 if (iz.ge.1 .or. iz.le.94) then 447 rhop(i) = rhop(i) + xdm_rhop(iz,r) 448 elseif (iz.lt.0 .or. iz.gt.94) then 449 call errquit('xc_xdm: wrong atomic number',j,iz) 450 endif 451 enddo 452 end do 453 454 do isp = 1, ipol 455 do i = 1, nq 456c calculate dipole moment 457 if (ipol .eq. 1) then 458 rhot = max(rho(i,1),eps) 459 rhos = max(rhot / 2d0,eps/2d0) 460 grho = dsqrt(delrho(i,1,1)**2 + delrho(i,2,1)**2 + 461 + delrho(i,3,1)**2) / 2d0 462 taus = ttau(i,1) 463 laps = lap(i,1) / 2d0 464 else 465 rhot = max(rho(i,2),eps) 466 if (isp .eq. 1) then 467 rhos = rho(i,2) - rho(i,3) 468 grho = dsqrt((delrho(i,1,1)-delrho(i,1,2))**2 + 469 + (delrho(i,2,1)-delrho(i,2,2))**2 + 470 + (delrho(i,3,1)-delrho(i,3,2))**2) 471 taus = (ttau(i,1)-ttau(i,2))*2d0 472 laps = (lap(i,1)-lap(i,2)) 473 else 474 rhos = rho(i,3) 475 grho = dsqrt(delrho(i,1,2)**2 + delrho(i,2,2)**2 + 476 + delrho(i,3,2)**2) 477 taus = ttau(i,2)*2d0 478 laps = lap(i,2) 479 endif 480 rhos = max(rhos,0.5*eps) 481 endif 482 483 if(abs(rhos).gt.eps) then 484 ds = taus - 0.25d0 * grho**2 / rhos 485 qs = 1d0/6d0 * (laps - 2d0 * ds) 486 rhs = 2d0/3d0*pi**(2d0/3d0)*(rhos)**(5d0/3d0)/qs 487 else 488 rhs= 0d0 489 endif 490 if (rhs > 0d0) then 491 xroot = 3d0 492 xshift = 1d0 493 do while ((xroot*exp(-2d0*xroot/3d0))/(xroot-2d0)<rhs) 494 xshift = xshift * 0.1d0 495 xroot = 2d0 + xshift 496 end do 497 else 498 xroot = 1d0 499 xshift = 1d0 500 do while ((xroot*exp(-2d0*xroot/3d0))/(xroot-2d0)>rhs) 501 xshift = xshift * 0.1d0 502 xroot = 2d0 - xshift 503 end do 504 end if 505 if(abs(xroot).gt.eps) then 506 xold = xroot + 1d0 507 do while (abs(xroot - xold) > thresh) 508 xold = xroot 509 expx = exp(-2d0 * xroot / 3d0) 510 gx = (xroot * expx) / (xroot - 2d0) 511 fx = gx - rhs 512 ffx = gx * (1d0 / xroot - 2d0/3d0 - 1d0 / (xroot - 2d0)) 513 xroot = xroot - fx / ffx 514 end do 515 endif 516 b(i) = xroot * (exp(-xroot) / (8d0*pi*rhos))**(1d0/3d0) 517 enddo 518 519c calculate contribution to atomic volumes, alphas and moments 520c skip ghost atoms 521 do iat = 1, natoms 522 it = iatype(iat) 523 iz = znuc_atom_type(it) 524 if (iz.eq.0) cycle 525 526 vsum = 0d0 527 do i = 1, nq 528 if (ipol .eq. 1) then 529 rhot = rho(i,1) 530 else 531 rhot = rho(i,2) 532 endif 533 534c atomic density of iat at the grid node 535 x = qxyz(1,i) - xyz(1,iat) 536 y = qxyz(2,i) - xyz(2,iat) 537 z = qxyz(3,i) - xyz(3,iat) 538 r = dsqrt(x*x + y*y + z*z) 539 rhoi = xdm_rhop(iz,r) 540 wei = rhoi / max(rhop(i),eps) * qwght(i) 541 542c contribution to atomic volume 543 vsum = vsum + r**3 * rhot * wei 544 545c contribution to atomic moments 546 db = max(r-b(i),0d0) 547 db2 = db * db 548 db3 = db2 * db 549 r2 = r * r 550 r3 = r2 * r 551 ml(1,it) = ml(1,it) + 552 + wei * rhot * (r - db)**2 553 ml(2,it) = ml(2,it) + 554 + wei * rhot * (r2 - db2)**2 555 ml(3,it) = ml(3,it) + 556 + wei * rhot * (r3 - db3)**2 557 enddo 558 v(it) = v(it) + vsum 559 enddo 560 enddo 561 end 562 563 subroutine xc_xdm_cleanup(rtdb) 564 565#include "mafdecls.fh" 566#include "rtdb.fh" 567c 568 integer rtdb 569c xdm common 570 logical savedvml 571 integer nxdm, lxdm 572 integer ixdm_v, ixdm_ml, lxdm_v, lxdm_ml 573 common /xdmd/ nxdm, ixdm_v, lxdm_v, ixdm_ml, lxdm_ml 574 575 if (.not. rtdb_get(rtdb, 'dft:xdm', mt_int, 1, 576 & lxdm)) lxdm = 0 577 if(lxdm.eq.0) return 578 579 580 if (.not. rtdb_get(rtdb,'dft:xdmsave', mt_log, 1, savedvml)) 581 & savedvml = .false. 582 if (.not.savedvml) then 583 if (.not.MA_free_heap(lxdm_ml)) 584 & call errquit('xc_xdm_cleanup: cannot clean ml',1,MA_ERR) 585 if (.not.MA_free_heap(lxdm_v)) 586 & call errquit('xc_xdm_cleanup: cannot clean v',1,MA_ERR) 587 endif 588 589 end 590 591 subroutine xc_xdm_setdefaults(a1,a2) 592#include "cdft.fh" 593c X C a1 a2 594c becke86b pbe96 0.82 1.16 595 596 double precision a1, a2 ! output 597 598 if (abs(xfac(55)).gt.1d-2 .and. abs(cfac(12)).gt.1d-2) then 599 write (*,*) "WARNING: B86b!" 600 a1 = 0.82d0 601 a2 = 1.16d0 602 else 603 write(6,*)'WARNING:' 604 write(6,*)'untested DF in XDM!!' 605 a1 = 0.82d0 606 a2 = 1.16d0 607 endif 608 609 end 610 611 function xdm_rhop(iz,r) 612 613 integer iz ! input atomic number 614 double precision r ! input distance to atom 615 double precision xdm_rhop ! output promolecular density 616 617c LDA atomic densities 618c rho(r) = sum_i C_i * exp(-r/zeta_i) 619 double precision rc1(94), rc2(94), rc3(94), rc4(94), 620 + rc5(94), rc6(94), rc7(94), zeta1(94), zeta2(94), 621 + zeta3(94), zeta4(94), zeta5(94), zeta6(94), zeta7(94) 622 DATA (rc1(I),I=1,94)/ 623 + 2.38247581d-01, 2.65077191d+00, 1.35064937d+01, 3.52085244d+01, 624 + 7.05975356d+01, 1.25179296d+02, 2.02122544d+02, 3.05256892d+02, 625 + 4.38665900d+02, 6.07145508d+02, 8.39844590d+02, 1.10364452d+03, 626 + 1.42251250d+03, 1.79552242d+03, 2.23091423d+03, 2.73461374d+03, 627 + 3.31278653d+03, 3.97189405d+03, 4.70704295d+03, 5.55728188d+03, 628 + 6.51013713d+03, 7.57599544d+03, 8.76508298d+03, 1.00881454d+04, 629 + 1.15553877d+04, 1.31780299d+04, 1.49681265d+04, 1.69392910d+04, 630 + 1.91364487d+04, 2.14808938d+04, 2.42414395d+04, 2.69764334d+04, 631 + 3.00224115d+04, 3.33619034d+04, 3.70058615d+04, 4.09716964d+04, 632 + 4.52399350d+04, 1.44691990d+04, 1.75811295d+04, 1.98707444d+04, 633 + 2.74722280d+04, 3.04390736d+04,-2.32888935d+00, 3.87816243d+04, 634 + 4.37829818d+04, 4.98625265d+04, 5.57415421d+04, 3.56643213d+03, 635 + 7.17222714d+04,-1.14919886d+01, 8.97267971d+04, 1.00937307d+05, 636 + 1.13620429d+05, 1.27857026d+05, 1.40502758d+05, 1.57775964d+05, 637 + 1.78064645d+05, 2.00731613d+05, 2.26118301d+05, 2.54550489d+05, 638 + 2.86238347d+05, 3.21681109d+05, 3.61267800d+05, 4.05477009d+05, 639 + 4.54870353d+05, 5.09923518d+05, 5.71460462d+05, 6.40210146d+05, 640 + 7.16876859d+05, 8.02519272d+05, 8.94824315d+05, 1.00284493d+06, 641 + 1.12164004d+06, 1.25799024d+06, 1.41707468d+06, 1.56657370d+06, 642 + 1.75110216d+06, 1.95685098d+06, 2.18876560d+06, 2.45050333d+06, 643 + 2.73584853d+06, 3.06874608d+06, 3.44077168d+06, 3.85778178d+06, 644 + 4.69176541d+06, 5.33909792d+06, 6.13639622d+06, 6.98851575d+06, 645 + 7.92638571d+06, 9.01616115d+06, 1.02412946d+07, 1.16535915d+07, 646 + 1.32723280d+07, 1.51679680d+07/ 647 DATA (zeta1(I),I=1,94)/ 648 + 5.94389750d-01, 3.63728528d-01, 1.72139700d-01, 1.21899504d-01, 649 + 1.02003256d-01, 8.49661099d-02, 7.32065807d-02, 6.45649696d-02, 650 + 5.78507920d-02, 5.22970367d-02, 4.14514521d-02, 3.73253761d-02, 651 + 3.46867228d-02, 3.16764946d-02, 2.92112390d-02, 2.71143558d-02, 652 + 2.52900857d-02, 2.36835198d-02, 2.20789393d-02, 2.08540286d-02, 653 + 1.97256033d-02, 1.87075897d-02, 1.77753759d-02, 1.69434569d-02, 654 + 1.61518812d-02, 1.54331276d-02, 1.47661483d-02, 1.41466847d-02, 655 + 1.36293029d-02, 1.30268084d-02, 1.27388256d-02, 1.20747615d-02, 656 + 1.15587091d-02, 1.10931067d-02, 1.06586168d-02, 1.02470835d-02, 657 + 9.84539854d-03, 6.81246963d-04, 8.60318829d-04, 7.16123706d-04, 658 + 1.82083380d-03, 1.49186880d-03, 8.37924126d-07, 1.20208187d-03, 659 + 1.09125302d-03, 1.05045931d-03, 9.15785226d-04, 4.45100103d-03, 660 + 8.42267282d-04, 4.55784641d-06, 6.49666842d-04, 6.03861716d-04, 661 + 5.65447958d-04, 5.33042747d-04, 3.92130359d-04, 3.61775974d-04, 662 + 3.55549092d-04, 3.48815329d-04, 3.42580721d-04, 3.36626058d-04, 663 + 3.30233653d-04, 3.24004067d-04, 3.17860234d-04, 3.12193286d-04, 664 + 3.05756593d-04, 3.00104329d-04, 2.94369895d-04, 2.88669044d-04, 665 + 2.83327379d-04, 2.78093140d-04, 3.59124379d-04, 3.01791198d-04, 666 + 2.88675155d-04, 2.57064717d-04, 1.32937491d-04, 2.70674086d-04, 667 + 2.64755263d-04, 2.60031593d-04, 2.52814234d-04, 2.44193852d-04, 668 + 2.43300956d-04, 2.32084917d-04, 2.23676582d-04, 2.16475971d-04, 669 + 9.83996587d-05, 8.86389763d-05, 7.29954688d-05, 6.85803806d-05, 670 + 6.78189840d-05, 6.58140892d-05, 6.50938107d-05, 6.38443503d-05, 671 + 6.26090037d-05, 6.04309983d-05/ 672 DATA (rc2(I),I=1,94)/ 673 + 0.00000000d+00, 0.00000000d+00, 4.68434853d-02, 2.52134142d-01, 674 + 2.62330732d-01, 5.80715241d-01, 1.09645788d+00, 1.86382178d+00, 675 + 2.97428667d+00, 4.62608071d+00, 2.21437182d+01, 3.56786858d+01, 676 + 4.36941392d+01, 6.48665951d+01, 9.00636335d+01, 1.20264312d+02, 677 + 1.56401427d+02, 1.99211797d+02, 2.67125067d+02, 3.21464903d+02, 678 + 3.85771242d+02, 4.57510139d+02, 5.37648439d+02, 6.20149697d+02, 679 + 7.18923103d+02, 8.22091519d+02, 9.34290727d+02, 1.05514023d+03, 680 + 1.14249171d+03, 1.32583496d+03, 1.28645677d+03, 1.60012205d+03, 681 + 1.84950772d+03, 2.09753208d+03, 2.35990684d+03, 2.64509680d+03, 682 + 2.97636124d+03, 4.49374019d+04, 4.79280798d+04, 5.31722423d+04, 683 + 5.00658603d+04, 5.59792506d+04, 6.72129157d+04, 6.73325347d+04, 684 + 7.36240653d+04, 7.95039519d+04, 8.76089515d+04, 1.02661441d+05, 685 + 1.02056673d+05, 1.20930700d+05, 1.23237713d+05, 1.33810482d+05, 686 + 1.44995277d+05, 1.56862063d+05, 1.77417216d+05, 1.90991139d+05, 687 + 2.05075073d+05, 2.20185428d+05, 2.36331156d+05, 2.53473962d+05, 688 + 2.71763871d+05, 2.91260169d+05, 3.12044059d+05, 3.34190817d+05, 689 + 3.57837847d+05, 3.82982425d+05, 4.09826224d+05, 4.38465778d+05, 690 + 4.68930154d+05, 5.01396788d+05, 4.85613275d+05, 5.51030851d+05, 691 + 5.93454095d+05, 6.54789149d+05, 5.99621174d+05, 7.21896735d+05, 692 + 7.70758213d+05, 8.21242604d+05, 8.79338995d+05, 9.44887717d+05, 693 + 9.99564233d+05, 1.08189801d+06, 1.16519392d+06, 1.25282932d+06, 694 + 1.57499234d+06, 1.79801905d+06, 2.16513590d+06, 2.44913677d+06, 695 + 2.69497434d+06, 2.99908623d+06, 3.30131969d+06, 3.65338787d+06, 696 + 4.04499005d+06, 4.52885470d+06/ 697 DATA (zeta2(I),I=1,94)/ 698 + 1.00000000d+00, 1.00000000d+00, 1.01441639d+00, 7.08517236d-01, 699 + 7.30003765d-01, 6.15383618d-01, 5.33576000d-01, 4.72577360d-01, 700 + 4.24378943d-01, 3.83754227d-01, 2.38826415d-01, 1.99477996d-01, 701 + 1.92829772d-01, 1.67225112d-01, 1.48741905d-01, 1.34281775d-01, 702 + 1.22429085d-01, 1.12511023d-01, 9.82839214d-02, 9.20130613d-02, 703 + 8.66687270d-02, 8.18541208d-02, 7.75184939d-02, 7.44554398d-02, 704 + 7.02827730d-02, 6.71643705d-02, 6.43177023d-02, 6.17266464d-02, 705 + 6.13921845d-02, 5.71521336d-02, 6.13478031d-02, 5.43903881d-02, 706 + 5.05712825d-02, 4.76791348d-02, 4.52299778d-02, 4.30400573d-02, 707 + 4.01557534d-02, 1.12290074d-02, 1.13416879d-02, 1.07424855d-02, 708 + 1.19482805d-02, 1.14350394d-02, 1.00142006d-02, 1.06774861d-02, 709 + 1.03201321d-02, 1.00700543d-02, 9.64166067d-03, 8.51583692d-03, 710 + 9.15978091d-03, 7.98513470d-03, 8.34889598d-03, 8.04927525d-03, 711 + 7.77194190d-03, 7.51012133d-03, 6.60822060d-03, 6.22807064d-03, 712 + 6.07149684d-03, 5.91742394d-03, 5.77498801d-03, 5.63905896d-03, 713 + 5.50038331d-03, 5.36796683d-03, 5.23954025d-03, 5.11839482d-03, 714 + 4.99401318d-03, 4.87751356d-03, 4.76426996d-03, 4.65444537d-03, 715 + 4.54930660d-03, 4.44815005d-03, 5.12909519d-03, 4.61664077d-03, 716 + 4.44800980d-03, 4.05938174d-03, 1.52242450d-03, 4.15624532d-03, 717 + 4.06185185d-03, 3.98192150d-03, 3.87525947d-03, 3.75252197d-03, 718 + 3.71912608d-03, 3.56493520d-03, 3.44266883d-03, 3.33335862d-03, 719 + 1.04842450d-03, 8.93405393d-04, 6.57894932d-04, 6.00392630d-04, 720 + 5.94698138d-04, 5.71622306d-04, 5.66927029d-04, 5.54520209d-04, 721 + 5.42218623d-04, 5.16230841d-04/ 722 DATA (rc3(I),I=1,94)/ 723 + 0.00000000d+00, 0.00000000d+00, 0.00000000d+00, 0.00000000d+00, 724 + 0.00000000d+00, 0.00000000d+00, 0.00000000d+00, 0.00000000d+00, 725 + 0.00000000d+00, 0.00000000d+00, 6.20076037d-02, 2.76051005d-01, 726 + 1.61014857d-01, 3.93461956d-01, 7.64466624d-01, 1.30215679d+00, 727 + 2.04180763d+00, 3.00984485d+00, 6.98610402d+00, 1.00834127d+01, 728 + 1.23076665d+01, 1.50188788d+01, 1.82395980d+01, 1.97075623d+01, 729 + 2.64041322d+01, 3.15174107d+01, 3.73934453d+01, 4.40633918d+01, 730 + 4.26128540d+01, 6.00190607d+01, 4.45592617d+01, 7.46129815d+01, 731 + 1.02455412d+02, 1.31470485d+02, 1.62925314d+02, 1.97827848d+02, 732 + 2.65304373d+02, 1.57117900d+03, 1.38597704d+03, 1.71450052d+03, 733 + 1.10431209d+03, 1.23423271d+03, 1.74411089d+03, 1.49332905d+03, 734 + 1.64723180d+03, 1.75852626d+03, 2.01523534d+03, 2.69723853d+03, 735 + 2.31628288d+03, 3.21014982d+03, 3.16820224d+03, 3.57493869d+03, 736 + 4.00748342d+03, 4.48040398d+03, 8.18987220d+03, 1.07287229d+04, 737 + 1.14521494d+04, 1.22608234d+04, 1.30255796d+04, 1.37766798d+04, 738 + 1.46837652d+04, 1.55931040d+04, 1.65429433d+04, 1.74661109d+04, 739 + 1.85712156d+04, 1.96379950d+04, 2.07455462d+04, 2.18891141d+04, 740 + 2.30389584d+04, 2.41952252d+04, 1.13020962d+04, 1.75777161d+04, 741 + 1.98740166d+04, 2.95731860d+04, 4.09107217d+05, 2.31029894d+04, 742 + 2.43065901d+04, 2.51685750d+04, 2.69880676d+04, 2.96363456d+04, 743 + 2.90462255d+04, 3.34752965d+04, 3.71679846d+04, 4.07226602d+04, 744 + 7.55624974d+05, 8.77788402d+05, 1.10863240d+06, 1.22508590d+06, 745 + 1.29264761d+06, 1.38835313d+06, 1.46321597d+06, 1.55554758d+06, 746 + 1.65439402d+06, 1.79235741d+06/ 747 DATA (zeta3(I),I=1,94)/ 748 + 1.00000000d+00, 1.00000000d+00, 1.00000000d+00, 1.00000000d+00, 749 + 1.00000000d+00, 1.00000000d+00, 1.00000000d+00, 1.00000000d+00, 750 + 1.00000000d+00, 1.00000000d+00, 1.00835672d+00, 7.69492311d-01, 751 + 8.95052085d-01, 7.49562876d-01, 6.50577105d-01, 5.78352607d-01, 752 + 5.22541213d-01, 4.77987125d-01, 3.65799244d-01, 3.09336153d-01, 753 + 2.97043149d-01, 2.87814545d-01, 2.79727989d-01, 2.91308753d-01, 754 + 2.64243561d-01, 2.56546792d-01, 2.48992596d-01, 2.41688705d-01, 755 + 2.55650963d-01, 2.27972392d-01, 2.57800251d-01, 2.17760648d-01, 756 + 1.95357993d-01, 1.79039271d-01, 1.65938044d-01, 1.54814449d-01, 757 + 1.34999942d-01, 6.51041574d-02, 7.21674364d-02, 6.50557950d-02, 758 + 8.90037912d-02, 8.60216892d-02, 7.34617487d-02, 8.14402358d-02, 759 + 7.89401490d-02, 7.80479431d-02, 7.36580206d-02, 6.49836533d-02, 760 + 7.11605517d-02, 6.14493259d-02, 6.10610277d-02, 5.77205998d-02, 761 + 5.47391222d-02, 5.19382050d-02, 3.30078991d-02, 2.52817223d-02, 762 + 2.45802867d-02, 2.40980268d-02, 2.39041263d-02, 2.36226590d-02, 763 + 2.32201481d-02, 2.28792384d-02, 2.25607988d-02, 2.23758534d-02, 764 + 2.19870231d-02, 2.17465623d-02, 2.15260243d-02, 2.13292973d-02, 765 + 2.11855810d-02, 2.10847684d-02, 4.09382448d-02, 3.03603604d-02, 766 + 2.82918344d-02, 2.03824377d-02, 5.62977448d-03, 2.72440875d-02, 767 + 2.68637757d-02, 2.67959682d-02, 2.60212616d-02, 2.48170625d-02, 768 + 2.57591136d-02, 2.36991792d-02, 2.23967941d-02, 2.13654035d-02, 769 + 4.55423037d-03, 4.25297993d-03, 3.72189295d-03, 3.53030502d-03, 770 + 3.47317821d-03, 3.36862712d-03, 3.31841195d-03, 3.24626475d-03, 771 + 3.17505028d-03, 3.05979289d-03/ 772 DATA (rc4(I),I=1,94)/ 773 + 0.00000000d+00, 0.00000000d+00, 0.00000000d+00, 0.00000000d+00, 774 + 0.00000000d+00, 0.00000000d+00, 0.00000000d+00, 0.00000000d+00, 775 + 0.00000000d+00, 0.00000000d+00, 0.00000000d+00, 0.00000000d+00, 776 + 0.00000000d+00, 0.00000000d+00, 0.00000000d+00, 0.00000000d+00, 777 + 0.00000000d+00, 0.00000000d+00, 6.56020424d-02, 2.77840007d-01, 778 + 3.64029769d-01, 4.23460753d-01, 4.63193641d-01, 2.89810882d-01, 779 + 5.14790980d-01, 5.33655411d-01, 5.50170763d-01, 5.65407164d-01, 780 + 2.84096529d-01, 5.93010550d-01, 1.19942620d-01, 3.60612540d-01, 781 + 7.40570005d-01, 1.28553223d+00, 2.02073112d+00, 2.97615952d+00, 782 + 8.39557269d+00, 6.02606011d+01, 4.87023540d+01, 6.80933520d+01, 783 + 1.62265279d+01, 1.72400236d+01, 3.36917884d+01, 1.83723888d+01, 784 + 1.98403203d+01, 1.85143190d+01, 2.48910274d+01, 4.01353587d+01, 785 + 2.64509614d+01, 4.78655062d+01, 6.01556846d+01, 7.90020599d+01, 786 + 1.01162437d+02, 1.28370348d+02, 7.11169654d+02, 1.64398006d+03, 787 + 1.80851978d+03, 1.91100792d+03, 1.94620437d+03, 2.02596238d+03, 788 + 2.12145362d+03, 2.21085640d+03, 2.29697962d+03, 2.34539613d+03, 789 + 2.45263808d+03, 2.51486491d+03, 2.56861396d+03, 2.61136834d+03, 790 + 2.63029593d+03, 2.62765250d+03, 1.95848638d+02, 6.52135512d+02, 791 + 8.28259407d+02, 2.73088861d+03, 1.26198542d+04, 8.86503684d+02, 792 + 9.16101499d+02, 9.11293836d+02, 9.93328324d+02, 1.14814275d+03, 793 + 1.00157533d+03, 1.30190843d+03, 1.54552793d+03, 1.77816217d+03, 794 + 2.04357126d+04, 2.35854277d+04, 3.42350892d+04, 4.03227664d+04, 795 + 4.16645314d+04, 4.53355366d+04, 4.65504326d+04, 4.90202580d+04, 796 + 5.16318071d+04, 5.74706407d+04/ 797 DATA (zeta4(I),I=1,94)/ 798 + 1.00000000d+00, 1.00000000d+00, 1.00000000d+00, 1.00000000d+00, 799 + 1.00000000d+00, 1.00000000d+00, 1.00000000d+00, 1.00000000d+00, 800 + 1.00000000d+00, 1.00000000d+00, 1.00000000d+00, 1.00000000d+00, 801 + 1.00000000d+00, 1.00000000d+00, 1.00000000d+00, 1.00000000d+00, 802 + 1.00000000d+00, 1.00000000d+00, 1.09950352d+00, 8.70037206d-01, 803 + 8.15822153d-01, 7.79528672d-01, 7.53236219d-01, 7.85765325d-01, 804 + 7.14317589d-01, 6.98660493d-01, 6.84606242d-01, 6.71753041d-01, 805 + 7.16027394d-01, 6.48884245d-01, 9.09499464d-01, 7.63611851d-01, 806 + 6.70634196d-01, 6.03533737d-01, 5.51786971d-01, 5.10051860d-01, 807 + 3.79712304d-01, 2.34159182d-01, 2.43754955d-01, 2.22148073d-01, 808 + 3.31676536d-01, 3.30775898d-01, 2.71976660d-01, 3.39485392d-01, 809 + 3.39067921d-01, 3.57675223d-01, 3.28498848d-01, 2.82667731d-01, 810 + 3.28954437d-01, 2.74280814d-01, 2.50819436d-01, 2.27953575d-01, 811 + 2.08767308d-01, 1.91777539d-01, 1.05361706d-01, 7.68835000d-02, 812 + 7.41638307d-02, 7.29762271d-02, 7.29383574d-02, 7.20265634d-02, 813 + 7.10019314d-02, 7.00798829d-02, 6.92361234d-02, 6.88348560d-02, 814 + 6.78372320d-02, 6.73390280d-02, 6.69450435d-02, 6.66732797d-02, 815 + 6.66395164d-02, 6.68229001d-02, 2.07796306d-01, 1.22939973d-01, 816 + 1.09882989d-01, 6.61756802d-02, 3.89737289d-02, 1.16012046d-01, 817 + 1.17207641d-01, 1.20249726d-01, 1.17166773d-01, 1.10793849d-01, 818 + 1.19240694d-01, 1.06979518d-01, 9.93514855d-02, 9.33089232d-02, 819 + 3.36075409d-02, 3.12070010d-02, 2.45901087d-02, 2.22122582d-02, 820 + 2.21170983d-02, 2.11671631d-02, 2.11828525d-02, 2.07857633d-02, 821 + 2.03904909d-02, 1.91739084d-02/ 822 DATA (rc5(I),I=1,94)/ 823 + 0.00000000d+00, 0.00000000d+00, 0.00000000d+00, 0.00000000d+00, 824 + 0.00000000d+00, 0.00000000d+00, 0.00000000d+00, 0.00000000d+00, 825 + 0.00000000d+00, 0.00000000d+00, 0.00000000d+00, 0.00000000d+00, 826 + 0.00000000d+00, 0.00000000d+00, 0.00000000d+00, 0.00000000d+00, 827 + 0.00000000d+00, 0.00000000d+00, 0.00000000d+00, 0.00000000d+00, 828 + 0.00000000d+00, 0.00000000d+00, 0.00000000d+00, 0.00000000d+00, 829 + 0.00000000d+00, 0.00000000d+00, 0.00000000d+00, 0.00000000d+00, 830 + 0.00000000d+00, 0.00000000d+00, 0.00000000d+00, 0.00000000d+00, 831 + 0.00000000d+00, 0.00000000d+00, 0.00000000d+00, 0.00000000d+00, 832 + 7.41352221d-02, 3.29043651d-01, 4.03366961d-01, 6.16991180d-01, 833 + 3.28706457d-01, 3.89407305d-01, 7.60824501d-01, 3.56433510d-01, 834 + 3.21154553d-01, 2.12638921d-01, 2.68825717d-01, 6.74790071d-01, 835 + 1.16080256d-01, 3.74840652d-01, 7.76154174d-01, 1.36358170d+00, 836 + 2.15461358d+00, 3.17719238d+00, 1.79344549d+01, 4.54639997d+01, 837 + 5.36279623d+01, 5.66894253d+01, 5.62578366d+01, 6.02089260d+01, 838 + 6.47699752d+01, 6.97534168d+01, 7.51227593d+01, 8.08249928d+01, 839 + 8.68693210d+01, 9.31064937d+01, 9.95858875d+01, 1.06259016d+02, 840 + 1.12891772d+02, 1.19452938d+02, 4.86622657d-01, 3.75386966d+01, 841 + 6.13420753d+01, 1.64863797d+02, 2.83371154d+02, 3.37760283d+01, 842 + 2.55288121d+01, 1.66310473d+01, 1.76333146d+01, 2.35863127d+01, 843 + 1.44663326d+01, 2.67219996d+01, 4.10151210d+01, 5.89740153d+01, 844 + 4.56772818d+02, 5.85287609d+02, 1.29543017d+03, 1.78567564d+03, 845 + 1.79552725d+03, 2.05105206d+03, 2.02552823d+03, 2.13209240d+03, 846 + 2.24557320d+03, 2.69991263d+03/ 847 DATA (zeta5(I),I=1,94)/ 848 + 1.00000000d+00, 1.00000000d+00, 1.00000000d+00, 1.00000000d+00, 849 + 1.00000000d+00, 1.00000000d+00, 1.00000000d+00, 1.00000000d+00, 850 + 1.00000000d+00, 1.00000000d+00, 1.00000000d+00, 1.00000000d+00, 851 + 1.00000000d+00, 1.00000000d+00, 1.00000000d+00, 1.00000000d+00, 852 + 1.00000000d+00, 1.00000000d+00, 1.00000000d+00, 1.00000000d+00, 853 + 1.00000000d+00, 1.00000000d+00, 1.00000000d+00, 1.00000000d+00, 854 + 1.00000000d+00, 1.00000000d+00, 1.00000000d+00, 1.00000000d+00, 855 + 1.00000000d+00, 1.00000000d+00, 1.00000000d+00, 1.00000000d+00, 856 + 1.00000000d+00, 1.00000000d+00, 1.00000000d+00, 1.00000000d+00, 857 + 1.11162276d+00, 8.92973194d-01, 8.48665299d-01, 7.81782861d-01, 858 + 8.20458147d-01, 7.79621642d-01, 7.10637387d-01, 7.49615537d-01, 859 + 7.43606586d-01, 7.03198459d-01, 7.34936372d-01, 6.65219859d-01, 860 + 9.33059092d-01, 7.92978117d-01, 7.04257390d-01, 6.39255341d-01, 861 + 5.88851087d-01, 5.47985182d-01, 3.64082401d-01, 2.77824918d-01, 862 + 2.62876282d-01, 2.61606275d-01, 2.68464970d-01, 2.65686159d-01, 863 + 2.62414808d-01, 2.58938476d-01, 2.55366523d-01, 2.49843912d-01, 864 + 2.48193305d-01, 2.44727449d-01, 2.41346773d-01, 2.38072521d-01, 865 + 2.34997746d-01, 2.32105807d-01, 8.03394712d-01, 2.85600858d-01, 866 + 2.53758662d-01, 2.08195570d-01, 1.84383408d-01, 2.84437965d-01, 867 + 3.08075946d-01, 3.63374585d-01, 3.64788118d-01, 3.33552690d-01, 868 + 3.95115818d-01, 3.33602667d-01, 2.93606072d-01, 2.62490124d-01, 869 + 1.63387243d-01, 1.50278726d-01, 1.11071589d-01, 9.65712781d-02, 870 + 9.73247043d-02, 9.18921858d-02, 9.40757894d-02, 9.28564165d-02, 871 + 9.15796780d-02, 8.47922939d-02/ 872 DATA (rc6(I),I=1,94)/ 873 + 0.00000000d+00, 0.00000000d+00, 0.00000000d+00, 0.00000000d+00, 874 + 0.00000000d+00, 0.00000000d+00, 0.00000000d+00, 0.00000000d+00, 875 + 0.00000000d+00, 0.00000000d+00, 0.00000000d+00, 0.00000000d+00, 876 + 0.00000000d+00, 0.00000000d+00, 0.00000000d+00, 0.00000000d+00, 877 + 0.00000000d+00, 0.00000000d+00, 0.00000000d+00, 0.00000000d+00, 878 + 0.00000000d+00, 0.00000000d+00, 0.00000000d+00, 0.00000000d+00, 879 + 0.00000000d+00, 0.00000000d+00, 0.00000000d+00, 0.00000000d+00, 880 + 0.00000000d+00, 0.00000000d+00, 0.00000000d+00, 0.00000000d+00, 881 + 0.00000000d+00, 0.00000000d+00, 0.00000000d+00, 0.00000000d+00, 882 + 0.00000000d+00, 0.00000000d+00, 0.00000000d+00, 0.00000000d+00, 883 + 0.00000000d+00, 0.00000000d+00, 0.00000000d+00, 0.00000000d+00, 884 + 0.00000000d+00, 0.00000000d+00, 0.00000000d+00, 0.00000000d+00, 885 + 0.00000000d+00, 0.00000000d+00, 0.00000000d+00, 0.00000000d+00, 886 + 0.00000000d+00, 0.00000000d+00, 8.41198181d-02, 3.30444084d-01, 887 + 4.91705562d-01, 4.94765259d-01, 3.69753458d-01, 3.69539339d-01, 888 + 3.69853666d-01, 3.70622471d-01, 3.71721718d-01, 4.34972852d-01, 889 + 3.75404514d-01, 3.77695171d-01, 3.80389931d-01, 3.83527706d-01, 890 + 3.86898212d-01, 3.90604654d-01, 2.97047551d-03, 3.99986047d-01, 891 + 7.20996101d-01, 1.06952842d+00, 1.37569576d+00, 1.25089939d+00, 892 + 1.23420443d+00, 6.29685108d-01, 5.66212693d-01, 1.04389432d+00, 893 + 7.78008312d-02, 2.99789533d-01, 6.88443909d-01, 1.27022912d+00, 894 + 2.43994094d+00, 3.56059743d+00, 1.77231239d+01, 4.17467686d+01, 895 + 3.83943193d+01, 5.43331103d+01, 3.96156778d+01, 3.92950589d+01, 896 + 3.94476243d+01, 5.61524915d+01/ 897 DATA (zeta6(I),I=1,94)/ 898 + 1.00000000d+00, 1.00000000d+00, 1.00000000d+00, 1.00000000d+00, 899 + 1.00000000d+00, 1.00000000d+00, 1.00000000d+00, 1.00000000d+00, 900 + 1.00000000d+00, 1.00000000d+00, 1.00000000d+00, 1.00000000d+00, 901 + 1.00000000d+00, 1.00000000d+00, 1.00000000d+00, 1.00000000d+00, 902 + 1.00000000d+00, 1.00000000d+00, 1.00000000d+00, 1.00000000d+00, 903 + 1.00000000d+00, 1.00000000d+00, 1.00000000d+00, 1.00000000d+00, 904 + 1.00000000d+00, 1.00000000d+00, 1.00000000d+00, 1.00000000d+00, 905 + 1.00000000d+00, 1.00000000d+00, 1.00000000d+00, 1.00000000d+00, 906 + 1.00000000d+00, 1.00000000d+00, 1.00000000d+00, 1.00000000d+00, 907 + 1.00000000d+00, 1.00000000d+00, 1.00000000d+00, 1.00000000d+00, 908 + 1.00000000d+00, 1.00000000d+00, 1.00000000d+00, 1.00000000d+00, 909 + 1.00000000d+00, 1.00000000d+00, 1.00000000d+00, 1.00000000d+00, 910 + 1.00000000d+00, 1.00000000d+00, 1.00000000d+00, 1.00000000d+00, 911 + 1.00000000d+00, 1.00000000d+00, 1.14652099d+00, 9.40419863d-01, 912 + 8.71830139d-01, 8.63645308d-01, 8.99781351d-01, 8.92847773d-01, 913 + 8.86263534d-01, 8.79891833d-01, 8.73693204d-01, 8.44849538d-01, 914 + 8.61492212d-01, 8.55512154d-01, 8.49557770d-01, 8.43601596d-01, 915 + 8.37682504d-01, 8.31763174d-01, 1.16006472d+00, 8.07528576d-01, 916 + 7.31082249d-01, 6.83425161d-01, 6.52505241d-01, 6.40350000d-01, 917 + 6.28116671d-01, 6.50985392d-01, 6.43444641d-01, 6.04768070d-01, 918 + 9.74028988d-01, 8.20486767d-01, 7.26711714d-01, 6.60125651d-01, 919 + 6.03396512d-01, 5.63782867d-01, 3.80991010d-01, 2.97211769d-01, 920 + 3.02765896d-01, 2.71958506d-01, 3.05459861d-01, 3.09950554d-01, 921 + 3.13551092d-01, 2.87960597d-01/ 922 DATA (rc7(I),I=1,94)/ 923 + 0.00000000d+00, 0.00000000d+00, 0.00000000d+00, 0.00000000d+00, 924 + 0.00000000d+00, 0.00000000d+00, 0.00000000d+00, 0.00000000d+00, 925 + 0.00000000d+00, 0.00000000d+00, 0.00000000d+00, 0.00000000d+00, 926 + 0.00000000d+00, 0.00000000d+00, 0.00000000d+00, 0.00000000d+00, 927 + 0.00000000d+00, 0.00000000d+00, 0.00000000d+00, 0.00000000d+00, 928 + 0.00000000d+00, 0.00000000d+00, 0.00000000d+00, 0.00000000d+00, 929 + 0.00000000d+00, 0.00000000d+00, 0.00000000d+00, 0.00000000d+00, 930 + 0.00000000d+00, 0.00000000d+00, 0.00000000d+00, 0.00000000d+00, 931 + 0.00000000d+00, 0.00000000d+00, 0.00000000d+00, 0.00000000d+00, 932 + 0.00000000d+00, 0.00000000d+00, 0.00000000d+00, 0.00000000d+00, 933 + 0.00000000d+00, 0.00000000d+00, 0.00000000d+00, 0.00000000d+00, 934 + 0.00000000d+00, 0.00000000d+00, 0.00000000d+00, 0.00000000d+00, 935 + 0.00000000d+00, 0.00000000d+00, 0.00000000d+00, 0.00000000d+00, 936 + 0.00000000d+00, 0.00000000d+00, 0.00000000d+00, 0.00000000d+00, 937 + 0.00000000d+00, 0.00000000d+00, 0.00000000d+00, 0.00000000d+00, 938 + 0.00000000d+00, 0.00000000d+00, 0.00000000d+00, 0.00000000d+00, 939 + 0.00000000d+00, 0.00000000d+00, 0.00000000d+00, 0.00000000d+00, 940 + 0.00000000d+00, 0.00000000d+00, 0.00000000d+00, 0.00000000d+00, 941 + 0.00000000d+00, 0.00000000d+00, 0.00000000d+00, 0.00000000d+00, 942 + 0.00000000d+00, 0.00000000d+00, 0.00000000d+00, 0.00000000d+00, 943 + 0.00000000d+00, 0.00000000d+00, 0.00000000d+00, 0.00000000d+00, 944 + 0.00000000d+00, 0.00000000d+00, 1.01580194d-01, 3.80390160d-01, 945 + 4.13217274d-01, 7.39534177d-01, 4.50487346d-01, 4.21044941d-01, 946 + 3.82783864d-01, 5.04264000d-01/ 947 DATA (zeta7(I),I=1,94)/ 948 + 1.00000000d+00, 1.00000000d+00, 1.00000000d+00, 1.00000000d+00, 949 + 1.00000000d+00, 1.00000000d+00, 1.00000000d+00, 1.00000000d+00, 950 + 1.00000000d+00, 1.00000000d+00, 1.00000000d+00, 1.00000000d+00, 951 + 1.00000000d+00, 1.00000000d+00, 1.00000000d+00, 1.00000000d+00, 952 + 1.00000000d+00, 1.00000000d+00, 1.00000000d+00, 1.00000000d+00, 953 + 1.00000000d+00, 1.00000000d+00, 1.00000000d+00, 1.00000000d+00, 954 + 1.00000000d+00, 1.00000000d+00, 1.00000000d+00, 1.00000000d+00, 955 + 1.00000000d+00, 1.00000000d+00, 1.00000000d+00, 1.00000000d+00, 956 + 1.00000000d+00, 1.00000000d+00, 1.00000000d+00, 1.00000000d+00, 957 + 1.00000000d+00, 1.00000000d+00, 1.00000000d+00, 1.00000000d+00, 958 + 1.00000000d+00, 1.00000000d+00, 1.00000000d+00, 1.00000000d+00, 959 + 1.00000000d+00, 1.00000000d+00, 1.00000000d+00, 1.00000000d+00, 960 + 1.00000000d+00, 1.00000000d+00, 1.00000000d+00, 1.00000000d+00, 961 + 1.00000000d+00, 1.00000000d+00, 1.00000000d+00, 1.00000000d+00, 962 + 1.00000000d+00, 1.00000000d+00, 1.00000000d+00, 1.00000000d+00, 963 + 1.00000000d+00, 1.00000000d+00, 1.00000000d+00, 1.00000000d+00, 964 + 1.00000000d+00, 1.00000000d+00, 1.00000000d+00, 1.00000000d+00, 965 + 1.00000000d+00, 1.00000000d+00, 1.00000000d+00, 1.00000000d+00, 966 + 1.00000000d+00, 1.00000000d+00, 1.00000000d+00, 1.00000000d+00, 967 + 1.00000000d+00, 1.00000000d+00, 1.00000000d+00, 1.00000000d+00, 968 + 1.00000000d+00, 1.00000000d+00, 1.00000000d+00, 1.00000000d+00, 969 + 1.00000000d+00, 1.00000000d+00, 1.10753948d+00, 9.23783778d-01, 970 + 9.01059856d-01, 8.16067064d-01, 8.69278115d-01, 8.68260354d-01, 971 + 8.71564352d-01, 8.34358766d-01/ 972 973 xdm_rhop = rc1(iz) * exp(-r/zeta1(iz)) 974 if (zeta2(iz)<1d-12) return 975 xdm_rhop = xdm_rhop + rc2(iz) * exp(-r/zeta2(iz)) 976 if (zeta3(iz)<1d-12) return 977 xdm_rhop = xdm_rhop + rc3(iz) * exp(-r/zeta3(iz)) 978 if (zeta4(iz)<1d-12) return 979 xdm_rhop = xdm_rhop + rc4(iz) * exp(-r/zeta4(iz)) 980 if (zeta5(iz)<1d-12) return 981 xdm_rhop = xdm_rhop + rc5(iz) * exp(-r/zeta5(iz)) 982 if (zeta6(iz)<1d-12) return 983 xdm_rhop = xdm_rhop + rc6(iz) * exp(-r/zeta6(iz)) 984 if (zeta7(iz)<1d-12) return 985 xdm_rhop = xdm_rhop + rc7(iz) * exp(-r/zeta7(iz)) 986 987 end 988 integer function xc_xdm_lxdm() 989 implicit none 990 991#include "cdft.fh" 992 xc_xdm_lxdm=lxdm 993 return 994 end 995