1 logical function raktest(rtdb) 2 implicit none 3#include "errquit.fh" 4c $Id$ 5#include "mafdecls.fh" 6#include "rtdb.fh" 7#include "context.fh" 8#include "global.fh" 9#include "stdio.fh" 10c::functions 11 logical task_hondo_deriv_check 12 external task_hondo_deriv_check 13 logical task_ecp_deriv_check, task_ecp_print_integrals 14 external task_ecp_deriv_check, task_ecp_print_integrals 15 logical task_computeSld, task_printsoints, task_pderiv 16 external task_computeSld, task_printsoints, task_pderiv 17 logical task_dddd, task_accy, raktask_intdd, raktask_ecppe 18 external task_dddd, task_accy, raktask_intdd, raktask_ecppe 19 logical raktask_intdd_3c 20 external raktask_intdd_3c 21 logical raktask_geomcalc 22 external raktask_geomcalc 23 logical raktask_fullsc 24 external raktask_fullsc 25 logical rak_justrunvib 26 external rak_justrunvib 27c::passed 28 integer rtdb ! rtdb handle 29c::local 30 integer me 31 integer raktask, rak_tmp 32c 33 call ga_sync() 34 call ga_sync() 35 raktask = 0 36 if (rtdb_get(rtdb,'raktask',MT_INT,1,rak_tmp)) 37 & raktask = rak_tmp 38c 39 call ga_sync() 40 call ga_sync() 41 me = ga_nodeid() 42c 43 if (raktask.eq.0) then !................................... 0 44 if (me.eq.0) then 45 write(luout,*)' default raktest task ' 46 write(luout,*)' test semi empirical interface ' 47 endif 48 call raktest_semi(rtdb) 49 raktest = .true. 50 else if (raktask.eq.1) then !................................. 1 51 if (me.eq.0)write(luout,*)' raktest task 1 stepper test' 52 call raktest_stpr(rtdb) 53 raktest = .true. 54 else if (raktask.eq.2) then !................................. 2 55 if (me.eq.0)write(luout,*)' raktest task 2 check int_init' 56 call raktest_init(rtdb) 57 raktest = .true. 58 else if (raktask.eq.3) then !................................. 3 59 if (me.eq.0)write(luout,*)' raktest task 3 check intd_init' 60 call raktest_initd(rtdb) 61 raktest = .true. 62 else if (raktask.eq.4) then !................................. 4 63 if (me.eq.0)write(luout,*)' raktest check 3ctr nai' 64 call raktest_3ctr(rtdb) 65 raktest = .true. 66 else if (raktask.eq.5) then !................................. 5 67 if (me.eq.0)write(luout,*)' test of general contraction code ' 68 call raktest_gc(rtdb) 69 raktest = .true. 70 else if (raktask.eq.6) then !................................. 6 71 if (me.eq.0)write(luout,*)' test of orbital printing code ' 72 call raktest_printorb(rtdb) 73 raktest = .true. 74 else if (raktask.eq.7) then !................................. 7 75 if (me.eq.0)write(luout,*)' test of writing geom objects out ' 76 call raktest_geomwrt(rtdb) 77 raktest = .true. 78 else if (raktask.eq.8) then !................................. 8 79 if (me.eq.0)write(luout,*)' test of spcart stuff ' 80 call raktest_spcart(rtdb) 81 raktest = .true. 82 else if (raktask.eq.9) then !................................. 9 83 if (me.eq.0)write(luout,*)' test of spcart stuff all in one' 84 call raktest_test9(rtdb) 85 raktest = .true. 86 else if (raktask.eq.10) then !................................. 10 87 if (me.eq.0)write(luout,*)' test of ecp stuff ' 88 call raktest_ecp(rtdb) 89 raktest = .true. 90 else if (raktask.eq.11) then !................................. 11 91 if (me.eq.0)write(luout,*)' bug in integrals test ' 92 call raktest_bug(rtdb) 93 raktest = .true. 94 else if (raktask.eq.12) then !................................. 12 95 if (me.eq.0)write(luout,*)' geometry printing routines ' 96 call raktest_geomprt(rtdb) 97 raktest = .true. 98 else if (raktask.eq.13) then !................................. 13 99 if (me.eq.0)write(luout,*)' test 3 center derivatives ' 100 call raktest_3cd(rtdb) 101 raktest = .true. 102 else if (raktask.eq.14) then !................................. 14 103 if (me.eq.0)write(luout,*)' test derivative overlap ' 104 call raktest_ovd(rtdb) 105 raktest = .true. 106 else if (raktask.eq.15) then !................................. 15 107 if (me.eq.0)write(luout,*) 108 & ' test blocking 2e derivative integral interface' 109 call raktest_bd2e(rtdb) 110 raktest = .true. 111 else if (raktask.eq.16) then !................................. 16 112 if (me.eq.0)write(luout,*)' raktest: disk test code ' 113 call raktest_diskspeed(rtdb) 114 raktest = .true. 115 else if (raktask.eq.17) then !................................. 17 116 if (me.eq.0)write(luout,*) 117 & ' raktest: compare 2e integrals from nwchem & texas' 118 call raktest_2ecompare(rtdb) 119 raktest = .true. 120 else if (raktask.eq.18) then !................................. 18 121 if (me.eq.0)write(luout,*) 122 & ' raktest: compute overlap and linear dependence ' 123 raktest = task_computeSld(rtdb) 124 else if (raktask.eq.19) then !................................. 19 125 if (me.eq.0)write(luout,*) 126 & ' raktest: task print SO integrals ' 127 raktest = task_printSOints(rtdb) 128 else if (raktask.eq.20) then !................................. 20 129 if (me.eq.0)write(luout,*) 130 & ' raktest: task test periodic deriv ' 131 raktest = task_pderiv(rtdb) 132 else if (raktask.eq.21) then !................................. 21 133 if (me.eq.0)write(luout,*) 134 & ' raktest: task test dddd bug' 135 raktest = task_dddd(rtdb) 136 else if (raktask.eq.22) then !................................. 22 137 if (me.eq.0)write(luout,*) 138 & ' raktest: accuracy test for ints/grads ' 139 raktest = task_accy(rtdb) 140 else if (raktask.eq.23) then !................................. 23 141 if (me.eq.0)write(luout,*) 142 & ' raktest: second derivative code ' 143 raktest = raktask_intdd(rtdb) 144 else if (raktask.eq.24) then !................................. 24 145 if (me.eq.0)write(luout,*) 146 & ' raktest: ecp PE code test' 147 raktest = raktask_ecppe(rtdb) 148 else if (raktask.eq.25) then !................................. 25 149 if (me.eq.0) write(luout,*) 150 & ' raktest: ecp deriv code test' 151 raktest = task_ecp_deriv_check(rtdb) 152 else if (raktask.eq.26) then !................................. 26 153 if (me.eq.0) write(luout,*) 154 & ' raktest: ecp print integrals' 155 raktest = task_ecp_print_integrals(rtdb) 156 else if (raktask.eq.27) then !................................. 27 157 if (me.eq.0) write(luout,*) 158 & ' raktest: hondo deriv code test' 159 raktest = task_hondo_deriv_check(rtdb) 160 else if (raktask.eq.28) then !................................. 28 161 raktest = raktask_geomcalc(rtdb) 162 else if (raktask.eq.29) then !................................. 29 163 raktest = raktask_fullsc(rtdb) 164 else if (raktask.eq.30) then !................................. 30 165 if (me.eq.0) write(luout,*) 166 & ' raktest: check 2e3c second derivatives' 167 raktest = raktask_intdd_3c(rtdb) 168 else if (raktask.eq.31) then !................................. 31 169 if (me.eq.0) write(luout,*) 170 & ' raktest: rerun vib' 171 raktest = rak_justrunvib(rtdb) 172 else 173 if (me.eq.0) then 174 write(luout,*)' unknown raktask number :',raktask 175 call errquit('raktest: fatal error',911, INPUT_ERR) 176 endif 177 endif 178c 179 end 180 subroutine raktest_geomwrt(rtdb) 181 implicit none 182#include "errquit.fh" 183c 184#include "stdio.fh" 185#include "geom.fh" 186c 187 integer rtdb 188c 189 character*40 new_geom_name 190 integer geom 191 integer igeom 192c 193 if (.not.geom_create(geom,'geometry')) 194 & call errquit('raktest_geomwrt: geom_create failed?',911, 195 & GEOM_ERR) 196c 197 if (.not. geom_rtdb_load(rtdb,geom,'geometry')) call errquit 198 & ('raktest_geomwrt: geom_rtdb_load -ref failed',911, 199 & RTDB_ERR) 200c 201 do igeom = 1,110 202 new_geom_name = ' ' 203 if (abs(igeom) .lt. 10) then 204 write(new_geom_name,'(''g-'',i1,''-step'')') 205 $ abs(igeom) 206 else if (abs(igeom) .lt. 100) then 207 write(new_geom_name,'(''g-'',i2,''-step'')') 208 $ abs(igeom) 209 else if (abs(igeom) .lt. 1000) then 210 write(new_geom_name,'(''g-'',i3,''-step'')') 211 $ abs(igeom) 212 else 213 write(new_geom_name,'(''g-'',i4,''-step'')') 214 $ abs(igeom) 215 endif 216 217 call sym_geom_project(geom, 1d-4) 218 219 if (.not.geom_rtdb_store(rtdb,geom,new_geom_name)) 220 & call errquit 221 & ('raktest_geomwrt: geom_rtdb_store (of copy) failed',911, 222 & RTDB_ERR) 223 write(luout,*)' stored geometry ',new_geom_name 224 enddo 225 end 226 subroutine raktest_printorb(rtdb) 227 implicit none 228c 229#include "rtdb.fh" 230#include "mafdecls.fh" 231#include "bas.fh" 232#include "geom.fh" 233c 234 integer rtdb ! [input] 235c 236 237c 238 end 239 subroutine gen_bf_tag(basis,i_bf,bf_tag) 240 implicit none 241#include "bas.fh" 242c 243c generate a string that tells all about the basis function 244c structure 245c bf_tag(1:16) = geom_tag() or bas_tag ! user symbol for basis function 246c bf_tag(17:17) = ' ' ! blank 247c bf_tag(18:18) = type ! l, s, p, d, .... 248c bf_tag(19:19) = ' ' ! blank 249c bf_tag(20:27) = xyz_tag() ! xyz's for bf 250c 251 integer basis 252 integer i_bf 253 character*27 bf_tag 254 integer cont 255 integer center 256 character*1 ch_type(-1:5) 257 data ch_type /'l','s','p','d','f','g','h'/ 258c 259c map bf -> cn 260c map bf -> ce -> tag 261c map cn -> type 262c map cn -> bfr -> ic 263 if (.not. bas_bf2cn(basis,i_bf,cont)) stop 'ceq' 264 if (.not. bas_bf2ce(basis,i_bf,center)) stop 'ceq' 265c 266 end 267 subroutine int_xyz_tag(lval,ic,xyz_tag,l_tag) 268 implicit none 269#include "errquit.fh" 270 integer lval ! [input] l value 271 integer ic ! [input] cartesean component 272 integer l_tag ! [input] length of xyz_tag character array 273 character*(*) xyz_tag ! [output] left justified 274c 275 integer nxyz(3) 276 character*1 pxyz(3) 277 integer ixyz, i, j 278 data pxyz /'x','y','z'/ 279 save pxyz 280c 281 if (lval.eq.0) then 282 xyz_tag(1:3) = ' s ' 283 ixyz = 4 284c 285 elseif (lval.eq.-1) then 286 if (ic.eq.1) then 287 xyz_tag(1:3) = ' s ' 288 elseif (ic.eq.2) then 289 xyz_tag(1:3) = ' x ' 290 elseif (ic.eq.2) then 291 xyz_tag(1:3) = ' y ' 292 elseif (ic.eq.2) then 293 xyz_tag(1:3) = ' z ' 294 else 295 call errquit('int_xyz_tag: error on lval=-1,ic=',ic, 296 & INPUT_ERR) 297 endif 298 ixyz = 4 299 elseif (lval.gt.0) then 300 call defNxyz(lval) 301 call getNxyz(lval,ic,nxyz) 302c 303 ixyz = 1 304c 305 do i=1,3 306 do j=1,nxyz(i) 307 xyz_tag(ixyz:ixyz) = pxyz(i) 308 ixyz = ixyz + 1 309 enddo 310 enddo 311 else 312 call errquit('int_xyz_tag: error on lval=',lval, INPUT_ERR) 313 endif 314 do i = ixyz, l_tag 315 xyz_tag(i:i) = ' ' 316 enddo 317 end 318c............................................................................... 319 subroutine raktest_spcart(rtdb) 320 implicit none 321#include "errquit.fh" 322#include "mafdecls.fh" 323#include "bas.fh" 324#include "geom.fh" 325#include "rtdb.fh" 326 integer rtdb ! [input] rtdb handle 327c 328 integer geom, basis 329 integer sp_basis, nshell_sp 330 integer nbf, nbfsq, nbf_sp, nbf_chk, nshell 331 integer max1e, max2e, mscr1, mscr2, m_scr, m_buf 332 integer h_cart_s, h_sph_s, h_scr, h_buf, h_2bfr 333 integer k_cart_s, k_sph_s, k_scr, k_buf, k_2bfr 334 integer h_cart_s2, h_sph_s2 335 integer k_cart_s2, k_sph_s2 336 integer h_eri_1, h_eri_2, h_eri_sp_1, h_eri_sp_2 337 integer k_eri_1, k_eri_2, k_eri_sp_1, k_eri_sp_2 338 double precision norm_cart, norm_sph 339 double precision ddot 340 external ddot 341 logical status 342c 343 logical int_normalize, int_norm_2c 344 external int_normalize, int_norm_2c 345c 346 if (.not.geom_create(geom,'geometry')) call errquit 347 & ('geom create failed',911, GEOM_ERR) 348 if (.not.geom_rtdb_load(rtdb,geom,'geometry')) call errquit 349 & ('geom_rtdb_load failed',911, RTDB_ERR) 350c 351 if (.not.bas_create(basis,'ao basis')) call errquit 352 & ('bas_create failed',911, BASIS_ERR) 353 if (.not.bas_rtdb_load(rtdb,geom,basis,'ao basis')) call errquit 354 & ('bas_rtdb_load failed',911, RTDB_ERR) 355c 356 write(6,*)' geom/basis loaded' 357c 358 write(6,*)' raw basis ' 359 if (.not. bas_print(basis)) 360 $ call errquit(' basis print failed', 0, BASIS_ERR) 361 write(6,*)' first normalize' 362 if (.not.int_normalize(rtdb,basis)) stop ' norm error 1' 363 if (.not. bas_print(basis)) 364 $ call errquit(' basis print failed', 0, BASIS_ERR) 365 write(6,*)' second normalize' 366 if (.not.int_normalize(rtdb,basis)) stop ' norm error 2' 367 if (.not. bas_print(basis)) 368 $ call errquit(' basis print failed', 0, BASIS_ERR) 369c 370 if (.not.bas_numbf(basis,nbf)) call errquit 371 & ('numbf failed',911, BASIS_ERR) 372c 373 nbfsq = nbf*nbf 374 if (.not.ma_push_get(mt_dbl,nbfsq,'square cart overlap', 375 & h_cart_s, k_cart_s)) call errquit 376 & (' cart overlap ma failed ',911, MA_ERR) 377 if (.not.ma_push_get(mt_dbl,nbfsq,'square spher overlap', 378 & h_sph_s, k_sph_s)) call errquit 379 & (' sph overlap ma failed ',911, MA_ERR) 380 if (.not.ma_push_get(mt_dbl,nbfsq,'square cart overlap 2', 381 & h_cart_s2, k_cart_s2)) call errquit 382 & (' cart2 overlap ma failed ',911, MA_ERR) 383 if (.not.ma_push_get(mt_dbl,nbfsq,'square spher overlap 2', 384 & h_sph_s2, k_sph_s2)) call errquit 385 & (' sph2 overlap ma failed ',911, MA_ERR) 386c 387 if (.not.bas_numcont(basis,nshell)) call errquit 388 & ('numcont error',911, BASIS_ERR) 389c 390 call int_init(rtdb,1,basis) 391 call int_mem(max1e,max2e,mscr1,mscr2) 392 m_buf = max(max1e*2,max2e*2) 393 m_scr = max(mscr1*2,mscr2) 394 m_buf = m_buf + (m_buf*110)/100 395 m_scr = m_scr + (m_scr*110)/100 396 397c 398 if (.not.ma_push_get(mt_dbl,m_scr,'scr for 1e',h_scr,k_scr)) 399 & call errquit('ma scr failed',911, MA_ERR) 400c 401 if (.not.ma_push_get(mt_dbl,m_buf,'buf for 1e',h_buf,k_buf)) 402 & call errquit('ma buf failed',911, MA_ERR) 403c 404 if (.not.ma_push_get(mt_int,2*nshell,'buf for sp cn2bfr', 405 & h_2bfr,k_2bfr)) 406 & call errquit('ma buf failed',911, MA_ERR) 407 call rak_tospbfr(basis,nshell,nbf_chk,nbf_sp,int_mb(k_2bfr)) 408c 409 if (nbf_chk.ne.nbf) then 410 write(6,*)' nbf not right ',nbf_chk, nbf 411 endif 412c 413 write(6,*)' nbf ',nbf 414 write(6,*)' nbf_sp ',nbf_sp 415c 416 if (.not.ma_push_get(mt_dbl,(nbf*nbf*nbf*nbf),'eri cart 1', 417 & h_eri_1,k_eri_1)) call errquit('ma failed',911, MA_ERR) 418 if (.not.ma_push_get(mt_dbl,(nbf_sp*nbf_sp*nbf_sp*nbf_sp), 419 & 'eri_sp 1', 420 & h_eri_sp_1,k_eri_sp_1)) call errquit('ma failed',911, 421 & MA_ERR) 422 if (.not.ma_push_get(mt_dbl,(nbf*nbf*nbf*nbf), 423 & 'eri cart 2', 424 & h_eri_2,k_eri_2)) call errquit('ma failed',922, MA_ERR) 425 if (.not.ma_push_get(mt_dbl,(nbf_sp*nbf_sp*nbf_sp*nbf_sp), 426 & 'eri_sp 2', 427 & h_eri_sp_2,k_eri_sp_2)) call errquit('ma failed',922, 428 & MA_ERR) 429c 430c 431 call rak_ovlap_test_sp(basis,nbf,nbf_sp,nshell, 432 & dbl_mb(k_scr),m_scr, 433 & dbl_mb(k_buf),m_buf, 434 & dbl_mb(k_cart_s),dbl_mb(k_sph_s), 435 & int_mb(k_2bfr)) 436c 437 call rak_2el_test_sp(basis,nbf,nbf_sp,nshell, 438 & dbl_mb(k_scr),m_scr, 439 & dbl_mb(k_buf),m_buf, 440 & int_mb(k_2bfr), 441 & dbl_mb(k_eri_sp_1),dbl_mb(k_eri_1)) 442c 443 call rak_ovlap(basis,nbf,nshell, 444 & dbl_mb(k_scr),m_scr, 445 & dbl_mb(k_buf),m_buf, 446 & dbl_mb(k_cart_s2),.false.,'cartcart') 447 448 call rak_2el(basis,nbf,nshell, 449 & dbl_mb(k_scr),m_scr, 450 & dbl_mb(k_buf),m_buf, 451 & dbl_mb(k_eri_2), 452 & .false., ' cartcart ') 453 if (.not.bas_create(sp_basis,'ao sp_basis')) call errquit 454 & ('bas_create failed',911, BASIS_ERR) 455 if (.not.bas_rtdb_load(rtdb,geom,sp_basis,'ao sp_basis')) 456 & call errquit 457 & ('bas_rtdb_load failed',911, RTDB_ERR) 458c 459 write(6,*)' geom/sp_basis loaded' 460c 461 write(6,*)' raw sp_basis ' 462 if (.not. bas_print(sp_basis)) 463 $ call errquit(' sp_basis print failed', 0, BASIS_ERR) 464 write(6,*)' first normalize' 465 if (.not.int_normalize(rtdb,sp_basis)) stop ' norm error 1' 466 if (.not. bas_print(sp_basis)) 467 $ call errquit(' sp_basis print failed', 0, BASIS_ERR) 468 if (.not.bas_numbf(sp_basis,nbf_sp)) call errquit 469 & ('numbf failed',911, BASIS_ERR) 470 if (.not.bas_numcont(sp_basis,nshell_sp)) call errquit 471 & ('numcont error',911, BASIS_ERR) 472 write(6,*)' sp_basis b4 rak_ovlap' 473 if (.not. bas_print(sp_basis)) 474 $ call errquit(' sp_basis print failed', 0, BASIS_ERR) 475 call rak_core(sp_basis,nbf_sp,nshell_sp, 476 & dbl_mb(k_scr),m_scr, 477 & dbl_mb(k_buf),m_buf, 478 & dbl_mb(k_sph_s2),.true.,'spsp') 479 call rak_ovlap(sp_basis,nbf_sp,nshell_sp, 480 & dbl_mb(k_scr),m_scr, 481 & dbl_mb(k_buf),m_buf, 482 & dbl_mb(k_sph_s2),.false.,'spsp') 483 call rak_2el(sp_basis,nbf_sp,nshell_sp, 484 & dbl_mb(k_scr),m_scr, 485 & dbl_mb(k_buf),m_buf, 486 & dbl_mb(k_eri_sp_2), 487 & .true.,' spsp ') 488c 489 call print_diff_vec((nbf*nbf), 490 & dbl_mb(k_cart_s), 491 & dbl_mb(k_cart_s2), 492 & 1.0d-05,' cart s ') 493 call print_diff_vec((nbf_sp*nbf_sp), 494 & dbl_mb(k_sph_s), 495 & dbl_mb(k_sph_s2), 496 & 1.0d-05,' spher s ') 497 call daxpy((nbf*nbf),-1.0d00, 498 & dbl_mb(k_cart_s2),1, 499 & dbl_mb(k_cart_s),1) 500 norm_cart = ddot((nbf*nbf),dbl_mb(k_cart_s),1,dbl_mb(k_cart_s),1) 501 call daxpy((nbf*nbf),-1.0d00, 502 & dbl_mb(k_sph_s2),1, 503 & dbl_mb(k_sph_s),1) 504 norm_sph = ddot((nbf*nbf),dbl_mb(k_sph_s),1,dbl_mb(k_sph_s),1) 505c 506 write(6,*)'1e diff norm_cart:',norm_cart 507 write(6,*)'1e diff norm_sph :',norm_sph 508c 509 call print_diff_vec((nbf*nbf*nbf*nbf), 510 & dbl_mb(k_eri_1), 511 & dbl_mb(k_eri_2), 512 & 1.0d-05,' eri cart ') 513 call print_diff_vec((nbf_sp*nbf_sp*nbf_sp*nbf_sp), 514 & dbl_mb(k_eri_sp_1), 515 & dbl_mb(k_eri_sp_2), 516 & 1.0d-05,' eri spherical ') 517 call daxpy((nbf*nbf*nbf*nbf),-1.0d00, 518 & dbl_mb(k_eri_2),1, 519 & dbl_mb(k_eri_1),1) 520 norm_cart = ddot((nbf*nbf*nbf*nbf), 521 & dbl_mb(k_eri_1),1,dbl_mb(k_eri_1),1) 522 call daxpy((nbf_sp*nbf_sp*nbf_sp*nbf_sp),-1.0d00, 523 & dbl_mb(k_eri_sp_2),1, 524 & dbl_mb(k_eri_sp_1),1) 525 norm_sph = ddot((nbf_sp*nbf_sp*nbf_sp*nbf_sp), 526 & dbl_mb(k_eri_sp_1),1,dbl_mb(k_eri_sp_1),1) 527c 528 write(6,*)'2e diff norm_cart:',norm_cart 529 write(6,*)'2e diff norm_sph :',norm_sph 530c 531 call int_terminate() 532c 533 status = .true. 534 status = status.and.ma_pop_stack(h_eri_sp_2) 535 status = status.and.ma_pop_stack(h_eri_2) 536 status = status.and.ma_pop_stack(h_eri_sp_1) 537 status = status.and.ma_pop_stack(h_eri_1) 538 status = status.and.ma_pop_stack(h_2bfr) 539 status = status.and.ma_pop_stack(h_buf) 540 status = status.and.ma_pop_stack(h_scr) 541 status = status.and.ma_pop_stack(h_sph_s2) 542 status = status.and.ma_pop_stack(h_cart_s2) 543 status = status.and.ma_pop_stack(h_sph_s) 544 status = status.and.ma_pop_stack(h_cart_s) 545c 546 if (.not.status) call errquit('ma pop fail',911, MA_ERR) 547c 548 if(.not.bas_destroy(basis)) call errquit 549 & ('basis bas_destroy failed',911, BASIS_ERR) 550 if(.not.bas_destroy(sp_basis)) call errquit 551 & ('sp_basis bas_destroy failed',911, BASIS_ERR) 552 if(.not.geom_destroy(geom)) call errquit 553 & ('geom_destroy failed',911, GEOM_ERR) 554c 555 end 556*....................................................................... 557 subroutine rak_2el(basis,nbf,nshell, 558 & scr,mscr,buf,mbuf,eri,print_int,msg) 559 implicit none 560#include "bas.fh" 561#include "stdio.fh" 562#include "util.fh" 563 integer basis,nbf,nshell,mscr,mbuf 564 double precision eri(*), scr(mscr), buf(mbuf) 565 logical print_int 566 character*(*) msg 567c 568 integer ish, jsh, ksh, lsh 569 integer ilo, jlo, klo, llo 570 integer ihi, jhi, khi, lhi 571 integer count, indx 572 logical stat_indx 573 integer ii,jj,kk,ll 574c 575 integer i,j,k,l,isym2,isym4 576 isym2(i,j)=max(i,j)*(max(i,j)-1)/2+min(i,j) 577 isym4(i,j,k,l)=max(isym2(i,j),isym2(k,l))* 578 & (max(isym2(i,j),isym2(k,l))-1)/2+ 579 & min(isym2(i,j),isym2(k,l)) 580c 581 write(6,*)' 2el ',msg 582c 583 call dfill((nbf*nbf*nbf*nbf),0.0d00,eri,1) 584c 585 do ish = 1,nshell 586 if (.not.bas_cn2bfr(basis,ish,ilo,ihi)) 587 & stop 'cn2bfr error i' 588 do jsh = 1,ish 589 if (.not.bas_cn2bfr(basis,jsh,jlo,jhi)) 590 & stop 'cn2bfr error j' 591 do ksh = 1,ish 592 if (.not.bas_cn2bfr(basis,ksh,klo,khi)) 593 & stop 'cn2bfr error k' 594 do lsh = 1,ksh 595 if (.not.bas_cn2bfr(basis,lsh,llo,lhi)) 596 & stop 'cn2bfr error l' 597 call int_2e4c 598 & (basis,ish,jsh,basis,ksh,lsh,mscr,scr,mbuf,buf) 599 count = 0 600 do i=ilo,ihi 601 do j=jlo,jhi 602 do k=klo,khi 603 do l=llo,lhi 604 count = count + 1 605 indx = isym4(i,j,k,l) 606 stat_indx = .false. 607 if (stat_indx) then 608 write(6,*)'indx:elel:shells',msg, 609 & indx,ish,jsh,ksh,lsh 610 write(6,*)'indx:elel:labels',msg, 611 & indx,i,j,k,l 612 endif 613 if (print_int.and.(abs(buf(count)).gt.0.0d00)) 614 & then 615 call int_canon(i,j,k,l,ii,jj,kk,ll) 616 write(69,*)ii,jj,kk,ll,buf(count),' 2el ',msg 617 endif 618 eri(indx) = buf(count) 619 enddo 620 enddo 621 enddo 622 enddo 623c.... end of shell loops 62499999 continue 625 enddo 626 enddo 627 enddo 628 enddo 629 630 end 631*....................................................................... 632 subroutine rak_2el_test_sp(basis,nbf, nbf_sp, nshell, 633 & scr,mscr,buf,mbuf,cn2bfr_sp,eri_sp,eri) 634 implicit none 635#include "bas.fh" 636#include "stdio.fh" 637#include "util.fh" 638 integer nbf, nbf_sp, nshell, mscr, mbuf, basis 639 integer cn2bfr_sp(2,nshell) 640 double precision scr(mscr), buf(mbuf) 641 double precision eri_sp(*), eri(*) 642c 643 integer ish, jsh, ksh, lsh 644 integer ilo, jlo, klo, llo 645 integer ihi, jhi, khi, lhi 646 integer ilosp, jlosp, klosp, llosp 647 integer ihisp, jhisp, khisp, lhisp 648 integer inbf, jnbf, knbf, lnbf 649 integer inbf_sp, jnbf_sp, knbf_sp, lnbf_sp 650 integer itype, jtype, ktype, ltype 651 integer igen, jgen, kgen, lgen 652 integer iatom, jatom, katom, latom 653 integer count, indx, lshtop 654 logical stat_indx 655* integer junk2 656 integer junk1,junk3 657 double precision ttrans_w, ttrans_c 658 double precision tcomp_w, tcomp_c 659 double precision tadd_w, tadd_c 660c 661* logical spcart_init, spcart_terminate 662* external spcart_init, spcart_terminate 663c 664 integer i,j,k,l,isym2,isym4 665 isym2(i,j)=max(i,j)*(max(i,j)-1)/2+min(i,j) 666 isym4(i,j,k,l)=max(isym2(i,j),isym2(k,l))* 667 & (max(isym2(i,j),isym2(k,l))-1)/2+ 668 & min(isym2(i,j),isym2(k,l)) 669c 670* if (.not.spcart_init(5,.true.,.false.)) 671* & stop ' spcart_init failed' 672c 673 call dfill((nbf_sp*nbf_sp*nbf_sp*nbf_sp),0.0d00,eri_sp,1) 674 call dfill((nbf*nbf*nbf*nbf),0.0d00,eri,1) 675 ttrans_w = 0.0d00 676 ttrans_c = 0.0d00 677 tcomp_w = 0.0d00 678 tcomp_c = 0.0d00 679 do ish = 1,nshell 680 if (.not.bas_cn2bfr(basis,ish,ilo,ihi)) 681 & stop 'cn2bfr error i' 682 if (.not.bas_continfo 683 & (basis,ish,itype,junk1,igen,junk3)) 684 & stop 'bas_continfo error i' 685 if (.not.bas_cn2ce(basis,ish,iatom)) 686 & stop 'bas_cn2ce error i' 687 ilosp = cn2bfr_sp(1,ish) 688 ihisp = cn2bfr_sp(2,ish) 689 inbf = ihi-ilo + 1 690 inbf_sp = ihisp-ilosp + 1 691 do jsh = 1,ish 692 if (.not.bas_cn2bfr(basis,jsh,jlo,jhi)) 693 & stop 'cn2bfr error j' 694 if (.not.bas_continfo 695 & (basis,jsh,jtype,junk1,jgen,junk3)) 696 & stop 'bas_continfo error j' 697 if (.not.bas_cn2ce(basis,jsh,jatom)) 698 & stop 'bas_cn2ce error j' 699 jlosp = cn2bfr_sp(1,jsh) 700 jhisp = cn2bfr_sp(2,jsh) 701 jnbf = jhi-jlo + 1 702 jnbf_sp = jhisp-jlosp + 1 703 do ksh = 1,ish 704 if (.not.bas_cn2bfr(basis,ksh,klo,khi)) 705 & stop 'cn2bfr error k' 706 if (.not.bas_continfo 707 & (basis,ksh,ktype,junk1,kgen,junk3)) 708 & stop 'bas_continfo error k' 709 if (.not.bas_cn2ce(basis,ksh,katom)) 710 & stop 'bas_cn2ce error k' 711 klosp = cn2bfr_sp(1,ksh) 712 khisp = cn2bfr_sp(2,ksh) 713 knbf = khi-klo + 1 714 knbf_sp = khisp-klosp + 1 715 lshtop = ksh 716 if (ksh.eq.ish) lshtop = jsh 717 do lsh = 1,lshtop 718 if (.not.bas_cn2bfr(basis,lsh,llo,lhi)) 719 & stop 'cn2bfr error l' 720 if (.not.bas_continfo 721 & (basis,lsh,ltype,junk1,lgen,junk3)) 722 & stop 'bas_continfo error l' 723 if (.not.bas_cn2ce(basis,lsh,latom)) 724 & stop 'bas_cn2ce error l' 725 llosp = cn2bfr_sp(1,lsh) 726 lhisp = cn2bfr_sp(2,lsh) 727 lnbf = lhi-llo + 1 728 lnbf_sp = lhisp-llosp + 1 729 tadd_c = util_cpusec() 730 tadd_w = util_wallsec() 731 call int_2e4c 732 & (basis,ish,jsh,basis,ksh,lsh,mscr,scr,mbuf,buf) 733c 734 tadd_c = util_cpusec() - tadd_c 735 tadd_w = util_wallsec() - tadd_w 736 if (tadd_c.gt.0.0d00) tcomp_c = tcomp_c + tadd_c 737 if (tadd_w.gt.0.0d00) tcomp_w = tcomp_w + tadd_w 738*rak: write(luout,10000)ish,jsh,ksh,lsh, 739*rak: & itype,jtype,ktype,ltype, 740*rak: & iatom,jatom,katom,latom 741 count = 0 742 do i=ilo,ihi 743 do j=jlo,jhi 744 do k=klo,khi 745 do l=llo,lhi 746 count = count + 1 747 indx = isym4(i,j,k,l) 748 stat_indx = indx.eq.5566 749 stat_indx = stat_indx.or. 750 & (indx.ge.5772.and.indx.le.5784) 751 stat_indx = stat_indx.or. 752 & (indx.ge.5887.and.indx.le.5901) 753 stat_indx = .false. 754 if (stat_indx) then 755 write(6,*)'indx:cart:shells', 756 & indx,ish,jsh,ksh,lsh 757 write(6,*)'indx:cart:labels', 758 & indx,i,j,k,l 759 endif 760 eri(indx) = buf(count) 761 enddo 762 enddo 763 enddo 764 enddo 765* 766 tadd_c = util_cpusec() 767 tadd_w = util_wallsec() 768 call spcart_bra2etran(buf,scr, 769 & jnbf,inbf,jnbf_sp,inbf_sp, 770 & jtype,itype,jgen,igen, 771 & (knbf*lnbf),.false.) 772 call spcart_ket2etran(buf,scr, 773 & lnbf,knbf,lnbf_sp,knbf_sp, 774 & ltype,ktype,lgen,kgen, 775 & (inbf_sp*jnbf_sp),.false.) 776c 777 tadd_c = util_cpusec() - tadd_c 778 tadd_w = util_wallsec() - tadd_w 779 if (tadd_c.gt.0.0d00) ttrans_c = ttrans_c + tadd_c 780 if (tadd_w.gt.0.0d00) ttrans_w = ttrans_w + tadd_w 781* 782*rak: write(luout,10000)ish,jsh,ksh,lsh, 783*rak: & itype,jtype,ktype,ltype, 784*rak: & iatom,jatom,katom,latom 785*rak: count = 0 786*rak: do i=ilosp,ihisp 787*rak: do j=jlosp,jhisp 788*rak: do k=klosp,khisp 789*rak: do l=llosp,lhisp 790*rak: count = count + 1 791*rak: if (abs(buf(count)).gt.1.0d-07) then 792*rak: write(luout,10002)i,j,k,l,buf(count),count 793*rak: endif 794*rak: enddo 795*rak: enddo 796*rak: enddo 797*rak: enddo 798 count = 0 799 do i=ilosp,ihisp 800 do j=jlosp,jhisp 801 do k=klosp,khisp 802 do l=llosp,lhisp 803 count = count + 1 804 indx = isym4(i,j,k,l) 805 stat_indx = indx.eq.5566 806 stat_indx = stat_indx.or. 807 & (indx.ge.5772.and.indx.le.5784) 808 stat_indx = stat_indx.or. 809 & (indx.ge.5887.and.indx.le.5901) 810 stat_indx = .false. 811 if (stat_indx) then 812 write(6,*)'indx:sp:shells', 813 & indx,ish,jsh,ksh,lsh 814 write(6,*)'indx:sp:labels', 815 & indx,i,j,k,l 816 endif 817 eri_sp(indx) = buf(count) 818 enddo 819 enddo 820 enddo 821 enddo 822c.... end of shell loops 82399999 continue 824 enddo 825 enddo 826 enddo 827 enddo 828c 829* if (.not.spcart_terminate()) stop 'term error' 830 write(luout,*)' total compute time ( cpu): ',tcomp_c 831 write(luout,*)' total compute time (wall): ',tcomp_w 832 write(luout,*)' total tranfrm time ( cpu): ',ttrans_c 833 write(luout,*)' total tranfrm time (wall): ',ttrans_w 834 if (tcomp_c.gt.1.0d-30) 835 & write(luout,'(1x,a,f10.2)') 836 & ' % overhead ( cpu): ', 837 & (ttrans_c/tcomp_c*100.0d00) 838 if (tcomp_w.gt.1.0d-30) 839 & write(luout,'(1x,a,f10.2)') 840 & ' % overhead (wall): ', 841 & (ttrans_w/tcomp_w*100.0d00) 842c 84310000 format( 844 & 'Shells <',i5,i5,'|',i5,i5,'>',5x, 845 & 'Types {',i5,i5,'|',i5,i5,'}',5x, 846 & 'Atoms (',i5,i5,'|',i5,i5,')') 84710001 format('(',i5,i5,'|',i5,i5,') =',1pd20.10,' cart',1x,i10) 84810002 format('[',i5,i5,'|',i5,i5,'] =',1pd20.10,' sphr',1x,i10) 849c 850 end 851*....................................................................... 852 subroutine rak_tospbfr(basis,nshell,nbf_chk,nbf_sp,cn2bfr) 853 implicit none 854#include "errquit.fh" 855#include "bas.fh" 856 integer nshell, nbf_chk, nbf_sp, basis 857 integer cn2bfr(2,nshell) 858 integer type, nprim, ngen, spsp, ish 859c 860 nbf_chk = 0 861 nbf_sp = 0 862 do ish = 1,nshell 863 if(.not.bas_continfo(basis,ish,type,nprim,ngen,spsp)) 864 & call errquit(' continfo failed ',911, BASIS_ERR) 865 cn2bfr(1,ish) = nbf_sp + 1 866 cn2bfr(2,ish) = nbf_sp + 2*type+1 867 nbf_chk = nbf_chk + (type+1)*(type+2)/2 868 nbf_sp = nbf_sp + 2*type+1 869* write(6,10000)ish,type,cn2bfr(1,ish),cn2bfr(2,ish) 870 enddo 871*10000 format('<ish:type><',i3,':',i3,'> range ',i3,' to ',i3) 872 end 873*------------------------------------------------------------------------------- 874 subroutine rak_ovlap(basis,nbf,nshell, 875 & scr,mscr,buf,mbuf,s,print_int,msg) 876 implicit none 877c 878#include "mafdecls.fh" 879#include "bas.fh" 880c 881 integer nbf, nshell, mscr, mbuf, basis 882 double precision scr(mscr), buf(mbuf) 883 double precision s(nbf,nbf) 884 logical print_int 885 character*(*) msg 886c 887 integer count 888 integer ibf, ish, ilow, ihi, nbfi 889 integer jbf, jsh, jlow, jhi, nbfj 890 double precision value 891 logical do_print 892c 893 call dfill((nbf*nbf),0.0d00,s,1) 894 call dfill(mscr,0.0d00,scr,1) 895 call dfill(mbuf,0.0d00,buf,1) 896c 897 do ish = 1,nshell 898 if (.not.bas_cn2bfr(basis,ish,ilow,ihi)) stop 'dead i' 899 nbfi = ihi - ilow + 1 900 do jsh = 1,ish 901 if (.not.bas_cn2bfr(basis,jsh,jlow,jhi)) stop 'dead j' 902 nbfj = jhi - jlow + 1 903c 904 call int_1eov(basis,ish,basis,jsh,mscr,scr,mbuf,buf) 905c 906 count = 0 907 do ibf = ilow,ihi 908 do jbf = jlow, jhi 909 count = count + 1 910 value = buf(count) 911 s(ibf,jbf) = value 912 s(jbf,ibf) = value 913 do_print = print_int 914 do_print = do_print .and. (ibf.ge.jbf) 915 do_print = do_print .and. (abs(value).gt.0.0d00) 916 if (do_print) then 917 write(69,*)ibf, jbf, value, ' ovlap ',msg 918 endif 919 enddo 920 enddo 921 922c. close jsh/ish loops 923 enddo 924 enddo 925 write(6,*)' generic overlap matrix (nbf=',nbf,') <',msg,'>' 926 call output(s,1,nbf,1,nbf,nbf,nbf,1) 927 end 928*------------------------------------------------------------------------------- 929 subroutine rak_core(basis,nbf,nshell, 930 & scr,mscr,buf,mbuf,s,print_int,msg) 931 implicit none 932c 933#include "mafdecls.fh" 934#include "bas.fh" 935c 936 integer nbf, nshell, mscr, mbuf, basis 937 double precision scr(mscr), buf(mbuf) 938 double precision s(nbf,nbf) 939 logical print_int 940 character*(*) msg 941c 942 integer count 943 integer ibf, ish, ilow, ihi, nbfi 944 integer jbf, jsh, jlow, jhi, nbfj 945 double precision value 946 logical do_print 947c 948 call dfill((nbf*nbf),0.0d00,s,1) 949 call dfill(mscr,0.0d00,scr,1) 950 call dfill(mbuf,0.0d00,buf,1) 951c 952 do ish = 1,nshell 953 if (.not.bas_cn2bfr(basis,ish,ilow,ihi)) stop 'dead i' 954 nbfi = ihi - ilow + 1 955 do jsh = 1,ish 956 if (.not.bas_cn2bfr(basis,jsh,jlow,jhi)) stop 'dead j' 957 nbfj = jhi - jlow + 1 958c 959 call int_1eh1(basis,ish,basis,jsh,mscr,scr,mbuf,buf) 960c 961 count = 0 962 do ibf = ilow,ihi 963 do jbf = jlow, jhi 964 count = count + 1 965 value = buf(count) 966 s(ibf,jbf) = value 967 s(jbf,ibf) = value 968 do_print = print_int 969 do_print = do_print .and. (ibf.ge.jbf) 970 do_print = do_print .and. (abs(value).gt.0.0d00) 971 if (do_print) then 972 write(69,*)ibf, jbf, value, ' h1 ',msg 973 endif 974 enddo 975 enddo 976 977c. close jsh/ish loops 978 enddo 979 enddo 980 write(6,*)' generic h1 matrix (nbf=',nbf,') <',msg,'>' 981 call output(s,1,nbf,1,nbf,nbf,nbf,1) 982 end 983*------------------------------------------------------------------------------- 984 subroutine rak_ovlap_test_sp(basis,nbf,nbf_sp,nshell, 985 & scr,mscr,buf,mbuf,s,s_sp,cn2bfr_sp) 986 implicit none 987#include "errquit.fh" 988c 989#include "mafdecls.fh" 990#include "bas.fh" 991c 992 integer nbf, nbf_sp, nshell, mscr, mbuf, basis 993 integer cn2bfr_sp(2,nshell) 994 double precision scr(mscr), buf(mbuf) 995 double precision s(nbf,nbf), s_sp(nbf_sp,nbf_sp) 996c 997 integer ibf, ish, ilow, ihi, ilow_sp, ihi_sp, nbfi, nbfi_sp 998 integer jbf, jsh, jlow, jhi, jlow_sp, jhi_sp, nbfj, nbfj_sp 999 integer typei, typej, nprim, igen, jgen, spsp 1000 integer nint, nint_sp, count, hi_ang, st_ang 1001 integer ii, jj 1002 double precision value 1003*rak: double precision pi, fact 1004c 1005* logical spcart_init, spcart_terminate 1006* external spcart_init, spcart_terminate 1007c 1008* write(6,*)'inside rak_ovlap' 1009 if (.not.bas_high_angular(basis,hi_ang)) stop ' dead ang' 1010 st_ang = hi_ang/2 1011* if (.not.spcart_init(st_ang,.true.,.false.)) stop ' dead sp' 1012* if (.not.spcart_init(hi_ang,.true.,.false.)) stop ' dead sp' 1013c 1014* write(6,*)' mscr = ',mscr 1015* write(6,*)' mbuf = ',mbuf 1016c 1017 call dfill((nbf_sp*nbf_sp),0.0d00,s_sp,1) 1018 call dfill((nbf*nbf),0.0d00,s,1) 1019 call dfill(mscr,0.0d00,scr,1) 1020 call dfill(mbuf,0.0d00,buf,1) 1021 1022 do ish = 1,nshell 1023 if (.not.bas_cn2bfr(basis,ish,ilow,ihi)) stop 'dead i' 1024 ilow_sp = cn2bfr_sp(1,ish) 1025 ihi_sp = cn2bfr_sp(2,ish) 1026 nbfi = ihi - ilow + 1 1027 nbfi_sp = ihi_sp - ilow_sp + 1 1028 if(.not.bas_continfo(basis,ish,typei,nprim,igen,spsp)) 1029 & call errquit(' continfo failed ',911, BASIS_ERR) 1030 do jsh = 1,ish 1031 if (.not.bas_cn2bfr(basis,jsh,jlow,jhi)) stop 'dead j' 1032 jlow_sp = cn2bfr_sp(1,jsh) 1033 jhi_sp = cn2bfr_sp(2,jsh) 1034 nbfj = jhi - jlow + 1 1035 nbfj_sp = jhi_sp - jlow_sp + 1 1036 if(.not.bas_continfo(basis,jsh,typej,nprim,jgen,spsp)) 1037 & call errquit(' continfo failed ',911, BASIS_ERR) 1038c 1039 nint = nbfi*nbfj 1040 nint_sp = nbfi_sp*nbfj_sp 1041c 1042*rak: write(6,*)' ' 1043*rak: write(6,*)'<ish,ilow,ihi,nbfi,typei>',ish,ilow,ihi,nbfi,typei 1044*rak: write(6,*)'<jsh,jlow,jhi,nbfj,typej>',jsh,jlow,jhi,nbfj,typej 1045*rak: write(6,*)' nint = ',nint,' nint(sp) = ',nint_sp 1046*rak: write(6,*)'<ish,ilowsp,ihisp,nbfisp,typei>', 1047*rak: & ish,ilow_sp,ihi_sp,nbfi_sp,typei 1048*rak: write(6,*)'<jsh,jlowsp,jhisp,nbfjsp,typej>', 1049*rak: & jsh,jlow_sp,jhi_sp,nbfj_sp,typej 1050*rak: write(6,*) ' ma broke 1' 1051*rak: if (.not.ma_verify_allocator_stuff()) stop ' ma broke 1' 1052*rak: write(6,*)' ' 1053c 1054 call int_1eov(basis,ish,basis,jsh,mscr,scr,mbuf,buf) 1055c 1056 count = 0 1057 do ibf = ilow,ihi 1058 do jbf = jlow, jhi 1059 count = count + 1 1060 value = buf(count) 1061 s(ibf,jbf) = value 1062 s(jbf,ibf) = value 1063 enddo 1064 enddo 1065#define NEW_WAY 1066#if defined(NEW_WAY) 1067 call spcart_tran1e(buf,scr, 1068 & nbfj,nbfi,typej,jgen, 1069 & nbfj_sp,nbfi_sp,typei,igen,.false.) 1070#else 1071*rak: write(6,*) ' ma broke 2' 1072*rak: if (.not.ma_verify_allocator_stuff()) stop ' ma broke 2' 1073c.... buf is now -- buf(nbfj,nbfi) 1074 write(6,*)' integral buffer cart,cart ' 1075 call output(buf,1,nbfj,1,nbfi,nbfj,nbfi,1) 1076 call spcart_a_s(buf,scr,nbfj,typei,.false.,.false.) 1077c.... scr is now -- scr(nbfj,nbfi_sp) 1078 write(6,*)' integral buffer cart,sph' 1079 call output(scr,1,nbfj,1,nbfi_sp,nbfj,nbfi_sp,1) 1080*rak: write(6,*) ' ma broke 3' 1081*rak: if (.not.ma_verify_allocator_stuff()) stop ' ma broke 3' 1082 call spcart_s_a(scr,buf,nbfi_sp,typej,.false.,.false.) 1083c.... buf is now -- buf(nbfj_sp,nbfi_sp) 1084 write(6,*)' integral buffer sph,sph' 1085 call output(buf,1,nbfj_sp,1,nbfi_sp,nbfj_sp,nbfi_sp,1) 1086*rak: write(6,*) ' ma broke 4' 1087*rak: if (.not.ma_verify_allocator_stuff()) stop ' ma broke 4' 1088#endif 1089 count = 0 1090 do ibf = ilow_sp,ihi_sp 1091 do jbf = jlow_sp, jhi_sp 1092 count = count + 1 1093 value = buf(count) 1094 s_sp(ibf,jbf) = value 1095 s_sp(jbf,ibf) = value 1096 enddo 1097 enddo 1098*rak: write(6,*) ' ma broke 5' 1099*rak: if (.not.ma_verify_allocator_stuff()) stop ' ma broke 5' 1100 enddo 1101 enddo 1102* write(6,*)' loops done ' 1103 write(6,*)' cartesian overlap matrix ' 1104 call output(s,1,nbf,1,nbf,nbf,nbf,1) 1105 write(6,*)' nbf ',nbf 1106 write(6,*)' nbf_sp ',nbf_sp 1107 write(6,*)' spherical overlap matrix ' 1108 call output(s_sp,1,nbf_sp,1,nbf_sp,nbf_sp,nbf_sp,1) 1109 write(6,*)' nbf ',nbf 1110 write(6,*)' nbf_sp ',nbf_sp 1111c 1112 count = 0 1113 do ii = 1,nbf_sp 1114 do jj = 1,(ii-1) 1115 if (abs(s_sp(ii,jj)).gt.1.0d-10) then 1116 count = count + 1 1117 write(6,'(a,i3,a,i3,a,1pd30.20,i6)') 1118 & 'non diag element.gt.1.0d-10 s_sp(', 1119 & ii,',',jj,') = ', 1120 & s_sp(ii,jj),count 1121 endif 1122 enddo 1123 enddo 1124 do ii = 1,nbf_sp 1125 if (abs(s_sp(ii,ii)-1.0d00).gt.1.0d-05) 1126 & write(6,'(a,i3,a,i3,a,1pd30.20)') 1127 & ' diagonal element.ne.1 s_sp(',ii,',',ii,') = ', 1128 & s_sp(ii,ii) 1129 enddo 1130*rak: do ii = 2, nbf_sp 1131*rak: do jj = 1,(ii-1) 1132*rak: write(6,'(a,2i5,2f15.8)') 1133*rak: & ' ratios ',ii,jj,(s_sp(ii,ii)/s_sp(jj,jj)), 1134*rak: & (s_sp(jj,jj)/s_sp(ii,ii)) 1135*rak: enddo 1136*rak: enddo 1137*rak: PI=2.0d00*acos(0.0d00) 1138*rak: write(6,*)pi,(pi-3.1415926535898D0) 1139*rak: do ii = 1, nbf_sp 1140*rak:*rak: write(6,*)' indx, val, val**2 ',ii,s_sp(ii,ii), 1141*rak:*rak: & (s_sp(ii,ii)*s_sp(ii,ii)) 1142*rak:*rak: write(6,*)' indx, val, 1/val ',ii,s_sp(ii,ii), 1143*rak:*rak: & (1.0d00/s_sp(ii,ii)) 1144*rak:*rak: write(6,*)' indx, val, sqrt ',ii,s_sp(ii,ii), 1145*rak:*rak: & (sqrt(s_sp(ii,ii))) 1146*rak:*rak: write(6,*)' indx, val, /sqrt(2) ',ii,s_sp(ii,ii), 1147*rak:*rak: & (s_sp(ii,ii)/sqrt(2.0d00)) 1148*rak:*rak: write(6,*)' indx, val, /sqrt(3) ',ii,s_sp(ii,ii), 1149*rak:*rak: & (s_sp(ii,ii)/sqrt(3.0d00)) 1150*rak: fact = 1.0d00 1151*rak: if (ii.ge.6.and.ii.le.10) fact = pi*8.0d00/5.0d00 1152*rak: if (ii.ge.11.and.ii.le.17) fact = pi*8.0d00/7.0d00 1153*rak: if (ii.ge.18.and.ii.le.26) fact = pi*8.0d00/9.0d00 1154*rak: write(6,*)' indx, val, *fact',ii,s_sp(ii,ii), 1155*rak: & (s_sp(ii,ii)*fact) 1156*rak: write(6,*)' indx, val, *pi',ii,s_sp(ii,ii), 1157*rak: & (s_sp(ii,ii)*pi) 1158*rak:*rak: write(6,*)' indx, val, /pi',ii,s_sp(ii,ii), 1159*rak:*rak: & (s_sp(ii,ii)/pi) 1160*rak: write(6,*)' ' 1161*rak: enddo 1162* if (.not.spcart_terminate()) stop ' dead sp term' 1163 end 1164 subroutine raktest_3ctr(rtdb) 1165 implicit none 1166#include "errquit.fh" 1167#include "mafdecls.fh" 1168#include "rtdb.fh" 1169#include "context.fh" 1170#include "geom.fh" 1171#include "bas.fh" 1172#include "nwc_const.fh" 1173#include "basP.fh" 1174#include "basdeclsP.fh" 1175#include "geomP.fh" 1176#include "geobasmapP.fh" 1177#include "bas_exndcf_dec.fh" 1178#include "bas_ibs_dec.fh" 1179c 1180 logical int_normalize 1181 external int_normalize 1182c 1183c test hf3 nai type routines 1184 integer rtdb 1185 integer geom,basis, basis_id 1186 integer nshell, memscr, membuf 1187 integer h_scr, k_scr, h_buf, k_buf 1188 integer ish, jsh, ucont 1189 integer li, i_prim, i_gen, i_iexp, i_icfp, i_cent, i_geom 1190 integer lj, j_prim, j_gen, j_iexp, j_icfp, j_cent, j_geom 1191 integer nint_out 1192 logical status 1193 character*255 mo_basis, geom_name 1194c 1195#include "bas_exndcf_sfn.fh" 1196#include "bas_ibs_sfn.fh" 1197c 1198 if (.not.context_rtdb_match(rtdb,'ao basis',mo_basis)) 1199 & mo_basis = 'ao basis' 1200 if (.not.context_rtdb_match(rtdb,'geometry',geom_name)) 1201 & geom_name = 'geometry' 1202c 1203 if(.not.geom_create(geom,geom_name))call errquit 1204 & ('raktest_3ctr: geom create error',911, GEOM_ERR) 1205 if(.not.bas_create(basis,mo_basis))call errquit 1206 & ('raktest_3ctr: basis create error',911, BASIS_ERR) 1207c 1208 if(.not.geom_rtdb_load(rtdb,geom,geom_name)) call errquit 1209 & ('raktest_3ctr: geom load ',911, RTDB_ERR) 1210 if(.not.bas_rtdb_load(rtdb,geom,basis,mo_basis)) call errquit 1211 & ('raktest_3ctr: basis load ',911, RTDB_ERR) 1212c 1213 basis_id = basis + BASIS_HANDLE_OFFSET 1214 nshell = ncont_tot_gb(basis_id) 1215 if (.not.int_normalize(rtdb,basis)) call errquit 1216 & ('raktest_3ctr: error normalizing ',911, INT_ERR) 1217c 1218 call int_init(rtdb,1,basis) 1219 memscr = 100 000 1220 membuf = 1000 1221 if (.not.ma_push_get(mt_dbl,memscr,' scratch ', 1222 & h_scr, k_scr)) call errquit 1223 & (' ma error 1',911, MA_ERR) 1224 if (.not.ma_push_get(mt_dbl,membuf,' buf ', 1225 & h_buf, k_buf)) call errquit 1226 & (' ma error 2',911, MA_ERR) 1227c 1228 do ish = 1,nshell 1229 do jsh = 1,ish 1230 write(6,*)' ============= shells <',ish,'|',jsh,'>', 1231 & '==================== start ==========' 1232 write(6,*)' ' 1233 1234 ucont = (sf_ibs_cn2ucn(ish,basis_id)) 1235 Li = infbs_cont(CONT_TYPE ,ucont,basis_id) 1236 i_prim = infbs_cont(CONT_NPRIM,ucont,basis_id) 1237 i_gen = infbs_cont(CONT_NGEN ,ucont,basis_id) 1238 i_iexp = infbs_cont(CONT_IEXP ,ucont,basis_id) 1239 i_icfp = infbs_cont(CONT_ICFP ,ucont,basis_id) 1240 i_cent = (sf_ibs_cn2ce(ish,basis_id)) 1241 i_geom = ibs_geom(basis_id) 1242c 1243 ucont = (sf_ibs_cn2ucn(jsh,basis_id)) 1244 Lj = infbs_cont(CONT_TYPE ,ucont,basis_id) 1245 j_prim = infbs_cont(CONT_NPRIM,ucont,basis_id) 1246 j_gen = infbs_cont(CONT_NGEN ,ucont,basis_id) 1247 j_iexp = infbs_cont(CONT_IEXP ,ucont,basis_id) 1248 j_icfp = infbs_cont(CONT_ICFP ,ucont,basis_id) 1249 j_cent = (sf_ibs_cn2ce(jsh,basis_id)) 1250 j_geom = ibs_geom(basis_id) 1251 1252 call hf1tmp( 1253 & coords(1,i_cent,i_geom), 1254 & dbl_mb(mb_exndcf(i_iexp,basis_id)), 1255 & dbl_mb(mb_exndcf(i_icfp,basis_id)), i_prim, i_gen, Li, 1256 & coords(1,j_cent,j_geom), 1257 & dbl_mb(mb_exndcf(j_iexp,basis_id)), 1258 & dbl_mb(mb_exndcf(j_icfp,basis_id)), j_prim, j_gen, Lj, 1259 & coords(1,1,i_geom),charge(1,i_geom),ncenter(i_geom), 1260 & dbl_mb(k_scr),dbl_mb(k_scr),dbl_mb(k_buf),membuf, 1261 & .false., .false., .true., .false., .false., 1262 & dbl_mb(k_scr), memscr) 1263 write(6,*)' i = c center ' 1264 call hf3pot( 1265 & coords(1,i_cent,i_geom), 1266 & dbl_mb(mb_exndcf(i_iexp,basis_id)), 1267 & dbl_mb(mb_exndcf(i_icfp,basis_id)), i_prim, i_gen, Li, 1268 & coords(1,j_cent,j_geom), 1269 & dbl_mb(mb_exndcf(j_iexp,basis_id)), 1270 & dbl_mb(mb_exndcf(j_icfp,basis_id)), j_prim, j_gen, Lj, 1271 & coords(1,i_cent,i_geom),0.0d00, 1.0d00, 1, 1, 0, 1272 & dbl_mb(k_buf), membuf, nint_out, .false., 1273 & dbl_mb(k_scr), memscr) 1274 write(6,*)' j = c center ' 1275 call hf3pot( 1276 & coords(1,i_cent,i_geom), 1277 & dbl_mb(mb_exndcf(i_iexp,basis_id)), 1278 & dbl_mb(mb_exndcf(i_icfp,basis_id)), i_prim, i_gen, Li, 1279 & coords(1,j_cent,j_geom), 1280 & dbl_mb(mb_exndcf(j_iexp,basis_id)), 1281 & dbl_mb(mb_exndcf(j_icfp,basis_id)), j_prim, j_gen, Lj, 1282 & coords(1,j_cent,j_geom),0.0d00, 1.0d00, 1, 1, 0, 1283 & dbl_mb(k_buf), membuf, nint_out, .false., 1284 & dbl_mb(k_scr), memscr) 1285 write(6,*)' i = c center swap' 1286 call hf3pot( 1287 & coords(1,i_cent,i_geom),0.0d00, 1.0d00, 1, 1, 0, 1288 & coords(1,j_cent,j_geom), 1289 & dbl_mb(mb_exndcf(j_iexp,basis_id)), 1290 & dbl_mb(mb_exndcf(j_icfp,basis_id)), j_prim, j_gen, Lj, 1291 & coords(1,i_cent,i_geom), 1292 & dbl_mb(mb_exndcf(i_iexp,basis_id)), 1293 & dbl_mb(mb_exndcf(i_icfp,basis_id)), i_prim, i_gen, Li, 1294 & dbl_mb(k_buf), membuf, nint_out, .false., 1295 & dbl_mb(k_scr), memscr) 1296 write(6,*)' j = c center swap' 1297 call hf3pot( 1298 & coords(1,i_cent,i_geom), 1299 & dbl_mb(mb_exndcf(i_iexp,basis_id)), 1300 & dbl_mb(mb_exndcf(i_icfp,basis_id)), i_prim, i_gen, Li, 1301 & coords(1,j_cent,j_geom),0.0d00, 1.0d00, 1, 1, 0, 1302 & coords(1,j_cent,j_geom), 1303 & dbl_mb(mb_exndcf(j_iexp,basis_id)), 1304 & dbl_mb(mb_exndcf(j_icfp,basis_id)), j_prim, j_gen, Lj, 1305 & dbl_mb(k_buf), membuf, nint_out, .false., 1306 & dbl_mb(k_scr), memscr) 1307 write(6,*)' ============= shells <',ish,'|',jsh,'>', 1308 & '==================== end ==========' 1309 write(6,*)' ' 1310 enddo 1311 enddo 1312c 1313 call int_terminate() 1314 status = ma_pop_stack(h_buf) 1315 status = status.and.ma_pop_stack(h_scr) 1316 if (.not.status) call errquit('pop failed',911, MA_ERR) 1317 status = bas_destroy(basis) 1318 status = status.and.geom_destroy(geom) 1319 if (.not.status) call errquit('b/g destroy failed',911, MA_ERR) 1320 end 1321 subroutine raktest_stpr(rtdb) 1322 implicit none 1323#include "errquit.fh" 1324#include "rtdb.fh" 1325#include "geom.fh" 1326#include "mafdecls.fh" 1327#include "global.fh" 1328 integer rtdb 1329 integer stpr_walk 1330 external stpr_walk 1331c 1332 integer geom 1333 integer nat 1334 integer k_grad, h_grad 1335 logical flag 1336 double precision energy 1337c 1338 if (ga_nodeid().eq.0) then 1339 if (.not. geom_create(geom, 'geometry')) 1340 & call errquit('raktest_stpr: geom_create?', 911, GEOM_ERR) 1341 if (.not. geom_rtdb_load(rtdb, geom, 'geometry')) 1342 & call errquit('raktest_stpr: no geometry ', 911, RTDB_ERR) 1343c get number of atoms = nat 1344 if (.not. geom_ncent(geom,nat)) 1345 & call errquit('raktest_stpr: geom_ncent?',911, GEOM_ERR) 1346 if (.not. geom_destroy(geom)) 1347 & call errquit('raktest_stpr: geom_destroy?',911, GEOM_ERR) 1348 if (.not. 1349 & MA_Push_Get(MT_DBL,(3*nat),'stpr fake gradient', 1350 & h_grad,k_grad)) 1351 & call errquit 1352 & ('raktest_stpr: allocation for gradient failed?',911, 1353 & MA_ERR) 1354 dbl_mb((k_grad+ 0)) = -0.008697D00 1355 dbl_mb((k_grad+ 1)) = -0.004076D00 1356 dbl_mb((k_grad+ 2)) = 0.004591D00 1357 dbl_mb((k_grad+ 3)) = 0.000926D00 1358 dbl_mb((k_grad+ 4)) = 0.000922D00 1359 dbl_mb((k_grad+ 5)) = -0.001206D00 1360 dbl_mb((k_grad+ 6)) = 0.000639D00 1361 dbl_mb((k_grad+ 7)) = 0.000263D00 1362 dbl_mb((k_grad+ 8)) = -0.003188D00 1363 dbl_mb((k_grad+ 9)) = 0.007578D00 1364 dbl_mb((k_grad+ 10)) = -0.002686D00 1365 dbl_mb((k_grad+ 11)) = 0.005190D00 1366 dbl_mb((k_grad+ 12)) = -0.000872D00 1367 dbl_mb((k_grad+ 13)) = 0.000692D00 1368 dbl_mb((k_grad+ 14)) = -0.001103D00 1369 dbl_mb((k_grad+ 15)) = -0.000782D00 1370 dbl_mb((k_grad+ 16)) = 0.000717D00 1371 dbl_mb((k_grad+ 17)) = -0.001393D00 1372 dbl_mb((k_grad+ 18)) = -0.008078D00 1373 dbl_mb((k_grad+ 19)) = 0.006782D00 1374 dbl_mb((k_grad+ 20)) = 0.000074D00 1375 dbl_mb((k_grad+ 21)) = 0.000549D00 1376 dbl_mb((k_grad+ 22)) = -0.000865D00 1377 dbl_mb((k_grad+ 23)) = -0.002881D00 1378 dbl_mb((k_grad+ 24)) = 0.006377D00 1379 dbl_mb((k_grad+ 25)) = 0.001920D00 1380 dbl_mb((k_grad+ 26)) = 0.000909D00 1381 dbl_mb((k_grad+ 27)) = -0.000695D00 1382 dbl_mb((k_grad+ 28)) = -0.000362D00 1383 dbl_mb((k_grad+ 29)) = -0.000411D00 1384 dbl_mb((k_grad+ 30)) = -0.001585D00 1385 dbl_mb((k_grad+ 31)) = -0.000732D00 1386 dbl_mb((k_grad+ 32)) = -0.000523D00 1387 dbl_mb((k_grad+ 33)) = -0.001062D00 1388 dbl_mb((k_grad+ 34)) = -0.000824D00 1389 dbl_mb((k_grad+ 35)) = -0.000303D00 1390 dbl_mb((k_grad+ 36)) = -0.022866D00 1391 dbl_mb((k_grad+ 37)) = 0.011419D00 1392 dbl_mb((k_grad+ 38)) = -0.010230D00 1393 dbl_mb((k_grad+ 39)) = 0.013343D00 1394 dbl_mb((k_grad+ 40)) = -0.017703D00 1395 dbl_mb((k_grad+ 41)) = -0.018748D00 1396 dbl_mb((k_grad+ 42)) = 0.002823D00 1397 dbl_mb((k_grad+ 43)) = 0.000704D00 1398 dbl_mb((k_grad+ 44)) = 0.000638D00 1399 dbl_mb((k_grad+ 45)) = 0.001783D00 1400 dbl_mb((k_grad+ 46)) = 0.000141D00 1401 dbl_mb((k_grad+ 47)) = 0.000499D00 1402 dbl_mb((k_grad+ 48)) = 0.003863D00 1403 dbl_mb((k_grad+ 49)) = 0.007955D00 1404 dbl_mb((k_grad+ 50)) = -0.003668D00 1405 dbl_mb((k_grad+ 51)) = -0.000253D00 1406 dbl_mb((k_grad+ 52)) = -0.000335D00 1407 dbl_mb((k_grad+ 53)) = 0.003427D00 1408 dbl_mb((k_grad+ 54)) = 0.003510D00 1409 dbl_mb((k_grad+ 55)) = -0.004396D00 1410 dbl_mb((k_grad+ 56)) = -0.000280D00 1411 dbl_mb((k_grad+ 57)) = -0.001050D00 1412 dbl_mb((k_grad+ 58)) = 0.000598D00 1413 dbl_mb((k_grad+ 59)) = 0.000407D00 1414 dbl_mb((k_grad+ 60)) = -0.001190D00 1415 dbl_mb((k_grad+ 61)) = 0.000896D00 1416 dbl_mb((k_grad+ 62)) = 0.000453D00 1417 dbl_mb((k_grad+ 63)) = -0.000536D00 1418 dbl_mb((k_grad+ 64)) = 0.000599D00 1419 dbl_mb((k_grad+ 65)) = 0.000487D00 1420 dbl_mb((k_grad+ 66)) = 0.008229D00 1421 dbl_mb((k_grad+ 67)) = 0.006460D00 1422 dbl_mb((k_grad+ 68)) = -0.000069D00 1423 dbl_mb((k_grad+ 69)) = -0.000490D00 1424 dbl_mb((k_grad+ 70)) = -0.000613D00 1425 dbl_mb((k_grad+ 71)) = -0.003029D00 1426 dbl_mb((k_grad+ 72)) = -0.006348D00 1427 dbl_mb((k_grad+ 73)) = 0.002016D00 1428 dbl_mb((k_grad+ 74)) = 0.000844D00 1429 dbl_mb((k_grad+ 75)) = 0.001549D00 1430 dbl_mb((k_grad+ 76)) = -0.000898D00 1431 dbl_mb((k_grad+ 77)) = -0.000570D00 1432 dbl_mb((k_grad+ 78)) = 0.000697D00 1433 dbl_mb((k_grad+ 79)) = -0.000361D00 1434 dbl_mb((k_grad+ 80)) = -0.000393D00 1435 dbl_mb((k_grad+ 81)) = 0.001078D00 1436 dbl_mb((k_grad+ 82)) = -0.000863D00 1437 dbl_mb((k_grad+ 83)) = -0.000271D00 1438 dbl_mb((k_grad+ 84)) = -0.007499D00 1439 dbl_mb((k_grad+ 85)) = -0.002478D00 1440 dbl_mb((k_grad+ 86)) = 0.005769D00 1441 dbl_mb((k_grad+ 87)) = 0.000808D00 1442 dbl_mb((k_grad+ 88)) = 0.000701D00 1443 dbl_mb((k_grad+ 89)) = -0.001448D00 1444 dbl_mb((k_grad+ 90)) = 0.000827D00 1445 dbl_mb((k_grad+ 91)) = 0.000750D00 1446 dbl_mb((k_grad+ 92)) = -0.001107D00 1447 dbl_mb((k_grad+ 93)) = 0.007929D00 1448 dbl_mb((k_grad+ 94)) = -0.004455D00 1449 dbl_mb((k_grad+ 95)) = 0.003384D00 1450 dbl_mb((k_grad+ 96)) = -0.000671D00 1451 dbl_mb((k_grad+ 97)) = 0.000441D00 1452 dbl_mb((k_grad+ 98)) = -0.003097D00 1453 dbl_mb((k_grad+ 99)) = -0.000863D00 1454 dbl_mb((k_grad+ 100)) = 0.000987D00 1455 dbl_mb((k_grad+ 101)) = -0.001353D00 1456 dbl_mb((k_grad+ 102)) = 0.007929D00 1457 dbl_mb((k_grad+ 103)) = 0.004455D00 1458 dbl_mb((k_grad+ 104)) = -0.003384D00 1459 dbl_mb((k_grad+ 105)) = -0.000671D00 1460 dbl_mb((k_grad+ 106)) = -0.000441D00 1461 dbl_mb((k_grad+ 107)) = 0.003097D00 1462 dbl_mb((k_grad+ 108)) = -0.000863D00 1463 dbl_mb((k_grad+ 109)) = -0.000987D00 1464 dbl_mb((k_grad+ 110)) = 0.001353D00 1465 dbl_mb((k_grad+ 111)) = -0.007499D00 1466 dbl_mb((k_grad+ 112)) = 0.002477D00 1467 dbl_mb((k_grad+ 113)) = -0.005768D00 1468 dbl_mb((k_grad+ 114)) = 0.000808D00 1469 dbl_mb((k_grad+ 115)) = -0.000701D00 1470 dbl_mb((k_grad+ 116)) = 0.001448D00 1471 dbl_mb((k_grad+ 117)) = 0.000827D00 1472 dbl_mb((k_grad+ 118)) = -0.000751D00 1473 dbl_mb((k_grad+ 119)) = 0.001107D00 1474 dbl_mb((k_grad+ 120)) = 0.008230D00 1475 dbl_mb((k_grad+ 121)) = -0.006460D00 1476 dbl_mb((k_grad+ 122)) = 0.000068D00 1477 dbl_mb((k_grad+ 123)) = -0.000491D00 1478 dbl_mb((k_grad+ 124)) = 0.000613D00 1479 dbl_mb((k_grad+ 125)) = 0.003029D00 1480 dbl_mb((k_grad+ 126)) = -0.006349D00 1481 dbl_mb((k_grad+ 127)) = -0.002015D00 1482 dbl_mb((k_grad+ 128)) = -0.000844D00 1483 dbl_mb((k_grad+ 129)) = 0.001078D00 1484 dbl_mb((k_grad+ 130)) = 0.000863D00 1485 dbl_mb((k_grad+ 131)) = 0.000270D00 1486 dbl_mb((k_grad+ 132)) = 0.001549D00 1487 dbl_mb((k_grad+ 133)) = 0.000898D00 1488 dbl_mb((k_grad+ 134)) = 0.000570D00 1489 dbl_mb((k_grad+ 135)) = 0.000697D00 1490 dbl_mb((k_grad+ 136)) = 0.000361D00 1491 dbl_mb((k_grad+ 137)) = 0.000393D00 1492 dbl_mb((k_grad+ 138)) = 0.003863D00 1493 dbl_mb((k_grad+ 139)) = -0.007955D00 1494 dbl_mb((k_grad+ 140)) = 0.003667D00 1495 dbl_mb((k_grad+ 141)) = -0.000253D00 1496 dbl_mb((k_grad+ 142)) = 0.000335D00 1497 dbl_mb((k_grad+ 143)) = -0.003427D00 1498 dbl_mb((k_grad+ 144)) = 0.003510D00 1499 dbl_mb((k_grad+ 145)) = 0.004397D00 1500 dbl_mb((k_grad+ 146)) = 0.000279D00 1501 dbl_mb((k_grad+ 147)) = -0.001190D00 1502 dbl_mb((k_grad+ 148)) = -0.000897D00 1503 dbl_mb((k_grad+ 149)) = -0.000453D00 1504 dbl_mb((k_grad+ 150)) = -0.000536D00 1505 dbl_mb((k_grad+ 151)) = -0.000599D00 1506 dbl_mb((k_grad+ 152)) = -0.000487D00 1507 dbl_mb((k_grad+ 153)) = -0.001050D00 1508 dbl_mb((k_grad+ 154)) = -0.000598D00 1509 dbl_mb((k_grad+ 155)) = -0.000407D00 1510 dbl_mb((k_grad+ 156)) = 0.013343D00 1511 dbl_mb((k_grad+ 157)) = 0.017702D00 1512 dbl_mb((k_grad+ 158)) = 0.018749D00 1513 dbl_mb((k_grad+ 159)) = 0.002823D00 1514 dbl_mb((k_grad+ 160)) = -0.000703D00 1515 dbl_mb((k_grad+ 161)) = -0.000639D00 1516 dbl_mb((k_grad+ 162)) = 0.001783D00 1517 dbl_mb((k_grad+ 163)) = -0.000141D00 1518 dbl_mb((k_grad+ 164)) = -0.000499D00 1519 dbl_mb((k_grad+ 165)) = -0.022865D00 1520 dbl_mb((k_grad+ 166)) = -0.011418D00 1521 dbl_mb((k_grad+ 167)) = 0.010231D00 1522 dbl_mb((k_grad+ 168)) = 0.000385D00 1523 dbl_mb((k_grad+ 169)) = 0.023722D00 1524 dbl_mb((k_grad+ 170)) = 0.016418D00 1525 dbl_mb((k_grad+ 171)) = 0.022263D00 1526 dbl_mb((k_grad+ 172)) = -0.013938D00 1527 dbl_mb((k_grad+ 173)) = -0.000893D00 1528 dbl_mb((k_grad+ 174)) = 0.000883D00 1529 dbl_mb((k_grad+ 175)) = -0.000287D00 1530 dbl_mb((k_grad+ 176)) = -0.002376D00 1531 dbl_mb((k_grad+ 177)) = -0.022844D00 1532 dbl_mb((k_grad+ 178)) = -0.010601D00 1533 dbl_mb((k_grad+ 179)) = -0.000400D00 1534 dbl_mb((k_grad+ 180)) = -0.000446D00 1535 dbl_mb((k_grad+ 181)) = -0.001743D00 1536 dbl_mb((k_grad+ 182)) = -0.001937D00 1537 dbl_mb((k_grad+ 183)) = 0.020689D00 1538 dbl_mb((k_grad+ 184)) = -0.010232D00 1539 dbl_mb((k_grad+ 185)) = 0.017869D00 1540 dbl_mb((k_grad+ 186)) = -0.022845D00 1541 dbl_mb((k_grad+ 187)) = 0.010601D00 1542 dbl_mb((k_grad+ 188)) = 0.000400D00 1543 dbl_mb((k_grad+ 189)) = -0.000446D00 1544 dbl_mb((k_grad+ 190)) = 0.001743D00 1545 dbl_mb((k_grad+ 191)) = 0.001937D00 1546 dbl_mb((k_grad+ 192)) = 0.022264D00 1547 dbl_mb((k_grad+ 193)) = 0.013938D00 1548 dbl_mb((k_grad+ 194)) = 0.000892D00 1549 dbl_mb((k_grad+ 195)) = 0.000883D00 1550 dbl_mb((k_grad+ 196)) = 0.000287D00 1551 dbl_mb((k_grad+ 197)) = 0.002376D00 1552 dbl_mb((k_grad+ 198)) = 0.020690D00 1553 dbl_mb((k_grad+ 199)) = 0.010232D00 1554 dbl_mb((k_grad+ 200)) = -0.017869D00 1555 dbl_mb((k_grad+ 201)) = 0.002430D00 1556 dbl_mb((k_grad+ 202)) = -0.022402D00 1557 dbl_mb((k_grad+ 203)) = 0.003835D00 1558 dbl_mb((k_grad+ 204)) = -0.001212D00 1559 dbl_mb((k_grad+ 205)) = -0.000669D00 1560 dbl_mb((k_grad+ 206)) = 0.002178D00 1561 dbl_mb((k_grad+ 207)) = -0.002379D00 1562 dbl_mb((k_grad+ 208)) = 0.025844D00 1563 dbl_mb((k_grad+ 209)) = -0.005825D00 1564 dbl_mb((k_grad+ 210)) = 0.001510D00 1565 dbl_mb((k_grad+ 211)) = 0.000491D00 1566 dbl_mb((k_grad+ 212)) = -0.001989D00 1567 dbl_mb((k_grad+ 213)) = 0.000386D00 1568 dbl_mb((k_grad+ 214)) = -0.023721D00 1569 dbl_mb((k_grad+ 215)) = -0.016418D00 1570 dbl_mb((k_grad+ 216)) = 0.022354D00 1571 dbl_mb((k_grad+ 217)) = -0.009889D00 1572 dbl_mb((k_grad+ 218)) = -0.001034D00 1573 dbl_mb((k_grad+ 219)) = 0.000445D00 1574 dbl_mb((k_grad+ 220)) = -0.001737D00 1575 dbl_mb((k_grad+ 221)) = -0.001918D00 1576 dbl_mb((k_grad+ 222)) = -0.023384D00 1577 dbl_mb((k_grad+ 223)) = -0.006755D00 1578 dbl_mb((k_grad+ 224)) = 0.020970D00 1579 dbl_mb((k_grad+ 225)) = -0.000504D00 1580 dbl_mb((k_grad+ 226)) = 0.001171D00 1581 dbl_mb((k_grad+ 227)) = 0.000301D00 1582 dbl_mb((k_grad+ 228)) = -0.019164D00 1583 dbl_mb((k_grad+ 229)) = 0.008952D00 1584 dbl_mb((k_grad+ 230)) = -0.018541D00 1585 dbl_mb((k_grad+ 231)) = -0.023385D00 1586 dbl_mb((k_grad+ 232)) = 0.006755D00 1587 dbl_mb((k_grad+ 233)) = -0.020970D00 1588 dbl_mb((k_grad+ 234)) = -0.000503D00 1589 dbl_mb((k_grad+ 235)) = -0.001171D00 1590 dbl_mb((k_grad+ 236)) = -0.000301D00 1591 dbl_mb((k_grad+ 237)) = 0.022354D00 1592 dbl_mb((k_grad+ 238)) = 0.009888D00 1593 dbl_mb((k_grad+ 239)) = 0.001034D00 1594 dbl_mb((k_grad+ 240)) = 0.000445D00 1595 dbl_mb((k_grad+ 241)) = 0.001737D00 1596 dbl_mb((k_grad+ 242)) = 0.001918D00 1597 dbl_mb((k_grad+ 243)) = -0.019164D00 1598 dbl_mb((k_grad+ 244)) = -0.008951D00 1599 dbl_mb((k_grad+ 245)) = 0.018541D00 1600 dbl_mb((k_grad+ 246)) = -0.002378D00 1601 dbl_mb((k_grad+ 247)) = -0.025843D00 1602 dbl_mb((k_grad+ 248)) = 0.005825D00 1603 dbl_mb((k_grad+ 249)) = 0.001510D00 1604 dbl_mb((k_grad+ 250)) = -0.000491D00 1605 dbl_mb((k_grad+ 251)) = 0.001990D00 1606 dbl_mb((k_grad+ 252)) = 0.002430D00 1607 dbl_mb((k_grad+ 253)) = 0.022402D00 1608 dbl_mb((k_grad+ 254)) = -0.003835D00 1609 dbl_mb((k_grad+ 255)) = -0.001213D00 1610 dbl_mb((k_grad+ 256)) = 0.000669D00 1611 dbl_mb((k_grad+ 257)) = -0.002178D00 1612 dbl_mb((k_grad+ 258)) = -0.008078D00 1613 dbl_mb((k_grad+ 259)) = -0.006781D00 1614 dbl_mb((k_grad+ 260)) = -0.000074D00 1615 dbl_mb((k_grad+ 261)) = 0.000549D00 1616 dbl_mb((k_grad+ 262)) = 0.000865D00 1617 dbl_mb((k_grad+ 263)) = 0.002881D00 1618 dbl_mb((k_grad+ 264)) = 0.006378D00 1619 dbl_mb((k_grad+ 265)) = -0.001920D00 1620 dbl_mb((k_grad+ 266)) = -0.000909D00 1621 dbl_mb((k_grad+ 267)) = -0.000695D00 1622 dbl_mb((k_grad+ 268)) = 0.000362D00 1623 dbl_mb((k_grad+ 269)) = 0.000411D00 1624 dbl_mb((k_grad+ 270)) = -0.001585D00 1625 dbl_mb((k_grad+ 271)) = 0.000732D00 1626 dbl_mb((k_grad+ 272)) = 0.000524D00 1627 dbl_mb((k_grad+ 273)) = -0.001062D00 1628 dbl_mb((k_grad+ 274)) = 0.000824D00 1629 dbl_mb((k_grad+ 275)) = 0.000303D00 1630 dbl_mb((k_grad+ 276)) = 0.007577D00 1631 dbl_mb((k_grad+ 277)) = 0.002686D00 1632 dbl_mb((k_grad+ 278)) = -0.005190D00 1633 dbl_mb((k_grad+ 279)) = -0.000782D00 1634 dbl_mb((k_grad+ 280)) = -0.000716D00 1635 dbl_mb((k_grad+ 281)) = 0.001394D00 1636 dbl_mb((k_grad+ 282)) = -0.000872D00 1637 dbl_mb((k_grad+ 283)) = -0.000692D00 1638 dbl_mb((k_grad+ 284)) = 0.001103D00 1639 dbl_mb((k_grad+ 285)) = -0.008696D00 1640 dbl_mb((k_grad+ 286)) = 0.004077D00 1641 dbl_mb((k_grad+ 287)) = -0.004592D00 1642 dbl_mb((k_grad+ 288)) = 0.000639D00 1643 dbl_mb((k_grad+ 289)) = -0.000263D00 1644 dbl_mb((k_grad+ 290)) = 0.003187D00 1645 dbl_mb((k_grad+ 291)) = 0.000926D00 1646 dbl_mb((k_grad+ 292)) = -0.000923D00 1647 dbl_mb((k_grad+ 293)) = 0.001207D00 1648 energy = -317.5656589D00 1649 flag = .true. 1650c put scf:converged = logical true 1651 if (.not. rtdb_put(rtdb,'scf:converged', MT_LOG, 1, flag)) 1652 & call errquit 1653 & ('raktest_stpr: failed to read converged in rtdb', 911, 1654 & RTDB_ERR) 1655c put scf:energy = real value 1656 if (.not. rtdb_put(rtdb,'scf:energy', MT_DBL, 1, energy)) 1657 & call errquit 1658 & ('raktest_stpr: failed to read energy in rtdb', 911, 1659 & RTDB_ERR) 1660c put scf:gradient = 3*nat reals 1661 if (.not. rtdb_put(rtdb, 'scf:gradient', MT_DBL, 1662 & (3*nat),dbl_mb(k_grad))) 1663 & call errquit 1664 & ('raktest_stpr: reading gradients failed',911, RTDB_ERR) 1665c free memory 1666 if (.not. ma_pop_stack(h_grad)) 1667 & call errquit('raktest_stpr: pop failed',911, MA_ERR) 1668 endif 1669 call ga_sync() 1670 call ga_sync() 1671 call ga_sync() 1672 if (stpr_walk(rtdb).eq.1) then 1673 write(6,*)' walk converged' 1674 else 1675 write(6,*)' walk NOT converged' 1676 endif 1677 end 1678 subroutine raktest_initd(rtdb) 1679c raktest = 4 1680 implicit none 1681#include "errquit.fh" 1682#include "rtdb.fh" 1683#include "geom.fh" 1684#include "bas.fh" 1685#include "nwc_const.fh" 1686#include "basP.fh" 1687#include "mafdecls.fh" 1688c 1689 integer rtdb 1690 integer geom 1691 integer mx1e, mxg, mxs1, mxs2 1692c 1693 logical status 1694c 1695 integer nbas, bases(6), mynbas 1696c 1697 if (.not.bas_rtdb_in(rtdb)) 1698 & call errquit('raktest4: error loading known basis sets',911, 1699 & BASIS_ERR) 1700c 1701 write(6,*)' number of basis sets in rtdb ',nbasis_rtdb 1702c 1703 do 00100 nbas = 1,nbasis_rtdb 1704 write(6,*)' basis ',nbas,' is ',bs_names_rtdb(nbas) 170500100 continue 1706c 1707 if (.not.geom_create(geom,'geometry')) 1708 & call errquit('raktest4: geom_create failed',911, GEOM_ERR) 1709 if (.not.geom_rtdb_load(rtdb,geom,'geometry')) 1710 & call errquit('raktest4: geom_create failed',911, RTDB_ERR) 1711c 1712 mynbas = 0 1713 do 00200 nbas = 1,nbasis_rtdb 1714 if (bs_names_rtdb(nbas).ne.'ecp basis') then 1715 mynbas = mynbas + 1 1716 if(.not.bas_create(bases(mynbas),bs_names_rtdb(mynbas))) 1717 & call errquit('raktest4: bas_create choked',911, BASIS_ERR) 1718 if(.not. 1719 & bas_rtdb_load 1720 & (rtdb,geom,bases(mynbas),bs_names_rtdb(mynbas))) 1721 & call errquit('raktest4: bas_rtdb_load choked',911, 1722 & RTDB_ERR) 1723 status = bas_print(bases(mynbas)) 1724 status = gbs_map_print(bases(mynbas)) 1725 endif 172600200 continue 1727 call intd_init(rtdb,mynbas,bases) 1728c 1729 call int_mem(mx1e, mxg, mxs1, mxs2) 1730 write(6,*)' one electron buffer size :',mx1e 1731 write(6,*)' two electron buffer size :',mxg 1732 write(6,*)' one electron scratch buffer size:',mxs1 1733 write(6,*)' two electron scratch buffer size:',mxs2 1734c 1735 call int_mem_print() 1736c 1737 do nbas = 1,mynbas 1738 if (.not.bas_destroy(bases(nbas))) call errquit 1739 & ('raktest_initd: bas_destroy failed',911, BASIS_ERR) 1740 enddo 1741 if (.not.geom_destroy(geom)) call errquit 1742 & ('raktest_initd: _destroy failed',911, GEOM_ERR) 1743c 1744 call intd_terminate() 1745c 1746 end 1747 subroutine raktest_init(rtdb) 1748c raktest = 5 1749 implicit none 1750#include "errquit.fh" 1751#include "rtdb.fh" 1752#include "geom.fh" 1753#include "bas.fh" 1754#include "nwc_const.fh" 1755#include "basP.fh" 1756#include "mafdecls.fh" 1757c 1758 integer rtdb 1759 integer geom 1760 integer mx1e, mxg, mxs1, mxs2 1761c 1762 logical status 1763c 1764 integer nbas, bases(6) 1765c 1766 if (.not.bas_rtdb_in(rtdb)) 1767 & call errquit('raktest4: error loading known basis sets',911, 1768 & BASIS_ERR) 1769c 1770 write(6,*)' number of basis sets in rtdb ',nbasis_rtdb 1771c 1772 do 00100 nbas = 1,nbasis_rtdb 1773 write(6,*)' basis ',nbas,' is ',bs_names_rtdb(nbas) 177400100 continue 1775c 1776 if (.not.geom_create(geom,'geometry')) 1777 & call errquit('raktest4: geom_create failed',911, GEOM_ERR) 1778 if (.not.geom_rtdb_load(rtdb,geom,'geometry')) 1779 & call errquit('raktest4: geom_create failed',911, RTDB_ERR) 1780c 1781 do 00200 nbas = 1,nbasis_rtdb 1782 if(.not.bas_create(bases(nbas),bs_names_rtdb(nbas))) 1783 & call errquit('raktest4: bas_create choked',911, BASIS_ERR) 1784 if(.not. 1785 & bas_rtdb_load 1786 & (rtdb,geom,bases(nbas),bs_names_rtdb(nbas))) 1787 & call errquit('raktest4: bas_rtdb_load choked',911, RTDB_ERR) 1788 status = bas_print(bases(nbas)) 1789 status = gbs_map_print(bases(nbas)) 179000200 continue 1791 call int_init(rtdb,nbasis_rtdb,bases) 1792c 1793 call int_mem(mx1e, mxg, mxs1, mxs2) 1794 write(6,*)' one electron buffer size :',mx1e 1795 write(6,*)' two electron buffer size :',mxg 1796 write(6,*)' one electron scratch buffer size:',mxs1 1797 write(6,*)' two electron scratch buffer size:',mxs2 1798c 1799 call int_mem_print() 1800c 1801 do nbas = 1,nbasis_rtdb 1802 if (.not.bas_destroy(bases(nbas))) call errquit 1803 & ('raktest_init: bas_destroy failed',911, BASIS_ERR) 1804 enddo 1805 if (.not.geom_destroy(geom)) call errquit 1806 & ('raktest_init: _destroy failed',911, GEOM_ERR) 1807c 1808 call int_terminate() 1809c 1810 end 1811 1812 Subroutine hf1mkrtmp(Axyz,Bxyz,Cxyz,zan,ncenters, 1813 & alpha,Pxyz,RS,PC,ff,R, 1814 & R0,R0C,IJK,NPP,Lp,Lp3,CENTER) 1815 1816 Implicit real*8 (a-h,o-z) 1817 Implicit integer (i-n) 1818 1819 Logical CENTER 1820 1821 Parameter (PI=3.1415926535898D0,PI4=4.D0/PI) 1822 1823c--> Cartesian Coordinates 1824 1825 Dimension Axyz(3),Bxyz(3) 1826 1827c--> Nuclear Cartesian Coordinates & Charges 1828 1829 Dimension Cxyz(3,ncenters),zan(ncenters) 1830 1831c--> Exponents 1832 1833 Dimension alpha(2,NPP) 1834 1835c--> Auxiliary Function Integrals & Index 1836 1837 Dimension R0(NPP,Lp3),R0C(ncenters,NPP,Lp3),IJK(0:Lp,0:Lp,0:Lp) 1838 1839c--> Scratch Space 1840 1841 Dimension Pxyz(3,NPP),RS(NPP),PC(NPP,3),ff(2,NPP),R(NPP,0:Lp,Lp3) 1842 1843c 1844c Define the auxiliary function integrals necessary to compute the nuclear 1845c attraction integrals (NAIs). These integrals are scaled by an appropriate 1846c factor, RS, defined as 1847c 1848c / a + b \ 1/2 1849c RS = | ------- | 1850c \ PI/4 / 1851c 1852c The scale factor for the Hermite expansion coefficients is assumed to be 1853c 1854c / PI \ 3/2 / a b __2 \ 1855c ES = | ------- | EXP| - ----- AB | 1856c \ a + b / \ a + b / 1857c 1858c Therefore, 1859c 1860c 2 PI / a b __2 \ 1861c ES RS = ------- EXP| - ----- AB | 1862c a + b \ a + b / 1863c 1864c****************************************************************************** 1865 1866 do 100 mp = 1,NPP 1867 do 100 j = 1,Lp3 1868 R0(mp,j) = 0.D0 1869 100 continue 1870 1871 do 110 mp = 1,NPP 1872 1873c Define the center "P". 1874 1875 a = alpha(1,mp) 1876 b = alpha(2,mp) 1877 1878 f1 = a/(a+b) 1879 f2 = b/(a+b) 1880 1881 Pxyz(1,mp) = f1*Axyz(1) + f2*Bxyz(1) 1882 Pxyz(2,mp) = f1*Axyz(2) + f2*Bxyz(2) 1883 Pxyz(3,mp) = f1*Axyz(3) + f2*Bxyz(3) 1884 1885c Define the scaling factor. 1886 1887 RS(mp) = sqrt((a+b)*PI4) 1888 1889 110 continue 1890 1891c Sum over all centers. 1892 1893 do 150 ic = 1,ncenters 1894 1895c Define factors necessary to compute incomplete gamma function and the 1896c auxiliary functions. 1897 1898 do 120 m = 1,NPP 1899 1900 alpha_t = alpha(1,m) + alpha(2,m) 1901 1902 ff(1,m) = RS(m) 1903 ff(2,m) = -2.D0*alpha_t 1904 1905 PCx = Pxyz(1,m) - Cxyz(1,ic) 1906 PCy = Pxyz(2,m) - Cxyz(2,ic) 1907 PCz = Pxyz(3,m) - Cxyz(3,ic) 1908 1909 R(m,0,1) = alpha_t*(PCx**2 + PCy**2 + PCz**2) 1910 1911 PC(m,1) = PCx 1912 PC(m,2) = PCy 1913 PC(m,3) = PCz 1914 1915 120 continue 1916 1917c Evaluate the incomplete gamma function. 1918 1919 call igamma(R,NPP,Lp) 1920 1921c Define the initial auxiliary functions (i.e., R000j, j=1,Lr). 1922 1923 do 135 j = 0,Lp 1924 do 130 m = 1,NPP 1925 R(m,j,1) = ff(1,m)*R(m,j,1) 1926 ff(1,m) = ff(1,m)*ff(2,m) 1927 130 continue 1928 135 continue 1929 1930c Recursively build the remaining auxiliary functions (i.e., RIJKj, j=0). 1931 1932 call hfmkr(R,IJK,PC,NPP,Lp,Lp3) 1933 1934c Transfer to R0 array. 1935 1936 if( CENTER )then 1937 do 141 n = 1,Lp3 1938 do 140 m = 1,NPP 1939 R0C(ic,m,n) = -zan(ic)*R(m,0,n) 1940 R0(m,n) = R0(m,n) + R0C(ic,m,n) 1941 140 continue 1942 141 continue 1943 else 1944 do 146 n = 1,Lp3 1945 do 145 m = 1,NPP 1946c* R0(m,n) = R0(m,n) - zan(ic)*R(m,0,n) 1947 R0C(ic,m,n) = R(m,0,n) 1948 145 continue 1949 146 continue 1950 end if 1951 1952 150 continue 1953 1954 end 1955 Subroutine hf1tmp(Axyz,Aprims,Acoefs,NPA,NCA,La, 1956 & Bxyz,Bprims,Bcoefs,NPB,NCB,Lb, 1957 & Cxyz,zan,ncenters, 1958 & bO2I,bKEI,bNAI,Nint,O2I,KEI,NAI,canAB, 1959 & DryRun,W0,maxW0) 1960 1961 Implicit real*8 (a-h,o-z) 1962 Implicit integer (i-n) 1963#include "errquit.fh" 1964 1965 Logical O2I,KEI,NAI,canAB 1966 1967 Logical GenCon,DryRun 1968 1969c--> Cartesian Coordinates, Primitives & Contraction Coefficients 1970 1971 Dimension Axyz(3),Aprims(NPA),Acoefs(NPA,NCA) 1972 Dimension Bxyz(3),Bprims(NPB),Bcoefs(NPB,NCB) 1973 1974c--> Nuclear Cartesian Coordinates & Charges 1975 1976 Dimension Cxyz(3,ncenters),zan(ncenters) 1977 1978c--> Blocks of Overlap, Kinetic Energy & Nuclear Attraction Integrals 1979 1980 Dimension bO2I(Nint),bKEI(Nint),bNAI(Nint) 1981 1982c--> Scratch Space. 1983 1984 Dimension W0(maxW0) 1985c 1986c Compute the overlap, kinetic energy, and nuclear attraction integrals for 1987c two shells of contracted Gaussians functions. This driver is NOT capable of 1988c evaluating integral derivatives. 1989c 1990c****************************************************************************** 1991#if defined(INTDEBUG) 1992 write(6,*)' inside hf1 ' 1993 write(6,*)' npa,nca,la = ',npa,nca,la 1994 write(6,*)' npb,ncb,lb = ',npb,ncb,lb 1995 write(6,*)' ncenters = ',ncenters 1996 write(6,*)' NINT = ',nint 1997 write(6,*)' maxW0 = ',maxw0 1998 write(6,*)' <canAB:DryRun>-<',canab,':',dryrun,'>' 1999 write(6,*)' <o2i:kei:nai>-<',o2i,':',kei,':',nai,'>' 2000 write(6,'(a,3(2x,1pd20.10))')' Axyz =',Axyz 2001 write(6,'(a,3(2x,1pd20.10))')' Bxyz =',Bxyz 2002 write(6,'(a,100(3(2x,1pd20.10/)))')' Cxyz =',Cxyz 2003 do iiii = 1,npa 2004 write(6,'(a,i3,a,2(2x,1pd20.10))') 2005 & 'Aprims:Acoeffs:(',iiii,') =',Aprims(iiii), 2006 & Acoefs(iiii,1) 2007 enddo 2008 do iiii = 1,npb 2009 write(6,'(a,i3,a,2(2x,1pd20.10))') 2010 & 'Bprims:Bcoeffs:(',iiii,') =',Bprims(iiii), 2011 & Bcoefs(iiii,1) 2012 enddo 2013#endif 2014* if (.not.dryrun) then 2015* call hf_print('hf1: a shell',axyz,aprims,acoefs,npa,nca,la) 2016* call hf_print('hf1: b shell',bxyz,bprims,bcoefs,npb,ncb,lb) 2017* endif 2018 MXD = 0 2019 if (KEI) call errquit('hf1tmp: only for pot ints',911, INT_ERR) 2020 if (O2I) call errquit('hf1tmp: only for pot ints',911, INT_ERR) 2021 2022c Determine whether general or segmented contraction is used. 2023 2024 NCP = NCA*NCB 2025 2026 GenCon = NCP.ne.1 2027 2028 if( GenCon )then 2029 write(*,*) 'HF1: Not yet ready for general contraction.' 2030 stop 2031 end if 2032 2033c To determine all the Hermite expansion coefficients required to evaluate 2034c the kinetic energy integrals, increment the angular momenta by one. 2035 2036 if( KEI )then 2037 Li = 1 2038 else 2039 Li = 0 2040 end if 2041 2042c Define the angular momentum of the overlap distribution. 2043 2044 Lp = La + Lb 2045 2046c Increment "Lp" to account for the order of differentiation. 2047 2048 Lp = Lp + MXD 2049 2050c Define the accumulated number of angular momentum functions <= Lp. 2051 2052 Lp3 = ((Lp+1)*(Lp+2)*(Lp+3))/6 2053 2054c Define the prefactor of the overlap distribution "P". 2055 2056c Assign pointers to scratch space. 2057 2058 i_ALPHAp = 1 2059 i_IPAIRp = i_ALPHAp + 2*(NPA*NPB) 2060 i_left = i_IPAIRp + 2*(NPA*NPB) - 1 2061 2062 i_ESp = (maxW0+1) - 3*(NPA*NPB) 2063 i_right = i_ESp 2064 2065 if( i_left.ge.i_right )then 2066 2067 write(*,*) 'HF1: Insufficient scratch space.' 2068 write(*,*) ' needed ',i_left + (maxW0-(i_right-1)) 2069 write(*,*) ' allocated ',maxW0 2070 2071 write(*,*) 'From the left ' 2072 write(*,*) 'ALPHAp: ',i_ALPHAp 2073 write(*,*) 'IPAIRp: ',i_IPAIRp 2074 write(*,*) 'From the right ' 2075 write(*,*) 'ESp : ',i_ESp 2076 2077 stop 2078 2079 end if 2080 2081 if( DryRun )then 2082 2083 MaxMem = i_left + (maxW0 - (i_right-1)) 2084 NPP = NPA*NPB 2085 2086 else 2087 2088 call hfset(Axyz,Aprims,Acoefs,NPA,NCA, 2089 & Bxyz,Bprims,Bcoefs,NPB,NCB, 2090 & GenCon,W0(i_ALPHAp),W0(i_IPAIRp),W0(i_ESp),NPP) 2091 2092 end if 2093 2094 if (npp.eq.0) then 2095 if (O2I) call dfill(Nint,0.0d00,bO2I,1) 2096 if (KEI) call dfill(Nint,0.0d00,bKEI,1) 2097 if (NAI) call dfill(Nint,0.0d00,bNAI,1) 2098 return 2099 endif 2100c Define the Hermite linear expansion coefficients. 2101 2102c Assign pointers to scratch space. 2103 2104 lprod = ((La+Li)+(Lb+Li)+1)*((La+Li)+1)*((Lb+Li)+1) 2105 2106 i_Ep = i_IPAIRp + 2*(NPA*NPB) 2107 i_pf = i_Ep + 3*NPP*(MXD+1)*lprod 2108 i_left = i_pf + 2*NPP - 1 2109 2110 if( i_left.ge.i_right )then 2111 2112 write(*,*) 'HF1: Insufficient scratch space.' 2113 write(*,*) ' needed ',i_left + (maxW0-(i_right-1)) 2114 write(*,*) ' allocated ',maxW0 2115 2116 write(*,*) 'From the right ' 2117 write(*,*) 'ALPHAp: ',i_ALPHAp 2118 write(*,*) 'IPAIRp: ',i_IPAIRp 2119 write(*,*) 'Ep : ',i_Ep 2120 write(*,*) 'pf : ',i_pf 2121 write(*,*) 'From the left ' 2122 write(*,*) 'ESp : ',i_ESp 2123 2124 stop 2125 2126 end if 2127 2128 if( DryRun )then 2129 2130 MaxMem = max( MaxMem, i_left + (maxW0 - (i_right-1)) ) 2131 2132 else 2133 2134 do 100 nd = 0,MXD 2135 call hfmke(Axyz,Bxyz,W0(i_ALPHAp),W0(i_ESp),W0(i_Ep),W0(i_pf), 2136 & nd,NPP,MXD,La+Li,Lb+Li) 2137 100 continue 2138 2139 end if 2140 2141c Compute the 2-center overlap integrals, <a|S|b>. 2142 2143 if( O2I )then 2144 if( .not. DryRun )then 2145 call hf2oi(W0(i_Ep),bO2I,Nint,NPP,La,Lb,Li,canAB) 2146 end if 2147 end if 2148 2149c Compute kinetic energy integrals, <a|T|b>. 2150 2151 if( KEI )then 2152 2153c Assign pointers to scratch space. 2154 2155 i_Ti = i_Ep + 3*NPP*(MXD+1)*lprod 2156 i_top = i_Ti + NPP - 1 2157 2158 if( i_top.gt.maxW0 )then 2159 2160 write(*,*) 'HF1: Insufficient scratch space.' 2161 write(*,*) ' needed ',i_top 2162 write(*,*) ' allocated ',maxW0 2163 2164 write(*,*) 'ALPHAp: ',i_ALPHAp 2165 write(*,*) 'IPAIRp: ',i_IPAIRp 2166 write(*,*) 'Ep : ',i_Ep 2167 write(*,*) 'Ti : ',i_Ti 2168 2169 stop 2170 2171 end if 2172 2173 if( DryRun )then 2174 2175 MaxMem = max( MaxMem, i_top ) 2176 2177 else 2178 2179 call hfkei(W0(i_ALPHAp),W0(i_Ep),bKEI,W0(i_Ti), 2180 & Nint,NPP,La,Lb,Li,canAB) 2181 end if 2182 2183 end if 2184 2185c Compute nuclear attraction integrals, <a|V|b>. 2186 2187 if( NAI )then 2188 2189c Define the auxiliary function integrals. 2190 2191c Assign scratch space. 2192 2193 i_R0 = i_Ep + 3*NPP*(MXD+1)*lprod 2194 i_IJK = i_R0 + NPP*Lp3*ncenters 2195 i_P = i_IJK + (Lp+1)**3 2196 i_RS = i_P + NPP*3 2197 i_PC = i_RS + NPP 2198 i_ff = i_PC + NPP*3 2199 i_Rj = i_ff + NPP*2 2200 i_top = i_Rj + NPP*(Lp+1)*Lp3 - 1 2201 2202 if( i_top.gt.maxW0 )then 2203 2204 write(*,*) 'HF1: Insufficient scratch space.' 2205 write(*,*) ' needed ',i_top 2206 write(*,*) ' allocated ',maxW0 2207 2208 write(*,*) 'ALPHAp: ',i_ALPHAp 2209 write(*,*) 'IPAIRp: ',i_IPAIRp 2210 write(*,*) 'Ep : ',i_Ep 2211 write(*,*) 'R0 : ',i_R0 2212 write(*,*) 'IJK : ',i_IJK 2213 write(*,*) 'P : ',i_P 2214 write(*,*) 'RS : ',i_RS 2215 write(*,*) 'PC : ',i_PC 2216 write(*,*) 'ff : ',i_ff 2217 write(*,*) 'Rj : ',i_Rj 2218 2219 stop 2220 2221 end if 2222 2223 if( DryRun )then 2224 2225 MaxMem = max( MaxMem, i_top ) 2226 2227 else 2228 2229 call hf1mkrtmp(Axyz,Bxyz,Cxyz,zan,ncenters, 2230 & W0(i_ALPHAp),W0(i_P),W0(i_RS),W0(i_PC),W0(i_ff), 2231 & W0(i_Rj),W0(i_R0),W0(i_R0),W0(i_IJK), 2232 & NPP,Lp,Lp3,.FALSE.) 2233 2234 call hfnaitmp(W0(i_Ep),W0(i_R0),W0(i_IJK),bNAI, 2235 & Nint,NPP,La,Lb,Li,Lp,Lp3,canAB,ncenters) 2236 2237 end if 2238 2239 end if 2240 2241c Return the maximum amount of scratch space required by a "dry run". 2242 2243 if( DryRun ) maxW0 = MaxMem 2244c 2245 end 2246 Subroutine hfnaitmp(E,R0,IJK,Vab,Nint,NPP, 2247 & La,Lb,Li,Lp,Lp3,canAB,ncenters) 2248 2249 Implicit real*8 (a-h,o-z) 2250 Implicit integer (i-n) 2251 2252 Logical canAB 2253 2254c--> Hermite Linear Expansion Coefficients 2255 2256 Dimension E(3,NPP,0:((La+Li)+(Lb+Li)),0:(La+Li),0:(Lb+Li)) 2257 2258c--> Auxiliary Function Integrals & Index 2259 2260 Dimension R0(ncenters,NPP,Lp3),IJK(0:Lp,0:Lp,0:Lp) 2261 2262c--> Nuclear Attraction Integrals 2263 2264 Dimension Vab(Nint) 2265 2266c--> Scratch Space 2267 2268 Dimension Nxyz(3) 2269c 2270c Compute the nuclear attraction integrals. 2271c 2272c formula: 2273c __ 2274c \ Ia,Ib Ja,Jb Ka,Kb 2275c Vab = / Ex * Ey * Ez * R 2276c -- Ip Jp Kp Ip,Jp,Kp 2277c Ip,Jp,Kp 2278c 2279c****************************************************************************** 2280 2281c Initialize the block of NAIs. 2282 2283 do 10 nn = 1,Nint 2284 Vab(nn) = 0.D0 2285 10 continue 2286 2287c Define the number of shell components on each center. 2288 2289 La2 = ((La+1)*(La+2))/2 2290 Lb2 = ((Lb+1)*(Lb+2))/2 2291 2292c loop over centers 2293 do 23456 icic = 1,ncenters 2294c Loop over shell components. 2295 2296 nn = 0 2297 2298 do 50 ma = 1,La2 2299 2300c Define the angular momentum indices for shell "A". 2301 2302 call getNxyz(La,ma,Nxyz) 2303 2304 Ia = Nxyz(1) 2305 Ja = Nxyz(2) 2306 Ka = Nxyz(3) 2307 2308 if( canAB )then 2309 mb_limit = ma 2310 else 2311 mb_limit = Lb2 2312 end if 2313 2314 do 40 mb = 1,mb_limit 2315 2316c Define the angular momentum indices for shell "B". 2317 2318 call getNxyz(Lb,mb,Nxyz) 2319 2320 Ib = Nxyz(1) 2321 Jb = Nxyz(2) 2322 Kb = Nxyz(3) 2323 2324 nn = nn + 1 2325 2326 do 30 Ip = 0,Ia+Ib 2327 do 30 Jp = 0,Ja+Jb 2328 do 30 Kp = 0,Ka+Kb 2329 2330 np = IJK(Ip,Jp,Kp) 2331 2332 do 20 mp = 1,NPP 2333 Vab(nn) = Vab(nn) + (E(1,mp,Ip,Ia,Ib)* 2334 & E(2,mp,Jp,Ja,Jb)* 2335 & E(3,mp,Kp,Ka,Kb))*R0(icic,mp,np) 2336 20 continue 2337 2338 30 continue 2339 2340 40 continue 2341 2342 50 continue 2343 2344 write(6,*)'==================================================', 2345 & '==============================' 2346 write(6,*)' for center ',icic,' modified V is ' 2347 do inn = 1,nn 2348 if (abs(Vab(inn)).gt.1.0d-07) then 2349 write(6,*)' Vab(',inn,' ) =',vab(inn),' of ',nn 2350 endif 2351 enddo 2352 write(6,*)'==================================================', 2353 & '==============================' 235423456 continue 2355 end 2356 subroutine raktest_gc(rtdb) 2357 implicit none 2358#include "errquit.fh" 2359#include "stdio.fh" 2360#include "bas.fh" 2361#include "rtdb.fh" 2362#include "geom.fh" 2363#include "mafdecls.fh" 2364#include "util.fh" 2365c 2366 integer rtdb 2367c 2368 integer geom, basisseg, basisns, bases(2) 2369 integer mx1e, mxg, mxs1, mxs2, membuf, memscr 2370 integer h_buf, k_buf, h_scr, k_scr 2371 integer h_s_seg, k_s_seg, h_s_ns, k_s_ns 2372 integer h_diff, k_diff 2373 integer nbf_seg, ncont_seg, nbf_ns, ncont_ns 2374 integer icount, niter 2375 double precision norm, tov_seg, tov_ns 2376 double precision thresh_norm 2377 logical status 2378 logical FF, FT 2379 parameter (FF=.false.,FT=.true.) 2380c 2381 logical int_normalize 2382 external int_normalize 2383c 2384 thresh_norm = 1.0d-06 2385c 2386 write(luout,*)' ============================================ ' 2387 write(luout,*)' ================ ================ ' 2388 write(luout,*)' ================ raktest_gc ================ ' 2389 write(luout,*)' ================ ================ ' 2390 write(luout,*)' ============================================ ' 2391c 2392 if (.not.geom_create(geom,'geometry')) 2393 & call errquit('raktest_gc: geom_create failed',911, GEOM_ERR) 2394 if (.not.geom_rtdb_load(rtdb,geom,'geometry')) 2395 & call errquit('raktest_gc: geom_rtdb_load failed',911, 2396 & RTDB_ERR) 2397c 2398 if(.not.bas_create(basisseg,'aobasisseg')) 2399 & call errquit('raktest_gc: bas_create choked',911, BASIS_ERR) 2400 if(.not. 2401 & bas_rtdb_load(rtdb,geom,basisseg,'aobasisseg')) 2402 & call errquit('raktest_gc: bas_rtdb_load choked',911, 2403 & RTDB_ERR) 2404c 2405 if(.not.bas_create(basisns,'aobasisns')) 2406 & call errquit('raktest_gc: bas_create choked',911, BASIS_ERR) 2407 if(.not. 2408 & bas_rtdb_load(rtdb,geom,basisns,'aobasisns')) 2409 & call errquit('raktest_gc: bas_rtdb_load choked',911, 2410 & RTDB_ERR) 2411c 2412 bases(1) = basisseg 2413 bases(2) = basisns 2414c 2415 status = geom_print(geom) 2416 status = bas_print(basisseg) 2417 status = gbs_map_print(basisseg) 2418 status = bas_print(basisns) 2419 status = gbs_map_print(basisns) 2420c 2421c... get memory requirements for segmented basis 2422 call int_init(rtdb,1,basisseg) 2423 call int_mem(mx1e, mxg, mxs1, mxs2) 2424 write(luout,*)' seg: one electron buffer size :',mx1e 2425 write(luout,*)' seg: two electron buffer size :',mxg 2426 write(luout,*)' seg: one electron scratch buffer size:',mxs1 2427 write(luout,*)' seg: two electron scratch buffer size:',mxs2 2428 call int_terminate() 2429c 2430c... get memory requirements for non-segmented basis 2431 call int_init(rtdb,1,basisns) 2432 call int_mem(mx1e, mxg, mxs1, mxs2) 2433 write(luout,*)' ns: one electron buffer size :',mx1e 2434 write(luout,*)' ns: two electron buffer size :',mxg 2435 write(luout,*)' ns: one electron scratch buffer size:',mxs1 2436 write(luout,*)' ns: two electron scratch buffer size:',mxs2 2437 call int_terminate() 2438 call int_init(rtdb,2,bases) 2439 if(.not.int_normalize(rtdb,basisseg)) call errquit 2440 & ('raktest_gc: int_normalize seg failed',911, INT_ERR) 2441 if(.not.int_normalize(rtdb,basisns)) call errquit 2442 & ('raktest_gc: int_normalize ns failed',911, INT_ERR) 2443 write(6,*)' after normalization' 2444 status = bas_print(basisseg) 2445 status = gbs_map_print(basisseg) 2446 status = bas_print(basisns) 2447 status = gbs_map_print(basisns) 2448c 2449 call int_mem(mx1e, mxg, mxs1, mxs2) 2450 write(luout,*)' both: one electron buffer size :',mx1e 2451 write(luout,*)' both: two electron buffer size :',mxg 2452 write(luout,*)' both: one electron scratch buffer size:',mxs1 2453 write(luout,*)' both: two electron scratch buffer size:',mxs2 2454c 2455 membuf = max(mx1e,mxg,10000) 2456 memscr = max(mxs1,mxs2,100000) 2457c 2458 if(.not.ma_push_get(mt_dbl,membuf,'int buff',h_buf,k_buf)) 2459 & call errquit('raktest_gc: ma failed 1',911, MA_ERR) 2460 if(.not.ma_push_get(mt_dbl,memscr,'int scr',h_scr,k_scr)) 2461 & call errquit('raktest_gc: ma failed 2',911, MA_ERR) 2462c 2463 if(.not.bas_numcont(basisseg,ncont_seg)) 2464 & call errquit('raktest_gc: bas_numcont failed',911, 2465 & BASIS_ERR) 2466 write(luout,*)' ao basis seg number of contractions :', 2467 & ncont_seg 2468 if(.not.bas_numbf(basisseg,nbf_seg)) 2469 & call errquit('raktest_gc: bas_numbf failed',911, 2470 & BASIS_ERR) 2471 write(luout,*)' ao basis seg number of basis functions :', 2472 & nbf_seg 2473 if(.not.bas_numcont(basisns,ncont_ns)) 2474 & call errquit('raktest_gc: bas_numcont failed',911, 2475 & BASIS_ERR) 2476 write(luout,*)' ao basis noseg number of contractions :', 2477 & ncont_ns 2478 if(.not.bas_numbf(basisns,nbf_ns)) 2479 & call errquit('raktest_gc: bas_numbf failed',911, 2480 & BASIS_ERR) 2481 write(luout,*)' ao basis noseg number of basis functions :', 2482 & nbf_ns 2483c 2484 if (nbf_ns.ne.nbf_seg) call errquit 2485 & ('raktest_gc: nbf_seg.ne.nbf_ns',(nbf_ns-nbf_seg), 2486 & UNKNOWN_ERR) 2487c 2488 if(.not.ma_push_get(mt_dbl,(nbf_seg*nbf_seg),'overlap seg', 2489 & h_s_seg,k_s_seg)) 2490 & call errquit('raktest_gc: ma failed 3',911, MA_ERR) 2491 if(.not.ma_push_get(mt_dbl,(nbf_ns*nbf_ns),'overlap noseg', 2492 & h_s_ns,k_s_ns)) 2493 & call errquit('raktest_gc: ma failed 4',911, MA_ERR) 2494 if(.not.ma_push_get(mt_dbl,(nbf_ns*nbf_ns),'difference matrix', 2495 & h_diff,k_diff)) 2496 & call errquit('raktest_gc: ma failed 4',911, MA_ERR) 2497 2498c 2499c set niter number of iterations 2500c 2501 niter = 20 2502 write(6,'(/,a,i5)') 2503 & 'number of instances computed for each seg/noseg set:',niter 2504c 2505c zero ma segments 2506 2507 write(6,*)' dfill k_diff',k_diff 2508 call dfill((nbf_ns*nbf_ns), 0.0d00,dbl_mb(k_diff),1) 2509 write(6,*)' dfill k_ns',k_s_ns 2510 call dfill((nbf_ns*nbf_ns), 0.0d00,dbl_mb(k_s_ns),1) 2511 write(6,*)' dfill k_seg',k_s_seg 2512 call dfill((nbf_seg*nbf_seg),0.0d00,dbl_mb(k_s_seg),1) 2513 write(6,*)' dfill k_scr',k_scr 2514 call dfill(memscr, 0.0d00,dbl_mb(k_scr),1) 2515 write(6,*)' dfill k_buf',k_buf 2516 call dfill(membuf, 0.0d00,dbl_mb(k_buf),1) 2517c 2518c overlap check 2519 tov_seg = util_cpusec() 2520 do icount = 1,niter 2521 call raktest_bs_gc(basisseg, 2522 & memscr,dbl_mb(k_scr), 2523 & membuf,dbl_mb(k_buf), 2524 & nbf_seg,ncont_seg,dbl_mb(k_s_seg), 2525 & FT,FF,FF,' segment overlap') 2526 enddo 2527 tov_seg = util_cpusec() - tov_seg 2528 tov_ns = util_cpusec() 2529 do icount = 1,niter 2530 call raktest_bs_gc(basisns, 2531 & memscr,dbl_mb(k_scr), 2532 & membuf,dbl_mb(k_buf), 2533 & nbf_ns,ncont_ns,dbl_mb(k_s_ns), 2534 & FT,FF,FF,'nosegment overlap') 2535 2536 enddo 2537 tov_ns = util_cpusec() - tov_ns 2538c 2539 call dcopy((nbf_ns*nbf_ns),dbl_mb(k_s_ns),1,dbl_mb(k_diff),1) 2540 call daxpy((nbf_seg*nbf_seg),(-1.0d00), 2541 & dbl_mb(k_s_seg),1,dbl_mb(k_diff),1) 2542 norm = ddot((nbf_ns*nbf_ns),dbl_mb(k_diff),1,dbl_mb(k_diff),1) 2543 if (norm.gt.thresh_norm) 2544 & call print_diff_gc(nbf_ns, 2545 & dbl_mb(k_diff), 2546 & dbl_mb(k_s_seg), 2547 & dbl_mb(k_s_ns), 2548 & thresh_norm, 2549 & 'overlap') 2550 write(luout,*)' ' 2551 write(luout,*)' time for segmented overlap :',tov_seg 2552 write(luout,*)' time for non-segmented overlap :',tov_ns 2553 write(luout,'(a,f10.2)') 2554 & ' % speedup :', 2555 & (tov_seg-tov_ns)/tov_seg*100.0d00 2556 write(luout,*)'raktest_gc: overlap difference norm :',norm 2557 write(luout,*)' ' 2558c 2559c zero ma segments 2560 call dfill((nbf_ns*nbf_ns), 0.0d00,dbl_mb(k_diff),1) 2561 call dfill((nbf_ns*nbf_ns), 0.0d00,dbl_mb(k_s_ns),1) 2562 call dfill((nbf_seg*nbf_seg),0.0d00,dbl_mb(k_s_seg),1) 2563 call dfill(memscr, 0.0d00,dbl_mb(k_scr),1) 2564 call dfill(membuf, 0.0d00,dbl_mb(k_buf),1) 2565c kinetic energy check 2566 tov_seg = util_cpusec() 2567 do icount = 1,niter 2568 call raktest_bs_gc(basisseg, 2569 & memscr,dbl_mb(k_scr), 2570 & membuf,dbl_mb(k_buf), 2571 & nbf_seg,ncont_seg,dbl_mb(k_s_seg), 2572 & FF,FT,FF,' segment kinetic') 2573 enddo 2574 tov_seg = util_cpusec() - tov_seg 2575 tov_ns = util_cpusec() 2576 do icount = 1,niter 2577 call raktest_bs_gc(basisns, 2578 & memscr,dbl_mb(k_scr), 2579 & membuf,dbl_mb(k_buf), 2580 & nbf_ns,ncont_ns,dbl_mb(k_s_ns), 2581 & FF,FT,FF,'nosegment kinetic') 2582 2583 enddo 2584 tov_ns = util_cpusec() - tov_ns 2585c 2586 call dcopy((nbf_ns*nbf_ns),dbl_mb(k_s_ns),1,dbl_mb(k_diff),1) 2587 call daxpy((nbf_seg*nbf_seg),(-1.0d00), 2588 & dbl_mb(k_s_seg),1,dbl_mb(k_diff),1) 2589 norm = ddot((nbf_ns*nbf_ns),dbl_mb(k_diff),1,dbl_mb(k_diff),1) 2590 if (norm.gt.thresh_norm) 2591 & call print_diff_gc(nbf_ns, 2592 & dbl_mb(k_diff), 2593 & dbl_mb(k_s_seg), 2594 & dbl_mb(k_s_ns), 2595 & thresh_norm, 2596 & 'kinetic') 2597 write(luout,*)' ' 2598 write(luout,*)' time for segmented kinetic :',tov_seg 2599 write(luout,*)' time for non-segmented kinetic :',tov_ns 2600 write(luout,'(a,f10.2)') 2601 & ' % speedup :', 2602 & (tov_seg-tov_ns)/tov_seg*100.0d00 2603 write(luout,*)'raktest_gc: kinetic difference norm :',norm 2604 write(luout,*)' ' 2605c 2606c zero ma segments 2607 call dfill((nbf_ns*nbf_ns), 0.0d00,dbl_mb(k_diff),1) 2608 call dfill((nbf_ns*nbf_ns), 0.0d00,dbl_mb(k_s_ns),1) 2609 call dfill((nbf_seg*nbf_seg),0.0d00,dbl_mb(k_s_seg),1) 2610 call dfill(memscr, 0.0d00,dbl_mb(k_scr),1) 2611 call dfill(membuf, 0.0d00,dbl_mb(k_buf),1) 2612c potential check 2613 tov_seg = util_cpusec() 2614 do icount = 1,niter 2615 call raktest_bs_gc(basisseg, 2616 & memscr,dbl_mb(k_scr), 2617 & membuf,dbl_mb(k_buf), 2618 & nbf_seg,ncont_seg,dbl_mb(k_s_seg), 2619 & FF,FF,FT,' segment potential') 2620 enddo 2621 tov_seg = util_cpusec() - tov_seg 2622 tov_ns = util_cpusec() 2623 do icount = 1,niter 2624 call raktest_bs_gc(basisns, 2625 & memscr,dbl_mb(k_scr), 2626 & membuf,dbl_mb(k_buf), 2627 & nbf_ns,ncont_ns,dbl_mb(k_s_ns), 2628 & FF,FF,FT,'nosegment potential') 2629 2630 enddo 2631 tov_ns = util_cpusec() - tov_ns 2632c 2633 call dcopy((nbf_ns*nbf_ns),dbl_mb(k_s_ns),1,dbl_mb(k_diff),1) 2634 call daxpy((nbf_seg*nbf_seg),(-1.0d00), 2635 & dbl_mb(k_s_seg),1,dbl_mb(k_diff),1) 2636 norm = ddot((nbf_ns*nbf_ns),dbl_mb(k_diff),1,dbl_mb(k_diff),1) 2637 if (norm.gt.thresh_norm) 2638 & call print_diff_gc(nbf_ns, 2639 & dbl_mb(k_diff), 2640 & dbl_mb(k_s_seg), 2641 & dbl_mb(k_s_ns), 2642 & thresh_norm, 2643 & 'potential') 2644 write(luout,*)' ' 2645 write(luout,*)' time for segmented potential :',tov_seg 2646 write(luout,*)' time for non-segmented potential :',tov_ns 2647 write(luout,'(a,f10.2)') 2648 & ' % speedup :', 2649 & (tov_seg-tov_ns)/tov_seg*100.0d00 2650 write(luout,*)'raktest_gc: potential difference norm :',norm 2651 write(luout,*)' ' 2652c 2653 status = ma_pop_stack(h_diff) 2654 status = status .and. ma_pop_stack(h_s_ns) 2655 status = status .and. ma_pop_stack(h_s_seg) 2656 status = status .and. ma_pop_stack(h_scr) 2657 status = status .and. ma_pop_stack(h_buf) 2658 if (.not.status) call errquit('raktest_gc: pop error',911, MA_ERR) 2659c 2660 if (.not.bas_destroy(basisseg)) call errquit 2661 & ('raktest_gc: bas_destroy failed ?',911, BASIS_ERR) 2662 if (.not.bas_destroy(basisns)) call errquit 2663 & ('raktest_gc: bas_destroy failed ?',911, BASIS_ERR) 2664 if (.not.geom_destroy(geom)) call errquit 2665 & ('raktest_gc: geom_destroy failed ?',911, GEOM_ERR) 2666c 2667 call int_terminate() 2668c 2669 end 2670 subroutine raktest_bs_gc(basis,nscr,scr,nbuf,buf,nbf,ncont,S, 2671 & OV,KE,PE,msg) 2672 implicit none 2673#include "errquit.fh" 2674c 2675#include "stdio.fh" 2676c 2677#include "bas.fh" 2678c 2679 integer basis,nscr,nbuf,nbf,ncont 2680 double precision scr(nscr), buf(nbuf) 2681 double precision S(nbf,nbf) 2682 logical OV,KE,PE 2683 character*(*) msg 2684c 2685 double precision val 2686 integer nint 2687 integer ish, ilo, ihi, ibf 2688 integer jsh, jlo, jhi, jbf 2689 integer icount 2690c 2691c check validity 2692 icount = 0 2693 if (OV) icount = icount + 1 2694 if (KE) icount = icount + 1 2695 if (PE) icount = icount + 1 2696 if (icount.eq.0) then 2697 write(luout,*)' no integral set defined ' 2698 call dfill ((nbf*nbf),0.0d00,S,1) 2699 return 2700 elseif(icount.ne.1) then 2701 write(luout,*)' OV =',OV 2702 write(luout,*)' KE =',KE 2703 write(luout,*)' PE =',PE 2704 stop 'error' 2705 endif 2706c 2707 do ish = 1,ncont 2708 do jsh = 1,ish 2709 if (.not.bas_cn2bfr(basis,ish,ilo,ihi)) call errquit 2710 & ('raktest_gc: cn2bfr failed',911, BASIS_ERR) 2711* write(luout,*)'ish = ',ish,' bfr =',ilo,ihi 2712 if (.not.bas_cn2bfr(basis,jsh,jlo,jhi)) call errquit 2713 & ('raktest_gc: cn2bfr failed',911, BASIS_ERR) 2714* write(luout,*)'jsh = ',jsh,' bfr =',jlo,jhi 2715 nint = (ihi-ilo+1)*(jhi-jlo+1) 2716* write(luout,*)'number of integrals = ',nint 2717 if (OV) then 2718 call int_1eov(basis,ish,basis,jsh, 2719 & nscr,scr,nbuf,buf) 2720 elseif (KE) then 2721 call int_1eke(basis,ish,basis,jsh, 2722 & nscr,scr,nbuf,buf) 2723 elseif (PE) then 2724 call int_1epe(basis,ish,basis,jsh, 2725 & nscr,scr,nbuf,buf) 2726 else 2727 write(luout,*)' no integral set defined ' 2728 call dfill ((nbf*nbf),0.0d00,S,1) 2729 return 2730 endif 2731 icount = 0 2732 nint = 0 2733 do ibf = ilo, ihi 2734 do jbf = jlo, jhi 2735 nint = nint+1 2736 val = buf(nint) 2737 s(ibf,jbf) = val 2738 s(jbf,ibf) = val 2739 if (abs(val).gt.1.0d-07) then 2740 icount = icount + 1 2741* write(69,10000)ibf,jbf,val,msg,' ints ',icount 2742* write(luout,10000)ibf,jbf,val,msg,' ints ',icount 2743 endif 2744 enddo 2745 enddo 2746 enddo 2747 enddo 2748c 2749* write(luout,*)' raktest_gc matrix ',msg 2750* call output(S,1,nbf,1,nbf,nbf,nbf,1) 2751c 275210000 format(i5,i5,1pd20.10,1x, a,a,i12) 2753c 2754 end 2755 subroutine print_diff_gc(nbf,diff,seg,ns,ths,msg) 2756 implicit none 2757#include "stdio.fh" 2758 integer nbf 2759 double precision diff(nbf,nbf) 2760 double precision seg(nbf,nbf) 2761 double precision ns(nbf,nbf) 2762 double precision ths 2763 character*(*) msg 2764c 2765 integer i,j 2766c 2767 write(luout,*)' differeces from ',msg 2768c 2769 do i=1,nbf 2770 do j=1,nbf 2771 if (abs(diff(i,j)).gt.ths) then 2772 write(6,10000)i,j,diff(i,j),seg(i,j),ns(i,j),msg 2773 endif 2774 enddo 2775 enddo 277610000 format('<',i3,'|',i3,'> <diff',1pd20.10,'> <seg',1pd20.10, 2777 & '> <noseg',1pd20.10,'>',1x,a) 2778 end 2779 subroutine print_diff_vec(n,a,b,ths,msg) 2780 implicit none 2781 integer n 2782 double precision a(n), b(n) 2783 double precision ths 2784 character*(*) msg 2785c 2786 integer i, count, nza, nzb 2787 double precision diff 2788c 2789 write(6,*)' print_diff_vec <<',msg,'>>' 2790 write(6,*)' print_diff_vec threshold:',ths 2791c 2792 nza = 0 2793 nzb = 0 2794 count = 0 2795 do i = 1,n 2796 if (a(i).ne.0.0d00) nza = nza + 1 2797 if (b(i).ne.0.0d00) nzb = nzb + 1 2798 diff = a(i) - b(i) 2799 if (abs(diff).gt.ths) then 2800 count = count + 1 2801 write(6,'(1x,i8,a,d12.6,a,d12.6,a,d12.6)') 2802 & i,' th element delta:',diff, 2803 & ' a:',a(i),' b:',b(i) 2804 endif 2805 enddo 2806 write(6,*)' print_diff_vec: number of different elements :', 2807 & count 2808 write(6,*)' print_diff_vec: number of nonzero elements in a :', 2809 & nza 2810 write(6,*)' print_diff_vec: number of nonzero elements in b :', 2811 & nzb 2812 end 2813 subroutine raktest_test9(rtdb) 2814c 2815 implicit none 2816#include "errquit.fh" 2817#include "mafdecls.fh" 2818#include "geom.fh" 2819#include "bas.fh" 2820#include "rtdb.fh" 2821c 2822 integer rtdb 2823c 2824 integer basis, geom, count 2825 integer i, j, k, l 2826 integer ish, jsh, ksh, lsh 2827 integer ilo, jlo, klo, llo 2828 integer ihi, jhi, khi, lhi 2829 integer k_buf, h_buf, k_scr, h_scr 2830 integer max1e, max2e, mscr1, mscr2, m_buf, m_scr 2831 integer inshell(4) 2832 logical status 2833 logical int_normalize 2834 external int_normalize 2835c 2836 if (.not.geom_create(geom,'geometry')) 2837 & call errquit('raktest_geomwrt: geom_create failed?',911, 2838 & GEOM_ERR) 2839 if (.not. geom_rtdb_load(rtdb,geom,'geometry')) call errquit 2840 & ('raktest_geomwrt: geom_rtdb_load -ref failed',911, RTDB_ERR) 2841c 2842 if (.not.bas_create(basis,'ao sp_basis')) call errquit 2843 & ('bas_create failed',911, BASIS_ERR) 2844 if (.not.bas_rtdb_load(rtdb,geom,basis,'ao sp_basis')) 2845 & call errquit 2846 & ('bas_rtdb_load failed',911, BASIS_ERR) 2847 if (.not.int_normalize(rtdb,basis)) stop ' norm error 1' 2848c 2849 if (.not.geom_print(geom)) stop ' print error' 2850 if (.not.bas_print(basis)) stop ' print error' 2851c 2852 call int_init(rtdb,1,basis) 2853 call int_mem(max1e,max2e,mscr1,mscr2) 2854 m_buf = max(max1e*2,max2e*2) 2855 m_scr = max(mscr1*2,mscr2) 2856 m_buf = m_buf + (m_buf*110)/100 2857 m_scr = m_scr + (m_scr*110)/100 2858c 2859 if (.not.ma_push_get(mt_dbl,m_scr,'scr for 1e',h_scr,k_scr)) 2860 & call errquit('ma scr failed',911, MA_ERR) 2861c 2862 if (.not.ma_push_get(mt_dbl,m_buf,'buf for 1e',h_buf,k_buf)) 2863 & call errquit('ma buf failed',911, MA_ERR) 2864c 2865 ish = 1 2866 jsh = 1 2867 ksh = 1 2868 lsh = 1 2869 if (rtdb_get(rtdb,'rak9:shells',mt_int,4,inshell)) then 2870 ish = inshell(1) 2871 jsh = inshell(2) 2872 ksh = inshell(3) 2873 lsh = inshell(4) 2874 else 2875 write(6,*)'rak9:shells not set on rtdb' 2876 endif 2877c 2878 write(6,*)' rak9 for shells ',ish,jsh,ksh,lsh 2879c 2880 call int_2e4c(basis,ish,jsh,basis,ksh,lsh, 2881 & m_scr,dbl_mb(k_scr),m_buf,dbl_mb(k_buf)) 2882c 2883 if (.not.bas_cn2bfr(basis,ish,ilo,ihi)) 2884 & stop 'cn2bfr error i' 2885 if (.not.bas_cn2bfr(basis,jsh,jlo,jhi)) 2886 & stop 'cn2bfr error j' 2887 if (.not.bas_cn2bfr(basis,ksh,klo,khi)) 2888 & stop 'cn2bfr error k' 2889 if (.not.bas_cn2bfr(basis,lsh,llo,lhi)) 2890 & stop 'cn2bfr error l' 2891 count = 0 2892 do i=ilo,ihi 2893 do j=jlo,jhi 2894 do k=klo,khi 2895 do l=llo,lhi 2896 write(6,*)i,j,k,l,dbl_mb(k_buf+count) 2897 write(69,*)i,j,k,l,dbl_mb(k_buf+count) 2898 count = count + 1 2899 enddo 2900 enddo 2901 enddo 2902 enddo 2903c 2904 status = .true. 2905 status = status.and.ma_pop_stack(h_buf) 2906 status = status.and.ma_pop_stack(h_scr) 2907c 2908 if (.not.status) call errquit('ma pop fail',911, MA_ERR) 2909 call int_terminate() 2910c 2911 if (.not.bas_destroy(basis)) call errquit 2912 & ('9: bas_destroy failed ?',911, BASIS_ERR) 2913 if (.not.geom_destroy(geom)) call errquit 2914 & ('9: geom_destroy failed ?',911, GEOM_ERR) 2915c 2916 end 2917 subroutine raktest_ecp(rtdb) 2918 implicit none 2919#include "errquit.fh" 2920#include "geom.fh" 2921#include "bas.fh" 2922#include "rtdb.fh" 2923 integer rtdb 2924c 2925 integer geom, basis, ecpid 2926c 2927 logical int_normalize, int_ecp_init 2928 external int_normalize, int_ecp_init 2929c 2930* if (.not.rtdb_print(rtdb,.true.)) stop ' rtdb_p 1?' 2931c 2932 if (.not.geom_create(geom,'geometry')) 2933 & call errquit('raktest_geomwrt: geom_create failed?',911, 2934 & GEOM_ERR) 2935 if (.not. geom_rtdb_load(rtdb,geom,'geometry')) call errquit 2936 & ('raktest_geomwrt: geom_rtdb_load -ref failed',911, 2937 & RTDB_ERR) 2938 if (.not.geom_print(geom)) stop ' print error' 2939c 2940 if (.not.bas_create(basis,'ao basis')) call errquit 2941 & ('bas_create failed',911, BASIS_ERR) 2942 if (.not.bas_rtdb_load(rtdb,geom,basis,'ao basis')) 2943 & call errquit 2944 & ('bas_rtdb_load failed',911, RTDB_ERR) 2945 if (.not.bas_print(basis)) stop ' print error' 2946* if (.not.rtdb_print(rtdb,.true.)) stop ' rtdb_p 2?' 2947c 2948 if (.not.bas_get_ecp_handle(basis,ecpid)) stop 'get/ecp/handle' 2949*ecp: if (.not.bas_create(ecpid,'ecp basis')) call errquit 2950*ecp: & ('bas_create failed',911, BASIS_ERR) 2951*ecp: if (.not.bas_rtdb_load(rtdb,geom,ecpid,'ecp basis')) 2952*ecp: & call errquit 2953*ecp: & ('bas_rtdb_load failed',911, RTDB_ERR) 2954 if (.not.bas_print(ecpid)) stop ' print error' 2955c 2956 if (.not.int_normalize(rtdb,basis)) stop ' norm error 1' 2957 if (.not.bas_print(basis)) stop ' print error' 2958 if (.not.gbs_map_print(basis)) stop ' gbs map print error' 2959 if (.not.bas_print(ecpid)) stop ' print error' 2960 if (.not.gbs_map_print(basis)) stop ' gbs map print error' 2961c 2962 if (.not.int_ecp_init(ecpid,0,0)) stop ' error in int_ecp_init' 2963 call int_ecp_terminate() 2964c 2965 if (.not.bas_destroy(basis)) stop ' bas_dest error 1' 2966 if (.not.bas_destroy(ecpid)) stop ' bas_dest error 2' 2967 if (.not.geom_destroy(geom)) stop ' geom_dest error 1' 2968 if (.not.bas_version()) stop ' bas_version error' 2969c 2970 end 2971 subroutine raktest_bug(rtdb) 2972 implicit none 2973#include "mafdecls.fh" 2974#include "geom.fh" 2975#include "bas.fh" 2976#include "nwc_const.fh" 2977#include "basP.fh" 2978#include "basdeclsP.fh" 2979#include "geomP.fh" 2980#include "geobasmapP.fh" 2981#include "bas_exndcf_dec.fh" 2982#include "bas_ibs_dec.fh" 2983c 2984 logical int_normalize 2985 external int_normalize 2986c 2987 integer rtdb 2988c 2989 integer basao, basecp 2990 integer indao, indecp 2991 integer geom 2992 integer maxg, maxs, maxg1, maxs1, maxg2, maxs2 2993 integer hao, kao, hscr, kscr, h3p, k3p 2994 logical ignore 2995c 2996 integer icont, ucont 2997 integer atype, aprim, ie_a, ic_a, a_cent 2998 integer btype, bprim, ie_b, ic_b, b_cent 2999 integer ctype, cprim, ie_c, ic_c, c_cent 3000 integer ac_prim, bc_prim 3001 integer ir_c 3002 integer i,ir,icnt,j 3003 integer l2a, l2b, nint 3004 logical FF, FT 3005 double precision xyza(3),xyzb(3),xyzc(3),xyzall(3,3) 3006 double precision zan(3) 3007c 3008 integer maxprim 3009 parameter (maxprim=100) 3010 integer rc(maxprim) 3011 double precision ea(maxprim),ca(maxprim) 3012 double precision eb(maxprim),cb(maxprim) 3013 double precision ec(maxprim),cc(maxprim) 3014 double precision eac(maxprim),cac(maxprim) 3015 double precision ebc(maxprim),cbc(maxprim) 3016c 3017#include "bas_exndcf_sfn.fh" 3018#include "bas_ibs_sfn.fh" 3019c 3020 ignore = geom_create(geom,'geometry') 3021 ignore = ignore.and.geom_rtdb_load(rtdb,geom,'geometry') 3022 ignore = ignore.and.bas_create(basao,'ao basis') 3023 ignore = ignore.and.bas_rtdb_load(rtdb,geom,basao,'ao basis') 3024 ignore = ignore.and.int_normalize(rtdb,basao) 3025 ignore = ignore.and.geom_print(geom) 3026 ignore = ignore.and.bas_print(basao) 3027 ignore = ignore.and.gbs_map_print(basao) 3028 ignore = ignore.and.bas_get_ecp_handle(basao,basecp) 3029 ignore = ignore.and.ecp_print(basecp) 3030 if (.not.ignore) stop ' something failed' 3031c 3032 call int_init(rtdb,1,basao) 3033 call int_mem_1e(maxg1,maxs1) 3034 call int_mem_2e4c(maxg2,maxs2) 3035 maxg = max(maxg1,maxg2) 3036 maxs = max(maxs1,maxs2) 3037c 3038 ignore = 3039 & ma_push_get(mt_dbl,maxg,'aobuf',hao,kao) 3040 ignore = 3041 & ma_push_get(mt_dbl,maxg,'3pbuf',h3p,k3p) 3042 ignore=ignore.and. 3043 & ma_push_get(mt_dbl,maxs,'scr',hscr,kscr) 3044 3045 indao = basao + basis_handle_offset 3046 indecp = basecp + basis_handle_offset 3047c 3048 call dfill(maxprim,0.0d00,ea,1) 3049 call dfill(maxprim,0.0d00,ca,1) 3050 call dfill(maxprim,0.0d00,eb,1) 3051 call dfill(maxprim,0.0d00,cb,1) 3052 call dfill(maxprim,0.0d00,ec,1) 3053 call dfill(maxprim,0.0d00,cc,1) 3054 call dfill(maxprim,0.0d00,eac,1) 3055 call dfill(maxprim,0.0d00,cac,1) 3056 call dfill(maxprim,0.0d00,ebc,1) 3057 call dfill(maxprim,0.0d00,cbc,1) 3058 call ifill(maxprim,0,rc,1) 3059 call dfill(3,0.0d00,xyza,1) 3060 call dfill(3,0.0d00,xyzb,1) 3061 call dfill(3,0.0d00,xyzc,1) 3062 call dfill(3*3,0.0d00,xyzall,1) 3063 call dfill(3,-1.0d00,zan,1) 3064c 3065c for cont1|ecp|cont1 it fails 3066c 3067c fill a 3068 icont = 1 3069 ucont = (sf_ibs_cn2ucn(icont,indao)) 3070 atype = infbs_cont(Cont_Type,ucont,indao) 3071 aprim = infbs_cont(Cont_Nprim,ucont,indao) 3072 if (infbs_cont(Cont_Ngen,ucont,indao).ne.1) 3073 & stop ' a ngen wrong ' 3074 ie_a = infbs_cont(Cont_Iexp,ucont,indao) 3075 ic_a = infbs_cont(Cont_Icfp,ucont,indao) 3076 a_cent = sf_ibs_cn2ce(icont,indao) 3077 call dcopy(3,coords(1,a_cent,geom),1,xyza,1) 3078 if (aprim.le.maxprim) then 3079 call dcopy(aprim,dbl_mb(mb_exndcf(ie_a,indao)),1,ea,1) 3080 call dcopy(aprim,dbl_mb(mb_exndcf(ic_a,indao)),1,ca,1) 3081 write(6,*)' ea ' 3082 call output(ea,1,aprim,1,1,aprim,1,1) 3083 write(6,*)' ca ' 3084 call output(ca,1,aprim,1,1,aprim,1,1) 3085 else 3086 stop ' aprim error ' 3087 endif 3088c fill b 3089 icont = 1 3090 ucont = (sf_ibs_cn2ucn(icont,indao)) 3091 btype = infbs_cont(Cont_Type,ucont,indao) 3092 bprim = infbs_cont(Cont_Nprim,ucont,indao) 3093 if (infbs_cont(Cont_Ngen,ucont,indao).ne.1) 3094 & stop ' b ngen wrong ' 3095 ie_b = infbs_cont(Cont_Iexp,ucont,indao) 3096 ic_b = infbs_cont(Cont_Icfp,ucont,indao) 3097 b_cent = sf_ibs_cn2ce(icont,indao) 3098 call dcopy(3,coords(1,b_cent,geom),1,xyzb,1) 3099 if (bprim.le.maxprim) then 3100 call dcopy(bprim,dbl_mb(mb_exndcf(ie_b,indao)),1,eb,1) 3101 call dcopy(bprim,dbl_mb(mb_exndcf(ic_b,indao)),1,cb,1) 3102 write(6,*)' eb ' 3103 call output(eb,1,bprim,1,1,bprim,1,1) 3104 write(6,*)' cb ' 3105 call output(cb,1,bprim,1,1,bprim,1,1) 3106 else 3107 stop ' bprim error ' 3108 endif 3109c fill c 3110 icont = 1 3111 ucont = sf_ibs_cn2ucn(icont,indecp) 3112 ctype = infbs_cont(Cont_Type,ucont,indecp) 3113 cprim = infbs_cont(Cont_Nprim,ucont,indecp) 3114 if (infbs_cont(Cont_Ngen,ucont,indecp).ne.1) 3115 & stop 'c ngen wrong' 3116 ie_c = infbs_cont(Cont_iexp, ucont,indecp) 3117 ic_c = infbs_cont(Cont_icfp, ucont,indecp) 3118 ir_c = infbs_cont(Cont_irexp,ucont,indecp) 3119 c_cent = sf_ibs_cn2ce(icont,indecp) 3120 call dcopy(3,coords(1,c_cent,geom),1,xyzc,1) 3121c 3122 icnt = 0 3123 do i = 1,cprim 3124 ir = int(sf_exndcf((ir_c-1+i),indecp) + 0.00001d0) 3125 if (ir.eq.1) then 3126 if ((icnt+1).le.maxprim) then 3127 icnt = icnt + 1 3128 ec(icnt) = sf_exndcf((ie_c-1+i),indecp) 3129 cc(icnt) = sf_exndcf((ic_c-1+i),indecp) 3130 rc(icnt) = ir 3131 else 3132 stop 'cprim > maxprim' 3133 endif 3134 endif 3135 enddo 3136 cprim = icnt 3137 write(6,*)' ec ' 3138 call output(ec,1,cprim,1,1,cprim,1,1) 3139 write(6,*)' cc ' 3140 call output(cc,1,cprim,1,1,cprim,1,1) 3141 write(6,*)' rc ' 3142 do i = 1,cprim 3143 write(6,*)i,rc(i) 3144 enddo 3145c 3146c copy all coords 3147 call dcopy(3,xyza,1,xyzall(1,1),1) 3148 call dcopy(3,xyzb,1,xyzall(1,2),1) 3149 call dcopy(3,xyzc,1,xyzall(1,3),1) 3150c print coords 3151 write(6,*)' xyza ',xyza 3152 write(6,*)' xyzb ',xyzb 3153 write(6,*)' xyzc ',xyzc 3154 write(6,*)' xyz all' 3155 call output(xyzall,1,3,1,3,3,3,1) 3156c 3157c 3158 if (a_cent.eq.b_cent.and.a_cent.eq.c_cent) then 3159 write(6,*)' one center case' 3160 else 3161 stop ' not one center case' 3162 endif 3163c 3164c form ac exponents 3165 icnt = 0 3166 do i = 1,aprim 3167 do j = 1,cprim 3168 icnt = icnt + 1 3169 eac(icnt) = ea(i)+ec(j) 3170 cac(icnt) = ca(i)*cc(j) 3171 enddo 3172 enddo 3173 ac_prim = icnt 3174 write(6,*)' eac' 3175 call output(eac,1,ac_prim,1,1,ac_prim,1,1) 3176 write(6,*)' cac' 3177 call output(cac,1,ac_prim,1,1,ac_prim,1,1) 3178c form bc exponents 3179 icnt = 0 3180 do i = 1,bprim 3181 do j = 1,cprim 3182 icnt = icnt + 1 3183 ebc(icnt) = eb(i)+ec(j) 3184 cbc(icnt) = cb(i)*cc(j) 3185 enddo 3186 enddo 3187 bc_prim = icnt 3188 write(6,*)' ebc' 3189 call output(ebc,1,bc_prim,1,1,bc_prim,1,1) 3190 write(6,*)' cbc' 3191 call output(cbc,1,bc_prim,1,1,bc_prim,1,1) 3192c 3193 FF = .false. 3194 FT = .true. 3195 call dfill(maxg,0.0d00,dbl_mb(kao),1) 3196 call dfill(maxs,0.0d00,dbl_mb(kscr),1) 3197 call hf1( 3198 & xyza,eac,cac,ac_prim,1,atype, 3199 & xyzb,eb,cb,bprim,1,btype, 3200 & xyza,zan,1, 3201 & dbl_mb(kao),dbl_mb(kao),dbl_mb(kao),maxg,FF,FF,FT,FF, 3202 & FF,dbl_mb(kscr),maxs) 3203c 3204 call dfill(maxg,0.0d00,dbl_mb(k3p),1) 3205 call dfill(maxs,0.0d00,dbl_mb(kscr),1) 3206 call hf3pot( 3207 & xyza,ea,ca,aprim,1,atype, 3208 & xyzb,eb,cb,bprim,1,btype, 3209 & xyzc,ec,cc,cprim,1,0, 3210 & dbl_mb(k3p),maxg,nint, 3211 & FF,dbl_mb(kscr),maxs) 3212c 3213 l2a = (atype+1)*(atype+2)/2 3214 l2b = (btype+1)*(btype+2)/2 3215 write(6,*)'from hf1' 3216 call intintp(dbl_mb(kao),l2a,l2b,'from hf1') 3217 write(6,*)'from hf3pot' 3218 call intintp(dbl_mb(k3p),l2a,l2b,'from hf3pot') 3219c 3220 call int_terminate() 3221 ignore = bas_destroy(basao) 3222c 3223 ignore = ma_pop_stack(hscr) 3224 ignore = ma_pop_stack(h3p) 3225 ignore = ma_pop_stack(hao) 3226c 3227 end 3228 subroutine intintp(z,r,c,msg) 3229 implicit none 3230 integer r,c 3231 double precision z(r,c) 3232 character*(*) msg 3233c 3234 integer ir, ic 3235 do ir = 1,r 3236 do ic = 1,c 3237 if (abs(z(ir,ic)).gt.1.0d-10) then 3238 write(6,10000)ir,ic,z(ir,ic),msg 3239 endif 3240 enddo 3241 enddo 324210000 format(1x,'(',i5,',',i5,')',1x,1pd20.10,a) 3243 end 3244 subroutine raktest_geomprt(rtdb) 3245 implicit none 3246#include "errquit.fh" 3247#include "rtdb.fh" 3248#include "geom.fh" 3249#include "stdio.fh" 3250 integer rtdb 3251c 3252 integer geom 3253c 3254 if (.not.geom_create(geom,'geometry')) 3255 & call errquit('raktest_geomprt:create error',911, GEOM_ERR) 3256 if (.not.geom_rtdb_load(rtdb,geom,'geometry')) 3257 & call errquit('raktest_geomprt:load error',911, RTDB_ERR) 3258 if (.not.geom_print_distances(geom)) 3259 & call errquit('raktest_geomprt:print_distance error',911, 3260 & GEOM_ERR) 3261 if (.not.geom_print_angles(geom)) 3262 & call errquit('raktest_geomprt:print_angles error',911, 3263 & GEOM_ERR) 3264 if (.not.geom_print_dihedrals(geom)) 3265 & call errquit('raktest_geomprt:print_dihedrals error',911, 3266 & GEOM_ERR) 3267 if (.not.geom_destroy(geom)) 3268 & call errquit('raktest_geomprt:destory',911, 3269 & GEOM_ERR) 3270 end 3271 subroutine raktest_3cd(rtdb) 3272 implicit none 3273#include "errquit.fh" 3274#include "mafdecls.fh" 3275#include "rtdb.fh" 3276#include "geom.fh" 3277#include "bas.fh" 3278 integer rtdb 3279c 3280 integer geom, basis 3281 integer nat, nbf, lmax, bufsz, dbufsz, nbftrip, szscr 3282 integer h_bovp,h_bovm,h_dbuf,h_fd3ov,h_3ov,h_scr,h_c,h_cp,h_cm 3283 integer k_bovp,k_bovm,k_dbuf,k_fd3ov,k_3ov,k_scr,k_c,k_cp,k_cm 3284 integer h_fdp3ov, h_fdm3ov 3285 integer k_fdp3ov, k_fdm3ov 3286 integer stackleft 3287c 3288 logical int_normalize 3289 external int_normalize 3290c 3291* if (.not.rtdb_print(rtdb,.true.)) stop ' rtdb_p 1?' 3292c 3293 if (.not.geom_create(geom,'geometry')) 3294 & call errquit('raktest_geomwrt: geom_create failed?',911, 3295 & GEOM_ERR) 3296 if (.not. geom_rtdb_load(rtdb,geom,'geometry')) call errquit 3297 & ('raktest_geomwrt: geom_rtdb_load -ref failed',911, 3298 & RTDB_ERR) 3299 if (.not.geom_print(geom)) stop ' print error' 3300c 3301 if (.not.bas_create(basis,'ao basis')) call errquit 3302 & ('bas_create failed',911, BASIS_ERR) 3303 if (.not.bas_rtdb_load(rtdb,geom,basis,'ao basis')) 3304 & call errquit 3305 & ('bas_rtdb_load failed',911, RTDB_ERR) 3306 if (.not.bas_print(basis)) stop ' print error' 3307 if (.not.rtdb_print(rtdb,.false.)) stop ' rtdb_p 2?' 3308 if (.not.int_normalize(rtdb,basis)) stop ' norm error 1' 3309 call int_init(rtdb,1,basis) 3310c 3311 if (.not.geom_ncent(geom,nat)) stop 'geom_ncent fe' 3312 if (.not.bas_numbf(basis,nbf)) stop 'bas_numbf fe' 3313 if (.not.bas_high_angular(basis,Lmax)) stop 'bas_ha fe' 3314c 3315 bufsz = (Lmax+1)*(Lmax+2)/2 3316 bufsz = bufsz*bufsz*bufsz 3317 dbufsz = 9*bufsz 3318 nbftrip = nbf*nbf*nbf 3319 szscr = 60000 3320 if (MA_verify_allocator_stuff()) write(6,*)' maok (0)' 3321 stackleft = ma_inquire_stack(mt_dbl) 3322 write(6,*)' <0> stack left ',stackleft, 3323 & ' next:',bufsz 3324 if (.not.ma_push_get(mt_dbl,bufsz,' buf ovlap +', 3325 & h_bovp,k_bovp)) stop ' ma_pg 1 fail ' 3326 if (MA_verify_allocator_stuff()) write(6,*)' maok (1)' 3327 stackleft = ma_inquire_stack(mt_dbl) 3328 write(6,*)' <1> stack left ',stackleft, 3329 & ' next:',bufsz 3330 if (.not.ma_push_get(mt_dbl,bufsz,' buf ovlap -', 3331 & h_bovm,k_bovm)) stop ' ma_pg 2 fail ' 3332 if (MA_verify_allocator_stuff()) write(6,*)' maok (2)' 3333 stackleft = ma_inquire_stack(mt_dbl) 3334 write(6,*)' <2> stack left ',stackleft, 3335 & ' next:',dbufsz 3336 if (.not.ma_push_get(mt_dbl,dbufsz,'deriv buf ', 3337 & h_dbuf,k_dbuf)) stop ' ma_pg 3 fail ' 3338 if (MA_verify_allocator_stuff()) write(6,*)' maok (3)' 3339 stackleft = ma_inquire_stack(mt_dbl) 3340 write(6,*)' <3> stack left ',stackleft, 3341 & ' next:',(nbftrip*nat*3) 3342 if (.not.ma_push_get(mt_dbl,(nbftrip*nat*3),' fd 3ov matrix ', 3343 & h_fd3ov, k_fd3ov)) stop ' ma_pg 4 fail ' 3344 if (MA_verify_allocator_stuff()) write(6,*)' maok (4)' 3345 stackleft = ma_inquire_stack(mt_dbl) 3346 write(6,*)' <4> stack left ',stackleft, 3347 & ' next:',(nbftrip*nat*3) 3348 if (.not.ma_push_get(mt_dbl,(nbftrip*nat*3),' fd + 3ov matrix ', 3349 & h_fdp3ov, k_fdp3ov)) stop ' ma_pg 5 fail ' 3350 if (MA_verify_allocator_stuff()) write(6,*)' maok (5)' 3351 stackleft = ma_inquire_stack(mt_dbl) 3352 write(6,*)' <5> stack left ',stackleft, 3353 & ' next:',(nbftrip*nat*3) 3354 if (.not.ma_push_get(mt_dbl,(nbftrip*nat*3),' fd - 3ov matrix ', 3355 & h_fdm3ov, k_fdm3ov)) stop ' ma_pg 6 fail ' 3356 if (MA_verify_allocator_stuff()) write(6,*)' maok (6)' 3357 stackleft = ma_inquire_stack(mt_dbl) 3358 write(6,*)' <6> stack left ',stackleft, 3359 & ' next:',(nbftrip*nat*3) 3360 if (.not.ma_push_get(mt_dbl,(nbftrip*nat*3),' 3ov matrix ', 3361 & h_3ov, k_3ov)) stop ' ma_pg 7 fail ' 3362 if (MA_verify_allocator_stuff()) write(6,*)' maok (7)' 3363 stackleft = ma_inquire_stack(mt_dbl) 3364 write(6,*)' <7> stack left ',stackleft, 3365 & ' next:',szscr 3366 if (.not.ma_push_get(mt_dbl,szscr,' scratch buffer', 3367 & h_scr, k_scr)) stop ' ma_pg 8 fail ' 3368 if (MA_verify_allocator_stuff()) write(6,*)' maok (8)' 3369 stackleft = ma_inquire_stack(mt_dbl) 3370 write(6,*)' <8> stack left ',stackleft, 3371 & ' next:',(3*nat) 3372 if (.not.ma_push_get(mt_dbl,(3*nat),' coords ', 3373 & h_c, k_c)) stop ' ma_pg 9 fail ' 3374 if (MA_verify_allocator_stuff()) write(6,*)' maok (9)' 3375 stackleft = ma_inquire_stack(mt_dbl) 3376 write(6,*)' <9> stack left ',stackleft, 3377 & ' next:',(3*nat) 3378 if (.not.ma_push_get(mt_dbl,(3*nat),' coords +', 3379 & h_cp, k_cp)) stop ' ma_pg 10 fail ' 3380 if (MA_verify_allocator_stuff()) write(6,*)' maok (10)' 3381 stackleft = ma_inquire_stack(mt_dbl) 3382 write(6,*)' <10> stack left ',stackleft, 3383 & ' next:',(3*nat) 3384 if (MA_verify_allocator_stuff()) write(6,*)' maok (11)' 3385 if (.not.ma_push_get(mt_dbl,(3*nat),' coords -', 3386 & h_cm, k_cm)) stop ' ma_pg 11 fail ' 3387 stackleft = ma_inquire_stack(mt_dbl) 3388 write(6,*)' <11> stack left ',stackleft 3389c 3390 call raktest_3cd_1(geom,basis,nbf,nat,szscr,bufsz,dbufsz, 3391 & dbl_mb(k_bovp), dbl_mb(k_bovm), dbl_mb(k_dbuf), 3392 & dbl_mb(k_fd3ov), dbl_mb(k_fdp3ov), dbl_mb(k_fdm3ov), 3393 & dbl_mb(k_3ov), dbl_mb(k_scr), 3394 & dbl_mb(k_c), dbl_mb(k_cp), dbl_mb(k_cm)) 3395c 3396 if (.not.ma_pop_stack(h_cm)) stop ' ma pop fail' 3397 if (.not.ma_pop_stack(h_cp)) stop ' ma pop fail' 3398 if (.not.ma_pop_stack(h_c)) stop ' ma pop fail' 3399 if (.not.ma_pop_stack(h_scr)) stop ' ma pop fail' 3400 if (.not.ma_pop_stack(h_3ov)) stop ' ma pop fail' 3401 if (.not.ma_pop_stack(h_fdm3ov)) stop ' ma pop fail' 3402 if (.not.ma_pop_stack(h_fdp3ov)) stop ' ma pop fail' 3403 if (.not.ma_pop_stack(h_fd3ov)) stop ' ma pop fail' 3404 if (.not.ma_pop_stack(h_dbuf)) stop ' ma pop fail' 3405 if (.not.ma_pop_stack(h_bovm)) stop ' ma pop fail' 3406 if (.not.ma_pop_stack(h_bovp)) stop ' ma pop fail' 3407c 3408 call int_terminate 3409 if (.not.bas_destroy(basis)) stop ' bas_dest error 1' 3410 if (.not.geom_destroy(geom)) stop ' geom_dest error 1' 3411 if (.not.bas_version()) stop ' bas_version error' 3412 end 3413 subroutine raktest_3cd_1(geom,basisin,nbf,nat,szscr,bufsz,dbufsz, 3414 & bufp, bufm, dbuf, fdo, fdp, fdm, o, scr, c, cp, cm) 3415 implicit none 3416#include "mafdecls.fh" 3417#include "bas.fh" 3418#include "nwc_const.fh" 3419#include "basP.fh" 3420#include "basdeclsP.fh" 3421#include "geomP.fh" 3422#include "geobasmapP.fh" 3423#include "bas_exndcf_dec.fh" 3424#include "bas_ibs_dec.fh" 3425#include "int_nbf.fh" 3426 integer geom 3427 integer basisin 3428 integer nbf 3429 integer nat 3430 integer szscr 3431 integer bufsz 3432 integer dbufsz 3433 double precision bufp(bufsz), bufm(bufsz), dbuf(dbufsz) 3434 double precision fdo(nbf,nbf,nbf,3,nat) 3435 double precision fdp(nbf,nbf,nbf,3,nat) 3436 double precision fdm(nbf,nbf,nbf,3,nat) 3437 double precision o(nbf,nbf,nbf,3,nat) 3438 double precision scr(szscr) 3439 double precision c(3,nat), cp(3,nat), cm(3,nat) 3440c 3441 double precision delta, delta2i, norm 3442 double precision ddot 3443 external ddot 3444 integer ucont 3445 integer ish, ilow, ihi, Li, i_prim, i_gen, i_iexp, i_icfp, i_cent 3446 integer jsh, jlow, jhi, Lj, j_prim, j_gen, j_iexp, j_icfp, j_cent 3447 integer ksh, klow, khi, Lk, k_prim, k_gen, k_iexp, k_icfp, k_cent 3448 integer ibf, jbf, kbf 3449 integer bs 3450 integer nshell, nint, count 3451 integer zcent(3), ixyz, zc 3452 integer pass 3453c 3454#include "bas_exndcf_sfn.fh" 3455#include "bas_ibs_sfn.fh" 3456c 3457 call dfill((3*nat*nbf*nbf*nbf),0.0d00,fdm,1) 3458 call dfill((3*nat*nbf*nbf*nbf),0.0d00,fdp,1) 3459 call dfill((3*nat*nbf*nbf*nbf),0.0d00,fdo,1) 3460 call dfill((3*nat*nbf*nbf*nbf),0.0d00, o,1) 3461 delta = 0.00001d00 3462 delta2i = 1.0d00/(2.0d00*delta) 3463 write(6,*)' nbf = ',nbf 3464 write(6,*)' nat = ',nat 3465 write(6,*)' delta,2inverse ',delta,delta2i 3466 bs = basisin + Basis_Handle_Offset 3467 nshell = ncont_tot_gb(bs) 3468 write(6,*)' nshell',nshell 3469c 3470 call dcopy (3*nat,coords(1,1,geom),1,c,1) 3471c 3472 pass = 0 3473 do ish = 1,nshell 3474 if (.not.bas_cn2bfr(basisin,ish,ilow,ihi)) stop 'cn2bfr i' 3475 ucont = (sf_ibs_cn2ucn(ish,bs)) 3476 Li = infbs_cont(CONT_TYPE ,ucont,bs) 3477 i_prim = infbs_cont(CONT_NPRIM,ucont,bs) 3478 i_gen = infbs_cont(CONT_NGEN ,ucont,bs) 3479 i_iexp = infbs_cont(CONT_IEXP ,ucont,bs) 3480 i_icfp = infbs_cont(CONT_ICFP ,ucont,bs) 3481 i_cent = (sf_ibs_cn2ce(ish,bs)) 3482 do jsh = 1,nshell 3483 if (.not.bas_cn2bfr(basisin,jsh,jlow,jhi)) stop 'cn2bfr j' 3484 ucont = (sf_ibs_cn2ucn(jsh,bs)) 3485 Lj = infbs_cont(CONT_TYPE ,ucont,bs) 3486 j_prim = infbs_cont(CONT_NPRIM,ucont,bs) 3487 j_gen = infbs_cont(CONT_NGEN ,ucont,bs) 3488 j_iexp = infbs_cont(CONT_IEXP ,ucont,bs) 3489 j_icfp = infbs_cont(CONT_ICFP ,ucont,bs) 3490 j_cent = (sf_ibs_cn2ce(jsh,bs)) 3491 do ksh = 1,nshell 3492 if (.not.bas_cn2bfr(basisin,ksh,klow,khi)) stop 'cn2bfr k' 3493 ucont = (sf_ibs_cn2ucn(ksh,bs)) 3494 Lk = infbs_cont(CONT_TYPE ,ucont,bs) 3495 k_prim = infbs_cont(CONT_NPRIM,ucont,bs) 3496 k_gen = infbs_cont(CONT_NGEN ,ucont,bs) 3497 k_iexp = infbs_cont(CONT_IEXP ,ucont,bs) 3498 k_icfp = infbs_cont(CONT_ICFP ,ucont,bs) 3499 k_cent = (sf_ibs_cn2ce(ksh,bs)) 3500 pass = pass + 1 3501 write(6,*)' pass ',pass 3502 if (i_cent.eq.j_cent.and.j_cent.eq.k_cent) goto 00100 3503c 3504 nint = int_nbf_x(Li)*int_nbf_x(Lj)*int_nbf_x(Lk) 3505 write(6,*)' nint = ',nint 3506c 3507c* icenter x + 3508 call dcopy(nat*3,c,1,cp,1) 3509 cp(1,i_cent) = cp(1,i_cent)+delta 3510 call hf3ois( 3511 & cp(1,i_cent), dbl_mb(mb_exndcf(i_iexp,bs)), 3512 & dbl_mb(mb_exndcf(i_icfp,bs)),i_prim, Li, 3513 & cp(1,j_cent), dbl_mb(mb_exndcf(j_iexp,bs)), 3514 & dbl_mb(mb_exndcf(j_icfp,bs)),j_prim, Lj, 3515 & cp(1,k_cent), dbl_mb(mb_exndcf(k_iexp,bs)), 3516 & dbl_mb(mb_exndcf(k_icfp,bs)),k_prim, Lk, 3517 & bufp,nint,.false.,.false.,scr,szscr) 3518c* icenter x - 3519 call dcopy(nat*3,c,1,cm,1) 3520 cm(1,i_cent) = cm(1,i_cent)-delta 3521 call hf3ois( 3522 & cm(1,i_cent), dbl_mb(mb_exndcf(i_iexp,bs)), 3523 & dbl_mb(mb_exndcf(i_icfp,bs)),i_prim, Li, 3524 & cm(1,j_cent), dbl_mb(mb_exndcf(j_iexp,bs)), 3525 & dbl_mb(mb_exndcf(j_icfp,bs)),j_prim, Lj, 3526 & cm(1,k_cent), dbl_mb(mb_exndcf(k_iexp,bs)), 3527 & dbl_mb(mb_exndcf(k_icfp,bs)),k_prim, Lk, 3528 & bufm,nint,.false.,.false.,scr,szscr) 3529 count = 0 3530 do ibf = ilow,ihi 3531 do jbf = jlow,jhi 3532 do kbf = klow,khi 3533 count = count + 1 3534 fdp(ibf,jbf,kbf,1,i_cent) = bufp(count) 3535 fdm(ibf,jbf,kbf,1,i_cent) = bufm(count) 3536 enddo 3537 enddo 3538 enddo 3539c 3540c* icenter y + 3541 call dcopy(nat*3,c,1,cp,1) 3542 cp(2,i_cent) = cp(2,i_cent)+delta 3543 call hf3ois( 3544 & cp(1,i_cent), dbl_mb(mb_exndcf(i_iexp,bs)), 3545 & dbl_mb(mb_exndcf(i_icfp,bs)),i_prim, Li, 3546 & cp(1,j_cent), dbl_mb(mb_exndcf(j_iexp,bs)), 3547 & dbl_mb(mb_exndcf(j_icfp,bs)),j_prim, Lj, 3548 & cp(1,k_cent), dbl_mb(mb_exndcf(k_iexp,bs)), 3549 & dbl_mb(mb_exndcf(k_icfp,bs)),k_prim, Lk, 3550 & bufp,nint,.false.,.false.,scr,szscr) 3551c* icenter y - 3552 call dcopy(nat*3,c,1,cm,1) 3553 cm(2,i_cent) = cm(2,i_cent)-delta 3554 call hf3ois( 3555 & cm(1,i_cent), dbl_mb(mb_exndcf(i_iexp,bs)), 3556 & dbl_mb(mb_exndcf(i_icfp,bs)),i_prim, Li, 3557 & cm(1,j_cent), dbl_mb(mb_exndcf(j_iexp,bs)), 3558 & dbl_mb(mb_exndcf(j_icfp,bs)),j_prim, Lj, 3559 & cm(1,k_cent), dbl_mb(mb_exndcf(k_iexp,bs)), 3560 & dbl_mb(mb_exndcf(k_icfp,bs)),k_prim, Lk, 3561 & bufm,nint,.false.,.false.,scr,szscr) 3562 count = 0 3563 do ibf = ilow,ihi 3564 do jbf = jlow,jhi 3565 do kbf = klow,khi 3566 count = count + 1 3567 fdp(ibf,jbf,kbf,2,i_cent) = bufp(count) 3568 fdm(ibf,jbf,kbf,2,i_cent) = bufm(count) 3569 enddo 3570 enddo 3571 enddo 3572c 3573c 3574c* icenter z + 3575 call dcopy(nat*3,c,1,cp,1) 3576 cp(3,i_cent) = cp(3,i_cent)+delta 3577 call hf3ois( 3578 & cp(1,i_cent), dbl_mb(mb_exndcf(i_iexp,bs)), 3579 & dbl_mb(mb_exndcf(i_icfp,bs)),i_prim, Li, 3580 & cp(1,j_cent), dbl_mb(mb_exndcf(j_iexp,bs)), 3581 & dbl_mb(mb_exndcf(j_icfp,bs)),j_prim, Lj, 3582 & cp(1,k_cent), dbl_mb(mb_exndcf(k_iexp,bs)), 3583 & dbl_mb(mb_exndcf(k_icfp,bs)),k_prim, Lk, 3584 & bufp,nint,.false.,.false.,scr,szscr) 3585c* icenter z - 3586 call dcopy(nat*3,c,1,cm,1) 3587 cm(3,i_cent) = cm(3,i_cent)-delta 3588 call hf3ois( 3589 & cm(1,i_cent), dbl_mb(mb_exndcf(i_iexp,bs)), 3590 & dbl_mb(mb_exndcf(i_icfp,bs)),i_prim, Li, 3591 & cm(1,j_cent), dbl_mb(mb_exndcf(j_iexp,bs)), 3592 & dbl_mb(mb_exndcf(j_icfp,bs)),j_prim, Lj, 3593 & cm(1,k_cent), dbl_mb(mb_exndcf(k_iexp,bs)), 3594 & dbl_mb(mb_exndcf(k_icfp,bs)),k_prim, Lk, 3595 & bufm,nint,.false.,.false.,scr,szscr) 3596 count = 0 3597 do ibf = ilow,ihi 3598 do jbf = jlow,jhi 3599 do kbf = klow,khi 3600 count = count + 1 3601 fdp(ibf,jbf,kbf,3,i_cent) = bufp(count) 3602 fdm(ibf,jbf,kbf,3,i_cent) = bufm(count) 3603 enddo 3604 enddo 3605 enddo 3606c 3607c 3608c* jcenter x + 3609 call dcopy(nat*3,c,1,cp,1) 3610 cp(1,j_cent) = cp(1,j_cent)+delta 3611 call hf3ois( 3612 & cp(1,i_cent), dbl_mb(mb_exndcf(i_iexp,bs)), 3613 & dbl_mb(mb_exndcf(i_icfp,bs)),i_prim, Li, 3614 & cp(1,j_cent), dbl_mb(mb_exndcf(j_iexp,bs)), 3615 & dbl_mb(mb_exndcf(j_icfp,bs)),j_prim, Lj, 3616 & cp(1,k_cent), dbl_mb(mb_exndcf(k_iexp,bs)), 3617 & dbl_mb(mb_exndcf(k_icfp,bs)),k_prim, Lk, 3618 & bufp,nint,.false.,.false.,scr,szscr) 3619c* jcenter x - 3620 call dcopy(nat*3,c,1,cm,1) 3621 cm(1,j_cent) = cm(1,j_cent)-delta 3622 call hf3ois( 3623 & cm(1,i_cent), dbl_mb(mb_exndcf(i_iexp,bs)), 3624 & dbl_mb(mb_exndcf(i_icfp,bs)),i_prim, Li, 3625 & cm(1,j_cent), dbl_mb(mb_exndcf(j_iexp,bs)), 3626 & dbl_mb(mb_exndcf(j_icfp,bs)),j_prim, Lj, 3627 & cm(1,k_cent), dbl_mb(mb_exndcf(k_iexp,bs)), 3628 & dbl_mb(mb_exndcf(k_icfp,bs)),k_prim, Lk, 3629 & bufm,nint,.false.,.false.,scr,szscr) 3630 count = 0 3631 do ibf = ilow,ihi 3632 do jbf = jlow,jhi 3633 do kbf = klow,khi 3634 count = count + 1 3635 fdp(ibf,jbf,kbf,1,j_cent) = bufp(count) 3636 fdm(ibf,jbf,kbf,1,j_cent) = bufm(count) 3637 enddo 3638 enddo 3639 enddo 3640c 3641c* jcenter y + 3642 call dcopy(nat*3,c,1,cp,1) 3643 cp(2,j_cent) = cp(2,j_cent)+delta 3644 call hf3ois( 3645 & cp(1,i_cent), dbl_mb(mb_exndcf(i_iexp,bs)), 3646 & dbl_mb(mb_exndcf(i_icfp,bs)),i_prim, Li, 3647 & cp(1,j_cent), dbl_mb(mb_exndcf(j_iexp,bs)), 3648 & dbl_mb(mb_exndcf(j_icfp,bs)),j_prim, Lj, 3649 & cp(1,k_cent), dbl_mb(mb_exndcf(k_iexp,bs)), 3650 & dbl_mb(mb_exndcf(k_icfp,bs)),k_prim, Lk, 3651 & bufp,nint,.false.,.false.,scr,szscr) 3652c* jcenter y - 3653 call dcopy(nat*3,c,1,cm,1) 3654 cm(2,j_cent) = cm(2,j_cent)-delta 3655 call hf3ois( 3656 & cm(1,i_cent), dbl_mb(mb_exndcf(i_iexp,bs)), 3657 & dbl_mb(mb_exndcf(i_icfp,bs)),i_prim, Li, 3658 & cm(1,j_cent), dbl_mb(mb_exndcf(j_iexp,bs)), 3659 & dbl_mb(mb_exndcf(j_icfp,bs)),j_prim, Lj, 3660 & cm(1,k_cent), dbl_mb(mb_exndcf(k_iexp,bs)), 3661 & dbl_mb(mb_exndcf(k_icfp,bs)),k_prim, Lk, 3662 & bufm,nint,.false.,.false.,scr,szscr) 3663 count = 0 3664 do ibf = ilow,ihi 3665 do jbf = jlow,jhi 3666 do kbf = klow,khi 3667 count = count + 1 3668 fdp(ibf,jbf,kbf,2,j_cent) = bufp(count) 3669 fdm(ibf,jbf,kbf,2,j_cent) = bufm(count) 3670 enddo 3671 enddo 3672 enddo 3673c 3674c 3675c* jcenter z + 3676 call dcopy(nat*3,c,1,cp,1) 3677 cp(3,j_cent) = cp(3,j_cent)+delta 3678 call hf3ois( 3679 & cp(1,i_cent), dbl_mb(mb_exndcf(i_iexp,bs)), 3680 & dbl_mb(mb_exndcf(i_icfp,bs)),i_prim, Li, 3681 & cp(1,j_cent), dbl_mb(mb_exndcf(j_iexp,bs)), 3682 & dbl_mb(mb_exndcf(j_icfp,bs)),j_prim, Lj, 3683 & cp(1,k_cent), dbl_mb(mb_exndcf(k_iexp,bs)), 3684 & dbl_mb(mb_exndcf(k_icfp,bs)),k_prim, Lk, 3685 & bufp,nint,.false.,.false.,scr,szscr) 3686c* jcenter z - 3687 call dcopy(nat*3,c,1,cm,1) 3688 cm(3,j_cent) = cm(3,j_cent)-delta 3689 call hf3ois( 3690 & cm(1,i_cent), dbl_mb(mb_exndcf(i_iexp,bs)), 3691 & dbl_mb(mb_exndcf(i_icfp,bs)),i_prim, Li, 3692 & cm(1,j_cent), dbl_mb(mb_exndcf(j_iexp,bs)), 3693 & dbl_mb(mb_exndcf(j_icfp,bs)),j_prim, Lj, 3694 & cm(1,k_cent), dbl_mb(mb_exndcf(k_iexp,bs)), 3695 & dbl_mb(mb_exndcf(k_icfp,bs)),k_prim, Lk, 3696 & bufm,nint,.false.,.false.,scr,szscr) 3697 count = 0 3698 do ibf = ilow,ihi 3699 do jbf = jlow,jhi 3700 do kbf = klow,khi 3701 count = count + 1 3702 fdp(ibf,jbf,kbf,3,j_cent) = bufp(count) 3703 fdm(ibf,jbf,kbf,3,j_cent) = bufm(count) 3704 enddo 3705 enddo 3706 enddo 3707c 3708c 3709c* kcenter x + 3710 call dcopy(nat*3,c,1,cp,1) 3711 cp(1,k_cent) = cp(1,k_cent)+delta 3712 call hf3ois( 3713 & cp(1,i_cent), dbl_mb(mb_exndcf(i_iexp,bs)), 3714 & dbl_mb(mb_exndcf(i_icfp,bs)),i_prim, Li, 3715 & cp(1,j_cent), dbl_mb(mb_exndcf(j_iexp,bs)), 3716 & dbl_mb(mb_exndcf(j_icfp,bs)),j_prim, Lj, 3717 & cp(1,k_cent), dbl_mb(mb_exndcf(k_iexp,bs)), 3718 & dbl_mb(mb_exndcf(k_icfp,bs)),k_prim, Lk, 3719 & bufp,nint,.false.,.false.,scr,szscr) 3720c* kcenter x - 3721 call dcopy(nat*3,c,1,cm,1) 3722 cm(1,k_cent) = cm(1,k_cent)-delta 3723 call hf3ois( 3724 & cm(1,i_cent), dbl_mb(mb_exndcf(i_iexp,bs)), 3725 & dbl_mb(mb_exndcf(i_icfp,bs)),i_prim, Li, 3726 & cm(1,j_cent), dbl_mb(mb_exndcf(j_iexp,bs)), 3727 & dbl_mb(mb_exndcf(j_icfp,bs)),j_prim, Lj, 3728 & cm(1,k_cent), dbl_mb(mb_exndcf(k_iexp,bs)), 3729 & dbl_mb(mb_exndcf(k_icfp,bs)),k_prim, Lk, 3730 & bufm,nint,.false.,.false.,scr,szscr) 3731 count = 0 3732 do ibf = ilow,ihi 3733 do jbf = jlow,jhi 3734 do kbf = klow,khi 3735 count = count + 1 3736 fdp(ibf,jbf,kbf,1,k_cent) = bufp(count) 3737 fdm(ibf,jbf,kbf,1,k_cent) = bufm(count) 3738 enddo 3739 enddo 3740 enddo 3741c 3742c* kcenter y + 3743 call dcopy(nat*3,c,1,cp,1) 3744 cp(2,k_cent) = cp(2,k_cent)+delta 3745 call hf3ois( 3746 & cp(1,i_cent), dbl_mb(mb_exndcf(i_iexp,bs)), 3747 & dbl_mb(mb_exndcf(i_icfp,bs)),i_prim, Li, 3748 & cp(1,j_cent), dbl_mb(mb_exndcf(j_iexp,bs)), 3749 & dbl_mb(mb_exndcf(j_icfp,bs)),j_prim, Lj, 3750 & cp(1,k_cent), dbl_mb(mb_exndcf(k_iexp,bs)), 3751 & dbl_mb(mb_exndcf(k_icfp,bs)),k_prim, Lk, 3752 & bufp,nint,.false.,.false.,scr,szscr) 3753c* kcenter y - 3754 call dcopy(nat*3,c,1,cm,1) 3755 cm(2,k_cent) = cm(2,k_cent)-delta 3756 call hf3ois( 3757 & cm(1,i_cent), dbl_mb(mb_exndcf(i_iexp,bs)), 3758 & dbl_mb(mb_exndcf(i_icfp,bs)),i_prim, Li, 3759 & cm(1,j_cent), dbl_mb(mb_exndcf(j_iexp,bs)), 3760 & dbl_mb(mb_exndcf(j_icfp,bs)),j_prim, Lj, 3761 & cm(1,k_cent), dbl_mb(mb_exndcf(k_iexp,bs)), 3762 & dbl_mb(mb_exndcf(k_icfp,bs)),k_prim, Lk, 3763 & bufm,nint,.false.,.false.,scr,szscr) 3764 count = 0 3765 do ibf = ilow,ihi 3766 do jbf = jlow,jhi 3767 do kbf = klow,khi 3768 count = count + 1 3769 fdp(ibf,jbf,kbf,2,k_cent) = bufp(count) 3770 fdm(ibf,jbf,kbf,2,k_cent) = bufm(count) 3771 enddo 3772 enddo 3773 enddo 3774c 3775c 3776c* kcenter z + 3777 call dcopy(nat*3,c,1,cp,1) 3778 cp(3,k_cent) = cp(3,k_cent)+delta 3779 call hf3ois( 3780 & cp(1,i_cent), dbl_mb(mb_exndcf(i_iexp,bs)), 3781 & dbl_mb(mb_exndcf(i_icfp,bs)),i_prim, Li, 3782 & cp(1,j_cent), dbl_mb(mb_exndcf(j_iexp,bs)), 3783 & dbl_mb(mb_exndcf(j_icfp,bs)),j_prim, Lj, 3784 & cp(1,k_cent), dbl_mb(mb_exndcf(k_iexp,bs)), 3785 & dbl_mb(mb_exndcf(k_icfp,bs)),k_prim, Lk, 3786 & bufp,nint,.false.,.false.,scr,szscr) 3787c* kcenter z - 3788 call dcopy(nat*3,c,1,cm,1) 3789 cm(3,k_cent) = cm(3,k_cent)-delta 3790 call hf3ois( 3791 & cm(1,i_cent), dbl_mb(mb_exndcf(i_iexp,bs)), 3792 & dbl_mb(mb_exndcf(i_icfp,bs)),i_prim, Li, 3793 & cm(1,j_cent), dbl_mb(mb_exndcf(j_iexp,bs)), 3794 & dbl_mb(mb_exndcf(j_icfp,bs)),j_prim, Lj, 3795 & cm(1,k_cent), dbl_mb(mb_exndcf(k_iexp,bs)), 3796 & dbl_mb(mb_exndcf(k_icfp,bs)),k_prim, Lk, 3797 & bufm,nint,.false.,.false.,scr,szscr) 3798 count = 0 3799 do ibf = ilow,ihi 3800 do jbf = jlow,jhi 3801 do kbf = klow,khi 3802 count = count + 1 3803 fdp(ibf,jbf,kbf,3,k_cent) = bufp(count) 3804 fdm(ibf,jbf,kbf,3,k_cent) = bufm(count) 3805 enddo 3806 enddo 3807 enddo 3808c 3809 call hfd3ois(c(1,i_cent),dbl_mb(mb_exndcf(i_iexp,bs)), 3810 & dbl_mb(mb_exndcf(i_icfp,bs)),i_prim, 1, Li, 3811 & c(1,j_cent),dbl_mb(mb_exndcf(j_iexp,bs)), 3812 & dbl_mb(mb_exndcf(j_icfp,bs)),j_prim, 1, Lj, 3813 & c(1,k_cent),dbl_mb(mb_exndcf(k_iexp,bs)), 3814 & dbl_mb(mb_exndcf(k_icfp,bs)),k_prim, 1, Lk, 3815 & dbuf,nint,.false.,scr,szscr) 3816 call raktest_3cd_2(i_cent,j_cent,k_cent,nint,dbuf,zcent) 3817c 3818 count = 0 3819 do zc = 1,3 3820 if (zcent(zc).ne.0) then 3821 do ixyz = 1,3 3822 do ibf = ilow,ihi 3823 do jbf = jlow,jhi 3824 do kbf = klow,khi 3825 count = count + 1 3826 o(ibf,jbf,kbf,ixyz,zcent(zc)) = dbuf(count) 3827 enddo 3828 enddo 3829 enddo 3830 enddo 3831 endif 3832 enddo 3833c 383400100 continue 3835c 3836 enddo 3837 enddo 3838 enddo 3839c 3840 call dcopy((nbf*nbf*nbf*3*nat),fdp,1,fdo,1) 3841 call daxpy((nbf*nbf*nbf*3*nat),-1.0d00,fdm,1,fdo,1) 3842 call dscal((nbf*nbf*nbf*3*nat),delta2i,fdo,1) 3843 call printboth_3cd(nbf,nat,fdo,o) 3844 call dcopy((nbf*nbf*nbf*3*nat),fdo,1,fdm,1) 3845 call daxpy((nbf*nbf*nbf*3*nat),-1.0d00,o,1,fdm,1) 3846 norm = ddot((nbf*nbf*nbf*3*nat),fdm,1,fdm,1) 3847c 3848 write(6,*)' difference norm: ',norm 3849c 3850 end 3851 subroutine raktest_3cd_2(ic,jc,kc,nint,dbuf,idcent) 3852 implicit none 3853 integer ic, jc, kc, nint 3854 integer idcent(3) 3855 double precision dbuf(nint,9) 3856c 3857 if ((ic.ne.jc).and.(ic.ne.kc).and.(jc.ne.kc)) then 3858 idcent(1) = ic 3859 idcent(2) = jc 3860 idcent(3) = kc 3861 goto 90000 3862 endif 3863 if (ic.eq.jc) then 3864 call daxpy(3*nint,1.0d00,dbuf(1,4),1,dbuf(1,1),1) 3865 call dcopy(3*nint,dbuf(1,7),1,dbuf(1,4),1) 3866 call dfill(3*nint,0.0d00,dbuf(1,7),1) 3867 write(6,*)' moving' 3868 idcent(1) = ic 3869 idcent(2) = kc 3870 idcent(3) = 0 3871 else if (ic.eq.kc) then 3872 call daxpy(3*nint,1.0d00,dbuf(1,7),1,dbuf(1,1),1) 3873 call dfill(3*nint,0.0d00,dbuf(1,7),1) 3874 idcent(1) = ic 3875 idcent(2) = jc 3876 idcent(3) = 0 3877 else if (jc.eq.kc) then 3878 call daxpy(3*nint,1.0d00,dbuf(1,7),1,dbuf(1,4),1) 3879 call dfill(3*nint,0.0d00,dbuf(1,7),1) 3880 idcent(1) = ic 3881 idcent(2) = jc 3882 idcent(3) = 0 3883 else 3884 write(6,*)ic,jc,kc 3885 stop ' how did I get here' 3886 endif 388790000 continue 3888 write(6,*)'idcent',idcent 3889 end 3890 subroutine printboth_3cd(nbf,nat,fdo,o) 3891 implicit none 3892 integer nbf 3893 integer nat 3894 double precision fdo(nbf,nbf,nbf,3,nat) 3895 double precision o(nbf,nbf,nbf,3,nat) 3896c 3897 double precision diff, thresh 3898 integer i,j,k,ixyz,ic, count 3899c 3900 thresh = 1.0d-06 3901* thresh = -1.0d00 3902c 3903 count = 0 3904 do ic = 1,nat 3905 do ixyz=1,3 3906 do k=1,nbf 3907 do j=1,nbf 3908 do i=1,nbf 3909 count = count + 1 3910 diff = fdo(i,j,k,ixyz,ic) - o(i,j,k,ixyz,ic) 3911 diff = abs(diff) 3912 if (diff.gt.thresh) then 3913 write(6,10000)count,i,j,k,ixyz,ic, 3914 & fdo(i,j,k,ixyz,ic), 3915 & o(i,j,k,ixyz,ic), diff 3916 endif 3917 enddo 3918 enddo 3919 enddo 3920 enddo 3921 enddo 392210000 format(1x,'<',i6,'>','(i=',i3,',j=', 3923 & i3,',k=',i3,',x=',i3,',at=',i3,') fd=', 3924 & 1pd12.5,' o=',1pd12.5,' diff=',1pd12.5) 3925 end 3926 subroutine raktest_ovd(rtdb) 3927 implicit none 3928#include "errquit.fh" 3929#include "mafdecls.fh" 3930#include "rtdb.fh" 3931#include "geom.fh" 3932#include "bas.fh" 3933 integer rtdb 3934c 3935 integer maxg1, maxs1 3936 integer h_buf, h_scr, h_cm, h_o, h_fdo, h_fdp, h_fdm 3937 integer k_buf, k_scr, k_cm, k_o, k_fdo, k_fdp, k_fdm 3938 integer geom, basis 3939 integer nat, nbf, bufsz, szscr, osz 3940c 3941 logical int_normalize 3942 external int_normalize 3943c 3944* if (.not.rtdb_print(rtdb,.true.)) stop ' rtdb_p 1?' 3945c 3946 if (.not.geom_create(geom,'geometry')) 3947 & call errquit('raktest_ovd: geom_create failed?',911, GEOM_ERR) 3948 if (.not. geom_rtdb_load(rtdb,geom,'geometry')) call errquit 3949 & ('raktest_ovd: geom_rtdb_load -ref failed',911, RTDB_ERR) 3950 if (.not.geom_print(geom)) stop ' print error' 3951c 3952 if (.not.bas_create(basis,'ao basis')) call errquit 3953 & ('raktest_ovd:bas_create failed',911, BASIS_ERR) 3954 if (.not.bas_rtdb_load(rtdb,geom,basis,'ao basis')) 3955 & call errquit 3956 & ('raktest_ovd:bas_rtdb_load failed',911, RTDB_ERR) 3957 if (.not.bas_print(basis)) stop ' print error' 3958 if (.not.gbs_map_print(basis)) stop ' print error' 3959 if (.not.rtdb_print(rtdb,.false.)) stop ' rtdb_p 2?' 3960 if (.not.int_normalize(rtdb,basis)) stop ' norm error 1' 3961 call intd_init(rtdb,1,basis) 3962 call int_mem_print() 3963c 3964 if (.not.geom_ncent(geom,nat)) stop 'geom_ncent fe' 3965 if (.not.bas_numbf(basis,nbf)) stop 'bas_numbf fe' 3966 write(6,*)' total number of basis functions:nbf:',nbf 3967 call int_mem_1e(maxg1,maxs1) 3968 bufsz = 2*maxg1 3969 szscr = 2*maxs1 3970 osz = nbf*nbf*3*nat 3971 if (.not.ma_push_get(mt_dbl,bufsz,'buf',h_buf,k_buf)) 3972 & stop ' ma buf failed' 3973 if (.not.ma_push_get(mt_dbl,szscr,'scr',h_scr,k_scr)) 3974 & stop ' ma scr failed' 3975 if (.not.ma_push_get(mt_dbl,(3*nat),'cm',h_cm,k_cm)) 3976 & stop ' ma cm failed' 3977 if (.not.ma_push_get(mt_dbl,osz,'o',h_o,k_o)) 3978 & stop ' ma o failed' 3979 if (.not.ma_push_get(mt_dbl,osz,'fdo',h_fdo,k_fdo)) 3980 & stop ' ma fdo failed' 3981 if (.not.ma_push_get(mt_dbl,osz,'fdp',h_fdp,k_fdp)) 3982 & stop ' ma fdp failed' 3983 if (.not.ma_push_get(mt_dbl,osz,'fdm',h_fdm,k_fdm)) 3984 & stop ' ma fdm failed' 3985 call raktest_ovd1(dbl_mb(k_buf),bufsz,dbl_mb(k_scr),szscr, 3986 & dbl_mb(k_fdo),dbl_mb(k_o),dbl_mb(k_fdp),dbl_mb(k_fdm), 3987 & nbf,nat,dbl_mb(k_cm),basis,geom) 3988 call intd_terminate() 3989 if (.not.ma_pop_stack(h_fdm)) stop 'h_fdm pop error' 3990 if (.not.ma_pop_stack(h_fdp)) stop 'h_fdp pop error' 3991 if (.not.ma_pop_stack(h_fdo)) stop 'h_fdo pop error' 3992 if (.not.ma_pop_stack(h_o)) stop 'h_o pop error' 3993 if (.not.ma_pop_stack(h_cm)) stop 'h_cm pop error' 3994 if (.not.ma_pop_stack(h_scr)) stop 'h_scr pop error' 3995 if (.not.ma_pop_stack(h_buf)) stop 'h_buf pop error' 3996 if (.not.bas_destroy(basis)) stop ' bas_destroy failed' 3997 if (.not.geom_destroy(geom)) stop ' geom_destroy failed' 3998 end 3999 subroutine raktest_ovd1(buf,lbuf,scr,lscr,fdo,o,fdp,fdm, 4000 & nbf,nat,cmaster,basis,geom) 4001 implicit none 4002#include "errquit.fh" 4003#include "mafdecls.fh" 4004#include "nwc_const.fh" 4005#include "geomP.fh" 4006#include "geom.fh" 4007#include "bas.fh" 4008#include "inp.fh" 4009#include "ecp_nwc.fh" 4010 double precision ddot 4011 external ddot 4012c 4013 integer lbuf, lscr 4014 double precision buf(lbuf),scr(lscr) 4015 integer nbf, nat 4016 double precision fdo(nbf,nbf,3,nat) 4017 double precision fdp(nbf,nbf,3,nat),fdm(nbf,nbf,3,nat) 4018 double precision o(nbf,nbf,3,nat) 4019 double precision cmaster(3,nat) 4020 integer basis, geom 4021* integer idatom(2) 4022c 4023 double precision delta, delta2i, norm 4024 integer ncont, count, zatom, ratom, ixyz 4025 integer ish,ilo,ihi,ibf 4026 integer jsh,jlo,jhi,jbf 4027 character*256 date_string 4028 integer lds 4029c 4030 if (.not.bas_numcont(basis,ncont)) stop ' bas_numcont' 4031 write(6,*)' buf size ', lbuf 4032 write(6,*)' lscr ',lscr 4033 call dfill(lbuf,0.0d00,buf,1) 4034 call dfill(lscr,0.0d00,scr,1) 4035 call dfill((nbf*nbf*3*nat),0.0d00,fdo,1) 4036 call dfill((nbf*nbf*3*nat),0.0d00,fdp,1) 4037 call dfill((nbf*nbf*3*nat),0.0d00,fdm,1) 4038 call dfill((nbf*nbf*3*nat),0.0d00,o,1) 4039 do ish = 1,ncont 4040 call util_date(date_string) 4041 lds = inp_strlen(date_string) 4042 write(6,*)' direct: ish: ',ish,' of ',ncont, 4043 & ' ',date_string(1:lds) 4044 call util_flush(6) 4045 if (.not.bas_cn2bfr(basis,ish,ilo,ihi)) stop 'cn2bfr i' 4046 do jsh = 1,ncont 4047 if (.not.bas_cn2bfr(basis,jsh,jlo,jhi)) stop 'cn2bfr i' 4048* call intd_1eov(basis,ish,basis,jsh,lscr,scr,lbuf,buf,idatom) 4049 call intd_1eh1(basis,ish,basis,jsh,lscr,scr,lbuf,buf) 4050 count = 0 4051* do zatom = 1,2 4052* if (idatom(zatom).gt.0) then 4053* ratom = idatom(zatom) 4054 do zatom = 1,nat 4055 ratom = zatom 4056 do ixyz = 1,3 4057 do ibf = ilo,ihi 4058 do jbf = jlo,jhi 4059 count = count + 1 4060 o(ibf,jbf,ixyz,ratom) = buf(count) 4061 enddo 4062 enddo 4063 enddo 4064* endif 4065 enddo 4066 enddo 4067 enddo 4068c 4069 call dcopy((3*nat),coords(1,1,geom),1,cmaster,1) 4070 delta = 0.00001d00 4071 delta2i = 1.0d00/(2.0d00*delta) 4072 do zatom = 1,nat 4073 do ixyz = 1,3 4074 call util_date(date_string) 4075 lds = inp_strlen(date_string) 4076 write(6,*)' finite diff: atom: ',zatom,' of ',nat, 4077 & ' coord(1:x,2:y,3:z) ',ixyz, 4078 & ' ',date_string(1:lds) 4079 call util_flush(6) 4080*-- 4081* + 4082 call dcopy((3*nat),cmaster,1,coords(1,1,geom),1) 4083 coords(ixyz,zatom,geom) = coords(ixyz,zatom,geom) + delta 4084c.. get coordinates for ecp centers. 4085 if (.not.geom_coords_ecp(geom,dbl_mb(k_xyzecp),n_ecp)) 4086 & call errquit('ecpdebug: geom_coords_ecp failed',911, 4087 & GEOM_ERR) 4088 do ish = 1,ncont 4089 call util_date(date_string) 4090 lds = inp_strlen(date_string) 4091 write(6,*)' finite diff: ish: ',ish,' of ',ncont, 4092 & ' ',date_string(1:lds) 4093 call util_flush(6) 4094 if (.not.bas_cn2bfr(basis,ish,ilo,ihi)) stop 'cn2bfr i' 4095 do jsh = 1,ncont 4096 if (.not.bas_cn2bfr(basis,jsh,jlo,jhi)) stop 'cn2bfr i' 4097* call int_1eov(basis,ish,basis,jsh,lscr,scr,lbuf,buf) 4098 call int_1eh1(basis,ish,basis,jsh,lscr,scr,lbuf,buf) 4099 count = 0 4100 do ibf = ilo,ihi 4101 do jbf = jlo,jhi 4102 count = count + 1 4103 fdp(ibf,jbf,ixyz,zatom) = buf(count) 4104 enddo 4105 enddo 4106 enddo 4107 enddo 4108* - 4109 call dcopy((3*nat),cmaster,1,coords(1,1,geom),1) 4110 coords(ixyz,zatom,geom) = coords(ixyz,zatom,geom) - delta 4111c.. get coordinates for ecp centers. 4112 if (.not.geom_coords_ecp(geom,dbl_mb(k_xyzecp),n_ecp)) 4113 & call errquit('ecpdebug: geom_coords_ecp failed',911, 4114 & GEOM_ERR) 4115 do ish = 1,ncont 4116 if (.not.bas_cn2bfr(basis,ish,ilo,ihi)) stop 'cn2bfr i' 4117 do jsh = 1,ncont 4118 if (.not.bas_cn2bfr(basis,jsh,jlo,jhi)) stop 'cn2bfr i' 4119* call int_1eov(basis,ish,basis,jsh,lscr,scr,lbuf,buf) 4120 call int_1eh1(basis,ish,basis,jsh,lscr,scr,lbuf,buf) 4121 count = 0 4122 do ibf = ilo,ihi 4123 do jbf = jlo,jhi 4124 count = count + 1 4125 fdm(ibf,jbf,ixyz,zatom) = buf(count) 4126 enddo 4127 enddo 4128 enddo 4129 enddo 4130*-- 4131 enddo 4132 enddo 4133 call dcopy((nbf*nbf*3*nat),fdp,1,fdo,1) 4134 call daxpy((nbf*nbf*3*nat),-1.0d00,fdm,1,fdo,1) 4135 call dscal((nbf*nbf*3*nat),delta2i,fdo,1) 4136 call printboth_ovd(nbf,nat,fdo,o) 4137 call dcopy((nbf*nbf*3*nat),fdo,1,fdm,1) 4138 call daxpy((nbf*nbf*3*nat),-1.0d00,o,1,fdm,1) 4139 norm = ddot((nbf*nbf*3*nat),fdm,1,fdm,1) 4140c 4141 write(6,*)' difference norm: ',norm 4142 end 4143 subroutine printboth_ovd(nbf,nat,fdo,o) 4144 implicit none 4145 integer nbf 4146 integer nat 4147 double precision fdo(nbf,nbf,3,nat) 4148 double precision o(nbf,nbf,3,nat) 4149c 4150 double precision diff, thresh 4151 integer i,j,ixyz,ic, count, nval_good, n_zero 4152c 4153 integer o_res, fd_res 4154 integer is_this_val_okay 4155 external is_this_val_okay 4156c 4157 thresh = 1.0d-06 4158c 4159 n_zero = 0 4160 nval_good = 0 4161 count = 0 4162 do ic = 1,nat 4163 do ixyz=1,3 4164 do j=1,nbf 4165 do i=1,nbf 4166 count = count + 1 4167 o_res = is_this_val_okay(o(i,j,ixyz,ic)) 4168* write(6,*)' o_res',o_res 4169 if (o_res .eq. 1) then 4170 write(6,10001)i,j,ixyz,ic 4171 elseif (o_res .eq. 2) then 4172 write(6,10003)i,j,ixyz,ic 4173 elseif (o_res.eq.0) then 4174 n_zero = n_zero + 1 4175 elseif (o_res.eq.3) then 4176 nval_good = nval_good + 1 4177 endif 4178 fd_res = is_this_val_okay(fdo(i,j,ixyz,ic)) 4179* write(6,*)' fd_res',fd_res 4180 if (fd_res .eq. 1) then 4181 write(6,10002)i,j,ixyz,ic 4182 elseif (fd_res .eq. 2) then 4183 write(6,10004)i,j,ixyz,ic 4184 endif 4185 diff = fdo(i,j,ixyz,ic) - o(i,j,ixyz,ic) 4186 diff = abs(diff) 4187 if (diff.gt.thresh) then 4188 write(6,10000)count,i,j,ixyz,ic, 4189 & fdo(i,j,ixyz,ic), 4190 & o(i,j,ixyz,ic), diff 4191* else 4192* if (abs(o(i,j,ixyz,ic)).gt.thresh) 4193* & nval_good = nval_good + 1 4194 endif 4195 enddo 4196 enddo 4197 enddo 4198 enddo 419910000 format(1x,'<',i6,'>','(i=',i3,',j=', 4200 & i3,',x=',i3,',at=',i3,') fd=', 4201 & 1pd12.5,' o=',1pd12.5,' diff=',1pd12.5) 4202 write(6,*)' zero values',n_zero 4203 write(6,*)' good values',nval_good 4204 write(6,*)' sum',(n_zero+nval_good) 4205 write(6,*)' count =',count 420610001 format(' calculated value (',i3,',',i3,',',i3,',', 4207 & i3,') is a NAN') 420810002 format(' finite diff value (',i3,',',i3,',',i3,',', 4209 & i3,') is a NAN') 421010003 format(' calculated value (',i3,',',i3,',',i3,',', 4211 & i3,') is an INFINITY') 421210004 format(' finite diff value (',i3,',',i3,',',i3,',', 4213 & i3,') is an INFINITY') 4214 end 4215 subroutine raktest_bd2e(rtdb) 4216 implicit none 4217#include "errquit.fh" 4218#include "geom.fh" 4219#include "bas.fh" 4220#include "mafdecls.fh" 4221#include "rtdb.fh" 4222#include "util.fh" 4223 integer rtdb 4224c 4225 integer basis 4226 integer geom 4227 integer nat 4228 integer nbf, nq_tot 4229 integer maxg, maxs, bufsz 4230 integer h_txs, h_nw, h_g, h_l, h_scr 4231 integer k_txs, k_nw, k_g, k_l, k_scr 4232 double precision norm 4233 logical int_normalize 4234 external int_normalize 4235c 4236 if (.not.geom_create(geom,'geometry')) 4237 & call errquit('raktest_geomwrt: geom_create failed?',911, 4238 & GEOM_ERR) 4239 if (.not. geom_rtdb_load(rtdb,geom,'geometry')) call errquit 4240 & ('raktest_geomwrt: geom_rtdb_load -ref failed',911, RTDB_ERR) 4241 if (.not.geom_print(geom)) stop ' print error' 4242c 4243 if (.not.bas_create(basis,'ao basis')) call errquit 4244 & ('bas_create failed',911, BASIS_ERR) 4245 if (.not.bas_rtdb_load(rtdb,geom,basis,'ao basis')) 4246 & call errquit 4247 & ('bas_rtdb_load failed',911, RTDB_ERR) 4248 if (.not.bas_print(basis)) stop ' print error' 4249 if (.not.rtdb_print(rtdb,.false.)) stop ' rtdb_p 2?' 4250 if (.not.int_normalize(rtdb,basis)) stop ' norm error 1' 4251c 4252 call intd_init(rtdb,1,basis) 4253 if (.not.geom_ncent(geom,nat)) stop 'geom_ncent fe' 4254 if (.not.bas_numbf(basis,nbf)) stop 'bas_numbf fe' 4255 if (.not.bas_numcont(basis,nq_tot)) stop 'bas_numcont fe' 4256 call int_mem_2e4c(maxg,maxs) 4257 maxg = maxg*2 4258 maxs = maxs*2 4259 bufsz = nbf*nbf*nbf*nbf*nat*3 4260c 4261 if (.not.ma_push_get(mt_dbl,bufsz,'texas',h_txs,k_txs)) 4262 & stop ' ma texas' 4263 call dfill(bufsz,0.0d00,dbl_mb(k_txs),1) 4264 if (.not.ma_push_get(mt_dbl,bufsz,'nwchem',h_nw,k_nw)) 4265 & stop ' ma nwchem' 4266 call dfill(bufsz,0.0d00,dbl_mb(k_nw),1) 4267 if (.not.ma_push_get(mt_dbl,12*maxg,'g',h_g,k_g)) 4268 & stop ' ma g' 4269 call dfill(12*maxg,0.0d00,dbl_mb(k_g),1) 4270 if (.not.ma_push_get(mt_int,4*maxg,'labels',h_l,k_l)) 4271 & stop ' ma labs' 4272 call ifill(4*maxg,0,int_mb(k_l),1) 4273 if (.not.ma_push_get(mt_dbl,maxs,'scratch',h_scr,k_scr)) 4274 & stop ' ma scratch' 4275 call dfill(maxs,0.0d00,dbl_mb(k_scr),1) 4276 write(6,*)' computing with texas ' 4277 call raktest_bd2e_getg(nq_tot, 4278 & nbf,nat,rtdb,basis,geom, 4279 & dbl_mb(k_txs), 4280 & maxg,dbl_mb(k_g),int_mb(k_l), 4281 & maxs,dbl_mb(k_scr)) 4282c 4283 call intd_terminate() 4284 call int_app_set_no_texas(rtdb) 4285 call intd_init(rtdb,1,basis) 4286 write(6,*)' computing with texas off' 4287 call raktest_bd2e_getg(nq_tot, 4288 & nbf,nat,rtdb,basis,geom, 4289 & dbl_mb(k_nw), 4290 & maxg,dbl_mb(k_g),int_mb(k_l), 4291 & maxs,dbl_mb(k_scr)) 4292 call intd_terminate() 4293 call int_app_unset_no_texas(rtdb) 4294 call raktest_bd2e_print(basis,nbf,nat, 4295 & dbl_mb(k_txs),dbl_mb(k_nw)) 4296 write(6,*)' computing norm' 4297 call daxpy(((nbf**4)*nat*3),-1.0d00, 4298 & dbl_mb(k_txs),1,dbl_mb(k_nw),1) 4299 norm = ddot(((nbf**4)*nat*3),dbl_mb(k_nw),1,dbl_mb(k_nw),1) 4300 write(6,*)' difference norm ',norm 4301c 4302 if (.not.ma_pop_stack(h_scr)) stop 'ma pop failed for h_scr' 4303 if (.not.ma_pop_stack(h_l)) stop 'ma pop failed for h_l' 4304 if (.not.ma_pop_stack(h_g)) stop 'ma pop failed for h_g' 4305 if (.not.ma_pop_stack(h_nw)) stop 'ma pop failed for h_nw' 4306 if (.not.ma_pop_stack(h_txs)) stop 'ma pop failed for h_txs' 4307 end 4308 subroutine raktest_bd2e_print(basis,nbf,nat,txs,nw) 4309 implicit none 4310#include "bas.fh" 4311 integer nbf,nat 4312 double precision txs(3,nat,nbf,nbf,nbf,nbf) 4313 double precision nw(3,nat,nbf,nbf,nbf,nbf) 4314 integer basis 4315c 4316 integer i,j,k,l,xyz,atom 4317 integer ia,ja,ka,la 4318 double precision tval 4319 double precision nval 4320 double precision dval 4321 double precision thresh 4322 logical printit 4323 thresh = 1.0d-08 4324 do i = 1,nbf 4325 do j = 1,nbf 4326 do k = 1,nbf 4327 do l = 1,nbf 4328 do atom = 1,nat 4329 do xyz = 1,3 4330 tval = txs(xyz,atom,i,j,k,l) 4331 nval = nw(xyz,atom,i,j,k,l) 4332 dval = abs((tval-nval)) 4333 printit = dval.gt.thresh 4334 if (printit) then 4335 write(6,10000)xyz,atom,i,j,k,l, 4336 & tval,nval,dval 4337 if (.not.bas_bf2ce(basis,i,ia)) stop 'bf2ce i' 4338 if (.not.bas_bf2ce(basis,j,ja)) stop 'bf2ce j' 4339 if (.not.bas_bf2ce(basis,k,ka)) stop 'bf2ce k' 4340 if (.not.bas_bf2ce(basis,l,la)) stop 'bf2ce l' 4341 write(6,10001)ia,ja,ka,la 4342 endif 4343 enddo 4344 enddo 4345 enddo 4346 enddo 4347 enddo 4348 enddo 434910000 format(1x,'deri [',6i5,'] ',1x,1pd20.10,1x,1pd20.10,1x,1pd20.10) 435010001 format(1x,6x,10x,4i5) 4351 end 4352 subroutine raktest_bd2e_getg(nq_total,nbf,nat,rtdb,basis,geom, 4353 & deri,lbuf,buf,labels,lscr,scr) 4354 implicit none 4355#include "errquit.fh" 4356#include "bas.fh" 4357 integer nq_total, nbf, nat, rtdb, basis,geom, lbuf, lscr 4358 double precision deri(3,nat,nbf,nbf,nbf,nbf) 4359 double precision buf(lbuf), scr(lscr) 4360 integer labels(lbuf,4) 4361c 4362 logical intbd_init4c, intbd_2e4c 4363 external intbd_init4c, intbd_2e4c 4364c 4365 integer nqs 4366 parameter (nqs = 20) 4367 double precision q4(nqs) 4368 integer iq(nqs),jq(nqs),kq(nqs),lq(nqs) 4369 integer qi,qj,qk,ql, q 4370 double precision dbl_dum 4371 integer ii,jj,kk,ll,iatom,jatom,katom,latom,xyz,int 4372 integer count, nintout 4373 logical more, lastq 4374c 4375 4376 q = 0 4377 do qi = 1,nq_total 4378 do qj = 1,qi 4379 do qk = 1,qj 4380 do ql = 1,qk 4381 lastq = qi.eq.nq_total 4382 lastq = lastq.and.qj.eq.qi 4383 lastq = lastq.and.qk.eq.qj 4384 lastq = lastq.and.ql.eq.qk 4385 if ((q+1).gt.nqs.or.lastq) then 4386 if (.not.intbd_init4c( 4387 & basis,iq,jq,basis,kq,lq,q,q4,.false., 4388 & lscr,scr,lbuf,dbl_dum)) call errquit 4389 & ('intbd_init4c failed',911, INT_ERR) 439000010 continue 4391 more = intbd_2e4c(basis,iq,jq,basis,kq,lq,q,q4,.false., 4392 & 1.0d-8,.false., 4393 & labels(1,1),labels(1,2),labels(1,3),labels(1,4), 4394 & buf,lbuf,nintout,lscr,scr) 4395 count = 0 4396 do int = 1,nintout 4397 ii = labels(int,1) 4398 jj = labels(int,2) 4399 kk = labels(int,3) 4400 ll = labels(int,4) 4401 if (.not.bas_bf2ce(basis,ii,iatom)) stop 'bf2ce i' 4402 if (.not.bas_bf2ce(basis,jj,jatom)) stop 'bf2ce j' 4403 if (.not.bas_bf2ce(basis,kk,katom)) stop 'bf2ce k' 4404 if (.not.bas_bf2ce(basis,ll,latom)) stop 'bf2ce l' 4405 do xyz = 1,3 4406 count = count + 1 4407 deri(xyz,iatom,ii,jj,kk,ll) = buf(count) 4408 enddo 4409 do xyz = 1,3 4410 count = count + 1 4411 deri(xyz,jatom,ii,jj,kk,ll) = buf(count) 4412 enddo 4413 do xyz = 1,3 4414 count = count + 1 4415 deri(xyz,katom,ii,jj,kk,ll) = buf(count) 4416 enddo 4417 do xyz = 1,3 4418 count = count + 1 4419 deri(xyz,latom,ii,jj,kk,ll) = buf(count) 4420 enddo 4421 if (iatom.eq.jatom) then 4422 do xyz = 1,3 4423 deri(xyz,iatom,ii,jj,kk,ll) = 4424 & deri(xyz,iatom,ii,jj,kk,ll) + 4425 & deri(xyz,jatom,ii,jj,kk,ll) 4426 deri(xyz,jatom,ii,jj,kk,ll) = 0.0d00 4427 enddo 4428 endif 4429 if (iatom.eq.katom) then 4430 do xyz = 1,3 4431 deri(xyz,iatom,ii,jj,kk,ll) = 4432 & deri(xyz,iatom,ii,jj,kk,ll) + 4433 & deri(xyz,katom,ii,jj,kk,ll) 4434 deri(xyz,katom,ii,jj,kk,ll) = 0.0d00 4435 enddo 4436 endif 4437 if (iatom.eq.latom) then 4438 do xyz = 1,3 4439 deri(xyz,iatom,ii,jj,kk,ll) = 4440 & deri(xyz,iatom,ii,jj,kk,ll) + 4441 & deri(xyz,latom,ii,jj,kk,ll) 4442 deri(xyz,latom,ii,jj,kk,ll) = 0.0d00 4443 enddo 4444 endif 4445 if (jatom.eq.katom) then 4446 do xyz = 1,3 4447 deri(xyz,jatom,ii,jj,kk,ll) = 4448 & deri(xyz,jatom,ii,jj,kk,ll) + 4449 & deri(xyz,katom,ii,jj,kk,ll) 4450 deri(xyz,katom,ii,jj,kk,ll) = 0.0d00 4451 enddo 4452 endif 4453 if (jatom.eq.latom) then 4454 do xyz = 1,3 4455 deri(xyz,jatom,ii,jj,kk,ll) = 4456 & deri(xyz,jatom,ii,jj,kk,ll) + 4457 & deri(xyz,latom,ii,jj,kk,ll) 4458 deri(xyz,latom,ii,jj,kk,ll) = 0.0d00 4459 enddo 4460 endif 4461 if (katom.eq.latom) then 4462 do xyz = 1,3 4463 deri(xyz,katom,ii,jj,kk,ll) = 4464 & deri(xyz,katom,ii,jj,kk,ll) + 4465 & deri(xyz,latom,ii,jj,kk,ll) 4466 deri(xyz,latom,ii,jj,kk,ll) = 0.0d00 4467 enddo 4468 endif 4469 4470 enddo 4471 if (more) goto 00010 4472 q = 1 4473 else 4474 q = q + 1 4475 endif 4476 iq(q) = qi 4477 jq(q) = qj 4478 kq(q) = qk 4479 lq(q) = ql 4480 enddo 4481 enddo 4482 enddo 4483 enddo 4484c 4485 end 4486 subroutine raktest_diskspeed(rtdb) 4487 implicit none 4488#include "stdio.fh" 4489#include "rtdb.fh" 4490#include "mafdecls.fh" 4491#include "util.fh" 4492#include "global.fh" 4493 integer rtdb 4494c 4495 integer default_size 4496 parameter (default_size = 512) 4497 integer default_count 4498 parameter (default_count = 10) 4499 integer default_iterations 4500 parameter (default_iterations = 10) 4501 integer size, count, iters 4502 integer h_io, k_io 4503 integer h_data, k_data 4504 integer data_size 4505c 4506 if (.not.rtdb_get(rtdb,'rak:disksize',mt_int,1,size)) 4507 & size = default_size 4508 if (.not.rtdb_get(rtdb,'rak:diskcount',mt_int,1,count)) 4509 & count = default_count 4510 if (.not.rtdb_get(rtdb,'rak:diskiters',mt_int,1,iters)) 4511 & iters = default_iterations 4512c 4513 if (ga_nodeid().eq.0) then 4514 write(luout,*) 4515 & ' size of buffer is : ',size, ' doubles' 4516 write(luout,*) 4517 & ' : ',(size*8), ' bytes' 4518 write(luout,*) 4519 & ' number of buffers put out is :',count 4520 write(luout,*) 4521 & ' number of iterations is :',iters 4522 write(luout,*) 4523 & ' output and reading of :', 4524 & 8*size*count*iters,' bytes' 4525 endif 4526c 4527 if (.not.ma_push_get(mt_dbl,size,'io buffer',h_io,k_io)) 4528 & stop ' ma_get of io buffer failed' 4529 data_size = 5*ga_nnodes() 4530 if (.not.ma_push_get(mt_dbl,data_size,'data',h_data,k_data)) 4531 & stop ' ma get of data buffer failed' 4532c 4533 call raktest_diskspeed_fill(size,dbl_mb(k_io)) 4534 call raktest_diskspeed_write(dbl_mb(k_io),size, 4535 & count,iters,dbl_mb(k_data),'rakdisk',.true.) 4536 call raktest_diskspeed_read(dbl_mb(k_io),size, 4537 & count,iters,dbl_mb(k_data),'rakdisk',.false.) 4538c 4539 if (.not.ma_pop_stack(h_data)) stop ' ma pop error ' 4540 if (.not.ma_pop_stack(h_io)) stop ' ma pop error ' 4541c 4542 end 4543 subroutine raktest_diskspeed_read(buf,size, 4544 & count,iters,adata,filename_stub,save_file) 4545 implicit none 4546#include "stdio.fh" 4547#include "eaf.fh" 4548#include "util.fh" 4549#include "global.fh" 4550#include "tcgmsg.fh" 4551#include "msgtypesf.h" 4552 integer size 4553 integer count 4554 integer iters 4555 double precision buf(size) 4556 double precision adata(5,*) 4557 character*(*)filename_stub 4558 logical save_file 4559c 4560 character*256 filename 4561 integer i, it 4562 integer fd 4563 integer eaf_rv 4564 double precision offset, bytes_tot 4565 double precision timec, timew, ratec, ratew 4566 integer inode 4567 integer adata_size 4568 integer mtype 4569 double precision timew_min, timew_max, timew_ave 4570 double precision ratew_min, ratew_max, ratew_ave 4571c 4572 call ga_sync() 4573 call util_file_name(filename_stub,.false.,.true.,filename) 4574c 4575 eaf_rv = eaf_open(filename,EAF_R,fd) 4576 if (eaf_rv.ne.0) then 4577 call eaf_errmsg(eaf_rv) 4578 stop 'eaf open error' 4579 endif 4580c 4581 timec = util_cpusec() 4582 timew = util_wallsec() 4583 do it = 1,iters 4584 offset = 0.0d00 4585 do i = 1,count 4586 eaf_rv = eaf_read(fd,offset,buf,8*size) 4587 if (eaf_rv.ne.0) then 4588 call eaf_errmsg(eaf_rv) 4589 stop ' eaf read error' 4590 endif 4591 offset = offset + dble(8*size) 4592 enddo 4593 enddo 4594 timec = util_cpusec() - timec 4595 timew = util_wallsec() - timew 4596 if (ga_nodeid().eq.0) then 4597 write(luout,*)' EAF info for node 0 only ' 4598 call eaf_print_stats(fd) 4599 endif 4600 call ga_sync() 4601c 4602 eaf_rv = eaf_close(fd) 4603 if (eaf_rv.ne.0) then 4604 call eaf_errmsg(eaf_rv) 4605 stop 'eaf close error' 4606 endif 4607c 4608 if (.not.save_file) call util_file_unlink(filename) 4609c 4610 bytes_tot = dble(count)*dble(8)*dble(size) 4611 bytes_tot = bytes_tot*dble(iters) 4612 bytes_tot = bytes_tot*1.0d-6 4613 ratec = bytes_tot/timec 4614 ratew = bytes_tot/timew 4615 adata_size = 5*ga_nnodes() 4616 call dfill(adata_size,0.0d00,adata,1) 4617 inode = ga_nodeid() + 1 4618 adata(1,inode) = bytes_tot 4619 adata(2,inode) = timec 4620 adata(3,inode) = ratec 4621 adata(4,inode) = timew 4622 adata(5,inode) = ratew 4623 mtype = 1234 + MSGDBL 4624 call ga_dgop(mtype,adata,adata_size,'+') 4625 if (ga_nodeid().eq.0) then 4626 do inode = 1,ga_nnodes() 4627 write(luout,*) 4628 & ' read statistics for node :',(inode-1) 4629 write(luout,10000) 4630 & ' read total Mbytes :',adata(1,inode) 4631 write(luout,10000) 4632 & ' read cpu time (s) :',adata(2,inode) 4633 write(luout,10000) 4634 & ' read cpu rate (mb/s) :',adata(3,inode) 4635 write(luout,10000) 4636 & ' read wall time (s) :',adata(4,inode) 4637 write(luout,10000) 4638 & ' read wall rate (mb/s) :',adata(5,inode) 4639 write(luout,*)' ' 4640 call util_flush(luout) 4641 enddo 4642 endif 4643 call ga_sync() 4644 call ga_sync() 4645 timew_min = timew 4646 call ga_dgop(5651,timew_min,1,'min') 4647 timew_max = timew 4648 call ga_dgop(5652,timew_max,1,'max') 4649 timew_ave = timew 4650 call ga_dgop(5653,timew_ave,1,'+') 4651 timew_ave = timew_ave/dble(ga_nnodes()) 4652 ratew_min = ratew 4653 call ga_dgop(5654,ratew_min,1,'min') 4654 ratew_max = ratew 4655 call ga_dgop(5655,ratew_max,1,'max') 4656 ratew_ave = ratew 4657 call ga_dgop(5656,ratew_ave,1,'+') 4658 ratew_ave = ratew_ave/dble(ga_nnodes()) 4659 call ga_sync() 4660 call ga_sync() 4661 call ga_sync() 4662 if (ga_nodeid().eq.0) then 4663 write(luout,10000)' read minimum wall time :',timew_min 4664 write(luout,10000)' read maximum wall time :',timew_max 4665 write(luout,10000)' read average wall time :',timew_ave 4666 write(luout,10000)' read minimum wall rate :',ratew_min 4667 write(luout,10000)' read maximum wall rate :',ratew_max 4668 write(luout,10000)' read average wall rate :',ratew_ave 4669 endif 4670c 467110000 format(1x,a,f10.5) 4672 end 4673 subroutine raktest_diskspeed_write(buf,size, 4674 & count,iters,adata,filename_stub,save_file) 4675 implicit none 4676#include "stdio.fh" 4677#include "eaf.fh" 4678#include "util.fh" 4679#include "global.fh" 4680#include "tcgmsg.fh" 4681#include "msgtypesf.h" 4682 integer size 4683 integer count 4684 integer iters 4685 double precision buf(size) 4686 double precision adata(5,*) 4687 character*(*)filename_stub 4688 logical save_file 4689c 4690 character*256 filename 4691 integer i, it 4692 integer fd 4693 integer eaf_rv 4694 double precision offset, bytes_tot 4695 double precision timec, timew, ratec, ratew 4696 integer inode 4697 integer adata_size 4698 integer mtype 4699 double precision timew_min, timew_max, timew_ave 4700 double precision ratew_min, ratew_max, ratew_ave 4701c 4702 call ga_sync() 4703 call util_file_name(filename_stub,.false.,.true.,filename) 4704c 4705 eaf_rv = eaf_open(filename,EAF_RW,fd) 4706 if (eaf_rv.ne.0) then 4707 call eaf_errmsg(eaf_rv) 4708 stop 'eaf open error' 4709 endif 4710c 4711 timec = util_cpusec() 4712 timew = util_wallsec() 4713 do it = 1,iters 4714 offset = 0.0d00 4715 do i = 1,count 4716 eaf_rv = eaf_write(fd,offset,buf,8*size) 4717 if (eaf_rv.ne.0) then 4718 call eaf_errmsg(eaf_rv) 4719 stop ' eaf write error' 4720 endif 4721 offset = offset + dble(8*size) 4722 enddo 4723 enddo 4724 timec = util_cpusec() - timec 4725 timew = util_wallsec() - timew 4726 if (ga_nodeid().eq.0) then 4727 write(luout,*)' EAF info for node 0 only ' 4728 call eaf_print_stats(fd) 4729 endif 4730 call ga_sync() 4731c 4732 eaf_rv = eaf_close(fd) 4733 if (eaf_rv.ne.0) then 4734 call eaf_errmsg(eaf_rv) 4735 stop 'eaf close error' 4736 endif 4737c 4738 if (.not.save_file) call util_file_unlink(filename) 4739c 4740 bytes_tot = dble(count)*dble(8)*dble(size) 4741 bytes_tot = bytes_tot*dble(iters) 4742 bytes_tot = bytes_tot*1.0d-6 4743 ratec = bytes_tot/timec 4744 ratew = bytes_tot/timew 4745 adata_size = 5*ga_nnodes() 4746 call dfill(adata_size,0.0d00,adata,1) 4747 inode = ga_nodeid() + 1 4748 adata(1,inode) = bytes_tot 4749 adata(2,inode) = timec 4750 adata(3,inode) = ratec 4751 adata(4,inode) = timew 4752 adata(5,inode) = ratew 4753 mtype = 2134 + MSGDBL 4754 call ga_dgop(mtype,adata,adata_size,'+') 4755 if (ga_nodeid().eq.0) then 4756 do inode = 1,ga_nnodes() 4757 write(luout,*) 4758 & ' write statistics for node :',(inode-1) 4759 write(luout,10000) 4760 & ' write total Mbytes :',adata(1,inode) 4761 write(luout,10000) 4762 & ' write cpu time (s) :',adata(2,inode) 4763 write(luout,10000) 4764 & ' write cpu rate (mb/s) :',adata(3,inode) 4765 write(luout,10000) 4766 & ' write wall time (s) :',adata(4,inode) 4767 write(luout,10000) 4768 & ' write wall rate (mb/s) :',adata(5,inode) 4769 write(luout,*)' ' 4770 call util_flush(luout) 4771 enddo 4772 endif 4773 call ga_sync() 4774 call ga_sync() 4775 timew_min = timew 4776 call ga_dgop(5657,timew_min,1,'min') 4777 timew_max = timew 4778 call ga_dgop(5658,timew_max,1,'max') 4779 timew_ave = timew 4780 call ga_dgop(5659,timew_ave,1,'+') 4781 timew_ave = timew_ave/dble(ga_nnodes()) 4782 ratew_min = ratew 4783 call ga_dgop(5660,ratew_min,1,'min') 4784 ratew_max = ratew 4785 call ga_dgop(5661,ratew_max,1,'max') 4786 ratew_ave = ratew 4787 call ga_dgop(5662,ratew_ave,1,'+') 4788 ratew_ave = ratew_ave/dble(ga_nnodes()) 4789 call ga_sync() 4790 call ga_sync() 4791 call ga_sync() 4792 if (ga_nodeid().eq.0) then 4793 write(luout,10000)' write minimum wall time :',timew_min 4794 write(luout,10000)' write maximum wall time :',timew_max 4795 write(luout,10000)' write average wall time :',timew_ave 4796 write(luout,10000)' write minimum wall rate :',ratew_min 4797 write(luout,10000)' write maximum wall rate :',ratew_max 4798 write(luout,10000)' write average wall rate :',ratew_ave 4799 call util_flush(luout) 4800 endif 4801c 480210000 format(1x,a,f10.5) 4803 end 4804 subroutine raktest_diskspeed_fill(size,buf) 4805 implicit none 4806 integer size 4807 double precision buf(size) 4808c 4809 integer i 4810c 4811 do i = 1,size 4812 buf(size) = sqrt(dble(i)) 4813 enddo 4814 end 4815 subroutine raktest_2ecompare(rtdb) 4816 implicit none 4817#include "errquit.fh" 4818#include "stdio.fh" 4819#include "mafdecls.fh" 4820#include "geom.fh" 4821#include "bas.fh" 4822#include "rtdb.fh" 4823#include "global.fh" 4824c:functions 4825 logical int_normalize 4826 external int_normalize 4827c:passed 4828 integer rtdb 4829c:local 4830 integer geom, basis 4831 integer szbuf, Lmax, nshell 4832 integer maxg, maxscr 4833 integer size_dbls, h_dbls, k_dbls 4834 integer size_ints, h_ints, k_ints 4835 integer k_bufnw, k_buftxs, k_scr, k_lab, k_eri 4836c 4837 if (.not.geom_create(geom,'geometry')) 4838 & call errquit('raktest_geomwrt: geom_create failed?',911, 4839 & GEOM_ERR) 4840 if (.not. geom_rtdb_load(rtdb,geom,'geometry')) call errquit 4841 & ('raktest_geomwrt: geom_rtdb_load -ref failed',911, RTDB_ERR) 4842c 4843 if (.not.bas_create(basis,'ao basis')) call errquit 4844 & ('bas_create failed',911, BASIS_ERR) 4845 if (.not.bas_rtdb_load(rtdb,geom,basis,'ao basis')) 4846 & call errquit 4847 & ('bas_rtdb_load failed',911, RTDB_ERR) 4848 if (ga_nodeid().eq.0) then 4849 if (.not.geom_print(geom)) stop ' print error' 4850 if (.not.bas_print(basis)) stop ' print error' 4851 if (.not.gbs_map_print(basis)) stop ' gbs_map_print 2?' 4852 endif 4853 if (.not.int_normalize(rtdb,basis)) stop ' norm error 1' 4854 if (.not.bas_numcont(basis,nshell)) stop 'bas_numcont fe' 4855 if (.not.bas_high_angular(basis,Lmax)) stop 'bas_ha fe' 4856c 4857 szbuf = (Lmax+1)*(Lmax+2)/2 4858 szbuf = szbuf**4 4859c 4860 call rak_init(rtdb,1,basis) 4861 call intb_mem_2e4c(maxg,maxscr) 4862 maxscr = maxscr + szbuf 4863 maxscr = maxscr + maxscr/5 + 1 4864 maxscr = max(51000,maxscr) 4865 call int_terminate() 4866c 4867 if (szbuf.lt.maxg) then 4868 write(luout,*)' szbuf =',szbuf 4869 write(luout,*)' maxg =',maxg 4870 call errquit('raktest_2ecompare:fatal error',911, UNKNOWN_ERR) 4871 endif 4872c 4873 size_dbls = 3*szbuf+maxscr 4874 size_ints = 4*szbuf 4875 if (ga_nodeid().eq.0) then 4876 write(luout,*)' raktest: maxscr :',maxscr 4877 write(luout,*)' raktest: szbuf :',szbuf 4878 write(luout,*)' raktest: size_dbls:',size_dbls 4879 write(luout,*)' raktest: size_ints:',size_ints 4880 call util_flush(luout) 4881 endif 4882 if (.not.ma_push_get(mt_dbl,size_dbls,'dbls scr',h_dbls,k_dbls)) 4883 & call errquit('raktest_2ecompare:fatal error:ma dbls',911, 4884 & MA_ERR) 4885 if (.not.ma_push_get(mt_int,size_ints,'ints scr',h_ints,k_ints)) 4886 & call errquit('raktest_2ecompare:fatal error:ma ints',911, 4887 & MA_ERR) 4888 k_bufnw = k_dbls 4889 k_buftxs = k_bufnw + szbuf 4890 k_eri = k_buftxs + szbuf 4891 k_scr = k_eri + szbuf 4892 k_lab = k_ints 4893c 4894 call ga_sync() 4895 4896 call raktest_2ecompare_a( 4897 & dbl_mb(k_bufnw), 4898 & dbl_mb(k_buftxs), 4899 & dbl_mb(k_scr), 4900 & dbl_mb(k_eri), 4901 & int_mb(k_lab), 4902 & szbuf,maxscr,nshell,basis,geom,rtdb) 4903c 4904 if (.not.ma_pop_stack(h_ints)) stop ' ma pop error' 4905 if (.not.ma_pop_stack(h_dbls)) stop ' ma pop error' 4906 if (.not.bas_destroy(basis)) stop 'bas_destroy fe' 4907 if (.not.geom_destroy(geom)) stop 'geom destroy fe' 4908c 4909 end 4910 subroutine raktest_2ecompare_a(bufnw, buftxs, scr, eri, 4911 & labs, szbuf, maxscr, nsh, basis, geom, rtdb) 4912 implicit none 4913#include "bas.fh" 4914#include "global.fh" 4915#include "util.fh" 4916#include "stdio.fh" 4917c:functions 4918 logical intb_init4c, intb_2e4c 4919 external intb_init4c, intb_2e4c 4920c:passed 4921 integer rtdb 4922 integer szbuf 4923 integer maxscr 4924 integer nsh 4925 integer basis 4926 integer geom 4927 double precision eri(szbuf) 4928 double precision bufnw(szbuf) 4929 double precision buftxs(szbuf) 4930 integer labs(szbuf,4) 4931 double precision scr(maxscr) 4932c:local 4933 integer me, nproc 4934 integer ish, jsh, ksh, lsh 4935 integer tish, tjsh, tksh, tlsh 4936 integer ilo, ihi, jlo, jhi, klo, khi, llo, lhi 4937 double precision q4, dummy, norm 4938 integer nint_t, nint_n 4939 logical more 4940 double precision zero,thresh 4941 double precision nwmax, txsmax 4942 double precision prev_time, now_time, delta_time 4943 integer count, wrong 4944 integer task_count 4945c 4946 me = ga_nodeid() 4947 nproc = ga_nnodes() 4948 zero = 1.0d-09 4949 thresh = zero*0.01d00 4950c 4951* call setdbg(1) 4952 call int_app_set_no_spint(rtdb) 4953 call rak_init(rtdb,1,basis) 4954c 4955 wrong = 0 4956 count = 0 4957 task_count = (me-1) 4958 do ish = nsh,1,-1 4959 if (ish.eq.nsh) then 4960 delta_time = 0.0d0 4961 prev_time = util_wallsec() 4962 else 4963 now_time = util_wallsec() 4964 delta_time = now_time - prev_time 4965 prev_time = now_time 4966 endif 4967 call ga_sync() 4968 if (me.eq.0) then 4969 write(luout,10001)ish,delta_time 4970 call util_flush(luout) 4971 endif 4972 if (.not.bas_cn2bfr(basis,ish,ilo,ihi)) stop 'cn2bfr i' 4973 do jsh = ish,1,-1 4974 if (.not.bas_cn2bfr(basis,jsh,jlo,jhi)) stop 'cn2bfr j' 4975 do ksh = jsh,1,-1 4976 if (.not.bas_cn2bfr(basis,ksh,klo,khi)) stop 'cn2bfr k' 4977 do lsh = ksh,1,-1 4978 if (.not.bas_cn2bfr(basis,lsh,llo,lhi)) stop 'cn2bfr l' 4979 count = count + 1 4980 task_count = task_count + 1 4981 if (mod(task_count,nproc).eq.0) then 4982* write(luout,*)ish,jsh,ksh,lsh,':', 4983* & task_count,count,me 4984* call util_flush(luout) 4985c... do default 4986 tish = ish 4987 tjsh = jsh 4988 tksh = ksh 4989 tlsh = lsh 4990 call ifill(4*szbuf,0,labs,1) 4991 call dfill(szbuf,0.0d00,buftxs,1) 4992* call rak_init(rtdb,1,basis) 4993 if (.not.intb_init4c( 4994 & basis,tish,tjsh,basis,tksh,tlsh,1, 4995 & q4,.false.,maxscr,scr,szbuf,dummy)) 4996 & stop 'intb_init failed 1' 499700100 continue 4998 call dfill (szbuf,0.0d00,eri,1) 4999* write(6,*)ish,jsh,ksh,lsh,':',' intb_2e4c 1',me 5000* call util_flush(luout) 5001 more = intb_2e4c( 5002 & basis,tish,tjsh,basis,tksh,tlsh,1, 5003 & q4,.false.,zero,.false., 5004 & labs(1,1),labs(1,2),labs(1,3),labs(1,4), 5005 & eri,szbuf,nint_t,maxscr,scr) 5006 call raktest_2ecompare_cp(buftxs,eri,nint_t, 5007 & labs(1,1),labs(1,2),labs(1,3),labs(1,4), 5008 & ilo,ihi,jlo,jhi,klo,khi,llo,lhi,nwmax) 5009 if (more) then 5010* write(luout,*)ish,jsh,ksh,lsh,':',' more 2 ?',me 5011* call util_flush(luout) 5012 goto 00100 5013 endif 5014* call int_terminate() 5015c... do no texas .eg., force nwchem integrals 5016 tish = ish 5017 tjsh = jsh 5018 tksh = ksh 5019 tlsh = lsh 5020 call ifill(4*szbuf,0,labs,1) 5021 call dfill(szbuf,0.0d00,bufnw,1) 5022* call int_app_set_no_texas(rtdb) 5023* call rak_init(rtdb,1,basis) 5024 call rak_notxs() 5025 if (.not.intb_init4c( 5026 & basis,tish,tjsh,basis,tksh,tlsh,1, 5027 & q4,.false.,maxscr,scr,szbuf,dummy)) 5028 & stop 'intb_init failed 2' 502900200 continue 5030 call dfill (szbuf,0.0d00,eri,1) 5031* write(6,*)ish,jsh,ksh,lsh,':',' intb_2e4c 2',me 5032* call util_flush(luout) 5033 more = intb_2e4c( 5034 & basis,tish,tjsh,basis,tksh,tlsh,1, 5035 & q4,.false.,zero,.false., 5036 & labs(1,1),labs(1,2),labs(1,3),labs(1,4), 5037 & eri,szbuf,nint_n,maxscr,scr) 5038 call raktest_2ecompare_cp(bufnw,eri,nint_n, 5039 & labs(1,1),labs(1,2),labs(1,3),labs(1,4), 5040 & ilo,ihi,jlo,jhi,klo,khi,llo,lhi,txsmax) 5041 if (more) then 5042* write(luout,*)ish,jsh,ksh,lsh,':',' more 2 ?',me 5043* call util_flush(luout) 5044 goto 00200 5045 endif 5046 call rak_nonotxs() 5047* call int_terminate() 5048* call int_app_unset_no_texas(rtdb) 5049c 5050* write(6,*)ish,jsh,ksh,lsh,':',' dcopy',me 5051* call util_flush(luout) 5052 call dcopy(szbuf,bufnw,1,eri,1) 5053 call daxpy(szbuf,-1.0d00,buftxs,1,eri,1) 5054 norm = ddot(szbuf,eri,1,eri,1) 5055 if (norm.gt.thresh)then 5056 wrong = wrong + 1 5057 call raktest_2ecompare_mdiff(szbuf,bufnw,buftxs,nwmax) 5058 write(luout,10000) 5059 & me,wrong,count,ish,jsh,ksh,lsh, 5060 & nint_t,nint_n,norm,nwmax 5061 call util_flush(luout) 5062 endif 5063 if (nint_t.ne.nint_n.and.norm.gt.thresh)then 5064 wrong = wrong + 1 5065 call raktest_2ecompare_mdiff(szbuf,bufnw,buftxs,nwmax) 5066 if (nwmax.gt.zero) then 5067 write(luout,10000) 5068 & me,wrong,count,ish,jsh,ksh,lsh, 5069 & nint_t,nint_n,norm, 5070 & nwmax 5071 call util_flush(luout) 5072 endif 5073 endif 5074c 5075 endif 5076 enddo 5077 enddo 5078 enddo 5079 enddo 5080c 5081* call ga_sync() 5082c 5083 delta_time = util_wallsec() - prev_time 5084 if (me.eq.0) then 5085 write(luout,10001) 0, delta_time 5086 call util_flush(luout) 5087 endif 5088 call ga_dgop(5662,wrong,1,'+') 5089 if (me.eq.0) write(luout,10002) wrong,count 5090c 5091 call int_terminate() 5092 call int_app_unset_no_spint(rtdb) 5093c 509410000 format(1x,i3,':',i5,' of',i8,'(',i4,i4,'|',i4,i4,') [t:',i5, 5095 & '|n:',i5,'] norm=',2(1pd15.6)) 509610001 format(1x,60('-'),'>',i5,1x,f8.2 ) 509710002 format(1x,i5,' of',i8,' quartets are bad') 5098c 5099 end 5100 subroutine rak_notxs() 5101 implicit none 5102* 5103* user_* variables determine .true. a user set some value for the 5104* specific integral code. 5105* def_* variables is the value that the user set. 5106* 5107* this means that if the user does not want to run the sp integral code 5108* he/she would set "int:cando_sp" false and the values of would be 5109* user_cando_sp = .true. and def_cando_sp = .false. 5110* 5111* to test then use: 5112* 5113* if(user_cando_sp.and.(.not.def_cando_sp) then 5114* do not do anything with sp code 5115* endif 5116* 5117* or 5118* 5119* if (.not.((user_cando_sp.and.(.not.def_cando_sp)))) call sp_code 5120* 5121* 5122* Ricky A. Kendall, HPCCG, EMSL, PNNL 5123* 5124 logical user_cando_sp ! did user set a value for sp 5125 logical user_cando_nw ! did user set a value for nw 5126 logical user_cando_txs ! did user set a value for txs 5127 logical def_cando_sp ! default user setable value for cando_sp 5128 logical def_cando_nw ! default user setable value for cando_nw 5129 logical def_cando_txs ! default user setable value for cando_txs 5130c 5131 logical app_stored_txs ! value stored in int_app_set_no_texas 5132 logical app_stored_spint ! value stored in int_app_set_no_spint 5133 integer rtdbIused 5134c 5135 common /clcando/ user_cando_sp, user_cando_nw, user_cando_txs, 5136 & def_cando_sp, def_cando_nw, def_cando_txs, 5137 & app_stored_txs, app_stored_spint, 5138 & rtdbIused 5139c 5140 user_cando_txs = .true. 5141 def_cando_txs = .false. 5142c 5143 end 5144 subroutine rak_nonotxs() 5145 implicit none 5146* 5147* user_* variables determine .true. a user set some value for the 5148* specific integral code. 5149* def_* variables is the value that the user set. 5150* 5151* this means that if the user does not want to run the sp integral code 5152* he/she would set "int:cando_sp" false and the values of would be 5153* user_cando_sp = .true. and def_cando_sp = .false. 5154* 5155* to test then use: 5156* 5157* if(user_cando_sp.and.(.not.def_cando_sp) then 5158* do not do anything with sp code 5159* endif 5160* 5161* or 5162* 5163* if (.not.((user_cando_sp.and.(.not.def_cando_sp)))) call sp_code 5164* 5165* 5166* Ricky A. Kendall, HPCCG, EMSL, PNNL 5167* 5168 logical user_cando_sp ! did user set a value for sp 5169 logical user_cando_nw ! did user set a value for nw 5170 logical user_cando_txs ! did user set a value for txs 5171 logical def_cando_sp ! default user setable value for cando_sp 5172 logical def_cando_nw ! default user setable value for cando_nw 5173 logical def_cando_txs ! default user setable value for cando_txs 5174c 5175 logical app_stored_txs ! value stored in int_app_set_no_texas 5176 logical app_stored_spint ! value stored in int_app_set_no_spint 5177 integer rtdbIused 5178c 5179 common /clcando/ user_cando_sp, user_cando_nw, user_cando_txs, 5180 & def_cando_sp, def_cando_nw, def_cando_txs, 5181 & app_stored_txs, app_stored_spint, 5182 & rtdbIused 5183c 5184 user_cando_txs = .false. 5185 def_cando_txs = .false. 5186c 5187 end 5188 subroutine raktest_2ecompare_mdiff(szbuf,bufnw,buftxs,valmax) 5189 implicit none 5190 integer szbuf 5191 double precision bufnw(szbuf),buftxs(szbuf),valmax 5192c 5193 double precision diff 5194 integer i 5195 valmax = 0.0d00 5196 do i = 1,szbuf 5197 diff = bufnw(i)-buftxs(i) 5198 valmax = max(abs(diff),valmax) 5199 enddo 5200 end 5201 subroutine raktest_2ecompare_cp(buf,eri,nint, 5202 & ilb,jlb,klb,llb, 5203 & ilo,ihi,jlo,jhi,klo,khi,llo,lhi,valmax) 5204 implicit none 5205 integer nint 5206 integer ilb(nint),jlb(nint),klb(nint),llb(nint) 5207 integer ilo,ihi,jlo,jhi,klo,khi,llo,lhi 5208 double precision eri(nint) 5209 double precision buf(llo:lhi,klo:khi,jlo:jhi,ilo:ihi) 5210 double precision valmax 5211c 5212 integer n,i,j,k,l 5213 valmax = 0.0d00 5214 do n = 1,nint 5215 i = ilb(n) 5216 j = jlb(n) 5217 k = klb(n) 5218 l = llb(n) 5219 buf(l,k,j,i) = eri(n) 5220 valmax = max(valmax,abs(eri(n))) 5221 enddo 5222 end 5223 subroutine rak_init(rtdb, nbas, bases) 5224c 5225c initializes integral code to data structers for a integral computation 5226c: no print 5227 implicit none 5228#include "errquit.fh" 5229#include "bas.fh" 5230#include "apiP.fh" 5231* 5232* user_* variables determine .true. a user set some value for the 5233* specific integral code. 5234* def_* variables is the value that the user set. 5235* 5236* this means that if the user does not want to run the sp integral code 5237* he/she would set "int:cando_sp" false and the values of would be 5238* user_cando_sp = .true. and def_cando_sp = .false. 5239* 5240* to test then use: 5241* 5242* if(user_cando_sp.and.(.not.def_cando_sp) then 5243* do not do anything with sp code 5244* endif 5245* 5246* or 5247* 5248* if (.not.((user_cando_sp.and.(.not.def_cando_sp)))) call sp_code 5249* 5250* 5251* Ricky A. Kendall, HPCCG, EMSL, PNNL 5252* 5253 logical user_cando_sp ! did user set a value for sp 5254 logical user_cando_nw ! did user set a value for nw 5255 logical user_cando_txs ! did user set a value for txs 5256 logical def_cando_sp ! default user setable value for cando_sp 5257 logical def_cando_nw ! default user setable value for cando_nw 5258 logical def_cando_txs ! default user setable value for cando_txs 5259c 5260 logical app_stored_txs ! value stored in int_app_set_no_texas 5261 logical app_stored_spint ! value stored in int_app_set_no_spint 5262 integer rtdbIused 5263c 5264 common /clcando/ user_cando_sp, user_cando_nw, user_cando_txs, 5265 & def_cando_sp, def_cando_nw, def_cando_txs, 5266 & app_stored_txs, app_stored_spint, 5267 & rtdbIused 5268 5269c 5270 5271#include "global.fh" 5272#include "mafdecls.fh" 5273#include "rtdb.fh" 5274#include "nwc_const.fh" 5275#include "int_nbf.fh" 5276c::functions 5277 logical spcart_init 5278 external spcart_init 5279 logical int_ecp_init 5280 external int_ecp_init 5281c::passed 5282c:tex-\begin{verbatim} 5283 integer rtdb ! [input] run time data base handle 5284 integer nbas ! [input] number of basis sets to be used 5285 integer bases(nbas) ! [input] basis set handles 5286c:tex-\end{verbatim} 5287c::local 5288 integer txs_mem_min ! memory from texas 5289 integer ibas, ang2use, angm, type 5290 logical status 5291 integer nqmax_texas ! maximum number of quartets in texas blocking interface 5292 parameter (nqmax_texas = 10000) 5293c 5294c block data api_data 5295c 5296c 5297c Block data structure to initialize the common block variables in the 5298c internal basis set object data structures 5299c 5300c 5301 call int_mem_zero() 5302c 5303 DCexp = 0.0D00 5304 DCcoeff = 1.0D00 5305 val_int_acc = 0.0d00 5306c 5307 intd_memthresh = 0 5308 numd_tot = 0 5309 numd_okay = 0 5310 numd_red = 0 5311c 5312 if(init_int.eq.1) then 5313 write(6,*)' warning nested int_inits' 5314 write(6,*)' int_init already called ' 5315 call util_flush(6) 5316 endif 5317c 5318c initialize type-> nbf maps 5319c 5320 int_nbf_x(-1) = 4 5321 int_nbf_s(-1) = 4 5322 do type = 0,int_nbf_max_ang 5323 int_nbf_x(type) = (type+1)*(type+2)/2 5324 int_nbf_s(type) = 2*type+1 5325 enddo 5326 5327c 5328c initialize cando information from rtdb 5329c 5330 user_cando_sp = .false. 5331 user_cando_nw = .false. 5332 user_cando_txs = .false. 5333 def_cando_sp = .false. 5334 def_cando_nw = .false. 5335 def_cando_txs = .false. 5336c 5337 if (rtdb_get(rtdb,'int:cando_sp',MT_LOG,1,status)) then 5338 user_cando_sp = .true. 5339 def_cando_sp = status 5340* if (ga_nodeid().eq.0) then 5341* write(6,*) 5342* & ' int_init: cando_sp set to always be ',def_cando_sp 5343* call util_flush(6) 5344* endif 5345 endif 5346c 5347 if (rtdb_get(rtdb,'int:cando_nw',MT_LOG,1,status)) then 5348 user_cando_nw = .true. 5349 def_cando_nw = status 5350*rak: if (ga_nodeid().eq.0) then 5351*rak: write(6,*) 5352*rak: & ' int_init: cando_nw set to always be ',def_cando_nw 5353*rak: call util_flush(6) 5354*rak: endif 5355 endif 5356c 5357 if (rtdb_get(rtdb,'int:cando_txs',MT_LOG,1,status)) then 5358 user_cando_txs = .true. 5359 def_cando_txs = status 5360*rak: if (ga_nodeid().eq.0) then 5361*rak: write(6,*) 5362*rak: & ' int_init: cando_txs set to always be ',def_cando_txs 5363*rak: call util_flush(6) 5364*rak: endif 5365 endif 5366* sanity checking: e.g., you only want to turn off a particular integral 5367* code never always turn it on. 5368* 5369 if (def_cando_sp.or.def_cando_nw.or.def_cando_txs) then 5370 if (ga_nodeid().eq.0) then 5371 write(6,*)' you are trying to turn an integral code on ? ' 5372 write(6,*)' sp ', def_cando_sp 5373 write(6,*)' nw ', def_cando_nw 5374 write(6,*)' txs ', def_cando_txs 5375 call util_flush(6) 5376 endif 5377 call errquit 5378 & ('int_init: logic error with user cando settings',911, 5379 & INT_ERR) 5380 endif 5381c 5382 status = .true. 5383 do 00100 ibas=1,nbas 5384 status = status .and. bas_check_handle(bases(ibas),'int_init') 538500100 continue 5386 5387 if (.not.status) then 5388 write(6,*)' at least one basis handle not valid' 5389 do 00200 ibas = 1,nbas 5390 write(6,'(a,i5)') 5391 & ' basis set handle ',bases(ibas) 539200200 continue 5393 call errquit('int_init: basis handles hosed ',nbas, INT_ERR) 5394 endif 5395* write(6,*)' int_init: basis set handles valid ' 5396c 5397c check for both sp and gc shells 5398c 5399 call int_bothsp_gc_check(bases,nbas,'int_init') 5400c 5401c initialize defnxyz routines 5402c 5403 ang2use = -1 5404 do 00300 ibas = 1,nbas 5405 if(.not.bas_high_angular(bases(ibas),angm)) 5406 & call errquit('int_init: angm error',angm, INPUT_ERR) 5407 ang2use = max(ang2use,angm) 540800300 continue 5409* 5410* test for higher than h functions 0123456 5411 if (ang2use.ge.6) call errquit 5412 & ('only basis sets with s through h functions are allowed', 5413 & 911, BASIS_ERR) 5414* 5415 call defNxyz(ang2use) 5416c 5417c initialize spcart stuff 5418c 5419 if (.not.(spcart_init(ang2use,.true.,.false.))) then 5420 call errquit('int_init: spcart_init failed',911, INT_ERR) 5421 endif 5422c 5423c... generate memory requirements and store in structures in apiP.fh 5424c 5425 call exact_mem(rtdb,bases,nbas) 5426 call sp_init(nbas,bases) 5427 call init70 ! To generate tables etc. 5428 call int_acc_std() 5429* def u=f d=f -> f.and.!f -> f -> e = t 5430* no txs u=t d=f -> t.and.!f -> t -> e = f 5431 if (.not.(user_cando_txs.and.(.not.def_cando_txs))) then 5432 call texas_init(rtdb,nbas,bases,nqmax_texas,txs_mem_min, 5433 * 'scfd_int') 5434 endif 5435c 5436c See if any basis has an attached ECP 5437c 5438 any_ecp = .false. 5439 ecp_bsh = 0 5440 do ibas = 1,nbas 5441 if (bas_get_ecp_handle(bases(ibas),ecp_bsh)) then 5442 any_ecp = .true. 5443 goto 00001 5444 endif 5445 enddo 544600001 continue 5447 if (any_ecp) then 5448 if (.not.ecp_check_handle(ecp_bsh,'int_init')) call errquit 5449 & ('int_init: ecp handle is invalid fatal error',911, 5450 & INT_ERR) 5451 if (.not.int_ecp_init(ecp_bsh,0,0)) call errquit 5452 & ('int_init: int_ecp_init failed ',911, INT_ERR) 5453 endif 5454 init_int = 1 5455 end 5456 5457