1c +++++++++++++++++++++++++++++++ 2c + calculate all Hyperfine AOs + 3c +++++++++++++++++++++++++++++++ 4c 1. (zpsox,zpsoy,zpsoz): 5c H^{ZPSO}_{mu nu,Aj}= \int dr K/r_A^3 6c \vec{r}_A x [chi_{mu}^* \nabla chi_{nu} - 7c chi_{nu}^* \nabla chi_{mu}^* ]_j 8c (Eq. 56 in J. Autschbach's write-up of 9c ZORA-NMR spin-spin coupling constants 10c Sept. 17, 2007's write-up) 11 12c 2. (fcsdxx,fcsdxy,fcsdxz, 13c fcsdyx,fcsdyy,fcsdyz, 14c fcsdzx,fcsdzy,fcsdzz): 15c H^{FC+SD}_{uv}=\int dr K U_{N,v} \nabla_{u} (chi_{mu}^* chi_{nu}) 16c where U_{N,v} is, 17c U_{N,v} c^{-2} r_{N,v}/r_N^3 (Eq. 10 in JA's draft of 18c 'Calculation of hyperfine tensor using zeroth-order regular approximation 19c and density functional theory: Expectation value versus linear response 20c approaches') 21 22 subroutine calc_zora_HFine_slow( 23 & ao_bas_han, ! in: AO basis handle 24 & geom, ! in: geometry handle 25 & ipol, ! in: nr. of polarizations 26 & g_dens, ! in: superposit. atomic densities 27 & chi_ao, ! in: basis functions 28 & delchi_ao, ! in: deriv. of basis functions 29 & qxyz, ! in: grid points 30 & qwght, ! in: weighting coeffs. 31 & nbf, ! in: nr. basis functions 32 & npts, ! in: nr. grid points 33 & natoms, ! in: nr. atoms 34 & ofinite, ! in: = .true. if Gaussian Nucl. Model of charges requested 35 & zetanuc_arr, ! in: zetanuc(i) i=1,natoms for Gaussian Nuclear Model 36 & atmass, ! in: atomic mass 37 & xyz_NMRcoords,! in : nuclear coordinates 38 & use_modelpotential, 39 & gexpo, 40 & gcoef, 41 & zpsox, ! out 42 & zpsoy, ! out 43 & zpsoz, ! out 44 & fcsdxx, ! out 45 & fcsdxy, ! out 46 & fcsdxz, ! out 47c & fcsdyx, ! out 48 & fcsdyy, ! out 49 & fcsdyz, ! out 50c & fcsdzx, ! out 51c & fcsdzy, ! out 52 & fcsdzz) ! out 53 implicit none 54#include "errquit.fh" 55#include "mafdecls.fh" 56#include "stdio.fh" 57#include "global.fh" 58#include "bas.fh" 59#include "zora.fh" 60 integer nbf,npts,ao_bas_han,natoms,geom 61 integer g_dens(2),ipol 62 double precision qwght(npts) 63 double precision qxyz(3,npts) 64 double precision chi_ao(npts,nbf) 65 double precision delchi_ao(npts,3,nbf) 66 double precision zpsox(nbf,nbf), 67 & zpsoy(nbf,nbf), 68 & zpsoz(nbf,nbf) 69 double precision fcsdxx(nbf,nbf), 70 & fcsdxy(nbf,nbf), 71 & fcsdxz(nbf,nbf), 72 & fcsdyy(nbf,nbf), 73 & fcsdyz(nbf,nbf), 74 & fcsdzz(nbf,nbf) 75 double precision ac_fcsd(3,3) 76 integer i,j,k,n 77 double precision amat_coul(npts,ipol) 78 double precision amat_nucl(npts),amat_NMRnucl(3,npts), 79 & amat_Pnucl(npts) 80 integer ipt,closegridpts(npts) 81 double precision clight_au2,tol 82 double precision amat_tot,Kzora 83 double precision fac1_arr(npts),fac2_arr(3,npts) 84 double precision ac_zpso(3) 85 double precision xyz_NMRcoords(3),atmass 86 double precision chi_cntr(3,nbf),threehalf 87 data threehalf /1.5d0/ 88 logical ofinite 89c ------- for Gaussian Nuclear Model --- START 90 double precision zetanuc_arr(natoms) 91c ------- for Gaussian Nuclear Model --- START 92c 93 logical use_modelpotential 94 double precision gexpo(natoms,50) 95 double precision gcoef(natoms,50) 96c 97c == preliminaries == 98 clight_au2 = clight_au*clight_au 99 do ipt = 1,npts 100 do i=1,ipol 101 amat_coul(ipt,i) = 0.d0 102 end do 103 amat_nucl(ipt) = 0.d0 104 amat_Pnucl(ipt) = 0.0d0 105 closegridpts(ipt) = 0 106 do i=1,3 107 amat_NMRnucl(i,ipt) = 0.d0 108 enddo 109 end do 110c 111c == calculate the total hartree potential on the grid == 112 call gridHartreePotential(use_modelpotential, 113 & ao_bas_han, geom, natoms, ipol, g_dens, npts, qxyz, qwght, 114 & closegridpts, gexpo, gcoef, amat_coul) 115c 116c == calculate the total nuclear potential on the grid == 117 if (ofinite) then 118c ------ Choosing Nuclear Model: erf(zetanuc^0.5 r_L) 119 call gridNuclearPotentialFinite(geom,natoms,npts,qxyz,qwght, 120 & zetanuc_arr, 121 & closegridpts, 122 & amat_nucl) 123c ------ Choosing Nuclear Model: P(1/2,zetanuc r_L^2) 124c call gridNuclearPotentialFinite2(geom,natoms,npts,qxyz,qwght, 125c & closegridpts,amat_nucl) 126 else ! default : point charge model for nuclei 127 call gridNuclearPotentialPoint(geom,natoms,npts,qxyz,qwght, 128 & closegridpts,amat_nucl) 129 endif 130 do k = 1,npts 131 if (k.eq.closegridpts(k)) qwght(k) = 0.d0 132 end do 133 call gridNMRPotential(amat_NMRnucl, ! out: NMR potential 134 & xyz_NMRcoords, 135 & npts,qxyz,closegridpts) 136 if (ofinite) then ! ====> GAUSSIAN charge nuclear model 137 call get_Pnucl(amat_Pnucl, ! out: P(3/2,r_N^2) 138 & atmass, ! in : atomic mass 139 & xyz_NMRcoords, ! in : EFG-nuclear coord. 140 & threehalf, 141 & npts,qxyz) 142c === define fac_arr 143 do k = 1,npts 144c == assemble hartree and nuclear contributions == 145 amat_tot = amat_nucl(k) + amat_coul(k,1) 146 Kzora=1.0d0/(1.0d0-0.5d0*amat_tot/clight_au2) 147c +++++++++++++++++++++++++++++++++++++++++++++++++++++++ 148 if (do_NonRel) then ! remove it after TEST 149 Kzora=1.0d0 ! remove it after TEST 150 endif ! remove it after TEST 151c +++++++++++++++++++++++++++++++++++++++++++++++++++++++ 152 fac1_arr(k)=Kzora*qwght(k)*amat_Pnucl(k) 153 do n=1,3 154 fac2_arr(n,k)=fac1_arr(k)*amat_NMRnucl(n,k) ! for ZPSO 155 enddo ! end-loop-n 156 enddo ! end-loop-k 157 else ! ====> POINT charge nuclear model (default)---START 158c === define fac_arr 159 do k = 1,npts 160c == assemble hartree and nuclear contributions == 161 amat_tot = amat_nucl(k) + amat_coul(k,1) 162 Kzora=1.0d0/(1.0d0-0.5d0*amat_tot/clight_au2) 163c +++++++++++++++++++++++++++++++++++++++++++++++++++++++ 164 if (do_NonRel) then ! remove it after TEST 165 Kzora=1.0d0 ! remove it after TEST 166 endif ! remove it after TEST 167c +++++++++++++++++++++++++++++++++++++++++++++++++++++++ 168 fac1_arr(k)=Kzora*qwght(k) 169 do n=1,3 170 fac2_arr(n,k)=fac1_arr(k)*amat_NMRnucl(n,k) ! for ZPSO 171 enddo ! end-loop-n 172 enddo ! end-loop-k 173 endif ! ====> POINT charge nuclear model (default)---END 174c == assemble zora correction == 175c ---- full matrix calc -------- START 176 do i = 1, nbf 177 do j = 1, nbf 178 call get_ints_zora_hfine_slow( 179 & nbf,npts,chi_ao,delchi_ao,i,j, 180 & fac2_arr, 181 & ac_zpso, ! out 182 & ac_fcsd) ! out 183 zpsox(i,j) = zpsox(i,j) + ac_zpso(1) 184 zpsoy(i,j) = zpsoy(i,j) + ac_zpso(2) 185 zpsoz(i,j) = zpsoz(i,j) + ac_zpso(3) 186 fcsdxx(i,j) = fcsdxx(i,j) + ac_fcsd(1,1) 187 fcsdxy(i,j) = fcsdxy(i,j) + ac_fcsd(1,2) 188 fcsdxz(i,j) = fcsdxz(i,j) + ac_fcsd(1,3) 189c fcsdyx(i,j) = fcsdyx(i,j) + ac_fcsd(2,1) 190 fcsdyy(i,j) = fcsdyy(i,j) + ac_fcsd(2,2) 191 fcsdyz(i,j) = fcsdyz(i,j) + ac_fcsd(2,3) 192c fcsdzx(i,j) = fcsdzx(i,j) + ac_fcsd(3,1) 193c fcsdzy(i,j) = fcsdzy(i,j) + ac_fcsd(3,2) 194 fcsdzz(i,j) = fcsdzz(i,j) + ac_fcsd(3,3) 195 enddo ! end-loop-j 196 enddo ! end-loop-i 197c ---- full matrix calc -------- END 198 return 199 end 200 201 subroutine get_ints_zora_hfine_slow( 202 & nbf, ! in: nr. basis functions 203 & npts, ! in: grid points 204 & chi_ao, ! in: basis functions 205 & delchi_ao, ! in: deriv. of basis functions 206 & i,j, ! in: (i,j) indices for delchi_ao 207 & fac2_arr, ! in 208 & ac_zpso, ! out : ZPSO term 209 & ac_fcsd) ! out : FC+SD term (n,m) component 210 implicit none 211#include "errquit.fh" 212#include "stdio.fh" 213#include "global.fh" 214 integer nbf,npts,i,j,k,m,n,a,b 215 double precision chi_ao(npts,nbf) 216 double precision delchi_ao(npts,3,nbf) 217 double precision fac2_arr(3,npts) 218 double precision ac_zpso(3), 219 & ac_fcsd(3,3) 220 double precision prod(3),prod1(3) 221 integer ind_nab(2,3) 222 data ind_nab / 2, 3, ! nab=123 223 & 3, 1, ! nab=231 224 & 1, 2 / ! nab=312 225 do n=1,3 ! reset 226 ac_zpso(n) = 0.0d0 227 do m=1,3 228 ac_fcsd(n,m) = 0.0d0 229 enddo 230 enddo 231 do k = 1, npts 232 do n=1,3 233 prod1(n) = chi_ao(k,i)*delchi_ao(k,n,j) 234 & +chi_ao(k,j)*delchi_ao(k,n,i) 235 enddo ! end-loop-n 236 do m=n,3 237 do n=1,3 238 ac_fcsd(n,m) = ac_fcsd(n,m) + 239 & fac2_arr(n,k)*prod1(m) 240 enddo ! end-loop-m 241 enddo ! end-loop-n 242 enddo ! end-loo-k 243 do k = 1, npts 244 do n=1,3 245 prod(n) = chi_ao(k,i)*delchi_ao(k,n,j) 246 & -chi_ao(k,j)*delchi_ao(k,n,i) 247 enddo ! end-loop-n 248 do n=1,3 249 a=ind_nab(1,n) 250 b=ind_nab(2,n) 251 ac_zpso(n) = ac_zpso(n) + 252 & fac2_arr(a,k)*prod(b)- 253 & fac2_arr(b,k)*prod(a) 254 enddo ! end-loop-n 255 enddo ! end-loo-k 256 return 257 end 258c +++++++++++++++++++++++++++++++++++++++ 259c +++++++++++++++++++++++++++++++++++++++ 260 subroutine calc_zora_HFine_fast( 261 & ao_bas_han, ! in: AO basis handle 262 & geom, ! in: geometry handle 263 & ipol, ! in: nr. of polarizations 264 & g_dens, ! in: superposit. atomic densities 265 & chi_ao, ! in: basis functions 266 & delchi_ao, ! in: deriv. of basis functions 267 & qxyz, ! in: grid points 268 & qwght, ! in: weighting coeffs. 269 & nbf, mbf,ibf, ! in: nr. basis functions 270 & npts, ! in: nr. grid points 271 & natoms, ! in: nr. atoms 272 & ofinite, ! in: = .true. if Gaussian Nucl. Model of charges requested 273 & zetanuc_arr, ! in: sqrt(zetanuc(i)) i=1,natoms for Gaussian Nuclear Model 274 & zetanuc_slc, ! in: zetanuc(i) 275 & Knucl, 276 & xyz_NMRcoords,! in: nuclear coordinates 277 & use_modelpotential, 278 & gexpo, 279 & gcoef, 280 & zpsoxy, ! out 281 & zpsoyz, ! out 282 & zpsozx, ! out 283 & fcsdxx, ! out 284 & fcsdxy, ! out 285 & fcsdxz, ! out 286 & fcsdyy, ! out 287 & fcsdyz, ! out 288 & fcsdzz) ! out 289 implicit none 290#include "errquit.fh" 291#include "mafdecls.fh" 292#include "stdio.fh" 293#include "global.fh" 294#include "bas.fh" 295#include "zora.fh" 296 integer nbf,npts,ao_bas_han,natoms,geom 297 integer mbf,ibf(*) 298 integer g_dens(2),ipol 299 double precision qwght(npts) 300 double precision qxyz(3,npts) 301 double precision chi_ao(npts,nbf) 302 double precision delchi_ao(npts,3,nbf) 303 double precision zpsoxy(nbf,nbf), 304 & zpsoxz(nbf,nbf), 305 & zpsoyx(nbf,nbf), 306 & zpsoyz(nbf,nbf), 307 & zpsozx(nbf,nbf), 308 & zpsozy(nbf,nbf) 309 double precision fcsdxx(nbf,nbf), 310 & fcsdxy(nbf,nbf), 311 & fcsdxz(nbf,nbf), 312c & fcsdyx(nbf,nbf), 313 & fcsdyy(nbf,nbf), 314 & fcsdyz(nbf,nbf), 315c & fcsdzx(nbf,nbf), 316c & fcsdzy(nbf,nbf), 317 & fcsdzz(nbf,nbf) 318 double precision ac_fcsd(3,3) 319 integer i,j,k,n 320 double precision amat_coul(npts,ipol) 321 double precision amat_nucl(npts),amat_NMRnucl(3,npts), 322 & amat_Pnucl(npts) 323 integer ipt,closegridpts(npts) 324 double precision clight_au2,tol 325 double precision amat_tot,Kzora 326 double precision fac1_arr(npts),fac2_arr(3,npts) 327 double precision ac_zpso(3,3) 328 double precision xyz_NMRcoords(3) 329 double precision chi_cntr(3,nbf),qxyz1(3) 330 double precision threehalf 331 data threehalf /1.5/ 332 logical ofinite,Knucl 333c 334 double precision zetanuc_arr(natoms),Pnucl 335 double precision zetanuc_slc 336 integer count_pt ! ONLY for checking get_Pnucl 337c 338 integer i0,j0 339 logical use_modelpotential 340 double precision gexpo(natoms,50) 341 double precision gcoef(natoms,50) 342c 343c == preliminaries == 344 clight_au2 = clight_au*clight_au 345 do ipt = 1,npts 346 do i=1,ipol 347 amat_coul(ipt,i) = 0.d0 348 end do 349 amat_nucl(ipt) = 0.d0 350 closegridpts(ipt) = 0 351 do i=1,3 352 amat_NMRnucl(i,ipt) = 0.d0 353 enddo 354 end do 355c 356c == calculate the total hartree potential on the grid == 357 call gridHartreePotential(use_modelpotential, 358 & ao_bas_han, geom, natoms, ipol, g_dens, npts, qxyz, qwght, 359 & closegridpts, gexpo, gcoef, amat_coul) 360c 361c == calculate the total nuclear potential on the grid == 362 if (ofinite) then 363c 364c ------ Choosing Nuclear Model: erf(zetanuc^0.5 r_L) 365 call gridNuclearPotentialFinite(geom,natoms,npts,qxyz,qwght, 366 & zetanuc_arr, 367 & closegridpts,amat_nucl) 368c ------ Choosing Nuclear Model: P(1/2,zetanuc r_L^2) 369c call gridNuclearPotentialFinite2(geom,natoms,npts,qxyz,qwght, 370c & closegridpts,amat_nucl) 371 else ! default : point charge model for nuclei 372c 373 call gridNuclearPotentialPoint(geom,natoms,npts,qxyz,qwght, 374 & closegridpts,amat_nucl) 375 endif 376c 377 do k = 1,npts 378 if (k.eq.closegridpts(k)) qwght(k) = 0.d0 379 end do 380c 381 call gridNMRPotential(amat_NMRnucl, ! out: NMR potential 382 & xyz_NMRcoords, 383 & npts,qxyz,closegridpts) 384 if (ofinite) then ! ====> GAUSSIAN charge nuclear model 385 if (Knucl) then !-- V=Vnucl (amat_tot) 386 count_pt=1 ! ONLY for checking get_Pnucl 387 do k = 1,npts 388 qxyz1(1)=qxyz(1,k) 389 qxyz1(2)=qxyz(2,k) 390 qxyz1(3)=qxyz(3,k) 391 call get_Pnucl1(Pnucl, ! out: P(3/2,r_N^2) 392 & zetanuc_slc, ! in : atomic mass 393 & xyz_NMRcoords, ! in : EFG-nuclear coord. 394 & threehalf, 395 & npts,qxyz1,count_pt) 396c == assemble hartree and nuclear contributions == 397 amat_tot = amat_nucl(k) ! V = Vnucl (ONLY) 398 Kzora=1.0d0/(1.0d0-0.5d0*amat_tot/clight_au2) 399c +++++++++++++++++++++++++++++++++++++++++++++++++++++++ 400 if (do_NonRel) then ! remove it after TEST 401 Kzora=1.0d0 ! remove it after TEST 402 endif ! remove it after TEST 403c +++++++++++++++++++++++++++++++++++++++++++++++++++++++ 404 fac1_arr(k)=Kzora*qwght(k)*Pnucl 405 do n=1,3 406 fac2_arr(n,k)=fac1_arr(k)*amat_NMRnucl(n,k) ! for ZPSO 407 enddo ! end-loop-n 408 enddo ! end-loop-k 409 else ! ------------ V=Vnucl+Vee (amat_tot) (default) 410c === define fac_arr 411 count_pt=1 ! ONLY for checking get_Pnucl 412 do k = 1,npts 413 qxyz1(1)=qxyz(1,k) 414 qxyz1(2)=qxyz(2,k) 415 qxyz1(3)=qxyz(3,k) 416 call get_Pnucl1(Pnucl, ! out: P(3/2,r_N^2) 417 & zetanuc_slc, ! in : atomic mass 418 & xyz_NMRcoords, ! in : EFG-nuclear coord. 419 & threehalf, 420 & npts,qxyz1,count_pt) 421c == assemble hartree and nuclear contributions == 422 amat_tot = amat_nucl(k) + amat_coul(k,1) ! V=Vnucl+Vee 423 Kzora=1.0d0/(1.0d0-0.5d0*amat_tot/clight_au2) 424c +++++++++++++++++++++++++++++++++++++++++++++++++++++++ 425 if (do_NonRel) then ! remove it after TEST 426 Kzora=1.0d0 ! remove it after TEST 427 endif ! remove it after TEST 428c +++++++++++++++++++++++++++++++++++++++++++++++++++++++ 429 fac1_arr(k)=Kzora*qwght(k)*Pnucl 430 do n=1,3 431 fac2_arr(n,k)=fac1_arr(k)*amat_NMRnucl(n,k) ! for ZPSO 432 enddo ! end-loop-n 433 enddo ! end-loop-k 434 endif ! end-if-Knucl 435 else ! ====> POINT charge nuclear model (default)---START 436 if (Knucl) then !-- V=Vnucl (amat_tot) 437c === define fac_arr 438 do k = 1,npts 439c == assemble hartree and nuclear contributions == 440 amat_tot = amat_nucl(k) ! V=Vnucl (ONLY) 441 Kzora=1.0d0/(1.0d0-0.5d0*amat_tot/clight_au2) 442c +++++++++++++++++++++++++++++++++++++++++++++++++++++++ 443 if (do_NonRel) then ! remove it after TEST 444 Kzora=1.0d0 ! remove it after TEST 445 endif ! remove it after TEST 446c +++++++++++++++++++++++++++++++++++++++++++++++++++++++ 447 fac1_arr(k)=Kzora*qwght(k) 448 do n=1,3 449 fac2_arr(n,k)=fac1_arr(k)*amat_NMRnucl(n,k) ! for ZPSO 450 enddo ! end-loop-n 451 enddo ! end-loop-k 452 else ! ------------ V=Vnucl+Vee (amat_tot) (default) 453c === define fac_arr 454 do k = 1,npts 455c == assemble hartree and nuclear contributions == 456 amat_tot = amat_nucl(k) + amat_coul(k,1) ! V=Vnucl+Vee 457 Kzora=1.0d0/(1.0d0-0.5d0*amat_tot/clight_au2) 458c +++++++++++++++++++++++++++++++++++++++++++++++++++++++ 459 if (do_NonRel) then ! remove it after TEST 460 Kzora=1.0d0 ! remove it after TEST 461 endif ! remove it after TEST 462c +++++++++++++++++++++++++++++++++++++++++++++++++++++++ 463 fac1_arr(k)=Kzora*qwght(k) 464 do n=1,3 465 fac2_arr(n,k)=fac1_arr(k)*amat_NMRnucl(n,k) ! for ZPSO 466 enddo ! end-loop-n 467 enddo ! end-loop-k 468 endif ! end-if-Knucl 469 endif ! ====> POINT charge nuclear model (default)---END 470c == assemble zora correction == 471c ---- full matrix calc -------- START 472 do i0 = 1, mbf 473 i=ibf(i0) 474 do j0 = i0, mbf 475 j=ibf(j0) 476 call get_ints_zora_hfine_fast( 477 & nbf,npts,chi_ao,delchi_ao,i0,j0, 478 & fac2_arr, 479 & ac_zpso, ! out 480 & ac_fcsd) ! out 481 zpsoxy(i,j) = zpsoxy(i,j) + ac_zpso(1,2) 482 zpsoxz(i,j) = zpsoxz(i,j) + ac_zpso(1,3) 483 zpsoyx(i,j) = zpsoyx(i,j) + ac_zpso(2,1) 484 zpsoyz(i,j) = zpsoyz(i,j) + ac_zpso(2,3) 485 zpsozx(i,j) = zpsozx(i,j) + ac_zpso(3,1) 486 zpsozy(i,j) = zpsozy(i,j) + ac_zpso(3,2) 487 fcsdxx(i,j) = fcsdxx(i,j) + ac_fcsd(1,1) 488 fcsdxy(i,j) = fcsdxy(i,j) + ac_fcsd(1,2) 489 fcsdxz(i,j) = fcsdxz(i,j) + ac_fcsd(1,3) 490c fcsdyx(i,j) = fcsdyx(i,j) + ac_fcsd(2,1) 491 fcsdyy(i,j) = fcsdyy(i,j) + ac_fcsd(2,2) 492 fcsdyz(i,j) = fcsdyz(i,j) + ac_fcsd(2,3) 493c fcsdzx(i,j) = fcsdzx(i,j) + ac_fcsd(3,1) 494c fcsdzy(i,j) = fcsdzy(i,j) + ac_fcsd(3,2) 495 fcsdzz(i,j) = fcsdzz(i,j) + ac_fcsd(3,3) 496 enddo ! end-loop-j 497 enddo ! end-loop-i 498crecover upper triangle 499 do i0 = 1, mbf 500 i=ibf(i0) 501 do j0 = i0+1, mbf 502 j=ibf(j0) 503 zpsoxy(j,i) = zpsoxy(i,j) 504 zpsoxz(j,i) = zpsoxz(i,j) 505 zpsoyx(j,i) = zpsoyx(i,j) 506 zpsoyz(j,i) = zpsoyz(i,j) 507 zpsozx(j,i) = zpsozx(i,j) 508 zpsozy(j,i) = zpsozy(i,j) 509 fcsdxx(j,i) = fcsdxx(i,j) 510 fcsdxy(j,i) = fcsdxy(i,j) 511 fcsdxz(j,i) = fcsdxz(i,j) 512 fcsdyy(j,i) = fcsdyy(i,j) 513 fcsdyz(j,i) = fcsdyz(i,j) 514 fcsdzz(j,i) = fcsdzz(i,j) 515 enddo ! end-loop-j 516 enddo ! end-loop-i 517c ---- full matrix calc -------- END 518 return 519 end 520 521 subroutine get_ints_zora_hfine_fast( 522 & nbf, ! in: nr. basis functions 523 & npts, ! in: grid points 524 & chi_ao, ! in: basis functions 525 & delchi_ao, ! in: deriv. of basis functions 526 & i,j, ! in: (i,j) indices for delchi_ao 527 & fac2_arr, ! in 528 & ac_zpso, ! out : ZPSO term 529 & ac_fcsd) ! out : FC+SD term (n,m) component 530 implicit none 531#include "errquit.fh" 532#include "stdio.fh" 533#include "global.fh" 534 integer nbf,npts,i,j,k,m,n,a,b 535 double precision chi_ao(npts,nbf) 536 double precision delchi_ao(npts,3,nbf) 537 double precision fac2_arr(3,npts) 538 double precision ac_zpso(3,3), 539 & ac_fcsd(3,3) 540 double precision prod(3),prod1(3),prod13 541 integer ind_nab(2,3) 542 data ind_nab / 2, 3, ! nab=123 543 & 3, 1, ! nab=231 544 & 1, 2 / ! nab=312 545 do n=1,3 ! reset 546 do m=1,3 547 ac_zpso(n,m) = 0.0d0 548 ac_fcsd(n,m) = 0.0d0 549 enddo 550 enddo 551 do k = 1, npts 552 prod13 = chi_ao(k,i)*delchi_ao(k,1,j) 553 & +chi_ao(k,j)*delchi_ao(k,1,i) 554 ac_fcsd(1,1) = ac_fcsd(1,1)+fac2_arr(1,k)*prod13 555 ac_fcsd(1,2) = ac_fcsd(1,2)+fac2_arr(2,k)*prod13 556 ac_fcsd(1,3) = ac_fcsd(1,3)+fac2_arr(3,k)*prod13 557 prod13 = chi_ao(k,i)*delchi_ao(k,2,j) 558 & +chi_ao(k,j)*delchi_ao(k,2,i) 559 ac_fcsd(2,2) = ac_fcsd(2,2)+fac2_arr(2,k)*prod13 560 ac_fcsd(2,3) = ac_fcsd(2,3)+fac2_arr(3,k)*prod13 561 562 prod13 = chi_ao(k,i)*delchi_ao(k,3,j) 563 & +chi_ao(k,j)*delchi_ao(k,3,i) 564 ac_fcsd(3,3) = ac_fcsd(3,3)+fac2_arr(3,k)*prod13 565 enddo ! end-loo-k 566#if 0 567 do k = 1, npts 568 do n=1,3 569 prod(n) = chi_ao(k,i)*delchi_ao(k,n,j) 570 enddo ! end-loop-n 571 do n=1,3 572 a=ind_nab(1,n) 573 b=ind_nab(2,n) 574 ac_zpso(a,b) = ac_zpso(a,b)+fac2_arr(a,k)*prod(b) 575 ac_zpso(b,a) = ac_zpso(b,a)+fac2_arr(b,k)*prod(a) 576 enddo ! end-loop-n 577 enddo ! end-loo-k 578#else 579 do n=1,3 580 a=ind_nab(1,n) 581 b=ind_nab(2,n) 582 do k = 1, npts 583 ac_zpso(a,b) = ac_zpso(a,b)+fac2_arr(a,k)* 584 B chi_ao(k,i)*delchi_ao(k,b,j) 585 ac_zpso(b,a) = ac_zpso(b,a)+fac2_arr(b,k)* 586 A chi_ao(k,i)*delchi_ao(k,a,j) 587 enddo ! end-loop-n 588 enddo ! end-loo-k 589#endif 590 return 591 end 592 593 subroutine calc_NMRHFine_F1ij( 594 & ao_bas_han, ! in: AO basis handle 595 & geom, ! in: geometry handle 596 & ipol, ! in: nr. of polarizations 597 & g_dens, ! in: superposit. atomic densities 598 & delchi_ao, ! in: deriv. of basis functions 599 & qxyz, ! in: grid points 600 & qwght, ! in: weighting coeffs. 601 & nbf, ! in: nr. basis functions 602 & npts, ! in: nr. grid points 603 & natoms, ! in: nr. atoms 604 & ofinite, ! in: = .true. if Gaussian Nucl. Model of charges requested 605 & zetanuc_arr, ! in: sqrt(zetanuc(i)) i=1,natoms for Gaussian Nuclear Model 606 & Knucl, 607 & use_modelpotential, 608 & gexpo, 609 & gcoef, 610 & zsrkineticx, ! out 611 & zsrkineticy, ! out 612 & zsrkineticz) ! out 613c Purpose : Evaluates AO matrix for operator 614c [\vec{p} K x \vec{p}]_v v=1,2,3=x,y,z 615 implicit none 616#include "errquit.fh" 617#include "mafdecls.fh" 618#include "stdio.fh" 619#include "global.fh" 620#include "bas.fh" 621#include "zora.fh" 622 integer nbf,npts,ao_bas_han,natoms,geom 623 integer g_dens(2),ipol 624 double precision qwght(npts) 625 double precision qxyz(3,npts) 626 double precision delchi_ao(npts,3,nbf) 627 double precision zsrkineticx(nbf,nbf), 628 & zsrkineticy(nbf,nbf), 629 & zsrkineticz(nbf,nbf) 630 integer i,j,k,n 631 double precision amat_coul(npts,ipol) 632 double precision amat_nucl(npts) 633 integer ipt,closegridpts(npts) 634 double precision clight_au2,tol 635 double precision amat_tot,Kzora 636 double precision fac1_arr(npts) 637 double precision ac_hfineF1ji(3) 638 double precision chi_cntr(3,nbf) 639 logical ofinite,Knucl 640 double precision zetanuc_arr(natoms) 641 external get_ints_zora_hfine_F1ji, 642 & gridNuclearPotentialFinite, 643 & gridNuclearPotentialPoint 644c 645 logical use_modelpotential 646 double precision gexpo(natoms,50) 647 double precision gcoef(natoms,50) 648c 649 clight_au2 = clight_au*clight_au 650c == preliminaries == 651 do ipt = 1,npts 652 do i=1,ipol 653 amat_coul(ipt,i) = 0.d0 654 end do 655 amat_nucl(ipt) = 0.d0 656 closegridpts(ipt) = 0 657 end do 658c 659c == calculate the total hartree potential on the grid == 660 call gridHartreePotential(use_modelpotential, 661 & ao_bas_han, geom, natoms, ipol, g_dens, npts, qxyz, qwght, 662 & closegridpts, gexpo, gcoef, amat_coul) 663c 664c == calculate the total nuclear potential on the grid == 665 if (ofinite) then 666c 667c ------ Choosing Nuclear Model: erf(zetanuc^0.5 r_L) 668 call gridNuclearPotentialFinite(geom,natoms,npts,qxyz,qwght, 669 & zetanuc_arr, 670 & closegridpts,amat_nucl) 671c ------ Choosing Nuclear Model: P(1/2,zetanuc r_L^2) 672c call gridNuclearPotentialFinite2(geom,natoms,npts,qxyz,qwght, 673c & closegridpts,amat_nucl) 674 else ! default : point charge model for nuclei 675 call gridNuclearPotentialPoint(geom,natoms,npts,qxyz,qwght, 676 & closegridpts,amat_nucl) 677 endif 678 do k = 1,npts 679 if (k.eq.closegridpts(k)) qwght(k) = 0.d0 680 end do 681c === define fac_arr 682 if (Knucl) then !-- V=Vnucl (amat_tot) 683 do k = 1,npts 684c == assemble hartree and nuclear contributions == 685 amat_tot = amat_nucl(k) 686 Kzora=1.0d0/(1.0d0-0.5d0*amat_tot/clight_au2)-1.0d0 ! Alternative expression 687 ! gives same value as K for this AO 688 ! but it is suppose to cancel noise 689c +++++++++++++++++++++++++++++++++++++++++++++++++++++++ 690 if (do_NonRel) then ! remove it after TEST 691 Kzora=1.0d0 ! remove it after TEST 692 endif ! remove it after TEST 693c +++++++++++++++++++++++++++++++++++++++++++++++++++++++ 694 fac1_arr(k)=Kzora*qwght(k) 695 enddo ! end-loop-k 696 else ! default V=Vnucl+Vhartee 697 do k = 1,npts 698c == assemble hartree and nuclear contributions == 699 amat_tot = amat_nucl(k) + amat_coul(k,1) 700 Kzora=1.0d0/(1.0d0-0.5d0*amat_tot/clight_au2)-1.0d0 ! Alternative expression 701 ! gives same value as K for this AO 702 ! but it is suppose to cancel noise 703c +++++++++++++++++++++++++++++++++++++++++++++++++++++++ 704 if (do_NonRel) then ! remove it after TEST 705 Kzora=1.0d0 ! remove it after TEST 706 endif ! remove it after TEST 707c +++++++++++++++++++++++++++++++++++++++++++++++++++++++ 708 fac1_arr(k)=Kzora*qwght(k) 709 enddo ! end-loop-k 710 endif 711c == assemble zora correction == 712c ---- main diagonal -------- START 713 do i = 1, nbf 714 j=i 715 call get_ints_zora_hfine_F1ji(nbf,npts,delchi_ao,i,j, 716 & fac1_arr, 717 & ac_hfineF1ji) ! out 718 zsrkineticx(i,j) = zsrkineticx(i,j) + ac_hfineF1ji(1) 719 zsrkineticy(i,j) = zsrkineticy(i,j) + ac_hfineF1ji(2) 720 zsrkineticz(i,j) = zsrkineticz(i,j) + ac_hfineF1ji(3) 721 enddo ! end-loop-i 722c ---- main diagonal -------- END 723c ---- off diagonal -------- START 724 do i = 1, nbf 725 do j = i+1, nbf 726 call get_ints_zora_hfine_F1ji(nbf,npts,delchi_ao,i,j, 727 & fac1_arr, 728 & ac_hfineF1ji) ! out 729 zsrkineticx(i,j) = zsrkineticx(i,j) + 2.0d0*ac_hfineF1ji(1) 730 zsrkineticy(i,j) = zsrkineticy(i,j) + 2.0d0*ac_hfineF1ji(2) 731 zsrkineticz(i,j) = zsrkineticz(i,j) + 2.0d0*ac_hfineF1ji(3) 732 enddo ! end-loop-j 733 enddo ! end-loop-i 734c ---- off diagonal -------- END 735 return 736 end 737 738 subroutine get_ints_zora_hfine_F1ji( 739 & nbf, ! in: nr. basis functions 740 & npts, ! in: grid points 741 & delchi_ao, ! in: deriv. of basis functions 742 & i,j, ! in: (i,j) indices for delchi_ao 743 & fac1_arr, ! in 744 & ac_hfineF1ji)! out 745 implicit none 746#include "errquit.fh" 747#include "stdio.fh" 748#include "global.fh" 749 integer nbf,npts,k,n,i,j 750 double precision delchi_ao(npts,3,nbf) 751 double precision fac1_arr(npts) 752 double precision ac_hfineF1ji(3) 753 double precision prod(3) 754 755 do n=1,3 ! reset 756 ac_hfineF1ji(n) = 0.0d0 757 enddo 758 do k = 1, npts 759 prod(1)= delchi_ao(k,2,i)*delchi_ao(k,3,j) 760 & -delchi_ao(k,3,i)*delchi_ao(k,2,j) 761 prod(2)= delchi_ao(k,3,i)*delchi_ao(k,1,j) 762 & -delchi_ao(k,1,i)*delchi_ao(k,3,j) 763 prod(3)= delchi_ao(k,1,i)*delchi_ao(k,2,j) 764 & -delchi_ao(k,2,i)*delchi_ao(k,1,j) 765 do n=1,3 766 ac_hfineF1ji(n) = ac_hfineF1ji(n) + fac1_arr(k)*prod(n) 767 enddo ! end-loop-n 768 enddo ! end-loo-k 769 return 770 end 771c 772 subroutine get_Pnucl(amat_Pnucl, ! out: potential 773 & atmass, ! in : atomic mass 774 & xyz_NMRcoords, ! in : EFG-nuclear coord. 775 & a_coeff, ! in : =3/2 for AOs =1/2 for Vnucl 776 & nqpts, ! in : nr. grid points 777 & qxyz) 778c About: a_coeff, 779c =3/2 when adding finite size charge Gaussian model in evaluation 780c of hyperfine AOs (calc_zora_HFine) 781c =1/2 for Vnucl (in gridNuclearPotential) 782 implicit none 783#include "geom.fh" 784#include "global.fh" 785#include "msgids.fh" 786#include "stdio.fh" 787 integer i,igrid,nqpts 788 double precision xyz_NMRcoords(3) 789 double precision qxyz(3,nqpts) 790 double precision rxyz(3),dist,dist2,ac_prod 791 double precision amat_Pnucl(nqpts) 792 character*16 element 793 character*2 symbol 794 character*16 tags 795 logical is_atom 796 double precision atmass,zetanuc 797 double precision rtemp,a_coeff 798c 799 double precision dgami 800 external dgami, 801 & get_znuc 802c 803 call get_znuc(atmass,zetanuc) 804 do igrid = 1,nqpts 805 ac_prod=0.0d0 806 do i=1,3 807 rxyz(i) = qxyz(i,igrid)-xyz_NMRcoords(i) 808 ac_prod=ac_prod+rxyz(i)*rxyz(i) 809 enddo 810 rtemp = zetanuc*ac_prod ! dist*dist 811 amat_Pnucl(igrid) = dgami(a_coeff,rtemp) ! P(3/2,\tilde{r}_N^2) 812 end do ! igrid 813c 814 return 815 end 816c 817c------- get_Pnucl1() ------------ START 818c Purpose : Get one single value of Pnucl 819 subroutine get_Pnucl1(Pnucl, ! out: potential 820 & zetanuc_slc, ! in : zetanuc 821 & xyz_NMRcoords,! in : EFG-nuclear coord. 822 & a_coeff, ! in : =3/2 for AOs =1/2 for Vnucl 823 & nqpts, ! in : nr. grid points 824 & qxyz, ! in : one single grid point 825 & count_pt) ! TO CHECK 826c About: a_coeff, 827c =3/2 when adding finite size charge Gaussian model in evaluation 828c of hyperfine AOs (calc_zora_HFine) 829c =1/2 for Vnucl (in gridNuclearPotential) 830 implicit none 831#include "msgids.fh" 832#include "stdio.fh" 833#include "global.fh" 834 835 integer count_pt ! to check 836 837 integer i,igrid,nqpts 838 double precision xyz_NMRcoords(3) 839 double precision qxyz(3),Pnucl 840 double precision rxyz(3),dist,dist2,ac_prod 841 double precision zetanuc_slc 842 double precision rtemp,a_coeff 843 double precision dgami 844 external dgami ! Incomplete Gamma 845c ---------- Defining values at hand to check --------- START 846 847c xyz_NMRcoords(1)= 0.07090063d0 848c xyz_NMRcoords(2)= -0.12532286d0 849c xyz_NMRcoords(3)= 0.00000000d0 850c qxyz(1)= 0.07090403d0 851c qxyz(2)=-0.12532425d0 852c qxyz(3)=-0.00000339d0 853c zetanuc_slc=140130060.38598028d0 854cNWChem: Grid coord=(0.07090403, -0.12532425, -0.00000339) 855c Nuclpos =(0.07090063, -0.12532286, 0.00000000) 856c (zetanuc,rtemp,gammap)=(140130060.38598028,0.00350104,0.00015551) 857cADF: Grid coord=(0.07090403, -0.12532425, -0.00000339) 858c Nuclpos =(0.07090063, -0.12532286, 0.00000000) 859c (zetanuc,rtemp,gammap)=(140130060.38598028,0.00349682,0.00015523) 860c ---------- Defining values at hand to check --------- ENd 861 ac_prod=0.0d0 862 do i=1,3 863 rxyz(i) = qxyz(i)-xyz_NMRcoords(i) 864 ac_prod=ac_prod+rxyz(i)*rxyz(i) 865 enddo 866 rtemp = zetanuc_slc*ac_prod ! dist*dist 867 Pnucl = dgami(a_coeff,rtemp) ! P(3/2,\tilde{r}_N^2) 868 return 869 end 870c 871c------- get_Pnucl1() ------------ END 872c Purpose: Evaluation of zetanuc 873c to be used in evaluation of Incomplete 874c Gamma Function [gratio(...)] 875c rtemp = zetanuc*ac_prod ! dist*dist 876c call dgratio(threehalf,rtemp,gammap,gammaq,0) 877c 878 subroutine get_znuc(atmass, zetanuc) 879c 880 implicit none 881#include "msgids.fh" 882#include "stdio.fh" 883c 884 double precision atmass 885 double precision parnuc1,parnuc2, 886 & one,two,three,threehalf, 887 & fm2bohr,zetanuc,rtemp 888 data parnuc1 /0.836d0/ 889 data parnuc2 /0.570d0/ 890 data one /1.0d0/ 891 data two /2.0d0/ 892 data three /3.0d0/ 893 data threehalf /1.5d0/ 894 data fm2bohr /52917.7249d0/ 895 rtemp = (parnuc1 * atmass**(one/three) + parnuc2) / fm2bohr 896 zetanuc = three / ( two * (rtemp**2)) 897 return 898 end 899c 900c Purpose: Evaluation of zetanuc arr 901c to be used in evaluation of Incomplete 902c Gamma Function [gratio(...)] 903c rtemp = zetanuc*ac_prod ! dist*dist 904c call dgratio(threehalf,rtemp,gammap,gammaq,0) 905c This routine is used in gridNuclearPotential() 906c 907 subroutine get_zetanuc_arr(geom, natoms, zetanuc_arr) 908c 909 implicit none 910#include "errquit.fh" 911#include "mafdecls.fh" 912#include "geom.fh" 913#include "global.fh" 914#include "msgids.fh" 915#include "stdio.fh" 916 integer geom ! handle for geometry 917 integer i,natoms 918 double precision zetanuc_arr(natoms), 919 & atmass,zetanuc 920 external get_znuc 921 do i=1,natoms 922 if(.not.geom_mass_get(geom,i,atmass)) call 923 & errquit(' mass_get failed ',i,GEOM_ERR) 924 call get_znuc(atmass,zetanuc) 925 zetanuc_arr(i)=zetanuc 926 enddo ! end-loop-i 927 return 928 end 929c $Id$ 930