1C> \ingroup task 2C> @{ 3 logical function task_ecp_deriv_check(rtdb) 4 implicit none 5#include "errquit.fh" 6* $Id$ 7c 8#include "stdio.fh" 9#include "mafdecls.fh" 10#include "geom.fh" 11#include "bas.fh" 12c 13 logical int_normalize, raktask25_a 14 external int_normalize, raktask25_a 15c 16 integer rtdb 17c 18 logical status 19 integer basis, geom 20 integer nbf, nat, nshell 21 integer maxg1, maxs1 22 integer hbuf, hbufp, hbufm, hscr, hg, hfd, hxyz 23 integer kbuf, kbufp, kbufm, kscr, kg, kfd, kxyz 24 integer hbufpp, hbufmm, hbuf3 25 integer kbufpp, kbufmm, kbuf3 26c 27 task_ecp_deriv_check = .false. 28c 29c 30 if (.not.geom_create(geom,'geometry')) call errquit 31 & ('geom create failed',911, GEOM_ERR) 32 if (.not.geom_rtdb_load(rtdb,geom,'geometry')) call errquit 33 & ('geom_rtdb_load failed',911, RTDB_ERR) 34c 35 if (.not.bas_create(basis,'ao basis')) call errquit 36 & ('bas_create failed',911, BASIS_ERR) 37 if (.not.bas_rtdb_load(rtdb,geom,basis,'ao basis')) call errquit 38 & ('bas_rtdb_load failed',911, RTDB_ERR) 39c 40 write(6,*)' geom/basis loaded' 41c 42 if (.not.int_normalize(rtdb,basis)) stop ' norm error 1' 43c 44 if (.not. bas_print(basis)) 45 $ call errquit(' basis print failed', 0, BASIS_ERR) 46c 47 if (.not.bas_numbf(basis,nbf)) call errquit 48 & ('numbf failed',911, BASIS_ERR) 49 if (.not.bas_numcont(basis,nshell)) call errquit 50 & ('numcont failed',911, BASIS_ERR) 51c 52 if (.not.geom_ncent(geom,nat)) stop 'geom_ncent fe' 53 write(6,*) 'number of atoms ', nat 54c 55 call intd_init(rtdb, 1, basis) 56 call int_mem_1e(maxg1, maxs1) 57 write(luout,*)' maxg1 = ',maxg1 58 write(luout,*)' maxs1 = ',maxs1 59 maxs1 = max(maxs1,(nbf*nbf*3*nat)) 60 write(luout,*)' maxs1 = ',maxs1, ' after max for copy ' 61 status = .true. 62 status = status .and. 63 & ma_alloc_get(mt_dbl,maxg1,'int buffer' ,hbuf,kbuf) 64 status = status .and. 65 & ma_alloc_get(mt_dbl,maxg1,'int buffer' ,hbuf3,kbuf3) 66 status = status .and. 67 & ma_alloc_get(mt_dbl,maxg1,'int buffer(+)' ,hbufp,kbufp) 68 status = status .and. 69 & ma_alloc_get(mt_dbl,maxg1,'int buffer(2+)' ,hbufpp,kbufpp) 70 status = status .and. 71 & ma_alloc_get(mt_dbl,maxg1,'int buffer(-)' ,hbufm,kbufm) 72 status = status .and. 73 & ma_alloc_get(mt_dbl,maxg1,'int buffer(2-)' ,hbufmm,kbufmm) 74 status = status .and. 75 & ma_alloc_get(mt_dbl,maxs1,'scr buffer' ,hscr,kscr) 76 status = status .and. 77 & ma_alloc_get(mt_dbl,(nbf*nbf*3*nat),'grad' ,hg,kg) 78 status = status .and. 79 & ma_alloc_get(mt_dbl,(nbf*nbf*3*nat), 80 & 'finite difference grad' ,hfd,kfd) 81 status = status .and. 82 & ma_alloc_get(mt_dbl,3*nat,'my coords',hxyz,kxyz) 83 if (.not.status) stop ' memory failed rak25 ' 84 call intd_terminate() 85 task_ecp_deriv_check = raktask25_a(rtdb, 86 & geom, basis, nbf, nat, nshell, maxg1, maxs1, 87 & dbl_mb(kg), dbl_mb(kfd), 88 & dbl_mb(kbuf), 89 & dbl_mb(kbuf3), 90 & dbl_mb(kscr), 91 & dbl_mb(kbufp), 92 & dbl_mb(kbufm), 93 & dbl_mb(kbufpp), 94 & dbl_mb(kbufmm), 95 & dbl_mb(kxyz)) 96 status = .true. 97 status = status.and.ma_free_heap(hg) 98 status = status.and.ma_free_heap(hfd) 99 status = status.and.ma_free_heap(hbuf) 100 status = status.and.ma_free_heap(hbuf3) 101 status = status.and.ma_free_heap(hscr) 102 status = status.and.ma_free_heap(hbufp) 103 status = status.and.ma_free_heap(hbufm) 104 status = status.and.ma_free_heap(hbufpp) 105 status = status.and.ma_free_heap(hbufmm) 106 status = status.and.ma_free_heap(hxyz) 107 status = status.and.bas_destroy(basis) 108 status = status.and.geom_destroy(geom) 109 task_ecp_deriv_check = task_ecp_deriv_check.and.status 110 end 111C> @} 112 logical function raktask25_a(rtdb, 113 & geom, basis, nbf, nat, nshell, 114 & sizeg, sizescr, 115 & grad, fdgrad, buf, buf3, scr, bufp, bufm, 116 & bufpp, bufmm, xyz) 117 implicit none 118#include "nwc_const.fh" 119#include "mafdecls.fh" 120#include "geom.fh" 121#include "geomP.fh" 122#include "basdeclsP.fh" 123#include "basP.fh" 124#include "bas.fh" 125#include "stdio.fh" 126#include "geobasmapP.fh" 127#include "bas_exndcf_dec.fh" 128#include "bas_ibs_dec.fh" 129c 130 logical int_ecp_init 131 external int_ecp_init 132 double precision ddot 133 external ddot 134c 135 logical ignore 136 integer rtdb, geom, basis, nbf, nat, nshell 137 integer sizeg, sizescr, ecpid, soid 138 double precision xyz(3,nat) 139 double precision grad(nbf,nbf,3,nat) 140 double precision fdgrad(nbf,nbf,3,nat) 141 double precision buf(sizeg), buf3(sizeg) 142 double precision bufp(sizeg), bufm(sizeg) 143 double precision bufpp(sizeg), bufmm(sizeg) 144 double precision scr(sizescr) 145c 146 double precision thresh, delta, factor, norm 147 double precision asum, fdsum 148 integer xbas, xsize 149 integer ii_np, ii_gen, ii_exp, ii_cf, ii_type, ii_atom 150 integer jj_np, jj_gen, jj_exp, jj_cf, jj_type, jj_atom 151 integer ishell, ilo, ihi, nbfshi 152 integer jshell, jlo, jhi, nbfshj 153 integer nbfsh, ucont 154 integer atom, ixyz, cnt, i, j 155 integer nat3, nshell_use 156c 157#include "bas_exndcf_sfn.fh" 158#include "bas_ibs_sfn.fh" 159c 160 call intd_init(rtdb,1,basis) 161c 162 nat3 = 3*nat 163 call dfill((3*nat),0.0d00,xyz,1) 164 call dfill((nbf*nbf*3*nat),0.0d00,grad,1) 165 call dfill((nbf*nbf*3*nat),0.0d00,fdgrad,1) 166 call dfill(sizeg,0.0d00,buf,1) 167 call dfill(sizeg,0.0d00,bufp,1) 168 call dfill(sizeg,0.0d00,bufm,1) 169 call dfill(sizeg,0.0d00,bufpp,1) 170 call dfill(sizeg,0.0d00,bufmm,1) 171 call dfill(sizescr,0.0d00,scr,1) 172* store original coordintates 173* write(6,*)' original coordinates ' 174* if (.not.geom_print(geom)) stop ' gp error' 175 call dcopy(nat3,coords(1,1,geom),1,xyz,1) 176c 177 ignore = bas_get_ecp_handle(basis,ecpid) 178 write(6,*)' ecp id ',ecpid 179 soid = 0 180 thresh = 1.0d-09 181 delta = 0.0025d00 182 write(6,*)' delta =',delta 183 xbas = basis + BASIS_HANDLE_OFFSET 184c 185 nshell_use = nshell 186* nshell_use = 1 187 do ishell = 1, nshell_use 188 write(6,*)' fd: ishell = ',ishell,' of ',nshell_use 189 call util_flush(6) 190 if (.not.bas_cn2bfr(basis,ishell,ilo,ihi)) 191 & stop 'cn2bfr error i' 192 nbfshi = ihi - ilo + 1 193 ucont = (sf_ibs_cn2ucn(ishell,xbas)) 194 ii_np = infbs_cont(CONT_NPRIM,ucont,xbas) 195 ii_gen = infbs_cont(CONT_NGEN,ucont,xbas) 196 ii_exp = infbs_cont(CONT_IEXP,ucont,xbas) 197 ii_cf = infbs_cont(CONT_ICFP,ucont,xbas) 198 ii_type = infbs_cont(CONT_TYPE,ucont,xbas) 199 ii_atom = (sf_ibs_cn2ce(ishell,xbas)) 200 do jshell = 1, ishell 201 if (.not.bas_cn2bfr(basis,jshell,jlo,jhi)) 202 & stop 'cn2bfr error j' 203 nbfshj = jhi - jlo + 1 204 nbfsh = nbfshi*nbfshj 205* write(6,*)' fd: jshell = ',jshell,' size =',nbfsh 206 ucont = (sf_ibs_cn2ucn(jshell,xbas)) 207 jj_np = infbs_cont(CONT_NPRIM,ucont,xbas) 208 jj_gen = infbs_cont(CONT_NGEN,ucont,xbas) 209 jj_exp = infbs_cont(CONT_IEXP,ucont,xbas) 210 jj_cf = infbs_cont(CONT_ICFP,ucont,xbas) 211 jj_type = infbs_cont(CONT_TYPE,ucont,xbas) 212 jj_atom = (sf_ibs_cn2ce(jshell,xbas)) 213 do atom = 1,nat 214 do ixyz = 1,3 215 call dcopy(nat3,xyz,1,coords(1,1,geom),1) 216 coords(ixyz,atom,geom) = 217 & coords(ixyz,atom,geom) + delta + delta 218 call dfill(sizescr,0.0d00,scr,1) 219 call dfill(sizeg,0.0d00,buf,1) 220 call dfill(sizeg,0.0d00,bufpp,1) 221c 222* write(6,*)' plus ', atom, ixyz 223* if (.not.geom_print(geom)) stop ' gp error' 224 call int_ecp_terminate() 225 ignore = int_ecp_init(ecpid, soid,1) 226 call int_ecp_hf1( 227 & coords(1,ii_atom,geom), 228 & dbl_mb(mb_exndcf(ii_exp,xbas)), 229 & dbl_mb(mb_exndcf(ii_cf,xbas)), 230 & ii_np, ii_gen, ii_type, ii_atom, 231 232 & coords(1,jj_atom,geom), 233 & dbl_mb(mb_exndcf(jj_exp,xbas)), 234 & dbl_mb(mb_exndcf(jj_cf,xbas)), 235 & jj_np, jj_gen, jj_type, jj_atom, 236 237 & bufpp,nbfsh,scr,sizescr,.false.) 238* if (nbfsh.eq.1) write(6,*) 239* & ' ',atom,' X=',ixyz,' +', bufp(1) 240* 241 call dcopy(nat3,xyz,1,coords(1,1,geom),1) 242 coords(ixyz,atom,geom) = coords(ixyz,atom,geom) + delta 243 call dfill(sizescr,0.0d00,scr,1) 244 call dfill(sizeg,0.0d00,bufp,1) 245c 246* write(6,*)' plus ', atom, ixyz 247* if (.not.geom_print(geom)) stop ' gp error' 248 call int_ecp_terminate() 249 ignore = int_ecp_init(ecpid, soid,1) 250 call int_ecp_hf1( 251 & coords(1,ii_atom,geom), 252 & dbl_mb(mb_exndcf(ii_exp,xbas)), 253 & dbl_mb(mb_exndcf(ii_cf,xbas)), 254 & ii_np, ii_gen, ii_type, ii_atom, 255 256 & coords(1,jj_atom,geom), 257 & dbl_mb(mb_exndcf(jj_exp,xbas)), 258 & dbl_mb(mb_exndcf(jj_cf,xbas)), 259 & jj_np, jj_gen, jj_type, jj_atom, 260 261 & bufp,nbfsh,scr,sizescr,.false.) 262* if (nbfsh.eq.1) write(6,*) 263* & ' ',atom,' X=',ixyz,' +', bufp(1) 264* 265 call dcopy(nat3,xyz,1,coords(1,1,geom),1) 266 coords(ixyz,atom,geom) = coords(ixyz,atom,geom) - delta 267 call dfill(sizescr,0.0d00,scr,1) 268 call dfill(sizeg,0.0d00,bufm,1) 269c 270* write(6,*)' minus ', atom, ixyz 271* if (.not.geom_print(geom)) stop ' gp error' 272 call int_ecp_terminate() 273 ignore = int_ecp_init(ecpid, soid,1) 274 call int_ecp_hf1( 275 & coords(1,ii_atom,geom), 276 & dbl_mb(mb_exndcf(ii_exp,xbas)), 277 & dbl_mb(mb_exndcf(ii_cf,xbas)), 278 & ii_np, ii_gen, ii_type, ii_atom, 279 280 & coords(1,jj_atom,geom), 281 & dbl_mb(mb_exndcf(jj_exp,xbas)), 282 & dbl_mb(mb_exndcf(jj_cf,xbas)), 283 & jj_np, jj_gen, jj_type, jj_atom, 284 285 & bufm,nbfsh,scr,sizescr,.false.) 286* 287 call dcopy(nat3,xyz,1,coords(1,1,geom),1) 288 coords(ixyz,atom,geom) = 289 & coords(ixyz,atom,geom) - delta - delta 290 call dfill(sizescr,0.0d00,scr,1) 291 call dfill(sizeg,0.0d00,bufmm,1) 292c 293* write(6,*)' minus ', atom, ixyz 294* if (.not.geom_print(geom)) stop ' gp error' 295 call int_ecp_terminate() 296 ignore = int_ecp_init(ecpid, soid,1) 297 call int_ecp_hf1( 298 & coords(1,ii_atom,geom), 299 & dbl_mb(mb_exndcf(ii_exp,xbas)), 300 & dbl_mb(mb_exndcf(ii_cf,xbas)), 301 & ii_np, ii_gen, ii_type, ii_atom, 302 303 & coords(1,jj_atom,geom), 304 & dbl_mb(mb_exndcf(jj_exp,xbas)), 305 & dbl_mb(mb_exndcf(jj_cf,xbas)), 306 & jj_np, jj_gen, jj_type, jj_atom, 307 308 & bufmm,nbfsh,scr,sizescr,.false.) 309* if (nbfsh.eq.1) write(6,*) 310* & ' ',atom,' X=',ixyz,' -', bufm(1) 311* 312*.............. 3 point 313 factor = 1.0d00/(2.0d00*delta) 314 call dcopy(nbfsh,bufp,1,buf3,1) 315 call daxpy(nbfsh,-1.0d00,bufm,1,buf3,1) 316 call dscal(nbfsh,factor,buf3,1) 317*.............. 5 point 318 factor = 2.0d0/(3.0d0*delta) 319 call dscal(nbfsh,factor,bufp,1) 320 call dscal(nbfsh,factor,bufm,1) 321 factor = 1.0d00/(12.0d00*delta) 322 call dscal(nbfsh,factor,bufpp,1) 323 call dscal(nbfsh,factor,bufmm,1) 324 call dcopy(nbfsh, bufp, 1,buf,1) 325 call daxpy(nbfsh,-1.0d00, bufm, 1,buf,1) 326 call daxpy(nbfsh,-1.0d00, bufpp, 1,buf,1) 327 call daxpy(nbfsh, 1.0d00, bufmm, 1,buf,1) 328 329* if (nbfsh.eq.1) write(6,*) 330* & ' ',atom,' X=',ixyz,' g', buf(1) 331 cnt = 1 332 do i = ilo,ihi 333 do j = jlo, jhi 334 fdgrad(i,j,ixyz,atom) = buf(cnt) 335 fdgrad(j,i,ixyz,atom) = buf(cnt) 336 grad(i,j,ixyz,atom) = buf3(cnt) ! grad has 3 point 337 grad(j,i,ixyz,atom) = buf3(cnt) 338 cnt = cnt + 1 339 enddo 340 enddo 341 if (.not.ma_verify_allocator_stuff()) 342 & stop ' ma broke 3' 343 enddo 344 enddo 345 enddo 346 enddo 347 write(6,*)' differences between 3 point and 5 point ' 348 call rak25_diff(fdgrad,grad,nbf,nat) 349c 350* call prt25(fdgrad,nbf,nat,'fd gradient') 351c 352 call dcopy(nat3,xyz,1,coords(1,1,geom),1) 353 call dfill((nbf*nbf*3*nat),0.0d00,grad,1) 354 call int_ecp_terminate() 355 ignore = int_ecp_init(ecpid, soid,1) 356* 357 nshell_use = nshell 358* nshell_use = 1 359 do ishell = 1, nshell_use 360 write(6,*)' a: ishell = ',ishell,' of ',nshell_use 361 call util_flush(6) 362 if (.not.bas_cn2bfr(basis,ishell,ilo,ihi)) 363 & stop 'cn2bfr error i' 364 nbfshi = ihi - ilo + 1 365 ucont = (sf_ibs_cn2ucn(ishell,xbas)) 366 ii_np = infbs_cont(CONT_NPRIM,ucont,xbas) 367 ii_gen = infbs_cont(CONT_NGEN,ucont,xbas) 368 ii_exp = infbs_cont(CONT_IEXP,ucont,xbas) 369 ii_cf = infbs_cont(CONT_ICFP,ucont,xbas) 370 ii_type = infbs_cont(CONT_TYPE,ucont,xbas) 371 ii_atom = (sf_ibs_cn2ce(ishell,xbas)) 372 do jshell = 1, ishell 373 if (.not.bas_cn2bfr(basis,jshell,jlo,jhi)) 374 & stop 'cn2bfr error j' 375 nbfshj = jhi - jlo + 1 376 nbfsh = nbfshi*nbfshj 377* write(6,*)' a: jshell = ',jshell,' size =',nbfsh 378 ucont = (sf_ibs_cn2ucn(jshell,xbas)) 379 jj_np = infbs_cont(CONT_NPRIM,ucont,xbas) 380 jj_gen = infbs_cont(CONT_NGEN,ucont,xbas) 381 jj_exp = infbs_cont(CONT_IEXP,ucont,xbas) 382 jj_cf = infbs_cont(CONT_ICFP,ucont,xbas) 383 jj_type = infbs_cont(CONT_TYPE,ucont,xbas) 384 jj_atom = (sf_ibs_cn2ce(jshell,xbas)) 385c 386 call dfill(sizescr,0.0d00,scr,1) 387 call dfill(sizeg,0.0d00,buf,1) 388c 389 call intd_ecp_hf1( 390 & coords(1,ii_atom,geom), 391 & dbl_mb(mb_exndcf(ii_exp,xbas)), 392 & dbl_mb(mb_exndcf(ii_cf,xbas)), 393 & ii_np, ii_gen, ii_type, ii_atom, 394 395 & coords(1,jj_atom,geom), 396 & dbl_mb(mb_exndcf(jj_exp,xbas)), 397 & dbl_mb(mb_exndcf(jj_cf,xbas)), 398 & jj_np, jj_gen, jj_type, jj_atom, 399 400 & buf,nbfsh,nat,scr,sizescr,.false.) 401 xsize = nbfsh*3*nat 402* if (nbfsh.eq.1) then 403* do i = 1, 3*nat 404* write(6,*)' grad buff ',i,' ',buf(i) 405* enddo 406* endif 407 cnt = 1 408 do atom=1,nat 409 do ixyz = 1,3 410 do i = ilo, ihi 411 do j = jlo, jhi 412 grad(i,j,ixyz,atom) = buf(cnt) 413 grad(j,i,ixyz,atom) = buf(cnt) 414 cnt = cnt + 1 415 enddo 416 enddo 417 enddo 418 enddo 419c 420 enddo 421 enddo 422 423c 424* call prt25(grad,nbf,nat,' gradient') 425c 426 xsize = nbf*nbf*3*nat 427 do i = 1,nbf 428 do j = 1,i 429 asum = 0.0d00 430 fdsum = 0.0d00 431 do atom = 1,nat 432 do ixyz = 1,3 433 asum = asum + grad(i,j,ixyz,atom) 434 fdsum = fdsum + fdgrad(i,j,ixyz,atom) 435 enddo 436 enddo 437 if ((abs(asum).gt.1.0d-06).or.(abs(fdsum).gt.1.0d-06)) 438 & write(6,*)i,j,asum,fdsum 439 enddo 440 enddo 441c 442 xsize = nbf*nbf*3*nat 443 call dcopy(xsize,grad,1,scr,1) 444 call daxpy(xsize,-1.0d00,fdgrad,1,scr,1) 445 norm = ddot(xsize,scr,1,scr,1) 446 write(luout,*)' difference norm for ecp derivs: ',norm 447 call util_flush(luout) 448c 449 call intd_terminate() 450 raktask25_a = norm .lt. thresh 451 if (raktask25_a) return 452c 453 call rak25_diff(fdgrad,grad,nbf,nat) 454c 455 end 456 subroutine prt25(g,nbf,nat,msg) 457 implicit none 458#include "stdio.fh" 459 integer nbf 460 integer nat 461 double precision g(nbf,nbf,3,nat) 462 character*(*) msg 463c 464 double precision thresh, val 465 integer i,j,x,a 466c 467 thresh = 1.0d-06 468 write(luout,*)'........................................ ',msg 469 do i = 1,nbf 470 do j = 1,i 471 do x = 1,3 472 do a = 1,nat 473 val = g(i,j,x,a) 474 if (abs(val).gt.thresh) then 475 write(luout,10000)i,j,x,a,val 476 endif 477 enddo 478 enddo 479 enddo 480 enddo 481 write(luout,*) 482 write(luout,*) 48310000 format(1x,'g(',4i5,') =',1pd20.10) 484 end 485 subroutine rak25_diff(fd,g,nbf,nat) 486 implicit none 487#include "stdio.fh" 488 integer nbf, nat 489 double precision fd(nbf,nbf,3,nat) 490 double precision g(nbf,nbf,3,nat) 491c 492 double precision aval, val, thresh 493 integer i, j, x, a, t 494 integer hgram(10) 495 double precision tt(9) 496 logical pheader 497c 498 call ifill(10,0,hgram,1) 499 do i = 1,9 500 tt(i) = 10.0d00**(-(12-i)) 501 enddo 502* write(6,*)' tt is ', tt 503 write(luout,*) 504 write(luout,*) 505 write(luout,*) 506 pheader = .false. 507 thresh = 1.0d-05 508 do i = 1,nbf 509 do j = 1,nbf 510 do x = 1,3 511 do a = 1,nat 512 val = g(i,j,x,a) - fd(i,j,x,a) 513 aval = abs(val) 514 if (aval.lt.tt(1)) then 515 hgram(1) = hgram(1) + 1 516 else if (aval.gt.tt(9)) then 517 hgram(10) = hgram(10) + 1 518 else 519 do t = 2,9 520 if (aval.le.tt(t).and.aval.gt.tt((t-1))) 521 & hgram(t) = hgram(t) + 1 522 enddo 523 endif 524 if (aval.gt.thresh)then 525 if (.not.pheader) then 526 write(luout,10000) 527 pheader = .true. 528 endif 529 write(luout,10001) 530 & i,j,x,a,g(i,j,x,a),fd(i,j,x,a),val 531 endif 532 enddo 533 enddo 534 enddo 535 enddo 536c 53710010 format(1x,'Difference count <',1x,1pd8.2,18x,i10,' values') 53810011 format(1x,'Difference count <',1x, 539 & 1pd8.2,' and >',1x,1pd8.2,3x,i10,' values') 54010012 format(1x,'Difference count ',15x,'>',1x,1pd8.2,3x,i10,' values') 541 write(luout,10010)tt(1),hgram(1) 542 do i = 2,9 543 write(luout,10011) tt(i),tt(i-1),hgram(i) 544 enddo 545 write(luout,10012)tt(9),hgram(10) 546* 12345 54710000 format(1x,' i ',' j ',' x ',' a ',3x, 548* 12345678901234567890123 549 & ' ---- analytic ------',3x, 550 & ' - FiniteDifference -',3x, 551 & ' --- difference -----') 55210001 format(1x,4i5,3(3x,1pd20.10)) 553 end 554 555