1C> \ingroup task 2C> @{ 3 logical function task_pderiv(rtdb) 4* 5* $Id$ 6* 7*********************************************************************** 8* This is about the nuclear coordinate derivatives of multicenter 9* integrals that we need. Notation: 10* i,j,p compound basis function (in practice, shell) indices, 11* including, specification of which center the function is on 12* 13* Xi, Xj orbital basis functions (chi), position of center not shown 14* explicitly here 15* 16* Gp density basis function, ditto 17* 18* R, R2 crystal lattice translation vectors 19* 20* rnuc position of some nucleus 21* 22* 1,2 electronic coordinates to be integrated (vol. elements d1,d2) 23* 24* I will just display the integrals whose derivatives we need, not the 25* differentiation itself. 26* 27* 1. orbital overlap Int d1 Xi(1) Xj(1-R) (i think we have this) eq. 78 28* 29* 2. Kinetic energy Int d1 Xi(1) [-del^2/2] Xj(1-R) (ditto) eq. 30 30* 31* 3. 1-center attraction -Sum(nuc)Z(nuc) Int d1 Gp(1-R) /|1-rnuc| eq. 60 32* 33* 4. 2-center repulsion Int d1 d2 Gp(1) (1/r12) Gp(2-R) [r12 obvious, 34* p = some other p index] eq. 52 35* 36* 5. 2-center attraction -Sum(nuc)Z(nuc) Int d1 Xi(1+R2)Xj(1+R2-R)/|1-rnuc| eq.42 37* 38* 6. 3-center repulsion Int d1 d2 Xi(1)Xj(1-R)(1/r12)Gp(2-R2) eq. 41 39* 40* All equations refer to our gross paper JCP 105, 10983 (1995), 41* doi:10.1063/1.472866 of which I believe you have a copy of, if not I will 42* bring one. 43* 44* Probably you will need some more details on how we want this, or we can 45* talk to you to understand how we get it, then we can reshuffle things as 46* necessary. For example do we get the derivative of any integral with 47* repect to any arbitrary nucleus that we specify, or with respect to all 48* nuclei in a single block of data, or only for those nuclei which are 49* relevant (nonzero deriv.) for a given integral? Call or message me as 50* needed and Ill come over.) 51* 52* Later 53* John 54*********************************************************************** 55* 56 57 implicit none 58#include "errquit.fh" 59#include "stdio.fh" 60#include "mafdecls.fh" 61#include "geom.fh" 62#include "bas.fh" 63#include "util.fh" 64c 65 logical int_normalize 66 external int_normalize 67 logical pderiv_compute_1e1cpe 68 external pderiv_compute_1e1cpe 69 logical pderiv_compute_1epe 70 external pderiv_compute_1epe 71 logical pderiv_compute_1eov 72 external pderiv_compute_1eov 73 logical pderiv_compute_1eke 74 external pderiv_compute_1eke 75 logical pderiv_compute_2e2c 76 external pderiv_compute_2e2c 77 logical pderiv_compute_2e3c 78 external pderiv_compute_2e3c 79 logical pderiv_compute_mpole 80 external pderiv_compute_mpole 81 logical pderiv_compute_3ov 82 external pderiv_compute_3ov 83 logical pderiv_compute_2e3ct, pderiv_compute_d2e3ct 84 external pderiv_compute_2e3ct, pderiv_compute_d2e3ct 85c 86 integer rtdb 87c 88 logical status 89 integer basis, geom 90 integer nat, nat3, size, maxg1, maxs1, maxg2, maxs2, maxg, maxs 91 integer nbf, nbfsq 92 integer hbuf, hbufp, hbufm, hscr, h1ec, h1efd, hxyz 93 integer kbuf, kbufp, kbufm, kscr, k1ec, k1efd, kxyz 94c 95 task_pderiv = .false. 96c 97 if (.not.geom_create(geom,'geometry')) call errquit 98 & ('geom create failed',911, GEOM_ERR) 99 if (.not.geom_rtdb_load(rtdb,geom,'geometry')) call errquit 100 & ('geom_rtdb_load failed',911, RTDB_ERR) 101c 102 if (.not.bas_create(basis,'ao basis')) call errquit 103 & ('bas_create failed',911, BASIS_ERR) 104 if (.not.bas_rtdb_load(rtdb,geom,basis,'ao basis')) call errquit 105 & ('bas_rtdb_load failed',911, RTDB_ERR) 106c 107 write(6,*)' geom/basis loaded' 108c 109 if (.not.int_normalize(rtdb,basis)) stop ' norm error 1' 110c 111 if (.not. bas_print(basis)) 112 $ call errquit(' basis print failed', 0, BASIS_ERR) 113c 114 if (.not.bas_numbf(basis,nbf)) call errquit 115 & ('numbf failed',911, BASIS_ERR) 116c 117 nbfsq = nbf*nbf 118 if (.not.geom_ncent(geom,nat)) stop 'geom_ncent fe' 119 write(6,*) 'number of atoms ', nat 120 nat3 = 3*nat 121 size = max(nbf,56) ! slmax = 35 for mpole 122 size = nat3*nbfsq*size 123 size = 2*size 124c 125 call intd_init(rtdb,1,basis) 126 call int_mem_print() 127 call int_mem_1e(maxg1,maxs1) 128 call int_mem_2e3c(maxg2,maxs2) 129 maxg2 = maxg2*3*4 130 maxg1 = maxg1*35 ! mpole 131 maxg = max(maxg1,maxg2) 132 maxs = max(maxs1,maxs2) 133 maxg = maxg + maxg/10 134 maxs = maxs + maxs/10 135 maxs = 2*maxs 136 maxg = 2*maxg 137 write(6,*)' normal maxg1 ',maxg1 138 write(6,*)' normal maxs1 ',maxs1 139 write(6,*)' normal maxg2 ',maxg2 140 write(6,*)' normal maxs2 ',maxs2 141 write(6,*)' normal maxg ',maxg 142 write(6,*)' normal maxs ',maxs 143 status = ma_alloc_get(mt_dbl,maxg,'int buffer' ,hbuf,kbuf) 144 status = status .and. 145 & ma_alloc_get(mt_dbl,maxg,'int buffer plus',hbufp,kbufp) 146 status = status .and. 147 & ma_alloc_get(mt_dbl,maxg,'int buffer minus',hbufm,kbufm) 148 status = status .and. 149 & ma_alloc_get(mt_dbl,maxs,'int scratch',hscr,kscr) 150 status = status .and. 151 & ma_alloc_get(mt_dbl,size,'block c',h1ec,k1ec) 152 status = status .and. 153 & ma_alloc_get(mt_dbl,size,'block fd',h1efd,k1efd) 154 status = status .and. 155 & ma_alloc_get(mt_dbl,3*nat,'my coords',hxyz,kxyz) 156 if (.not.status) stop ' memory alloc failed rak20 (1)' 157c 158 task_pderiv = pderiv_compute_1e1cpe( 159 & geom,basis,nbf,nat, 160 & maxg,maxs, 161 & dbl_mb(kbuf),dbl_mb(kbufp),dbl_mb(kbufm), 162 & dbl_mb(kscr),dbl_mb(k1ec),dbl_mb(k1efd), 163 & dbl_mb(kxyz)) 164c 165 task_pderiv = task_pderiv .and. 166 & pderiv_compute_1eov( 167 & geom,basis,nbf,nat, 168 & maxg,maxs, 169 & dbl_mb(kbuf),dbl_mb(kbufp),dbl_mb(kbufm), 170 & dbl_mb(kscr),dbl_mb(k1ec),dbl_mb(k1efd), 171 & dbl_mb(kxyz)) 172c 173 task_pderiv = task_pderiv .and. 174 & pderiv_compute_1eke( 175 & geom,basis,nbf,nat, 176 & maxg,maxs, 177 & dbl_mb(kbuf),dbl_mb(kbufp),dbl_mb(kbufm), 178 & dbl_mb(kscr),dbl_mb(k1ec),dbl_mb(k1efd), 179 & dbl_mb(kxyz)) 180c 181 task_pderiv = task_pderiv .and. 182 & pderiv_compute_1epe( 183 & geom,basis,nbf,nat, 184 & maxg,maxs, 185 & dbl_mb(kbuf),dbl_mb(kbufp),dbl_mb(kbufm), 186 & dbl_mb(kscr),dbl_mb(k1ec),dbl_mb(k1efd), 187 & dbl_mb(kxyz)) 188c 189 task_pderiv = task_pderiv .and. 190 & pderiv_compute_2e2c( 191 & geom,basis,nbf,nat, 192 & maxg,maxs, 193 & dbl_mb(kbuf),dbl_mb(kbufp),dbl_mb(kbufm), 194 & dbl_mb(kscr),dbl_mb(k1ec),dbl_mb(k1efd), 195 & dbl_mb(kxyz)) 196c 197 task_pderiv = task_pderiv .and. 198 & pderiv_compute_2e3c( 199 & geom,basis,nbf,nat, 200 & maxg,maxs, 201 & dbl_mb(kbuf),dbl_mb(kbufp),dbl_mb(kbufm), 202 & dbl_mb(kscr),dbl_mb(k1ec),dbl_mb(k1efd), 203 & dbl_mb(kxyz)) 204 205 task_pderiv = task_pderiv .and. 206 & pderiv_compute_3ov( 207 & geom,basis,nbf,nat, 208 & maxg,maxs, 209 & dbl_mb(kbuf),dbl_mb(kbufp),dbl_mb(kbufm), 210 & dbl_mb(kscr),dbl_mb(k1ec),dbl_mb(k1efd), 211 & dbl_mb(kxyz)) 212 213 task_pderiv = task_pderiv .and. 214 & pderiv_compute_mpole( 215 & geom,basis,nbf,nat, 216 & maxg,maxs, 217 & dbl_mb(kbuf),dbl_mb(kbufp),dbl_mb(kbufm), 218 & dbl_mb(kscr),dbl_mb(k1ec),dbl_mb(k1efd), 219 & dbl_mb(kxyz)) 220 221 task_pderiv = task_pderiv .and. 222 & pderiv_compute_2e3ct( 223 & geom,basis,nbf,nat, 224 & maxg,maxs, 225 & dbl_mb(kbuf),dbl_mb(kbufp),dbl_mb(kbufm), 226 & dbl_mb(kscr),dbl_mb(k1ec),dbl_mb(k1efd), 227 & dbl_mb(kxyz)) 228 229 task_pderiv = task_pderiv .and. 230 & pderiv_compute_d2e3ct( 231 & geom,basis,nbf,nat, 232 & maxg,maxs, 233 & dbl_mb(kbuf),dbl_mb(kbufp),dbl_mb(kbufm), 234 & dbl_mb(kscr),dbl_mb(k1ec),dbl_mb(k1efd), 235 & dbl_mb(kxyz)) 236 237 238 239 call intd_terminate() 240 if (.not.bas_destroy(basis)) stop ' bas destroy fail' 241 if (.not.geom_destroy(geom)) stop ' geom destroy fail' 242 status = .true. 243 status = status.and.ma_free_heap(hbuf) 244 status = status.and.ma_free_heap(hbufp) 245 status = status.and.ma_free_heap(hbufm) 246 status = status.and.ma_free_heap(hscr) 247 status = status.and.ma_free_heap(h1ec) 248 status = status.and.ma_free_heap(h1efd) 249 status = status.and.ma_free_heap(hxyz) 250 task_pderiv = status.and.task_pderiv 251c 252 end 253C> @} 254 logical function pderiv_compute_1e1cpe( 255 & geom,basis,nbf,nat,lbuf,lscr, 256 & buf,bufp,bufm,scr,p1c,p1fd,xyz) 257 implicit none 258#include "mafdecls.fh" 259#include "nwc_const.fh" 260#include "geomP.fh" 261#include "bas.fh" 262#include "stdio.fh" 263c 264 double precision ddot 265 external ddot 266c 267 integer geom 268 integer basis 269 integer nbf 270 integer nat 271 integer lbuf 272 integer lscr 273 double precision buf(lbuf), bufp(lbuf), bufm(lbuf) 274 double precision scr(lscr) 275 double precision p1c(nbf,3,nat) 276 double precision p1fd(nbf,3,nat) 277 double precision xyz(3,nat) 278c 279 integer nat3 280 integer atom, ixyz 281 double precision delta, factor, thresh, norm 282 double precision R(3) 283 integer nzero1, nzero2 284 integer ishell, nshell, ilo, ihi, nbfsh, cnt, i 285 integer IR 286c 287 call dfill(3,0.0d00,R,1) 288 call dfill((3*nat),0.0d00,xyz,1) 289 call dfill((nbf*3*nat),0.0d00,p1c,1) 290 call dfill((nbf*3*nat),0.0d00,p1fd,1) 291 call dfill(lbuf,0.0d00,buf,1) 292 call dfill(lbuf,0.0d00,bufp,1) 293 call dfill(lbuf,0.0d00,bufm,1) 294 call dfill(lscr,0.0d00,scr,1) 295 do IR = 1,2 296 if (IR.eq.1) call dfill(3,0.0d00,R,1) 297 if (IR.eq.2) then 298 R(1) = 1.0d00 299 R(2) = 2.0d00 300 R(3) = 3.0d00 301 endif 302 if (.not.ma_verify_allocator_stuff()) stop ' ma broke 1' 303 thresh = 1.0d-10 304 delta = 0.0001d00 305 factor = 1.0d00/(2.0d00*delta) 306 nat3 = 3*nat 307* store original coordintates 308 call dcopy(nat3,coords(1,1,geom),1,xyz,1) 309* 310 if (.not.bas_numcont(basis,nshell)) stop 'bas_numcont' 311 312 do ishell = 1, nshell 313 if (.not.bas_cn2bfr(basis,ishell,ilo,ihi)) 314 & stop 'cn2bfr error i' 315 nbfsh = ihi - ilo + 1 316 do atom = 1,nat 317 do ixyz = 1,3 318 call dcopy(nat3,xyz,1,coords(1,1,geom),1) 319 coords(ixyz,atom,geom) = coords(ixyz,atom,geom) + delta 320 call dfill(lscr,0.0d00,scr,1) 321 call dfill(lbuf,0.0d00,buf,1) 322 call dfill(lbuf,0.0d00,bufp,1) 323 call intp_1e1cpe(basis,ishell,R,lscr,scr,lbuf,bufp) 324* 325 call dcopy(nat3,xyz,1,coords(1,1,geom),1) 326 coords(ixyz,atom,geom) = coords(ixyz,atom,geom) - delta 327 call dfill(lscr,0.0d00,scr,1) 328 call dfill(lbuf,0.0d00,bufm,1) 329 call intp_1e1cpe(basis,ishell,R,lscr,scr,lbuf,bufm) 330* 331 call dcopy(nbfsh,bufp,1,buf,1) 332 call daxpy(nbfsh,-1.0d00,bufm,1,buf,1) 333 call dscal(nbfsh,factor,buf,1) 334 cnt = 1 335 do i = ilo,ihi 336 p1fd(i,ixyz,atom) = buf(cnt) 337 cnt = cnt + 1 338 enddo 339 enddo 340 enddo 341 enddo 342c 343 write(6,*)' fd: full list ' 344 nzero1 = 0 345 do atom = 1,nat 346 do ixyz = 1,3 347 do i = 1,nbf 348 if (abs(p1fd(i,ixyz,atom)).le.thresh) nzero1 = nzero1 + 1 349* write(6,10000)i,ixyz,atom,p1fd(i,ixyz,atom) 350 enddo 351 enddo 352 enddo 353 write(6,*)' fd: non-zero list ' 354 nzero2 = 0 355 do atom = 1,nat 356 do ixyz = 1,3 357 do i = 1,nbf 358 if (abs(p1fd(i,ixyz,atom)).gt.thresh) then 359* write(6,10000)i,ixyz,atom,p1fd(i,ixyz,atom) 360 continue 361 else 362 nzero2 = nzero2 + 1 363 endif 364 enddo 365 enddo 366 enddo 367 write(6,*)' fd: num zeros :1:', nzero1 368 write(6,*)' fd: num zeros :2:', nzero2, ':delta:', 369 & (nzero2-nzero1) 370 write(6,*)' fd: non-zero :1:', ((nbf*3*nat)-nzero1) 371 write(6,*)' fd: non-zero :2:', ((nbf*3*nat)-nzero2) 372 write(6,*)' fd: possible : :', (nbf*3*nat) 373c 374 call dcopy(nat3,xyz,1,coords(1,1,geom),1) 375 do ishell = 1,nshell 376 if (.not.bas_cn2bfr(basis,ishell,ilo,ihi)) 377 & stop 'cn2bfr error i' 378 call dfill(lscr,0.0d00,scr,1) 379 call dfill(lbuf,0.0d00,buf,1) 380 call intpd_1e1cpe(basis,ishell,R,lscr,scr,lbuf,buf) 381 cnt = 1 382 do atom=1,nat 383 do ixyz = 1,3 384 do i = ilo, ihi 385 p1c(i,ixyz,atom) = buf(cnt) 386 cnt = cnt + 1 387 enddo 388 enddo 389 enddo 390 enddo 391 write(6,*)' c: full list ' 392 nzero1 = 0 393 do atom = 1,nat 394 do ixyz = 1,3 395 do i = 1,nbf 396 if (abs(p1c(i,ixyz,atom)).le.thresh) nzero1 = nzero1 + 1 397* write(6,10001)i,ixyz,atom,p1c(i,ixyz,atom) 398 enddo 399 enddo 400 enddo 401 write(6,*)' c: non-zero list ' 402 nzero2 = 0 403 do atom = 1,nat 404 do ixyz = 1,3 405 do i = 1,nbf 406 if (abs(p1c(i,ixyz,atom)).gt.thresh) then 407* write(6,10001)i,ixyz,atom,p1c(i,ixyz,atom) 408 continue 409 else 410 nzero2 = nzero2 + 1 411 endif 412 enddo 413 enddo 414 enddo 415 write(6,*)' c: num zeros :1:', nzero1 416 write(6,*)' c: num zeros :2:', nzero2, ':delta:', 417 & (nzero2-nzero1) 418 write(6,*)' c: non-zero :1:', ((nbf*3*nat)-nzero1) 419 write(6,*)' c: non-zero :2:', ((nbf*3*nat)-nzero2) 420 write(6,*)' c: possible : :', (nbf*3*nat) 421c 42210000 format(1x,'p1fd(',i3,',',i2,',',i3,') = ',1pd20.10) 42310001 format(1x,' p1c(',i3,',',i2,',',i3,') = ',1pd20.10) 424c 425 call daxpy((nbf*3*nat),-1.0d00,p1fd,1,p1c,1) 426 norm = ddot((nbf*3*nat),p1c,1,p1c,1) 427 write(luout,*)' 1e1cpe difference norm:',ir,' ',norm 428c 429 pderiv_compute_1e1cpe = norm.lt.thresh 430 enddo 431 end 432 logical function pderiv_compute_1eov( 433 & geom,basis,nbf,nat,lbuf,lscr, 434 & buf,bufp,bufm,scr,ovc,ovfd,xyz) 435 implicit none 436#include "mafdecls.fh" 437#include "nwc_const.fh" 438#include "geomP.fh" 439#include "bas.fh" 440#include "stdio.fh" 441c 442 double precision ddot 443 external ddot 444c 445 integer geom 446 integer basis 447 integer nbf 448 integer nat 449 integer lbuf 450 integer lscr 451 double precision buf(lbuf), bufp(lbuf), bufm(lbuf) 452 double precision scr(lscr) 453 double precision ovc(nbf,nbf,3,nat) 454 double precision ovfd(nbf,nbf,3,nat) 455 double precision xyz(3,nat) 456c 457 integer nat3 458 integer atom, ixyz, IR 459 double precision delta, factor, thresh, norm 460 double precision R(3) 461 integer nzero1, nzero2 462 integer nshell, cnt, nbfsh, watom(2) 463 integer ishell, ilo, ihi, nbfshi, i 464 integer jshell, jlo, jhi, nbfshj, j 465c 466 call dfill(3,0.0d00,R,1) 467 call dfill((3*nat),0.0d00,xyz,1) 468 call dfill((nbf*nbf*3*nat),0.0d00,ovc,1) 469 call dfill((nbf*nbf*3*nat),0.0d00,ovfd,1) 470 call dfill(lbuf,0.0d00,buf,1) 471 call dfill(lbuf,0.0d00,bufp,1) 472 call dfill(lbuf,0.0d00,bufm,1) 473 call dfill(lscr,0.0d00,scr,1) 474 do IR = 1,2 475 if (IR.eq.1) call dfill(3,0.0d00,R,1) 476 if (IR.eq.2) then 477 R(1) = 1.0d00 478 R(2) = 2.0d00 479 R(3) = 3.0d00 480 endif 481 if (.not.ma_verify_allocator_stuff()) stop ' ma broke 2' 482 thresh = 1.0d-10 483 delta = 0.0001d00 484 factor = 1.0d00/(2.0d00*delta) 485 nat3 = 3*nat 486* store original coordintates 487 call dcopy(nat3,coords(1,1,geom),1,xyz,1) 488* 489 if (.not.bas_numcont(basis,nshell)) stop 'bas_numcont' 490 491 do ishell = 1, nshell 492 if (.not.bas_cn2bfr(basis,ishell,ilo,ihi)) 493 & stop 'cn2bfr error i' 494 nbfshi = ihi - ilo + 1 495 do jshell = 1, ishell 496 if (.not.bas_cn2bfr(basis,jshell,jlo,jhi)) 497 & stop 'cn2bfr error j' 498 nbfshj = jhi - jlo + 1 499 nbfsh = nbfshi*nbfshj 500 do atom = 1,nat 501 do ixyz = 1,3 502 call dcopy(nat3,xyz,1,coords(1,1,geom),1) 503 coords(ixyz,atom,geom) = coords(ixyz,atom,geom) + delta 504 call dfill(lscr,0.0d00,scr,1) 505 call dfill(lbuf,0.0d00,buf,1) 506 call dfill(lbuf,0.0d00,bufp,1) 507 call intp_1eov(basis,ishell,basis,jshell, 508 & R,lscr,scr,lbuf,bufp) 509* 510 call dcopy(nat3,xyz,1,coords(1,1,geom),1) 511 coords(ixyz,atom,geom) = coords(ixyz,atom,geom) - delta 512 call dfill(lscr,0.0d00,scr,1) 513 call dfill(lbuf,0.0d00,bufm,1) 514 call intp_1eov(basis,ishell,basis,jshell, 515 & R,lscr,scr,lbuf,bufm) 516* 517 call dcopy(nbfsh,bufp,1,buf,1) 518 call daxpy(nbfsh,-1.0d00,bufm,1,buf,1) 519 call dscal(nbfsh,factor,buf,1) 520 cnt = 1 521 do i = ilo,ihi 522 do j = jlo, jhi 523 ovfd(i,j,ixyz,atom) = buf(cnt) 524 ovfd(j,i,ixyz,atom) = buf(cnt) 525 cnt = cnt + 1 526 enddo 527 enddo 528 if (.not.ma_verify_allocator_stuff()) 529 & stop ' ma broke 3' 530 enddo 531 enddo 532 enddo 533 enddo 534c 535 write(6,*)' fd: full list i,j' 536 nzero1 = 0 537 do atom = 1,nat 538 do ixyz = 1,3 539 do i = 1,nbf 540 do j = 1,nbf 541 if (abs(ovfd(i,j,ixyz,atom)).le.thresh) 542 & nzero1 = nzero1 + 1 543* write(6,10000)i,j,ixyz,atom,ovfd(i,j,ixyz,atom) 544 enddo 545 enddo 546 enddo 547 enddo 548 write(6,*)' fd: non-zero list ' 549 nzero2 = 0 550 do atom = 1,nat 551 do ixyz = 1,3 552 do i = 1,nbf 553 do j = 1,nbf 554 if (abs(ovfd(i,j,ixyz,atom)).gt.thresh) then 555* write(6,10000)i,j,ixyz,atom,ovfd(i,j,ixyz,atom) 556 continue 557 else 558 nzero2 = nzero2 + 1 559 endif 560 enddo 561 enddo 562 enddo 563 enddo 564 write(6,*)' fd: num zeros :1:', nzero1 565 write(6,*)' fd: num zeros :2:', nzero2, ':delta:', 566 & (nzero2-nzero1) 567 write(6,*)' fd: non-zero :1:', ((nbf*nbf*3*nat)-nzero1) 568 write(6,*)' fd: non-zero :2:', ((nbf*nbf*3*nat)-nzero2) 569 write(6,*)' fd: possible : :', (nbf*nbf*3*nat) 570 nzero1 = 0 571 nzero2 = 0 572c 573 call dcopy(nat3,xyz,1,coords(1,1,geom),1) 574 do ishell = 1,nshell 575 if (.not.bas_cn2bfr(basis,ishell,ilo,ihi)) 576 & stop 'cn2bfr error i' 577 do jshell = 1,ishell 578 if (.not.bas_cn2bfr(basis,jshell,jlo,jhi)) 579 & stop 'cn2bfr error j' 580 call dfill(lscr,0.0d00,scr,1) 581 call dfill(lbuf,0.0d00,buf,1) 582 call intpd_1eov(basis,ishell,basis,jshell,R,lscr,scr, 583 & lbuf,buf,watom) 584* write(6,*)'watom 1eov',watom 585 cnt = 1 586 if (.not.ma_verify_allocator_stuff()) stop ' ma broke 4' 587 do atom=1,2 588 if (watom(atom).gt.0) then 589 do ixyz = 1,3 590 do i = ilo, ihi 591 do j = jlo, jhi 592*buffer is <jlo:jhi, ilo:ihi, 1:3, 2> 593 ovc(i,j,ixyz,watom(atom)) = buf(cnt) 594 ovc(j,i,ixyz,watom(atom)) = buf(cnt) 595 cnt = cnt + 1 596 enddo 597 enddo 598 enddo 599 endif 600 enddo 601 if (.not.ma_verify_allocator_stuff()) stop ' ma broke 5' 602 enddo 603 enddo 604 write(6,*)' c: full list i,j' 605 nzero1 = 0 606 do atom = 1,nat 607 do ixyz = 1,3 608 do i = 1,nbf 609 do j = 1,nbf 610 if (abs(ovc(i,j,ixyz,atom)).le.thresh) 611 & nzero1 = nzero1 + 1 612* write(6,10001)i,j,ixyz,atom,ovc(i,j,ixyz,atom) 613 enddo 614 enddo 615 enddo 616 enddo 617 write(6,*)' c: non-zero list ' 618 nzero2 = 0 619 do atom = 1,nat 620 do ixyz = 1,3 621 do i = 1,nbf 622 do j = 1,nbf 623 if (abs(ovc(i,j,ixyz,atom)).gt.thresh) then 624* write(6,10001)i,j,ixyz,atom,ovc(i,j,ixyz,atom) 625 continue 626 else 627 nzero2 = nzero2 + 1 628 endif 629 enddo 630 enddo 631 enddo 632 enddo 633 write(6,*)' c: num zeros :1:', nzero1 634 write(6,*)' c: num zeros :2:', nzero2, ':delta:', 635 & (nzero2-nzero1) 636 write(6,*)' c: non-zero :1:', ((nbf*nbf*3*nat)-nzero1) 637 write(6,*)' c: non-zero :2:', ((nbf*nbf*3*nat)-nzero2) 638 write(6,*)' c: possible : :', (nbf*nbf*3*nat) 639c 64010000 format(1x,'ovfd(',i3,',',i3,',',i2,',',i3,') = ',1pd20.10) 64110001 format(1x,' ovc(',i3,',',i3,',',i2,',',i3,') = ',1pd20.10) 642c 643 call daxpy((nbf*nbf*3*nat),-1.0d00,ovfd,1,ovc,1) 644 norm = ddot((nbf*nbf*3*nat),ovc,1,ovc,1) 645 write(luout,*)' 1eov difference norm:',ir,' ',norm 646c 647 pderiv_compute_1eov = norm.lt.thresh 648 enddo 649 end 650 logical function pderiv_compute_1eke( 651 & geom,basis,nbf,nat,lbuf,lscr, 652 & buf,bufp,bufm,scr,kec,kefd,xyz) 653 implicit none 654#include "mafdecls.fh" 655#include "nwc_const.fh" 656#include "geomP.fh" 657#include "bas.fh" 658#include "stdio.fh" 659c 660 double precision ddot 661 external ddot 662c 663 integer geom 664 integer basis 665 integer nbf 666 integer nat 667 integer lbuf 668 integer lscr 669 double precision buf(lbuf), bufp(lbuf), bufm(lbuf) 670 double precision scr(lscr) 671 double precision kec(nbf,nbf,3,nat) 672 double precision kefd(nbf,nbf,3,nat) 673 double precision xyz(3,nat) 674c 675 integer nat3 676 integer atom, ixyz, IR 677 double precision delta, factor, thresh, norm 678 double precision R(3) 679 integer nzero1, nzero2 680 integer nshell, cnt, nbfsh, watom(2) 681 integer ishell, ilo, ihi, nbfshi, i 682 integer jshell, jlo, jhi, nbfshj, j 683c 684 call dfill(3,0.0d00,R,1) 685 call dfill((3*nat),0.0d00,xyz,1) 686 call dfill((nbf*nbf*3*nat),0.0d00,kec,1) 687 call dfill((nbf*nbf*3*nat),0.0d00,kefd,1) 688 call dfill(lbuf,0.0d00,buf,1) 689 call dfill(lbuf,0.0d00,bufp,1) 690 call dfill(lbuf,0.0d00,bufm,1) 691 call dfill(lscr,0.0d00,scr,1) 692 do IR = 1,2 693 if (IR.eq.1) call dfill(3,0.0d00,R,1) 694 if (IR.eq.2) then 695 R(1) = 1.0d00 696 R(2) = 2.0d00 697 R(3) = 3.0d00 698 endif 699 if (.not.ma_verify_allocator_stuff()) stop ' ma broke 6' 700 thresh = 1.0d-10 701 delta = 0.0001d00 702 factor = 1.0d00/(2.0d00*delta) 703 nat3 = 3*nat 704* store original coordintates 705 call dcopy(nat3,coords(1,1,geom),1,xyz,1) 706* 707 if (.not.bas_numcont(basis,nshell)) stop 'bas_numcont' 708 709 do ishell = 1, nshell 710 if (.not.bas_cn2bfr(basis,ishell,ilo,ihi)) 711 & stop 'cn2bfr error i' 712 nbfshi = ihi - ilo + 1 713 do jshell = 1, ishell 714 if (.not.bas_cn2bfr(basis,jshell,jlo,jhi)) 715 & stop 'cn2bfr error j' 716 nbfshj = jhi - jlo + 1 717 nbfsh = nbfshi*nbfshj 718 do atom = 1,nat 719 do ixyz = 1,3 720 call dcopy(nat3,xyz,1,coords(1,1,geom),1) 721 coords(ixyz,atom,geom) = coords(ixyz,atom,geom) + delta 722 call dfill(lscr,0.0d00,scr,1) 723 call dfill(lbuf,0.0d00,buf,1) 724 call dfill(lbuf,0.0d00,bufp,1) 725 call intp_1eke(basis,ishell,basis,jshell, 726 & R,lscr,scr,lbuf,bufp) 727* 728 call dcopy(nat3,xyz,1,coords(1,1,geom),1) 729 coords(ixyz,atom,geom) = coords(ixyz,atom,geom) - delta 730 call dfill(lscr,0.0d00,scr,1) 731 call dfill(lbuf,0.0d00,bufm,1) 732 call intp_1eke(basis,ishell,basis,jshell, 733 & R,lscr,scr,lbuf,bufm) 734* 735 call dcopy(nbfsh,bufp,1,buf,1) 736 call daxpy(nbfsh,-1.0d00,bufm,1,buf,1) 737 call dscal(nbfsh,factor,buf,1) 738 cnt = 1 739 do i = ilo,ihi 740 do j = jlo, jhi 741 kefd(i,j,ixyz,atom) = buf(cnt) 742 kefd(j,i,ixyz,atom) = buf(cnt) 743 cnt = cnt + 1 744 enddo 745 enddo 746 if (.not.ma_verify_allocator_stuff()) 747 & stop ' ma broke 7' 748 enddo 749 enddo 750 enddo 751 enddo 752c 753 write(6,*)' fd: full list i,j' 754 nzero1 = 0 755 do atom = 1,nat 756 do ixyz = 1,3 757 do i = 1,nbf 758 do j = 1,nbf 759 if (abs(kefd(i,j,ixyz,atom)).le.thresh) 760 & nzero1 = nzero1 + 1 761* write(6,10000)i,j,ixyz,atom,kefd(i,j,ixyz,atom) 762 enddo 763 enddo 764 enddo 765 enddo 766 write(6,*)' fd: non-zero list ' 767 nzero2 = 0 768 do atom = 1,nat 769 do ixyz = 1,3 770 do i = 1,nbf 771 do j = 1,nbf 772 if (abs(kefd(i,j,ixyz,atom)).gt.thresh) then 773* write(6,10000)i,j,ixyz,atom,kefd(i,j,ixyz,atom) 774 continue 775 else 776 nzero2 = nzero2 + 1 777 endif 778 enddo 779 enddo 780 enddo 781 enddo 782 write(6,*)' fd: num zeros :1:', nzero1 783 write(6,*)' fd: num zeros :2:', nzero2, ':delta:', 784 & (nzero2-nzero1) 785 write(6,*)' fd: non-zero :1:', ((nbf*nbf*3*nat)-nzero1) 786 write(6,*)' fd: non-zero :2:', ((nbf*nbf*3*nat)-nzero2) 787 write(6,*)' fd: possible : :', (nbf*nbf*3*nat) 788 nzero1 = 0 789 nzero2 = 0 790c 791 call dcopy(nat3,xyz,1,coords(1,1,geom),1) 792 do ishell = 1,nshell 793 if (.not.bas_cn2bfr(basis,ishell,ilo,ihi)) 794 & stop 'cn2bfr error i' 795 do jshell = 1,ishell 796 if (.not.bas_cn2bfr(basis,jshell,jlo,jhi)) 797 & stop 'cn2bfr error j' 798 call dfill(lscr,0.0d00,scr,1) 799 call dfill(lbuf,0.0d00,buf,1) 800 call intpd_1eke(basis,ishell,basis,jshell,R,lscr,scr, 801 & lbuf,buf,watom) 802* write(6,*)'watom 1eke',watom 803 cnt = 1 804 if (.not.ma_verify_allocator_stuff()) stop ' ma broke 8' 805 do atom=1,2 806 if (watom(atom).gt.0) then 807 do ixyz = 1,3 808 do i = ilo, ihi 809 do j = jlo, jhi 810 kec(i,j,ixyz,watom(atom)) = buf(cnt) 811 kec(j,i,ixyz,watom(atom)) = buf(cnt) 812 cnt = cnt + 1 813 enddo 814 enddo 815 enddo 816 endif 817 enddo 818 if (.not.ma_verify_allocator_stuff()) stop ' ma broke 9' 819 enddo 820 enddo 821 write(6,*)' c: full list i,j' 822 nzero1 = 0 823 do atom = 1,nat 824 do ixyz = 1,3 825 do i = 1,nbf 826 do j = 1,nbf 827 if (abs(kec(i,j,ixyz,atom)).le.thresh) 828 & nzero1 = nzero1 + 1 829* write(6,10001)i,j,ixyz,atom,kec(i,j,ixyz,atom) 830 enddo 831 enddo 832 enddo 833 enddo 834 write(6,*)' c: non-zero list ' 835 nzero2 = 0 836 do atom = 1,nat 837 do ixyz = 1,3 838 do i = 1,nbf 839 do j = 1,nbf 840 if (abs(kec(i,j,ixyz,atom)).gt.thresh) then 841* write(6,10001)i,j,ixyz,atom,kec(i,j,ixyz,atom) 842 continue 843 else 844 nzero2 = nzero2 + 1 845 endif 846 enddo 847 enddo 848 enddo 849 enddo 850 write(6,*)' c: num zeros :1:', nzero1 851 write(6,*)' c: num zeros :2:', nzero2, ':delta:', 852 & (nzero2-nzero1) 853 write(6,*)' c: non-zero :1:', ((nbf*nbf*3*nat)-nzero1) 854 write(6,*)' c: non-zero :2:', ((nbf*nbf*3*nat)-nzero2) 855 write(6,*)' c: possible : :', (nbf*nbf*3*nat) 856c 85710000 format(1x,'kefd(',i3,',',i3,',',i2,',',i3,') = ',1pd20.10) 85810001 format(1x,' kec(',i3,',',i3,',',i2,',',i3,') = ',1pd20.10) 859c 860 call daxpy((nbf*nbf*3*nat),-1.0d00,kefd,1,kec,1) 861 norm = ddot((nbf*nbf*3*nat),kec,1,kec,1) 862 write(luout,*)' 1eke difference norm:',ir,' ',norm 863c 864 pderiv_compute_1eke = norm.lt.thresh 865 enddo 866 end 867 logical function pderiv_compute_1epe( 868 & geom,basis,nbf,nat,lbuf,lscr, 869 & buf,bufp,bufm,scr,pec,pefd,xyz) 870 implicit none 871#include "mafdecls.fh" 872#include "nwc_const.fh" 873#include "geomP.fh" 874#include "bas.fh" 875#include "stdio.fh" 876c 877 double precision ddot 878 external ddot 879c 880 integer geom 881 integer basis 882 integer nbf 883 integer nat 884 integer lbuf 885 integer lscr 886 double precision buf(lbuf), bufp(lbuf), bufm(lbuf) 887 double precision scr(lscr) 888 double precision pec(nbf,nbf,3,nat) 889 double precision pefd(nbf,nbf,3,nat) 890 double precision xyz(3,nat) 891c 892 integer nat3 893 integer atom, ixyz, IR 894 double precision delta, factor, thresh, norm 895 double precision R1(3), R2(3) 896 integer nzero1, nzero2 897 integer nshell, cnt, nbfsh 898 integer ishell, ilo, ihi, nbfshi, i 899 integer jshell, jlo, jhi, nbfshj, j 900c 901 call dfill(3,0.0d00,R1,1) 902 call dfill(3,0.0d00,R2,1) 903 call dfill((3*nat),0.0d00,xyz,1) 904 call dfill((nbf*nbf*3*nat),0.0d00,pec,1) 905 call dfill((nbf*nbf*3*nat),0.0d00,pefd,1) 906 call dfill(lbuf,0.0d00,buf,1) 907 call dfill(lbuf,0.0d00,bufp,1) 908 call dfill(lbuf,0.0d00,bufm,1) 909 call dfill(lscr,0.0d00,scr,1) 910 do IR = 1,3 911 if (IR.eq.1) then 912 call dfill(3,0.0d00,R1,1) 913 call dfill(3,0.0d00,R2,1) 914 elseif (IR.eq.2) then 915 call dfill(3,0.0d00,R2,1) 916 R1(1) = 1.0d00 917 R1(2) = 2.0d00 918 R1(3) = 3.0d00 919 elseif (IR.eq.3) then 920 R1(1) = 1.0d00 921 R1(2) = 2.0d00 922 R1(3) = 3.0d00 923 R2(1) = 3.0d00 924 R2(2) = 4.0d00 925 R2(3) = 5.0d00 926 else 927 stop ' how did IR get here' 928 endif 929 if (.not.ma_verify_allocator_stuff()) stop ' ma broke 10' 930 thresh = 1.0d-10 931 delta = 0.0001d00 932 factor = 1.0d00/(2.0d00*delta) 933 nat3 = 3*nat 934* store original coordintates 935 call dcopy(nat3,coords(1,1,geom),1,xyz,1) 936* 937 if (.not.bas_numcont(basis,nshell)) stop 'bas_numcont' 938 939 do ishell = 1, nshell 940 if (.not.bas_cn2bfr(basis,ishell,ilo,ihi)) 941 & stop 'cn2bfr error i' 942 nbfshi = ihi - ilo + 1 943 do jshell = 1, ishell 944 if (.not.bas_cn2bfr(basis,jshell,jlo,jhi)) 945 & stop 'cn2bfr error j' 946 nbfshj = jhi - jlo + 1 947 nbfsh = nbfshi*nbfshj 948 do atom = 1,nat 949 do ixyz = 1,3 950 call dcopy(nat3,xyz,1,coords(1,1,geom),1) 951 coords(ixyz,atom,geom) = coords(ixyz,atom,geom) + delta 952 call dfill(lscr,0.0d00,scr,1) 953 call dfill(lbuf,0.0d00,buf,1) 954 call dfill(lbuf,0.0d00,bufp,1) 955 call intp_1epe(basis,ishell,R1,basis,jshell, 956 & R2,lscr,scr,lbuf,bufp) 957* 958 call dcopy(nat3,xyz,1,coords(1,1,geom),1) 959 coords(ixyz,atom,geom) = coords(ixyz,atom,geom) - delta 960 call dfill(lscr,0.0d00,scr,1) 961 call dfill(lbuf,0.0d00,bufm,1) 962 call intp_1epe(basis,ishell,R1,basis,jshell, 963 & R2,lscr,scr,lbuf,bufm) 964* 965 call dcopy(nbfsh,bufp,1,buf,1) 966 call daxpy(nbfsh,-1.0d00,bufm,1,buf,1) 967 call dscal(nbfsh,factor,buf,1) 968 cnt = 1 969 do i = ilo,ihi 970 do j = jlo, jhi 971 pefd(i,j,ixyz,atom) = buf(cnt) 972 pefd(j,i,ixyz,atom) = buf(cnt) 973 cnt = cnt + 1 974 enddo 975 enddo 976 if (.not.ma_verify_allocator_stuff()) 977 & stop ' ma broke 11' 978 enddo 979 enddo 980 enddo 981 enddo 982c 983 write(6,*)' fd: full list i,j' 984 nzero1 = 0 985 do atom = 1,nat 986 do ixyz = 1,3 987 do i = 1,nbf 988 do j = 1,nbf 989 if (abs(pefd(i,j,ixyz,atom)).le.thresh) 990 & nzero1 = nzero1 + 1 991* write(6,10000)i,j,ixyz,atom,pefd(i,j,ixyz,atom) 992 enddo 993 enddo 994 enddo 995 enddo 996 write(6,*)' fd: non-zero list ' 997 nzero2 = 0 998 do atom = 1,nat 999 do ixyz = 1,3 1000 do i = 1,nbf 1001 do j = 1,nbf 1002 if (abs(pefd(i,j,ixyz,atom)).gt.thresh) then 1003* write(6,10000)i,j,ixyz,atom,pefd(i,j,ixyz,atom) 1004 continue 1005 else 1006 nzero2 = nzero2 + 1 1007 endif 1008 enddo 1009 enddo 1010 enddo 1011 enddo 1012 write(6,*)' fd: num zeros :1:', nzero1 1013 write(6,*)' fd: num zeros :2:', nzero2, ':delta:', 1014 & (nzero2-nzero1) 1015 write(6,*)' fd: non-zero :1:', ((nbf*nbf*3*nat)-nzero1) 1016 write(6,*)' fd: non-zero :2:', ((nbf*nbf*3*nat)-nzero2) 1017 write(6,*)' fd: possible : :', (nbf*nbf*3*nat) 1018 nzero1 = 0 1019 nzero2 = 0 1020c 1021 call dcopy(nat3,xyz,1,coords(1,1,geom),1) 1022 do ishell = 1,nshell 1023 if (.not.bas_cn2bfr(basis,ishell,ilo,ihi)) 1024 & stop 'cn2bfr error i' 1025 do jshell = 1,ishell 1026 if (.not.bas_cn2bfr(basis,jshell,jlo,jhi)) 1027 & stop 'cn2bfr error j' 1028 call dfill(lscr,0.0d00,scr,1) 1029 call dfill(lbuf,0.0d00,buf,1) 1030 call intpd_1epe(basis,ishell,R1,basis,jshell,R2, 1031 & lscr,scr,lbuf,buf) 1032 cnt = 1 1033 if (.not.ma_verify_allocator_stuff()) stop ' ma broke 12' 1034 do atom=1,nat 1035 do ixyz = 1,3 1036 do i = ilo, ihi 1037 do j = jlo, jhi 1038 pec(i,j,ixyz,atom) = buf(cnt) 1039 pec(j,i,ixyz,atom) = buf(cnt) 1040 cnt = cnt + 1 1041 enddo 1042 enddo 1043 enddo 1044 enddo 1045 if (.not.ma_verify_allocator_stuff()) stop ' ma broke 13' 1046 enddo 1047 enddo 1048 write(6,*)' c: full list i,j' 1049 nzero1 = 0 1050 do atom = 1,nat 1051 do ixyz = 1,3 1052 do i = 1,nbf 1053 do j = 1,nbf 1054 if (abs(pec(i,j,ixyz,atom)).le.thresh) 1055 & nzero1 = nzero1 + 1 1056* write(6,10001)i,j,ixyz,atom,pec(i,j,ixyz,atom) 1057 enddo 1058 enddo 1059 enddo 1060 enddo 1061 write(6,*)' c: non-zero list ' 1062 nzero2 = 0 1063 do atom = 1,nat 1064 do ixyz = 1,3 1065 do i = 1,nbf 1066 do j = 1,nbf 1067 if (abs(pec(i,j,ixyz,atom)).gt.thresh) then 1068* write(6,10001)i,j,ixyz,atom,pec(i,j,ixyz,atom) 1069 continue 1070 else 1071 nzero2 = nzero2 + 1 1072 endif 1073 enddo 1074 enddo 1075 enddo 1076 enddo 1077 write(6,*)' c: num zeros :1:', nzero1 1078 write(6,*)' c: num zeros :2:', nzero2, ':delta:', 1079 & (nzero2-nzero1) 1080 write(6,*)' c: non-zero :1:', ((nbf*nbf*3*nat)-nzero1) 1081 write(6,*)' c: non-zero :2:', ((nbf*nbf*3*nat)-nzero2) 1082 write(6,*)' c: possible : :', (nbf*nbf*3*nat) 1083c 108410000 format(1x,'pefd(',i3,',',i3,',',i2,',',i3,') = ',1pd20.10) 108510001 format(1x,' pec(',i3,',',i3,',',i2,',',i3,') = ',1pd20.10) 1086c 1087 call daxpy((nbf*nbf*3*nat),-1.0d00,pefd,1,pec,1) 1088 norm = ddot((nbf*nbf*3*nat),pec,1,pec,1) 1089 write(luout,*)' 1epe difference norm:',IR,' ',norm 1090c 1091 pderiv_compute_1epe = norm.lt.thresh 1092 enddo 1093 end 1094 logical function pderiv_compute_2e2c( 1095 & geom,basis,nbf,nat,lbuf,lscr, 1096 & buf,bufp,bufm,scr,eri2c,eri2fd,xyz) 1097 implicit none 1098#include "mafdecls.fh" 1099#include "nwc_const.fh" 1100#include "geomP.fh" 1101#include "bas.fh" 1102#include "stdio.fh" 1103c 1104 double precision ddot 1105 external ddot 1106c 1107 integer geom 1108 integer basis 1109 integer nbf 1110 integer nat 1111 integer lbuf 1112 integer lscr 1113 double precision buf(lbuf), bufp(lbuf), bufm(lbuf) 1114 double precision scr(lscr) 1115 double precision eri2c(nbf,nbf,3,nat) 1116 double precision eri2fd(nbf,nbf,3,nat) 1117 double precision xyz(3,nat) 1118c 1119 integer nat3 1120 integer atom, ixyz, IR 1121 double precision delta, factor, thresh, norm 1122 double precision R(3) 1123 integer nzero1, nzero2 1124 integer nshell, cnt, nbfsh, watom(2) 1125 integer ishell, ilo, ihi, nbfshi, i 1126 integer jshell, jlo, jhi, nbfshj, j 1127c 1128 call dfill(3,0.0d00,R,1) 1129 call dfill((3*nat),0.0d00,xyz,1) 1130 call dfill((nbf*nbf*3*nat),0.0d00,eri2c,1) 1131 call dfill((nbf*nbf*3*nat),0.0d00,eri2fd,1) 1132 call dfill(lbuf,0.0d00,buf,1) 1133 call dfill(lbuf,0.0d00,bufp,1) 1134 call dfill(lbuf,0.0d00,bufm,1) 1135 call dfill(lscr,0.0d00,scr,1) 1136 do IR = 1,2 1137 if (IR.eq.1) call dfill(3,0.0d00,R,1) 1138 if (IR.eq.2) then 1139 R(1) = 1.0d00 1140 R(2) = 2.0d00 1141 R(3) = 3.0d00 1142 endif 1143 if (.not.ma_verify_allocator_stuff()) stop ' ma broke 14' 1144 thresh = 1.0d-10 1145 delta = 0.0001d00 1146 factor = 1.0d00/(2.0d00*delta) 1147 nat3 = 3*nat 1148* store original coordintates 1149 call dcopy(nat3,coords(1,1,geom),1,xyz,1) 1150* 1151 if (.not.bas_numcont(basis,nshell)) stop 'bas_numcont' 1152 1153 do ishell = 1, nshell 1154 if (.not.bas_cn2bfr(basis,ishell,ilo,ihi)) 1155 & stop 'cn2bfr error i' 1156 nbfshi = ihi - ilo + 1 1157 do jshell = 1, ishell 1158 if (.not.bas_cn2bfr(basis,jshell,jlo,jhi)) 1159 & stop 'cn2bfr error j' 1160 nbfshj = jhi - jlo + 1 1161 nbfsh = nbfshi*nbfshj 1162 do atom = 1,nat 1163 do ixyz = 1,3 1164 call dcopy(nat3,xyz,1,coords(1,1,geom),1) 1165 coords(ixyz,atom,geom) = coords(ixyz,atom,geom) + delta 1166 call dfill(lscr,0.0d00,scr,1) 1167 call dfill(lbuf,0.0d00,buf,1) 1168 call dfill(lbuf,0.0d00,bufp,1) 1169 call intp_2e2c(basis,ishell,basis,jshell, 1170 & R,lscr,scr,lbuf,bufp) 1171* 1172 call dcopy(nat3,xyz,1,coords(1,1,geom),1) 1173 coords(ixyz,atom,geom) = coords(ixyz,atom,geom) - delta 1174 call dfill(lscr,0.0d00,scr,1) 1175 call dfill(lbuf,0.0d00,bufm,1) 1176 call intp_2e2c(basis,ishell,basis,jshell, 1177 & R,lscr,scr,lbuf,bufm) 1178* 1179 call dcopy(nbfsh,bufp,1,buf,1) 1180 call daxpy(nbfsh,-1.0d00,bufm,1,buf,1) 1181 call dscal(nbfsh,factor,buf,1) 1182 cnt = 1 1183 do i = ilo,ihi 1184 do j = jlo, jhi 1185 eri2fd(i,j,ixyz,atom) = buf(cnt) 1186 eri2fd(j,i,ixyz,atom) = buf(cnt) 1187 cnt = cnt + 1 1188 enddo 1189 enddo 1190 if (.not.ma_verify_allocator_stuff()) 1191 & stop ' ma broke 15' 1192 enddo 1193 enddo 1194 enddo 1195 enddo 1196c 1197 write(6,*)' fd: full list i,j' 1198 nzero1 = 0 1199 do atom = 1,nat 1200 do ixyz = 1,3 1201 do i = 1,nbf 1202 do j = 1,nbf 1203 if (abs(eri2fd(i,j,ixyz,atom)).le.thresh) 1204 & nzero1 = nzero1 + 1 1205* write(6,10000)i,j,ixyz,atom,eri2fd(i,j,ixyz,atom) 1206 enddo 1207 enddo 1208 enddo 1209 enddo 1210 write(6,*)' fd: non-zero list ' 1211 nzero2 = 0 1212 do atom = 1,nat 1213 do ixyz = 1,3 1214 do i = 1,nbf 1215 do j = 1,nbf 1216 if (abs(eri2fd(i,j,ixyz,atom)).gt.thresh) then 1217* write(6,10000)i,j,ixyz,atom,eri2fd(i,j,ixyz,atom) 1218 continue 1219 else 1220 nzero2 = nzero2 + 1 1221 endif 1222 enddo 1223 enddo 1224 enddo 1225 enddo 1226 write(6,*)' fd: num zeros :1:', nzero1 1227 write(6,*)' fd: num zeros :2:', nzero2, ':delta:', 1228 & (nzero2-nzero1) 1229 write(6,*)' fd: non-zero :1:', ((nbf*nbf*3*nat)-nzero1) 1230 write(6,*)' fd: non-zero :2:', ((nbf*nbf*3*nat)-nzero2) 1231 write(6,*)' fd: possible : :', (nbf*nbf*3*nat) 1232 nzero1 = 0 1233 nzero2 = 0 1234c 1235 call dcopy(nat3,xyz,1,coords(1,1,geom),1) 1236 do ishell = 1,nshell 1237 if (.not.bas_cn2bfr(basis,ishell,ilo,ihi)) 1238 & stop 'cn2bfr error i' 1239 do jshell = 1,ishell 1240 if (.not.bas_cn2bfr(basis,jshell,jlo,jhi)) 1241 & stop 'cn2bfr error j' 1242 call dfill(lscr,0.0d00,scr,1) 1243 call dfill(lbuf,0.0d00,buf,1) 1244* write(6,*)'rak20: ishell, jshell ',ishell,jshell 1245 call intpd_2e2c(basis,ishell,basis,jshell,R,lscr,scr, 1246 & lbuf,buf,watom) 1247* write(6,*)'watom 2e2c',watom 1248 cnt = 1 1249 if (.not.ma_verify_allocator_stuff()) stop ' ma broke 16' 1250 do atom=1,2 1251 if (watom(atom).gt.0) then 1252 do ixyz = 1,3 1253 do i = ilo, ihi 1254 do j = jlo, jhi 1255 eri2c(i,j,ixyz,watom(atom)) = buf(cnt) 1256 eri2c(j,i,ixyz,watom(atom)) = buf(cnt) 1257 cnt = cnt + 1 1258 enddo 1259 enddo 1260 enddo 1261 endif 1262 enddo 1263 if (.not.ma_verify_allocator_stuff()) stop ' ma broke 17' 1264 enddo 1265 enddo 1266 write(6,*)' c: full list i,j' 1267 nzero1 = 0 1268 do atom = 1,nat 1269 do ixyz = 1,3 1270 do i = 1,nbf 1271 do j = 1,nbf 1272 if (abs(eri2c(i,j,ixyz,atom)).le.thresh) 1273 & nzero1 = nzero1 + 1 1274* write(6,10001)i,j,ixyz,atom,eri2c(i,j,ixyz,atom) 1275 enddo 1276 enddo 1277 enddo 1278 enddo 1279 write(6,*)' c: non-zero list ' 1280 nzero2 = 0 1281 do atom = 1,nat 1282 do ixyz = 1,3 1283 do i = 1,nbf 1284 do j = 1,nbf 1285 if (abs(eri2c(i,j,ixyz,atom)).gt.thresh) then 1286* write(6,10001)i,j,ixyz,atom,eri2c(i,j,ixyz,atom) 1287 continue 1288 else 1289 nzero2 = nzero2 + 1 1290 endif 1291 enddo 1292 enddo 1293 enddo 1294 enddo 1295 write(6,*)' c: num zeros :1:', nzero1 1296 write(6,*)' c: num zeros :2:', nzero2, ':delta:', 1297 & (nzero2-nzero1) 1298 write(6,*)' c: non-zero :1:', ((nbf*nbf*3*nat)-nzero1) 1299 write(6,*)' c: non-zero :2:', ((nbf*nbf*3*nat)-nzero2) 1300 write(6,*)' c: possible : :', (nbf*nbf*3*nat) 1301c 130210000 format(1x,'eri2fd(',i3,',',i3,',',i2,',',i3,') = ',1pd20.10) 130310001 format(1x,' eri2c(',i3,',',i3,',',i2,',',i3,') = ',1pd20.10) 1304c 1305 call daxpy((nbf*nbf*3*nat),-1.0d00,eri2fd,1,eri2c,1) 1306 norm = ddot((nbf*nbf*3*nat),eri2c,1,eri2c,1) 1307 write(luout,*)' 2e2c difference norm:',ir,' ',norm 1308c 1309 pderiv_compute_2e2c = norm.lt.thresh 1310 enddo 1311 end 1312 logical function pderiv_compute_2e3c( 1313 & geom,basis,nbf,nat,lbuf,lscr, 1314 & buf,bufp,bufm,scr,eri3c,eri3fd,xyz) 1315 implicit none 1316#include "mafdecls.fh" 1317#include "nwc_const.fh" 1318#include "geomP.fh" 1319#include "bas.fh" 1320#include "stdio.fh" 1321c 1322 double precision ddot 1323 external ddot 1324c 1325 integer geom 1326 integer basis 1327 integer nbf 1328 integer nat 1329 integer lbuf 1330 integer lscr 1331 double precision buf(lbuf), bufp(lbuf), bufm(lbuf) 1332 double precision scr(lscr) 1333 double precision eri3c(nbf,nbf,nbf,3,nat) 1334 double precision eri3fd(nbf,nbf,nbf,3,nat) 1335 double precision xyz(3,nat) 1336c 1337 integer nat3, nintz 1338 integer atom, ixyz, IR 1339 double precision delta, factor, thresh, norm 1340 double precision R1(3), R2(3) 1341 integer nzero1, nzero2 1342 integer nshell, cnt, nbfsh, watom(4) 1343 integer ishell, ilo, ihi, nbfshi, i 1344 integer jshell, jlo, jhi, nbfshj, j 1345 integer kshell, klo, khi, nbfshk, k 1346c 1347 call dfill(3,0.0d00,R1,1) 1348 call dfill(3,0.0d00,R2,1) 1349 call dfill((3*nat),0.0d00,xyz,1) 1350 call dfill((nbf*nbf*nbf*3*nat),0.0d00,eri3c,1) 1351 call dfill((nbf*nbf*nbf*3*nat),0.0d00,eri3fd,1) 1352 call dfill(lbuf,0.0d00,buf,1) 1353 call dfill(lbuf,0.0d00,bufp,1) 1354 call dfill(lbuf,0.0d00,bufm,1) 1355 call dfill(lscr,0.0d00,scr,1) 1356 do IR = 1,3 1357 if (IR.eq.1) then 1358 call dfill(3,0.0d00,R1,1) 1359 call dfill(3,0.0d00,R2,1) 1360 elseif (IR.eq.2) then 1361 call dfill(3,0.0d00,R2,1) 1362 R1(1) = 1.0d00 1363 R1(2) = 2.0d00 1364 R1(3) = 3.0d00 1365 elseif (IR.eq.3) then 1366 R1(1) = 1.0d00 1367 R1(2) = 2.0d00 1368 R1(3) = 3.0d00 1369 R2(1) = 3.0d00 1370 R2(2) = 4.0d00 1371 R2(3) = 5.0d00 1372 else 1373 stop ' how did IR get here' 1374 endif 1375 if (.not.ma_verify_allocator_stuff()) stop ' ma broke 18' 1376 thresh = 1.0d-10 1377 delta = 0.0001d00 1378 factor = 1.0d00/(2.0d00*delta) 1379 nat3 = 3*nat 1380* store original coordintates 1381 call dcopy(nat3,coords(1,1,geom),1,xyz,1) 1382* 1383 if (.not.bas_numcont(basis,nshell)) stop 'bas_numcont' 1384 1385 do ishell = 1, nshell 1386 if (.not.bas_cn2bfr(basis,ishell,ilo,ihi)) 1387 & stop 'cn2bfr error i' 1388 nbfshi = ihi - ilo + 1 1389 do jshell = 1, nshell 1390 if (.not.bas_cn2bfr(basis,jshell,jlo,jhi)) 1391 & stop 'cn2bfr error j' 1392 nbfshj = jhi - jlo + 1 1393 do kshell = 1, jshell 1394 if (.not.bas_cn2bfr(basis,kshell,klo,khi)) 1395 & stop 'cn2bfr error k' 1396 nbfshk = khi - klo + 1 1397 nbfsh = nbfshi*nbfshj*nbfshk 1398 do atom = 1,nat 1399 do ixyz = 1,3 1400 call dcopy(nat3,xyz,1,coords(1,1,geom),1) 1401 coords(ixyz,atom,geom) = 1402 & coords(ixyz,atom,geom) + delta 1403 call dfill(lscr,0.0d00,scr,1) 1404 call dfill(lbuf,0.0d00,buf,1) 1405 call dfill(lbuf,0.0d00,bufp,1) 1406 call intp_2e3c(basis,ishell,basis,jshell,kshell, 1407 & R1,R2,lscr,scr,lbuf,bufp) 1408* 1409 call dcopy(nat3,xyz,1,coords(1,1,geom),1) 1410 coords(ixyz,atom,geom) = 1411 & coords(ixyz,atom,geom) - delta 1412 call dfill(lscr,0.0d00,scr,1) 1413 call dfill(lbuf,0.0d00,bufm,1) 1414 call intp_2e3c(basis,ishell,basis,jshell,kshell, 1415 & R1,R2,lscr,scr,lbuf,bufm) 1416* 1417 call dcopy(nbfsh,bufp,1,buf,1) 1418 call daxpy(nbfsh,-1.0d00,bufm,1,buf,1) 1419 call dscal(nbfsh,factor,buf,1) 1420 cnt = 1 1421 do i = ilo,ihi 1422 do j = jlo, jhi 1423 do k = klo, khi 1424 eri3fd(i,j,k,ixyz,atom) = buf(cnt) 1425 eri3fd(i,k,j,ixyz,atom) = buf(cnt) 1426 cnt = cnt + 1 1427 enddo 1428 enddo 1429 enddo 1430 if (.not.ma_verify_allocator_stuff()) 1431 & stop ' ma broke 19' 1432 enddo 1433 enddo 1434 enddo 1435 enddo 1436 enddo 1437c 1438 write(6,*)' fd: full list i,j' 1439 nzero1 = 0 1440 do atom = 1,nat 1441 do ixyz = 1,3 1442 do i = 1,nbf 1443 do j = 1,nbf 1444 do k = 1,nbf 1445 if (abs(eri3fd(i,j,k,ixyz,atom)).le.thresh) 1446 & nzero1 = nzero1 + 1 1447* write(6,10000)i,j,k,ixyz,atom,eri3fd(i,j,k,ixyz,atom) 1448 enddo 1449 enddo 1450 enddo 1451 enddo 1452 enddo 1453 write(6,*)' fd: non-zero list ' 1454 nzero2 = 0 1455 do atom = 1,nat 1456 do ixyz = 1,3 1457 do i = 1,nbf 1458 do j = 1,nbf 1459 do k = 1,nbf 1460 if (abs(eri3fd(i,j,k,ixyz,atom)).gt.thresh) then 1461* write(6,10000)i,j,k,ixyz,atom, 1462* & eri3fd(i,j,k,ixyz,atom) 1463 continue 1464 else 1465 nzero2 = nzero2 + 1 1466 endif 1467 enddo 1468 enddo 1469 enddo 1470 enddo 1471 enddo 1472 write(6,*)' fd: num zeros :1:', nzero1 1473 write(6,*)' fd: num zeros :2:', nzero2, ':delta:', 1474 & (nzero2-nzero1) 1475 write(6,*)' fd: non-zero :1:', ((nbf*nbf*nbf*3*nat)-nzero1) 1476 write(6,*)' fd: non-zero :2:', ((nbf*nbf*nbf*3*nat)-nzero2) 1477 write(6,*)' fd: possible : :', (nbf*nbf*nbf*3*nat) 1478 nzero1 = 0 1479 nzero2 = 0 1480c 1481 call dcopy(nat3,xyz,1,coords(1,1,geom),1) 1482 do ishell = 1,nshell 1483 if (.not.bas_cn2bfr(basis,ishell,ilo,ihi)) 1484 & stop 'cn2bfr error i' 1485 do jshell = 1,nshell 1486 if (.not.bas_cn2bfr(basis,jshell,jlo,jhi)) 1487 & stop 'cn2bfr error j' 1488 do kshell = 1, jshell 1489 if (.not.bas_cn2bfr(basis,kshell,klo,khi)) 1490 & stop 'cn2bfr error k' 1491 call dfill(lscr,0.0d00,scr,1) 1492 call dfill(lbuf,0.0d00,buf,1) 1493* write(6,*)'rak20: ishell, jshell, kshell ', 1494* & ishell,jshell,kshell 1495 call intpd_2e3c(basis,ishell,basis,jshell,kshell,R1,R2, 1496 & lscr,scr,lbuf,buf,watom) 1497* write(6,*)'watom 2e3c',watom 1498 if (.not.ma_verify_allocator_stuff()) stop ' ma broke 20' 1499 nintz = (ihi-ilo+1)*(jhi-jlo+1)*(khi-klo+1) 1500 do atom=1,4 1501 if (watom(atom).gt.0) then 1502 cnt = ((atom-1)*nintz*3) + 1 1503 do ixyz = 1,3 1504 do i = ilo, ihi 1505 do j = jlo, jhi 1506 do k = klo, khi 1507 eri3c(i,j,k,ixyz,watom(atom)) = buf(cnt) 1508 eri3c(i,k,j,ixyz,watom(atom)) = buf(cnt) 1509 cnt = cnt + 1 1510 enddo 1511 enddo 1512 enddo 1513 enddo 1514 endif 1515 enddo 1516 if (.not.ma_verify_allocator_stuff()) stop ' ma broke 21' 1517 enddo 1518 enddo 1519 enddo 1520 write(6,*)' c: full list i,j' 1521 nzero1 = 0 1522 do atom = 1,nat 1523 do ixyz = 1,3 1524 do i = 1,nbf 1525 do j = 1,nbf 1526 do k = 1,nbf 1527 if (abs(eri3c(i,j,k,ixyz,atom)).le.thresh) 1528 & nzero1 = nzero1 + 1 1529* write(6,10001)i,j,k,ixyz,atom,eri3c(i,j,k,ixyz,atom) 1530 enddo 1531 enddo 1532 enddo 1533 enddo 1534 enddo 1535 write(6,*)' c: non-zero list ' 1536 nzero2 = 0 1537 do atom = 1,nat 1538 do ixyz = 1,3 1539 do i = 1,nbf 1540 do j = 1,nbf 1541 do k = 1,nbf 1542 if (abs(eri3c(i,j,k,ixyz,atom)).gt.thresh) then 1543* write(6,10001)i,j,k,ixyz,atom, 1544* & eri3c(i,j,k,ixyz,atom) 1545 continue 1546 else 1547 nzero2 = nzero2 + 1 1548 endif 1549 enddo 1550 enddo 1551 enddo 1552 enddo 1553 enddo 1554 write(6,*)' c: num zeros :1:', nzero1 1555 write(6,*)' c: num zeros :2:', nzero2, ':delta:', 1556 & (nzero2-nzero1) 1557 write(6,*)' c: non-zero :1:', ((nbf*nbf*nbf*3*nat)-nzero1) 1558 write(6,*)' c: non-zero :2:', ((nbf*nbf*nbf*3*nat)-nzero2) 1559 write(6,*)' c: possible : :', (nbf*nbf*nbf*3*nat) 1560c 156110000 format(1x,'eri3fd(',i3,',',i3,',',i3,',',i2,',',i3, 1562 & ') = ',1pd20.10) 156310001 format(1x,' eri3c(',i3,',',i3,',',i3,',',i2,',',i3, 1564 & ') = ',1pd20.10) 1565c 1566 call daxpy((nbf*nbf*nbf*3*nat),-1.0d00,eri3fd,1,eri3c,1) 1567 norm = ddot((nbf*nbf*nbf*3*nat),eri3c,1,eri3c,1) 1568 write(luout,*)' 2e3c difference norm:',ir,' ',norm 1569c 1570 pderiv_compute_2e3c = norm.lt.thresh 1571 enddo 1572 end 1573 logical function pderiv_compute_3ov( 1574 & geom,basis,nbf,nat,lbuf,lscr, 1575 & buf,bufp,bufm,scr,ov3c,ov3fd,xyz) 1576 implicit none 1577#include "mafdecls.fh" 1578#include "nwc_const.fh" 1579#include "geomP.fh" 1580#include "bas.fh" 1581#include "stdio.fh" 1582c 1583 double precision ddot 1584 external ddot 1585c 1586 integer geom 1587 integer basis 1588 integer nbf 1589 integer nat 1590 integer lbuf 1591 integer lscr 1592 double precision buf(lbuf), bufp(lbuf), bufm(lbuf) 1593 double precision scr(lscr) 1594 double precision ov3c(nbf,nbf,nbf,3,nat) 1595 double precision ov3fd(nbf,nbf,nbf,3,nat) 1596 double precision xyz(3,nat) 1597c 1598 integer nat3, nintz 1599 integer atom, ixyz 1600 double precision delta, factor, thresh, norm 1601*rak: double precision zeta(50) 1602*rak: integer ztype, znp, zng, zsph 1603 integer nzero1, nzero2 1604 integer nshell, cnt, nbfsh, watom(3) 1605 integer ishell, ilo, ihi, nbfshi, i 1606 integer jshell, jlo, jhi, nbfshj, j 1607 integer kshell, klo, khi, nbfshk, k 1608c 1609 call dfill((3*nat),0.0d00,xyz,1) 1610 call dfill((nbf*nbf*nbf*3*nat),0.0d00,ov3c,1) 1611 call dfill((nbf*nbf*nbf*3*nat),0.0d00,ov3fd,1) 1612 call dfill(lbuf,0.0d00,buf,1) 1613 call dfill(lbuf,0.0d00,bufp,1) 1614 call dfill(lbuf,0.0d00,bufm,1) 1615 call dfill(lscr,0.0d00,scr,1) 1616 if (.not.ma_verify_allocator_stuff()) stop ' ma broke 18' 1617 thresh = 1.0d-10 1618 delta = 0.0001d00 1619 factor = 1.0d00/(2.0d00*delta) 1620 nat3 = 3*nat 1621* store original coordintates 1622 call dcopy(nat3,coords(1,1,geom),1,xyz,1) 1623* 1624 if (.not.bas_numcont(basis,nshell)) stop 'bas_numcont' 1625* 1626*rak: if (.not.bas_continfo(basis,1,ztype,znp,zng,zsph)) stop 'e1' 1627*rak: if (znp.gt.50) stop ' zeta dimension too small ' 1628*rak: call dfill(znp,0.0d00,zeta,1) 1629*rak: if (.not.bas_get_exponent(basis,1,zeta)) stop 'bas_get_e failed' 1630*rak: zeta(1) = 1.0d-50 1631*rak: if (.not.bas_set_exponent(basis,1,zeta,znp)) 1632*rak: & stop 'bas_set_e failed' 1633* 1634 do ishell = 1, nshell 1635 if (.not.bas_cn2bfr(basis,ishell,ilo,ihi)) 1636 & stop 'cn2bfr error i' 1637 nbfshi = ihi - ilo + 1 1638 do jshell = 1, nshell 1639 if (.not.bas_cn2bfr(basis,jshell,jlo,jhi)) 1640 & stop 'cn2bfr error j' 1641 nbfshj = jhi - jlo + 1 1642 do kshell = 1, nshell 1643 if (.not.bas_cn2bfr(basis,kshell,klo,khi)) 1644 & stop 'cn2bfr error k' 1645 nbfshk = khi - klo + 1 1646 nbfsh = nbfshi*nbfshj*nbfshk 1647 do atom = 1,nat 1648 do ixyz = 1,3 1649 call dcopy(nat3,xyz,1,coords(1,1,geom),1) 1650 coords(ixyz,atom,geom) = 1651 & coords(ixyz,atom,geom) + delta 1652 call dfill(lscr,0.0d00,scr,1) 1653 call dfill(lbuf,0.0d00,buf,1) 1654 call dfill(lbuf,0.0d00,bufp,1) 1655 1656 call int_1e3ov( 1657 & basis,ishell,basis,jshell,basis,kshell, 1658 & lscr,scr,lbuf,bufp) 1659* 1660 call dcopy(nat3,xyz,1,coords(1,1,geom),1) 1661 coords(ixyz,atom,geom) = 1662 & coords(ixyz,atom,geom) - delta 1663 call dfill(lscr,0.0d00,scr,1) 1664 call dfill(lbuf,0.0d00,bufm,1) 1665 call int_1e3ov( 1666 & basis,ishell,basis,jshell,basis,kshell, 1667 & lscr,scr,lbuf,bufm) 1668* 1669 call dcopy(nbfsh,bufp,1,buf,1) 1670 call daxpy(nbfsh,-1.0d00,bufm,1,buf,1) 1671 call dscal(nbfsh,factor,buf,1) 1672 cnt = 1 1673 do i = ilo,ihi 1674 do j = jlo, jhi 1675 do k = klo, khi 1676 ov3fd(i,j,k,ixyz,atom) = buf(cnt) 1677 cnt = cnt + 1 1678 enddo 1679 enddo 1680 enddo 1681 if (.not.ma_verify_allocator_stuff()) 1682 & stop ' ma broke 19' 1683 enddo 1684 enddo 1685 enddo 1686 enddo 1687 enddo 1688c 1689 write(6,*)' fd: full list i,j' 1690 nzero1 = 0 1691 do atom = 1,nat 1692 do ixyz = 1,3 1693 do i = 1,nbf 1694 do j = 1,nbf 1695 do k = 1,nbf 1696 if (abs(ov3fd(i,j,k,ixyz,atom)).le.thresh) 1697 & nzero1 = nzero1 + 1 1698* write(6,10000)i,j,k,ixyz,atom,ov3fd(i,j,k,ixyz,atom) 1699 enddo 1700 enddo 1701 enddo 1702 enddo 1703 enddo 1704 write(6,*)' fd: non-zero list ' 1705 nzero2 = 0 1706 do atom = 1,nat 1707 do ixyz = 1,3 1708 do i = 1,nbf 1709 do j = 1,nbf 1710 do k = 1,nbf 1711 if (abs(ov3fd(i,j,k,ixyz,atom)).gt.thresh) then 1712* write(6,10000)i,j,k,ixyz,atom, 1713* & ov3fd(i,j,k,ixyz,atom) 1714 continue 1715 else 1716 nzero2 = nzero2 + 1 1717 endif 1718 enddo 1719 enddo 1720 enddo 1721 enddo 1722 enddo 1723 write(6,*)' fd: num zeros :1:', nzero1 1724 write(6,*)' fd: num zeros :2:', nzero2, ':delta:', 1725 & (nzero2-nzero1) 1726 write(6,*)' fd: non-zero :1:', ((nbf*nbf*nbf*3*nat)-nzero1) 1727 write(6,*)' fd: non-zero :2:', ((nbf*nbf*nbf*3*nat)-nzero2) 1728 write(6,*)' fd: possible : :', (nbf*nbf*nbf*3*nat) 1729 nzero1 = 0 1730 nzero2 = 0 1731c 1732 call dcopy(nat3,xyz,1,coords(1,1,geom),1) 1733 do ishell = 1,nshell 1734 if (.not.bas_cn2bfr(basis,ishell,ilo,ihi)) 1735 & stop 'cn2bfr error i' 1736 do jshell = 1,nshell 1737 if (.not.bas_cn2bfr(basis,jshell,jlo,jhi)) 1738 & stop 'cn2bfr error j' 1739 do kshell = 1, jshell 1740 if (.not.bas_cn2bfr(basis,kshell,klo,khi)) 1741 & stop 'cn2bfr error k' 1742 call dfill(lscr,0.0d00,scr,1) 1743 call dfill(lbuf,0.0d00,buf,1) 1744* write(6,*)'rak20: ishell, jshell, kshell ', 1745* & ishell,jshell,kshell 1746 call intd_1e3ov( 1747 & basis,ishell,basis,jshell,basis,kshell, 1748 & lscr,scr,lbuf,buf,watom) 1749* write(6,*)'watom 3ov',watom 1750 if (.not.ma_verify_allocator_stuff()) stop ' ma broke 20' 1751 nintz = (ihi-ilo+1)*(jhi-jlo+1)*(khi-klo+1) 1752 do atom=1,3 1753 if (watom(atom).gt.0) then 1754 cnt = ((atom-1)*nintz*3) + 1 1755 do ixyz = 1,3 1756 do i = ilo, ihi 1757 do j = jlo, jhi 1758 do k = klo, khi 1759 ov3c(i,j,k,ixyz,watom(atom)) = buf(cnt) 1760 ov3c(i,k,j,ixyz,watom(atom)) = buf(cnt) 1761 cnt = cnt + 1 1762 enddo 1763 enddo 1764 enddo 1765 enddo 1766 endif 1767 enddo 1768 if (.not.ma_verify_allocator_stuff()) stop ' ma broke 21' 1769 enddo 1770 enddo 1771 enddo 1772 write(6,*)' c: full list i,j' 1773 nzero1 = 0 1774 do atom = 1,nat 1775 do ixyz = 1,3 1776 do i = 1,nbf 1777 do j = 1,nbf 1778 do k = 1,nbf 1779 if (abs(ov3c(i,j,k,ixyz,atom)).le.thresh) 1780 & nzero1 = nzero1 + 1 1781* write(6,10001)i,j,k,ixyz,atom,ov3c(i,j,k,ixyz,atom) 1782 enddo 1783 enddo 1784 enddo 1785 enddo 1786 enddo 1787 write(6,*)' c: non-zero list ' 1788 nzero2 = 0 1789 do atom = 1,nat 1790 do ixyz = 1,3 1791 do i = 1,nbf 1792 do j = 1,nbf 1793 do k = 1,nbf 1794 if (abs(ov3c(i,j,k,ixyz,atom)).gt.thresh) then 1795* write(6,10001)i,j,k,ixyz,atom, 1796* & ov3c(i,j,k,ixyz,atom) 1797 continue 1798 else 1799 nzero2 = nzero2 + 1 1800 endif 1801 enddo 1802 enddo 1803 enddo 1804 enddo 1805 enddo 1806 write(6,*)' c: num zeros :1:', nzero1 1807 write(6,*)' c: num zeros :2:', nzero2, ':delta:', 1808 & (nzero2-nzero1) 1809 write(6,*)' c: non-zero :1:', ((nbf*nbf*nbf*3*nat)-nzero1) 1810 write(6,*)' c: non-zero :2:', ((nbf*nbf*nbf*3*nat)-nzero2) 1811 write(6,*)' c: possible : :', (nbf*nbf*nbf*3*nat) 1812c 181310000 format(1x,'ov3fd(',i3,',',i3,',',i3,',',i2,',',i3, 1814 & ') = ',1pd20.10) 181510001 format(1x,' ov3c(',i3,',',i3,',',i3,',',i2,',',i3, 1816 & ') = ',1pd20.10) 1817c 1818 call daxpy((nbf*nbf*nbf*3*nat),-1.0d00,ov3fd,1,ov3c,1) 1819 norm = ddot((nbf*nbf*nbf*3*nat),ov3c,1,ov3c,1) 1820 write(luout,*)' 3ov difference norm: ',norm 1821c 1822 pderiv_compute_3ov = norm.lt.thresh 1823 end 1824 logical function pderiv_compute_mpole( 1825 & geom,basis,nbf,nat,lbuf,lscr, 1826 & buf,bufp,bufm,scr,mp3c,mp3fd,xyz) 1827 implicit none 1828#include "mafdecls.fh" 1829#include "nwc_const.fh" 1830#include "geomP.fh" 1831#include "bas.fh" 1832#include "stdio.fh" 1833c 1834 double precision ddot 1835 external ddot 1836c 1837 integer geom 1838 integer basis 1839 integer nbf 1840 integer nat 1841 integer lbuf 1842 integer lscr 1843 integer lval_max 1844 parameter (lval_max = 5) 1845 integer slmax 1846* parameter (slmax = 1) ! for lval_max = 0 1847* parameter (slmax = 4) ! for lval_max = 1 1848* parameter (slmax = 35) ! ! for lval_max = 4 1849 parameter (slmax = 56) ! 1 +3+6+10+15+21 = 35+21 = 56 1850 double precision buf(lbuf), bufp(lbuf), bufm(lbuf) 1851 double precision scr(lscr) 1852 double precision mp3c(slmax,nbf,nbf,3,(nat+1)) ! add one for the multipole center 1853 double precision mp3fd(slmax,nbf,nbf,3,(nat+1)) ! ditto 1854 double precision xyz(3,nat) 1855c 1856 integer nat3, nintz 1857 integer atom, ixyz, IR, lval, nbflval 1858 double precision delta, factor, thresh, norm 1859 double precision R1(3), l_center(3), l_orgcenter(3) 1860 integer nzero1, nzero2 1861 integer nshell, cnt, nbfsh, watom(3), nshell_bas 1862 integer ishell, ilo, ihi, nbfshi, i 1863 integer jshell, jlo, jhi, nbfshj, j 1864 integer k, klo, khi, num_ret 1865 integer h_diff, k_diff 1866 integer imp_cent 1867c 1868 integer pole_index 1869 pole_index(lval) = lval*((lval+3)*(lval+3)+2)/6 + 1 1870c 1871 call dfill(3,0.0d00,R1,1) 1872 call dfill((3*nat),0.0d00,xyz,1) 1873 call dfill((slmax*nbf*nbf*3*(nat+1)),0.0d00,mp3c,1) 1874 call dfill((slmax*nbf*nbf*3*(nat+1)),0.0d00,mp3fd,1) 1875 call dfill(lbuf,0.0d00,buf,1) 1876 call dfill(lbuf,0.0d00,bufp,1) 1877 call dfill(lbuf,0.0d00,bufm,1) 1878 call dfill(lscr,0.0d00,scr,1) 1879 do imp_cent = 1, 2 1880 if (imp_cent.eq.1) then 1881 call dfill(3,0.0d00,l_center,1) ! start with multipole center at origin 1882 else if (imp_cent.eq.2) then 1883 call dcopy(3,coords(1,1,geom),1,l_center,1) 1884 endif 1885 call dcopy(3,l_center,1,l_orgcenter,1) 1886 1887 do IR = 1, 3 1888 if (IR.eq.1) then 1889 call dfill(3,0.0d00,R1,1) 1890 elseif (IR.eq.2) then 1891 R1(1) = 1.0d00 1892 R1(2) = 2.0d00 1893 R1(3) = 3.0d00 1894 elseif (IR.eq.3) then 1895 R1(1) = 1.0d00 1896 R1(2) = 1.0d00 1897 R1(3) = 1.0d00 1898 else 1899 stop ' how did IR get here' 1900 endif 1901 if (.not.ma_verify_allocator_stuff()) stop ' ma broke 18' 1902 thresh = 1.0d-10 1903 delta = 0.0001d00 1904 factor = 1.0d00/(2.0d00*delta) 1905*rak: write(6,*)' factor = ',factor 1906 nat3 = 3*nat 1907* store original coordintates 1908 call dcopy(nat3,coords(1,1,geom),1,xyz,1) 1909*rak: write(6,*)' original coords ' 1910*rak: call output(coords(1,1,geom),1,3,1,nat,3,nat,1) 1911*rak: write(6,*)' l_center',l_center 1912*rak: write(6,*)' l_orgcenter',l_orgcenter 1913* 1914 if (.not.bas_numcont(basis,nshell_bas)) stop 'bas_numcont' 1915 nshell = nshell_bas 1916 nshell = 1 ! tmp change 1917*rak: write(6,*)' nshell = ',nshell 1918 1919 do ishell = 1, nshell 1920 if (.not.bas_cn2bfr(basis,ishell,ilo,ihi)) 1921 & stop 'cn2bfr error i' 1922 nbfshi = ihi - ilo + 1 1923 do jshell = 1, nshell 1924 if (.not.bas_cn2bfr(basis,jshell,jlo,jhi)) 1925 & stop 'cn2bfr error j' 1926 nbfshj = jhi - jlo + 1 1927 do lval = 0,lval_max 1928*rak: write(6,*)'................................1', 1929*rak: & ishell,jshell,lval 1930 nbflval = (lval+1)*(lval+2)/2 1931 nbfsh = nbfshj*nbfshi*nbflval 1932 klo = pole_index((lval-1)) + 1 1933 klo = max(klo,1) ! minumum value of klo is 1 1934 khi = pole_index(lval) 1935*rak: write(6,*)' klo, khi',klo,khi 1936 do atom = 1,nat 1937 do ixyz = 1,3 1938 call dcopy(nat3,xyz,1,coords(1,1,geom),1) 1939 coords(ixyz,atom,geom) = 1940 & coords(ixyz,atom,geom) + delta 1941 call dfill(lscr,0.0d00,scr,1) 1942 call dfill(lbuf,0.0d00,buf,1) 1943 call dfill(lbuf,0.0d00,bufp,1) 1944*rak: write(6,*)' coords intp_mpolel(1)',atom,ixyz 1945*rak: call output(coords(1,1,geom),1,3,1,nat,3,nat,1) 1946*rak: write(6,*)' l_center',l_center 1947*rak: write(6,*)' l_orgcenter',l_orgcenter,' -*' 1948 call intp_mpolel(basis,ishell,basis,jshell, 1949 & R1,lval,l_orgcenter,lscr,scr,lbuf,bufp,num_ret) 1950* call intp_1eov(basis,ishell,basis,jshell, 1951* & R1,lscr,scr,lbuf,buf) 1952* cnt = 0 1953* write(6,*) 1954* & ' i j cnt mpole ovl diff + atom ixyz' 1955* do i = ilo,ihi 1956* do j = jlo,jhi 1957* cnt = cnt + 1 1958* norm = abs(bufp(cnt)-buf(cnt)) 1959* write(6,*)i,j,cnt,bufp(cnt),buf(cnt),norm, 1960* & atom, ixyz 1961* enddo 1962* enddo 1963* norm = 0.0d00 1964* call dfill(lbuf,0.0d00,buf,1) 1965* 1966 call dcopy(nat3,xyz,1,coords(1,1,geom),1) 1967 coords(ixyz,atom,geom) = 1968 & coords(ixyz,atom,geom) - delta 1969 call dfill(lscr,0.0d00,scr,1) 1970 call dfill(lbuf,0.0d00,bufm,1) 1971*rak: write(6,*)' coords intp_mpolel(2)',atom,ixyz 1972*rak: call output(coords(1,1,geom),1,3,1,nat,3,nat,1) 1973*rak: write(6,*)' l_center',l_center 1974*rak: write(6,*)' l_orgcenter',l_orgcenter,' -*' 1975 call intp_mpolel(basis,ishell,basis,jshell, 1976 & R1,lval,l_orgcenter,lscr,scr,lbuf,bufm,num_ret) 1977*rak: call intp_1eov(basis,ishell,basis,jshell, 1978*rak: & R1,lscr,scr,lbuf,buf) 1979*rak: cnt = 0 1980*rak: write(6,*) 1981*rak: & ' i j cnt mpole ovl diff - atom ixyz' 1982*rak: do i = ilo,ihi 1983*rak: do j = jlo,jhi 1984*rak: cnt = cnt + 1 1985*rak: norm = abs(bufp(cnt)-buf(cnt)) 1986*rak: write(6,*)i,j,cnt,bufp(cnt),buf(cnt),norm, 1987*rak: & atom,ixyz 1988*rak: enddo 1989*rak: enddo 1990*rak: norm = 0.0d00 1991*rak: call dfill(lbuf,0.0d00,buf,1) 1992* 1993*rak: call prak3mat(bufp,nbfsh,nbfshi,nbfshj,nbflval, 1994*rak: & ilo,ihi,jlo,jhi,klo,khi, 1995*rak: & 'bufp in atom move',atom,ixyz) 1996*rak: call prak3mat(bufm,nbfsh,nbfshi,nbfshj,nbflval, 1997*rak: & ilo,ihi,jlo,jhi,klo,khi, 1998*rak: & 'bufm in atom move',atom,ixyz) 1999 call dcopy(nbfsh,bufp,1,buf,1) 2000 call daxpy(nbfsh,-1.0d00,bufm,1,buf,1) 2001*rak: call prak3mat(buf,nbfsh,nbfshi,nbfshj,nbflval, 2002*rak: & ilo,ihi,jlo,jhi,klo,khi, 2003*rak: & 'buf before scale in atom move',atom,ixyz) 2004 call dscal(nbfsh,factor,buf,1) 2005*rak: call prak3mat(buf,nbfsh,nbfshi,nbfshj,nbflval, 2006*rak: & ilo,ihi,jlo,jhi,klo,khi, 2007*rak: & 'buf in atom move',atom,ixyz) 2008 cnt = 1 2009 do i = ilo,ihi 2010 do j = jlo, jhi 2011 do k = klo, khi 2012 mp3fd(k,j,i,ixyz,atom) = buf(cnt) 2013*rak: write(70,*) 2014*rak: & 'fd<',k,j,i,ixyz,atom,'>=',buf(cnt) 2015 cnt = cnt + 1 2016 enddo 2017 enddo 2018 enddo 2019*rak: write(6,*)' nbfsh ',nbfsh,' cnt ',cnt 2020 if (.not.ma_verify_allocator_stuff()) 2021 & stop ' ma broke 19' 2022 enddo 2023 enddo 2024* now do multipole center move 2025 call dcopy(nat3,xyz,1,coords(1,1,geom),1) ! restore original atom coords 2026 do ixyz = 1,3 2027 call dcopy(3,l_orgcenter,1,l_center,1) 2028 l_center(ixyz) = l_center(ixyz) + delta 2029 call dfill(lscr,0.0d00,scr,1) 2030 call dfill(lbuf,0.0d00,buf,1) 2031 call dfill(lbuf,0.0d00,bufp,1) 2032*rak: write(6,*)' coords intp_mpolel(3)',atom,ixyz 2033*rak: call output(coords(1,1,geom),1,3,1,nat,3,nat,1) 2034*rak: write(6,*)' l_center',l_center,' -*' 2035*rak: write(6,*)' l_orgcenter',l_orgcenter 2036 call intp_mpolel(basis,ishell,basis,jshell, 2037 & R1,lval,l_center,lscr,scr,lbuf,bufp,num_ret) 2038* 2039 call dcopy(3,l_orgcenter,1,l_center,1) 2040 l_center(ixyz) = l_center(ixyz) - delta 2041 call dfill(lscr,0.0d00,scr,1) 2042 call dfill(lbuf,0.0d00,bufm,1) 2043*rak: write(6,*)' coords intp_mpolel(4)',atom,ixyz 2044*rak: call output(coords(1,1,geom),1,3,1,nat,3,nat,1) 2045*rak: write(6,*)' l_center',l_center,' -*' 2046*rak: write(6,*)' l_orgcenter',l_orgcenter 2047 call intp_mpolel(basis,ishell,basis,jshell, 2048 & R1,lval,l_center,lscr,scr,lbuf,bufm,num_ret) 2049* 2050*rak: call prak3mat(bufp,nbfsh,nbfshi,nbfshj,nbflval, 2051*rak: & ilo,ihi,jlo,jhi,klo,khi, 2052*rak: & 'bufp',(nat+1),ixyz) 2053*rak: call prak3mat(bufm,nbfsh,nbfshi,nbfshj,nbflval, 2054*rak: & ilo,ihi,jlo,jhi,klo,khi, 2055*rak: & 'bufm',(nat+1),ixyz) 2056 call dcopy(nbfsh,bufp,1,buf,1) 2057 call daxpy(nbfsh,-1.0d00,bufm,1,buf,1) 2058*rak: call prak3mat(buf,nbfsh,nbfshi,nbfshj,nbflval, 2059*rak: & ilo,ihi,jlo,jhi,klo,khi, 2060*rak: & 'buf b4 scale in mp center move',(nat+1),ixyz) 2061 call dscal(nbfsh,factor,buf,1) 2062*rak: call prak3mat(buf,nbfsh,nbfshi,nbfshj,nbflval, 2063*rak: & ilo,ihi,jlo,jhi,klo,khi, 2064*rak: & 'buf in mp center move',(nat+1),ixyz) 2065 cnt = 1 2066 do i = ilo,ihi 2067 do j = jlo, jhi 2068 do k = klo, khi 2069 mp3fd(k,j,i,ixyz,nat+1) = buf(cnt) 2070*rak: write(70,*) 2071*rak: & 'fd<',k,j,i,ixyz,(nat+1),'>=',buf(cnt) 2072 cnt = cnt + 1 2073 enddo 2074 enddo 2075 enddo 2076 if (.not.ma_verify_allocator_stuff()) 2077 & stop ' ma broke 19' 2078 2079 enddo 2080 enddo 2081 enddo 2082 enddo 2083c 2084 write(6,*)' fd: full list i,j' 2085 nzero1 = 0 2086 do atom = 1,nat+1 2087 do ixyz = 1,3 2088 do i = 1,nbf 2089 do j = 1,nbf 2090 do k = 1,slmax 2091 if (abs(mp3fd(k,j,i,ixyz,atom)).le.thresh) 2092 & nzero1 = nzero1 + 1 2093* write(6,10000)k,j,i,ixyz,atom,mp3fd(k,j,i,ixyz,atom) 2094 enddo 2095 enddo 2096 enddo 2097 enddo 2098 enddo 2099 write(6,*)' fd: non-zero list ' 2100 nzero2 = 0 2101 do atom = 1,nat+1 2102 do ixyz = 1,3 2103 do i = 1,nbf 2104 do j = 1,nbf 2105 do k = 1,slmax 2106 if (abs(mp3fd(k,j,i,ixyz,atom)).gt.thresh) then 2107* write(6,10000)k,j,i,ixyz,atom, 2108* & mp3fd(k,j,i,ixyz,atom) 2109 continue 2110 else 2111 nzero2 = nzero2 + 1 2112 endif 2113 enddo 2114 enddo 2115 enddo 2116 enddo 2117 enddo 2118 write(6,*)' fd: num zeros :1:', nzero1 2119 write(6,*)' fd: num zeros :2:', nzero2, ':delta:', 2120 & (nzero2-nzero1) 2121 write(6,*)' fd: non-zero :1:', 2122 & ((nbf*nbf*slmax*3*(nat+1))-nzero1) 2123 write(6,*)' fd: non-zero :2:', 2124 & ((nbf*nbf*slmax*3*(nat+1))-nzero2) 2125 write(6,*)' fd: possible : :', 2126 & (nbf*nbf*slmax*3*(nat+1)) 2127 nzero1 = 0 2128 nzero2 = 0 2129c 2130 call dcopy(nat3,xyz,1,coords(1,1,geom),1) 2131 do ishell = 1,nshell 2132 if (.not.bas_cn2bfr(basis,ishell,ilo,ihi)) 2133 & stop 'cn2bfr error i' 2134 do jshell = 1,nshell 2135 if (.not.bas_cn2bfr(basis,jshell,jlo,jhi)) 2136 & stop 'cn2bfr error j' 2137 do lval = 0, lval_max 2138*rak: write(6,*)'................................2', 2139*rak: & ishell,jshell,lval 2140 klo = pole_index((lval-1)) + 1 2141 klo = max(klo,1) ! minumum value of klo is 1 2142 khi = pole_index(lval) 2143*rak: write(6,*)' klo, khi',klo,khi 2144c 2145 call dfill(lscr,0.0d00,scr,1) 2146 call dfill(lbuf,0.0d00,buf,1) 2147 call dfill(lbuf,0.0d00,bufp,1) 2148*rak: write(6,*)'rak20: ishell, jshell,lval ', 2149*rak: & ishell,jshell,lval 2150*rak: write(6,*)' lscr ', lscr 2151*rak: call intpd_1eov(basis,ishell,basis,jshell,R1,lscr,scr, 2152*rak: & lbuf,bufp,watom) 2153*rak: write(6,*)' coords intpd_mpolel' 2154*rak: call output(coords(1,1,geom),1,3,1,nat,3,nat,1) 2155*rak: write(6,*)' l_center',l_center 2156*rak: write(6,*)' l_orgcenter',l_orgcenter,' -*' 2157 call intpd_mpolel(basis,ishell,basis,jshell,R1, 2158 & lval,l_orgcenter, 2159 & lscr,scr,lbuf,buf,num_ret,watom) 2160*rak: cnt = 0 2161*rak: write(6,*) ' atom ixyz i j cnt mpole ovl diff ' 2162*rak: do atom = 1,nat+1 2163*rak: do ixyz = 1,3 2164*rak: do i = ilo,ihi 2165*rak: do j = jlo,jhi 2166*rak: cnt = cnt + 1 2167*rak: norm = abs(bufp(cnt)-buf(cnt)) 2168*rak: write(6,*) 2169*rak: & atom,ixyz,i,j,cnt,buf(cnt),bufp(cnt),norm 2170*rak: enddo 2171*rak: enddo 2172*rak: enddo 2173*rak: enddo 2174*rak: call dfill(lbuf,0.0d00,bufp,1) 2175*rak: write(6,*)'watom mpole',watom 2176 if (.not.ma_verify_allocator_stuff()) stop ' ma broke 20' 2177 nintz = (ihi-ilo+1)*(jhi-jlo+1)*((lval+1)*(lval+2)/2) 2178*rak: write(6,*)'nintz ',nintz 2179*rak: write(6,*)'ilo/hi',ilo,ihi 2180*rak: write(6,*)'jlo/hi',jlo,jhi 2181*rak: write(6,*)'klo/hi',klo,khi 2182 do atom=1,3 2183*rak: write(6,*)' atom,watom(atom)',atom,watom(atom) 2184 if (watom(atom).gt.0) then 2185 cnt = ((atom-1)*nintz*3) + 1 2186 do ixyz = 1,3 2187*rak: k_diff = cnt + (ixyz-1)*nintz 2188*rak: call prak3mat(buf(k_diff),nbfsh, 2189*rak: & nbfshi,nbfshj,nbflval, 2190*rak: & ilo,ihi,jlo,jhi,klo,khi, 2191*rak: & 'buf in derivs call(1)',watom(atom),ixyz) 2192 do i = ilo, ihi 2193 do j = jlo, jhi 2194 do k = klo, khi 2195 mp3c(k,j,i,ixyz,watom(atom)) = buf(cnt) 2196*rak: write(71,*) 2197*rak: & ' c<',k,j,i,ixyz,watom(atom),'>=',buf(cnt) 2198 cnt = cnt + 1 2199 enddo 2200 enddo 2201 enddo 2202 enddo 2203 endif 2204 if (watom(atom).eq.-3) then 2205 cnt = ((atom-1)*nintz*3) + 1 2206 do ixyz = 1,3 2207*rak: k_diff = cnt + (ixyz-1)*nintz 2208*rak: call prak3mat(buf(k_diff),nbfsh, 2209*rak: & nbfshi,nbfshj,nbflval, 2210*rak: & ilo,ihi,jlo,jhi,klo,khi, 2211*rak: & 'buf in derivs call(2)',watom(atom),ixyz) 2212 do i = ilo, ihi 2213 do j = jlo, jhi 2214 do k = klo, khi 2215 mp3c(k,j,i,ixyz,nat+1) = buf(cnt) 2216*rak: write(71,*) 2217*rak: & ' c<',k,j,i,ixyz,(nat+1),'>=',buf(cnt) 2218 cnt = cnt + 1 2219 enddo 2220 enddo 2221 enddo 2222 enddo 2223 endif 2224 enddo 2225 if (.not.ma_verify_allocator_stuff()) stop ' ma broke 21' 2226 enddo 2227 enddo 2228 enddo 2229 write(6,*)' c: full list i,j' 2230 nzero1 = 0 2231 do atom = 1,nat+1 2232 do ixyz = 1,3 2233 do i = 1,nbf 2234 do j = 1,nbf 2235 do k = 1,slmax 2236 if (abs(mp3c(k,j,i,ixyz,atom)).le.thresh) 2237 & nzero1 = nzero1 + 1 2238* write(6,10001)k,j,i,ixyz,atom,mp3c(k,j,i,ixyz,atom) 2239 enddo 2240 enddo 2241 enddo 2242 enddo 2243 enddo 2244 write(6,*)' c: non-zero list ' 2245 nzero2 = 0 2246 do atom = 1,nat+1 2247 do ixyz = 1,3 2248 do i = 1,nbf 2249 do j = 1,nbf 2250 do k = 1,slmax 2251 if (abs(mp3c(k,j,i,ixyz,atom)).gt.thresh) then 2252* write(6,10001)k,j,i,ixyz,atom, 2253* & mp3c(k,j,i,ixyz,atom) 2254 continue 2255 else 2256 nzero2 = nzero2 + 1 2257 endif 2258 enddo 2259 enddo 2260 enddo 2261 enddo 2262 enddo 2263 write(6,*)' c: num zeros :1:', nzero1 2264 write(6,*)' c: num zeros :2:', nzero2, ':delta:', 2265 & (nzero2-nzero1) 2266 write(6,*)' c: non-zero :1:', 2267 & ((nbf*nbf*slmax*3*(nat+1))-nzero1) 2268 write(6,*)' c: non-zero :2:', 2269 & ((nbf*nbf*slmax*3*(nat+1))-nzero2) 2270 write(6,*)' c: possible : :', 2271 & (nbf*nbf*slmax*3*(nat+1)) 2272c 227310000 format(1x,'mp3fd(',i3,',',i3,',',i3,',',i2,',',i3, 2274 & ') = ',1pd20.10) 227510001 format(1x,' mp3c(',i3,',',i3,',',i3,',',i2,',',i3, 2276 & ') = ',1pd20.10) 2277c 2278 if (.not.ma_alloc_get(mt_dbl, 2279 & (nbf*nbf*slmax*3*(nat+1)), 2280 & 'diff buffer', 2281 & h_diff, 2282 & k_diff)) stop 'mpole: ma alloc failed' 2283 call dcopy ((nbf*nbf*slmax*3*(nat+1)),mp3c,1,dbl_mb(k_diff),1) 2284 call daxpy((nbf*nbf*slmax*3*(nat+1)), 2285 & -1.0d00,mp3fd,1,dbl_mb(k_diff),1) 2286 norm = ddot((nbf*nbf*slmax*3*(nat+1)), 2287 & dbl_mb(k_diff),1,dbl_mb(k_diff),1) 2288 write(luout,*)' mpole difference norm:ir=',ir, 2289 & ':imp_cent= ',imp_cent,': ',norm 2290c 2291 pderiv_compute_mpole = norm.lt.thresh 2292 if (.not.pderiv_compute_mpole) then 2293 do atom = 1,(nat+1) 2294 do ixyz = 1,3 2295 nintz = (atom-1)*3*slmax*nbf*nbf 2296 nintz = nintz + (ixyz-1)*slmax*nbf*nbf 2297 call prak3mat(dbl_mb(k_diff + nintz), 2298 & (nbf*nbf*slmax), 2299 & nbf,nbf,slmax, 2300 & 1,nbf,1,nbf,1,slmax, 2301 & 'diff matrix',atom,ixyz) 2302 enddo 2303 enddo 2304 endif 2305 if (.not.ma_free_heap(h_diff)) stop 'mpole: ma free failed ' 2306 enddo 2307 enddo 2308 end 2309 subroutine prak3mat(buf,nbf,nbfi,nbfj,nbfk, 2310 & ilo,ihi,jlo,jhi,klo,khi,msg,atom,ixyz) 2311 implicit none 2312 integer atom 2313 integer ixyz 2314 integer nbf, nbfi, nbfj, nbfk 2315 integer ilo,ihi,jlo,jhi,klo,khi 2316 character*(*) msg 2317 double precision buf(nbfk,nbfj,nbfi) 2318c 2319 integer i, j, k 2320c 2321 write(6,*)' ' 2322 write(6,*)' prak3mat<',msg,'> atom=',atom,' ixyz=',ixyz 2323c 2324 if (nbf.ne.(nbfk*nbfj*nbfi)) then 2325 write(6,*)' nbf error ' 2326 write(6,*)' nbf = ',nbf 2327 write(6,*)' nbfi = ',nbfi 2328 write(6,*)' nbfj = ',nbfj 2329 write(6,*)' nbfk = ',nbfk 2330 endif 2331 do i = ilo,ihi 2332 do j = jlo,jhi 2333 do k = klo,khi 2334 if (abs(buf(k,j,i)).gt.1.0d-07) then 2335 write(6,10000)k,j,i,buf(k,j,i) 2336 endif 2337 enddo 2338 enddo 2339 enddo 2340 write(6,*)' ' 234110000 format(1x,'buf(',i3,',',i3,',',i3,')=',1pd20.10) 2342 end 2343 logical function pderiv_compute_2e3ct( 2344 & geom,basis,nbf,nat,lbuf,lscr, 2345 & buf,bufp,bufm,scr,eri3,eri3t,xyz) 2346 implicit none 2347#include "mafdecls.fh" 2348#include "nwc_const.fh" 2349#include "geomP.fh" 2350#include "bas.fh" 2351#include "stdio.fh" 2352c 2353 double precision ddot 2354 external ddot 2355c 2356 integer geom 2357 integer basis 2358 integer nbf 2359 integer nat 2360 integer lbuf 2361 integer lscr 2362 double precision buf(lbuf), bufp(lbuf), bufm(lbuf) 2363 double precision scr(lscr) 2364 double precision eri3t(nbf,nbf,nbf) 2365 double precision eri3(nbf,nbf,nbf) 2366 double precision xyz(3,nat) 2367c 2368 integer IR 2369 double precision thresh, norm 2370 double precision RJ(3), RK(3) 2371 integer nshell, cnt, nbfsh 2372 integer ishell, ilo, ihi, nbfshi, i 2373 integer jshell, jlo, jhi, nbfshj, j 2374 integer kshell, klo, khi, nbfshk, k 2375 integer nzero, nzerot, npos, npost 2376 integer h_diff, k_diff 2377c 2378 call dfill(3,0.0d00,RJ,1) 2379 call dfill(3,0.0d00,RK,1) 2380 call dfill((3*nat),0.0d00,xyz,1) 2381 call dfill((nbf*nbf*nbf),0.0d00,eri3t,1) 2382 call dfill((nbf*nbf*nbf),0.0d00,eri3,1) 2383 call dfill(lbuf,0.0d00,buf,1) 2384 call dfill(lbuf,0.0d00,bufp,1) 2385 call dfill(lbuf,0.0d00,bufm,1) 2386 call dfill(lscr,0.0d00,scr,1) 2387 do IR = 1,3 2388 if (IR.eq.1) then 2389 call dfill(3,0.0d00,RJ,1) 2390 call dfill(3,0.0d00,RK,1) 2391 elseif (IR.eq.2) then 2392 call dfill(3,0.0d00,RK,1) 2393 RJ(1) = 1.0d00 2394 RJ(2) = 2.0d00 2395 RJ(3) = 3.0d00 2396 elseif (IR.eq.3) then 2397 RJ(1) = 1.0d00 2398 RJ(2) = 2.0d00 2399 RJ(3) = 3.0d00 2400 RK(1) = 3.0d00 2401 RK(2) = 4.0d00 2402 RK(3) = 5.0d00 2403 else 2404 stop ' how did IR get here' 2405 endif 2406 if (.not.ma_verify_allocator_stuff()) stop ' ma broke 18' 2407 thresh = 1.0d-10 2408* 2409 if (.not.bas_numcont(basis,nshell)) stop 'bas_numcont' 2410 2411 npos = 0 2412 nzero = 0 2413 nzerot = 0 2414 npost = 0 2415 do ishell = 1, nshell ! basis functon NOT translated 2416 if (.not.bas_cn2bfr(basis,ishell,ilo,ihi)) 2417 & stop 'cn2bfr error i' 2418 nbfshi = ihi - ilo + 1 2419 do jshell = 1, nshell ! basis functon translated 2420 if (.not.bas_cn2bfr(basis,jshell,jlo,jhi)) 2421 & stop 'cn2bfr error j' 2422 nbfshj = jhi - jlo + 1 2423 do kshell = 1, nshell ! fitting function 2424 if (.not.bas_cn2bfr(basis,kshell,klo,khi)) 2425 & stop 'cn2bfr error k' 2426 nbfshk = khi - klo + 1 2427 nbfsh = nbfshi*nbfshj*nbfshk 2428 call dfill(lscr,0.0d00,scr,1) 2429 call dfill(lbuf,0.0d00,buf,1) 2430 call intp_2e3c(basis,kshell,basis,jshell,ishell, 2431 & Rk,Rj,lscr,scr,lbuf,buf) 2432 2433 cnt = 1 2434 do k = klo, khi 2435 do j = jlo, jhi 2436 do i = ilo, ihi 2437* write(70,*)i,j,k,buf(cnt) 2438 npos = npos + 1 2439 if (abs(buf(cnt)).lt.thresh) nzero = nzero+1 2440 eri3(i,j,k) = buf(cnt) 2441 cnt = cnt + 1 2442 enddo 2443 enddo 2444 enddo 2445 2446 call dfill(lscr,0.0d00,scr,1) 2447 call dfill(lbuf,0.0d00,buf,1) 2448 call intp_2e3ct(basis,ishell,jshell,basis,kshell, 2449 & RJ,RK,lscr,scr,lbuf,buf) 2450* 2451 cnt = 1 2452 do i = ilo, ihi 2453 do j = jlo, jhi 2454 do k = klo, khi 2455* write(71,*)i,j,k,buf(cnt) 2456 npost = npost + 1 2457 if (abs(buf(cnt)).lt.thresh) nzerot=nzerot+1 2458 eri3t(i,j,k) = buf(cnt) 2459 cnt = cnt + 1 2460 enddo 2461 enddo 2462 enddo 2463c 2464 enddo 2465 enddo 2466 enddo 2467c 2468 if (.not.ma_alloc_get(mt_dbl, 2469 & (nbf*nbf*nbf), 2470 & 'diff buffer', 2471 & h_diff, 2472 & k_diff)) stop '2e3ct: ma alloc failed' 2473 call dcopy((nbf*nbf*nbf),eri3,1,dbl_mb(k_diff),1) 2474 call daxpy((nbf*nbf*nbf),-1.0d00,eri3t,1,dbl_mb(k_diff),1) 2475 norm = ddot((nbf*nbf*nbf),dbl_mb(k_diff),1,dbl_mb(k_diff),1) 2476 write(luout,*)' 2e3ct info:' 2477 write(luout,*)' number possible in regular:',npos 2478 write(luout,*)' number zero in regular:',nzero 2479 write(luout,*)' number non-zero in regular:',(npos-nzero) 2480 write(luout,*)' number possible in transpo:',npost 2481 write(luout,*)' number zero in transpo:',nzerot 2482 write(luout,*)' number non-zero in transpo:',(npost-nzerot) 2483 write(luout,*)' 2e3ct difference norm:',ir,' ',norm 2484c 2485 pderiv_compute_2e3ct = norm.lt.thresh 2486 if (.not.pderiv_compute_2e3ct) then 2487 call prak2e3c(nbf,eri3,eri3t,dbl_mb(k_diff)) 2488 endif 2489 if (.not.ma_free_heap(h_diff)) stop '2e3ct: ma free failed ' 2490 enddo 2491 end 2492 subroutine prak2e3c(nbf,eri3,eri3t,diff) 2493 implicit none 2494 integer nbf 2495 double precision eri3 (nbf,nbf,nbf) 2496 double precision eri3t(nbf,nbf,nbf) 2497 double precision diff (nbf,nbf,nbf) 2498c 2499 integer i,j,k 2500c 2501 do i = 1,nbf 2502 do j = 1,nbf 2503 do k = 1,nbf 2504 if (abs(diff(i,j,k)).gt.1.0d-07) then 2505 write(6,10000)i,j,k, 2506 & eri3(i,j,k), 2507 & eri3t(i,j,k), 2508 & diff(i,j,k) 2509 endif 2510 enddo 2511 enddo 2512 enddo 251310000 format(1x,'<',i3,',',i3,',',i3,'> =',3(1pd20.10)) 2514 end 2515 logical function pderiv_compute_d2e3ct( 2516 & geom,basis,nbf,nat,lbuf,lscr, 2517 & buf,bufp,bufm,scr,eri3,eri3t,xyz) 2518 implicit none 2519#include "mafdecls.fh" 2520#include "nwc_const.fh" 2521#include "geomP.fh" 2522#include "bas.fh" 2523#include "stdio.fh" 2524c 2525 double precision ddot 2526 external ddot 2527c 2528 integer geom 2529 integer basis 2530 integer nbf 2531 integer nat 2532 integer lbuf 2533 integer lscr 2534 double precision buf(lbuf), bufp(lbuf), bufm(lbuf) 2535 double precision scr(lscr) 2536 double precision eri3t(nbf,nbf,nbf,3,nat) 2537 double precision eri3(nbf,nbf,nbf,3,nat) 2538 double precision xyz(3,nat) 2539c 2540 integer nat3, nintz 2541 integer atom, ixyz, IR 2542 double precision thresh, norm 2543 double precision RJ(3), RK(3) 2544 integer nshell, cnt, nbfsh, watom(4) 2545 integer ishell, ilo, ihi, nbfshi, i 2546 integer jshell, jlo, jhi, nbfshj, j 2547 integer kshell, klo, khi, nbfshk, k 2548 integer nzero, nzerot, npos, npost 2549 integer h_diff, k_diff 2550c 2551 nat3 = 3*nat 2552 call dfill(3,0.0d00,RJ,1) 2553 call dfill(3,0.0d00,RK,1) 2554 call dfill((nat3),0.0d00,xyz,1) 2555 call dfill((nbf*nbf*nbf*nat3),0.0d00,eri3t,1) 2556 call dfill((nbf*nbf*nbf*nat3),0.0d00,eri3,1) 2557 call dfill(lbuf,0.0d00,buf,1) 2558 call dfill(lbuf,0.0d00,bufp,1) 2559 call dfill(lbuf,0.0d00,bufm,1) 2560 call dfill(lscr,0.0d00,scr,1) 2561 do IR = 1,3 2562 if (IR.eq.1) then 2563 call dfill(3,0.0d00,RJ,1) 2564 call dfill(3,0.0d00,RK,1) 2565 elseif (IR.eq.2) then 2566 call dfill(3,0.0d00,RK,1) 2567 RJ(1) = 1.0d00 2568 RJ(2) = 2.0d00 2569 RJ(3) = 3.0d00 2570 elseif (IR.eq.3) then 2571 RJ(1) = 1.0d00 2572 RJ(2) = 2.0d00 2573 RJ(3) = 3.0d00 2574 RK(1) = 3.0d00 2575 RK(2) = 4.0d00 2576 RK(3) = 5.0d00 2577 else 2578 stop ' how did IR get here' 2579 endif 2580 if (.not.ma_verify_allocator_stuff()) stop ' ma broke 18' 2581 thresh = 1.0d-10 2582* 2583 if (.not.bas_numcont(basis,nshell)) stop 'bas_numcont' 2584 2585 npos = 0 2586 nzero = 0 2587 nzerot = 0 2588 npost = 0 2589 do ishell = 1, nshell ! basis functon NOT translated 2590 if (.not.bas_cn2bfr(basis,ishell,ilo,ihi)) 2591 & stop 'cn2bfr error i' 2592 nbfshi = ihi - ilo + 1 2593 do jshell = 1, nshell ! basis functon translated 2594 if (.not.bas_cn2bfr(basis,jshell,jlo,jhi)) 2595 & stop 'cn2bfr error j' 2596 nbfshj = jhi - jlo + 1 2597 do kshell = 1, nshell ! fitting function 2598 if (.not.bas_cn2bfr(basis,kshell,klo,khi)) 2599 & stop 'cn2bfr error k' 2600 nbfshk = khi - klo + 1 2601 nbfsh = nbfshi*nbfshj*nbfshk 2602 call dfill(lscr,0.0d00,scr,1) 2603 call dfill(lbuf,0.0d00,buf,1) 2604 call intpd_2e3c(basis,kshell,basis,jshell,ishell, 2605 & Rk,Rj,lscr,scr,lbuf,buf,watom) 2606 2607 nintz = nbfsh 2608 do atom = 1,4 2609 if (watom(atom).gt.0) then 2610 cnt = ((atom-1)*nintz*3) + 1 2611 do ixyz = 1,3 2612 do k = klo, khi 2613 do j = jlo, jhi 2614 do i = ilo, ihi 2615 npos = npos + 1 2616 if (abs(buf(cnt)).lt.thresh) nzero = nzero+1 2617* write(70,*)i,j,k,ixyz,watom(atom),buf(cnt) 2618 eri3(i,j,k,ixyz,watom(atom)) = buf(cnt) 2619 cnt = cnt + 1 2620 enddo 2621 enddo 2622 enddo 2623 enddo 2624 endif 2625 enddo 2626 2627 call dfill(lscr,0.0d00,scr,1) 2628 call dfill(lbuf,0.0d00,buf,1) 2629 call intpd_2e3ct(basis,ishell,jshell,basis,kshell, 2630 & RJ,RK,lscr,scr,lbuf,buf,watom) 2631* 2632 do atom = 1,4 2633 if (watom(atom).gt.0) then 2634 cnt = ((atom-1)*nintz*3) + 1 2635 do ixyz = 1,3 2636 do i = ilo, ihi 2637 do j = jlo, jhi 2638 do k = klo, khi 2639 npost = npost + 1 2640 if (abs(buf(cnt)).lt.thresh) nzerot=nzerot+1 2641* write(71,*)i,j,k,ixyz,watom(atom),buf(cnt) 2642 eri3t(i,j,k,ixyz,watom(atom)) = buf(cnt) 2643 cnt = cnt + 1 2644 enddo 2645 enddo 2646 enddo 2647 enddo 2648 endif 2649 enddo 2650c 2651 enddo 2652 enddo 2653 enddo 2654c 2655 if (.not.ma_alloc_get(mt_dbl, 2656 & (nbf*nbf*nbf*nat3), 2657 & 'diff buffer', 2658 & h_diff, 2659 & k_diff)) stop 'd2e3ct: ma alloc failed' 2660 call dcopy((nbf*nbf*nbf*nat3),eri3,1,dbl_mb(k_diff),1) 2661 call daxpy((nbf*nbf*nbf*nat3),-1.0d00,eri3t,1,dbl_mb(k_diff),1) 2662 norm = 2663 & ddot((nbf*nbf*nbf*nat3),dbl_mb(k_diff),1,dbl_mb(k_diff),1) 2664 write(luout,*)' d2e3ct info:' 2665 write(luout,*)' number possible in regular:',npos 2666 write(luout,*)' number zero in regular:',nzero 2667 write(luout,*)' number non-zero in regular:',(npos-nzero) 2668 write(luout,*)' number possible in transpo:',npost 2669 write(luout,*)' number zero in transpo:',nzerot 2670 write(luout,*)' number non-zero in transpo:',(npost-nzerot) 2671 write(luout,*)' d2e3ct difference norm:',ir,' ',norm 2672c 2673 pderiv_compute_d2e3ct = norm.lt.thresh 2674 if (.not.pderiv_compute_d2e3ct) then 2675 call prakd2e3c(nbf,nat,eri3,eri3t,dbl_mb(k_diff)) 2676 endif 2677 if (.not.ma_free_heap(h_diff)) stop 'd2e3ct: ma free failed ' 2678 enddo 2679 end 2680 subroutine prakd2e3c(nbf,nat,eri3,eri3t,diff) 2681 implicit none 2682 integer nbf 2683 integer nat 2684 double precision eri3(nbf,nbf,nbf,3,nat) 2685 double precision eri3t(nbf,nbf,nbf,3,nat) 2686 double precision diff(nbf,nbf,nbf,3,nat) 2687c 2688 integer atom,ixyz,i,j,k 2689c 2690 do atom = 1,nat 2691 do ixyz = 1,3 2692 do i = 1,nbf 2693 do j = 1,nbf 2694 do k = 1,nbf 2695 if (abs(diff(i,j,k,ixyz,atom)).gt.1.0d-07) then 2696 write(6,10000)i,j,k,ixyz,atom, 2697 & eri3(i,j,k,ixyz,atom), 2698 & eri3t(i,j,k,ixyz,atom), 2699 & diff(i,j,k,ixyz,atom) 2700 endif 2701 enddo 2702 enddo 2703 enddo 2704 enddo 2705 enddo 270610000 format(1x,'<',i3,',',i3,',',i3,',',i3,',',i3,'> =',3(1pd20.10)) 2707 end 2708