1c $Id$ 2 logical function kgdtest (rtdb) 3 implicit none 4#include "mafdecls.fh" 5#include "rtdb.fh" 6#include "context.fh" 7#include "global.fh" 8#include "stdio.fh" 9c::passed 10 integer rtdb ! rtdb handle 11c::local 12 integer me 13 integer kgdtask, kgd_tmp 14c 15 call ga_sync() 16 call ga_sync() 17 kgdtask = 0 18 if (rtdb_get(rtdb,'kgdtask',MT_INT,1,kgd_tmp)) 19 & kgdtask = kgd_tmp 20c 21 call ga_sync() 22 call ga_sync() 23 me = ga_nodeid() 24 write (luout,*) 'in kgdtest' 25 call util_flush(luout) 26c 27 if (kgdtask.eq.0) then !................................... 0 28 if (me.eq.0) then 29 write(luout,*)' default kgdtest task ' 30 write(luout,*)' test use of kgdtest! ' 31 endif 32 kgdtest = .true. 33 else if (kgdtask.eq.1) then !................................. 1 34 if (me.eq.0) write(luout,*) 35 & ' kgdtest task 1: relativistic integral test' 36 call kgdtest_rel1e(rtdb) 37 kgdtest = .true. 38 else if (kgdtask.eq.2) then !................................. 2 39 if (me.eq.0) write(luout,*) 40 & ' kgdtest task 2: general contraction test' 41 call kgdtest_gencon(rtdb) 42 kgdtest = .true. 43 else if (kgdtask.eq.3) then !................................. 3 44 if (me.eq.0) write(luout,*) 45 & ' kgdtest task 3: so-ecp integral test' 46 call kgdtest_soecp(rtdb) 47 kgdtest = .true. 48 else if (kgdtask.eq.4) then !................................. 4 49 if (me.eq.0) write(luout,*) 50 & ' kgdtest task 4: rel2e integral test' 51 call kgdtest_rel2e(rtdb) 52 kgdtest = .true. 53 else if (kgdtask.eq.5) then !................................. 5 54 if (me.eq.0) write(luout,*) 55 & ' kgdtest task 5: ecp memory test' 56 call kgdtest_ecpmem(rtdb) 57 kgdtest = .true. 58 else if (kgdtask.eq.6) then !................................. 6 59 if (me.eq.0) write(luout,*) 60 & ' kgdtest task 6: rel general contraction integral test' 61 call kgdtest_relgen(rtdb) 62 kgdtest = .true. 63 end if 64 end 65************************************************************************ 66 subroutine kgdtest_rel1e(rtdb) 67 implicit none 68#include "errquit.fh" 69#include "mafdecls.fh" 70#include "rtdb.fh" 71#include "context.fh" 72#include "geom.fh" 73#include "bas.fh" 74#include "nwc_const.fh" 75#include "basP.fh" 76#include "basdeclsP.fh" 77#include "geomP.fh" 78#include "geobasmapP.fh" 79#include "bas_exndcf_dec.fh" 80#include "bas_ibs_dec.fh" 81c 82 logical int_normalize 83 external int_normalize 84c 85c test relativistic one-electron integrals 86c 87 integer rtdb 88 integer geom,basis, basis_id 89 integer nshell, memscr, membuf 90 integer h_scr, k_scr, h_buf, k_buf 91 integer ish, jsh, ucont 92 integer li, i_prim, i_gen, i_iexp, i_icfp, i_cent, i_geom 93 integer lj, j_prim, j_gen, j_iexp, j_icfp, j_cent, j_geom 94C integer nint_out 95 integer ihi,jhi 96 logical status 97 character*255 mo_basis, geom_name 98c 99#include "bas_exndcf_sfn.fh" 100#include "bas_ibs_sfn.fh" 101c 102 mo_basis = 'ao basis' 103 geom_name = 'geometry' 104c 105 if(.not.geom_create(geom,geom_name))call errquit 106 & ('kgdtest_rel1e: geom create error',911, GEOM_ERR) 107 if(.not.bas_create(basis,mo_basis))call errquit 108 & ('kgdtest_rel1e: basis create error',911, BASIS_ERR) 109c 110 if(.not.geom_rtdb_load(rtdb,geom,geom_name)) call errquit 111 & ('kgdtest_rel1e: geom load ',911, RTDB_ERR) 112 if(.not.bas_rtdb_load(rtdb,geom,basis,mo_basis)) call errquit 113 & ('kgdtest_rel1e: basis load ',911, RTDB_ERR) 114c 115 basis_id = basis + BASIS_HANDLE_OFFSET 116 nshell = ncont_tot_gb(basis_id) 117 if (.not.int_normalize(rtdb,basis)) call errquit 118 & ('kgdtest_rel1e: error normalizing ',911, INT_ERR) 119c 120 call int_init(rtdb,1,basis) 121 memscr = 100 000 122 membuf = 1000 123 if (.not.ma_push_get(mt_dbl,memscr,' scratch ', 124 & h_scr, k_scr)) call errquit 125 & (' ma error 1',911, MA_ERR) 126 if (.not.ma_push_get(mt_dbl,membuf*3,' buf ', 127 & h_buf, k_buf)) call errquit 128 & (' ma error 2',911, MA_ERR) 129c 130 do ish = 1,nshell 131 do jsh = 1,ish 132 write(6,*)' ============= shells <',ish,'|',jsh,'>', 133 & '==================== start ==========' 134 write(6,*)' ' 135 136 ucont = (sf_ibs_cn2ucn(ish,basis_id)) 137 Li = infbs_cont(CONT_TYPE ,ucont,basis_id) 138 i_prim = infbs_cont(CONT_NPRIM,ucont,basis_id) 139 i_gen = infbs_cont(CONT_NGEN ,ucont,basis_id) 140 i_iexp = infbs_cont(CONT_IEXP ,ucont,basis_id) 141 i_icfp = infbs_cont(CONT_ICFP ,ucont,basis_id) 142 i_cent = (sf_ibs_cn2ce(ish,basis_id)) 143 i_geom = ibs_geom(basis_id) 144 ihi = i_gen*(Li+1)*(Li+2)/2 145c 146 ucont = (sf_ibs_cn2ucn(jsh,basis_id)) 147 Lj = infbs_cont(CONT_TYPE ,ucont,basis_id) 148 j_prim = infbs_cont(CONT_NPRIM,ucont,basis_id) 149 j_gen = infbs_cont(CONT_NGEN ,ucont,basis_id) 150 j_iexp = infbs_cont(CONT_IEXP ,ucont,basis_id) 151 j_icfp = infbs_cont(CONT_ICFP ,ucont,basis_id) 152 j_cent = (sf_ibs_cn2ce(jsh,basis_id)) 153 j_geom = ibs_geom(basis_id) 154 jhi = j_gen*(Lj+1)*(Lj+2)/2 155* 156* Calculate overlap and kinetic energy integrals 157* 158 call int_hf1sp( 159 & coords(1,i_cent,i_geom), 160 & dbl_mb(mb_exndcf(i_iexp,basis_id)), 161 & dbl_mb(mb_exndcf(i_icfp,basis_id)),i_prim,i_gen,Li,i_cent, 162 & coords(1,j_cent,j_geom), 163 & dbl_mb(mb_exndcf(j_iexp,basis_id)), 164 & dbl_mb(mb_exndcf(j_icfp,basis_id)),j_prim,j_gen,Lj,j_cent, 165 & coords(1,1,i_geom),charge(1,i_geom),ncenter(i_geom), 166 & dbl_mb(k_buf),dbl_mb(k_buf+1000),dbl_mb(k_buf+2000), 167 & .true.,.true.,.false.,.false., 168 & .false.,dbl_mb(k_scr),memscr,'kgdtest_rel1e') 169 call ecp_matpr (dbl_mb(k_buf),1,jhi,1,ihi,1,jhi,1,ihi, 170 & 'overlap integrals','E',120,8) 171 call ecp_matpr (dbl_mb(k_buf+1000),1,jhi,1,ihi,1,jhi,1,ihi, 172 & 'kinetic integrals','E',120,8) 173C modified metric 174 call rel_onel ( 175 & coords(1,i_cent,i_geom), 176 & dbl_mb(mb_exndcf(i_iexp,basis_id)), 177 & dbl_mb(mb_exndcf(i_icfp,basis_id)), 178 & dbl_mb(mb_exndcf(i_icfp,basis_id)),i_prim,i_gen,Li, 179 & coords(1,j_cent,j_geom), 180 & dbl_mb(mb_exndcf(j_iexp,basis_id)), 181 & dbl_mb(mb_exndcf(j_icfp,basis_id)), 182 & dbl_mb(mb_exndcf(j_icfp,basis_id)),j_prim,j_gen,Lj, 183 & coords(1,1,i_geom),charge(1,i_geom), 184 & geom_invnucexp(1,i_geom),ncenter(i_geom), 185 & dbl_mb(k_buf),dbl_mb(k_buf+1000),dbl_mb(k_buf+2000), 186 & membuf,.true.,.false.,.false.,.false.,.true.,.false., 187 & .false.,.false.,dbl_mb(k_scr),memscr,0,1) 188 call ecp_matpr (dbl_mb(k_buf),1,jhi,1,ihi,1,jhi,1,ihi, 189 & 'modified overlap integrals','E',120,8) 190C modified kinetic energy 191 call rel_onel ( 192 & coords(1,i_cent,i_geom), 193 & dbl_mb(mb_exndcf(i_iexp,basis_id)), 194 & dbl_mb(mb_exndcf(i_icfp,basis_id)), 195 & dbl_mb(mb_exndcf(i_icfp,basis_id)),i_prim,i_gen,Li, 196 & coords(1,j_cent,j_geom), 197 & dbl_mb(mb_exndcf(j_iexp,basis_id)), 198 & dbl_mb(mb_exndcf(j_icfp,basis_id)), 199 & dbl_mb(mb_exndcf(j_icfp,basis_id)),j_prim,j_gen,Lj, 200 & coords(1,1,i_geom),charge(1,i_geom), 201 & geom_invnucexp(1,i_geom),ncenter(i_geom), 202 & dbl_mb(k_buf),dbl_mb(k_buf+1000),dbl_mb(k_buf+2000), 203 & membuf,.false.,.true.,.false.,.false.,.true.,.false., 204 & .false.,.false.,dbl_mb(k_scr),memscr,0,1) 205 call ecp_matpr (dbl_mb(k_buf+1000),1,jhi,1,ihi,1,jhi,1,ihi, 206 & 'modified kinetic integrals','E',120,8) 207C modified potential energy 208 call rel_onel ( 209 & coords(1,i_cent,i_geom), 210 & dbl_mb(mb_exndcf(i_iexp,basis_id)), 211 & dbl_mb(mb_exndcf(i_icfp,basis_id)), 212 & dbl_mb(mb_exndcf(i_icfp,basis_id)),i_prim,i_gen,Li, 213 & coords(1,j_cent,j_geom), 214 & dbl_mb(mb_exndcf(j_iexp,basis_id)), 215 & dbl_mb(mb_exndcf(j_icfp,basis_id)), 216 & dbl_mb(mb_exndcf(j_icfp,basis_id)),j_prim,j_gen,Lj, 217 & coords(1,1,i_geom),charge(1,i_geom), 218 & geom_invnucexp(1,i_geom),ncenter(i_geom), 219 & dbl_mb(k_buf),dbl_mb(k_buf+1000),dbl_mb(k_buf+2000), 220 & membuf,.false.,.false.,.true.,.false.,.true.,.false., 221 & .false.,.false.,dbl_mb(k_scr),memscr,0,1) 222 call ecp_matpr (dbl_mb(k_buf+2000),1,jhi,1,ihi,1,jhi,1,ihi, 223 & 'modified potential integrals','E',120,8) 224 write(6,*)' ============= shells <',ish,'|',jsh,'>', 225 & '==================== end ==========' 226 write(6,*)' ' 227 enddo 228 enddo 229c 230 call int_terminate() 231 status = ma_pop_stack(h_buf) 232 status = status.and.ma_pop_stack(h_scr) 233 if (.not.status) call errquit('pop failed',911, MA_ERR) 234 status = bas_destroy(basis) 235 status = status.and.geom_destroy(geom) 236 if (.not.status) call errquit('b/g destroy failed',911, GEOM_ERR) 237 return 238 end 239************************************************************************ 240 subroutine kgdtest_gencon(rtdb) 241 implicit none 242#include "errquit.fh" 243#include "mafdecls.fh" 244#include "rtdb.fh" 245#include "context.fh" 246#include "geom.fh" 247#include "bas.fh" 248#include "nwc_const.fh" 249#include "basP.fh" 250#include "basdeclsP.fh" 251#include "geomP.fh" 252#include "geobasmapP.fh" 253#include "bas_exndcf_dec.fh" 254#include "bas_ibs_dec.fh" 255#include "stdio.fh" 256c 257 logical int_normalize 258 external int_normalize 259 integer idamax 260 external idamax 261c 262c test general contraction in McMD 2e integrals. 263c 264 integer rtdb 265 integer geom,basis, basis_id 266 integer nshell, memscr, membuf 267 integer h_scr, k_scr, h_buf, k_buf 268 integer ish, jsh, ksh,lsh, ucont 269 integer li, i_prim, i_gen, i_iexp, i_icfp, i_cent, i_geom 270 integer lj, j_prim, j_gen, j_iexp, j_icfp, j_cent, j_geom 271 integer lk, k_prim, k_gen, k_iexp, k_icfp, k_cent, k_geom 272 integer ll, l_prim, l_gen, l_iexp, l_icfp, l_cent, l_geom 273 integer Nints 274 integer i_eri,i_c,j_c,k_c,l_c,ic,jc,kc,lc,i_seg 275 integer ihi,jhi,khi,lhi,i2,j2,k2,l2,n_cart,n_cont 276 character*255 mo_basis, geom_name 277 double precision difmax 278c 279#include "bas_exndcf_sfn.fh" 280#include "bas_ibs_sfn.fh" 281c 282 mo_basis = 'ao basis' 283 geom_name = 'geometry' 284c 285 if(.not.geom_create(geom,geom_name))call errquit 286 & ('kgdtest_gencon: geom create error',911, GEOM_ERR) 287 if(.not.bas_create(basis,mo_basis))call errquit 288 & ('kgdtest_gencon: basis create error',911, BASIS_ERR) 289c 290 if(.not.geom_rtdb_load(rtdb,geom,geom_name)) call errquit 291 & ('kgdtest_gencon: geom load ',911, RTDB_ERR) 292 if(.not.bas_rtdb_load(rtdb,geom,basis,mo_basis)) call errquit 293 & ('kgdtest_gencon: basis load ',911, RTDB_ERR) 294c 295 basis_id = basis + BASIS_HANDLE_OFFSET 296 nshell = ncont_tot_gb(basis_id) 297 if (.not.int_normalize(rtdb,basis)) call errquit 298 & ('kgdtest_gencon: error normalizing ',911, INT_ERR) 299c 300 call int_init(rtdb,1,basis) 301 memscr = 1 000 000 302 membuf = 10 000 303 if (.not.ma_push_get(mt_dbl,memscr,' scratch ', 304 & h_scr, k_scr)) call errquit 305 & (' ma error 1',911, MA_ERR) 306 if (.not.ma_push_get(mt_dbl,membuf,' buf ', 307 & h_buf, k_buf)) call errquit 308 & (' ma error 2',911, MA_ERR) 309c 310 difmax = 0.0d0 311 do ish = 1,nshell 312 ucont = (sf_ibs_cn2ucn(ish,basis_id)) 313 Li = infbs_cont(CONT_TYPE ,ucont,basis_id) 314 i_prim = infbs_cont(CONT_NPRIM,ucont,basis_id) 315 i_gen = infbs_cont(CONT_NGEN ,ucont,basis_id) 316 i_iexp = infbs_cont(CONT_IEXP ,ucont,basis_id) 317 i_icfp = infbs_cont(CONT_ICFP ,ucont,basis_id) 318 i_cent = (sf_ibs_cn2ce(ish,basis_id)) 319 i_geom = ibs_geom(basis_id) 320 i2 = (Li+1)*(Li+2)/2 321 ihi = i_gen*i2 322c 323 do jsh = 1,ish 324 ucont = (sf_ibs_cn2ucn(jsh,basis_id)) 325 Lj = infbs_cont(CONT_TYPE ,ucont,basis_id) 326 j_prim = infbs_cont(CONT_NPRIM,ucont,basis_id) 327 j_gen = infbs_cont(CONT_NGEN ,ucont,basis_id) 328 j_iexp = infbs_cont(CONT_IEXP ,ucont,basis_id) 329 j_icfp = infbs_cont(CONT_ICFP ,ucont,basis_id) 330 j_cent = (sf_ibs_cn2ce(jsh,basis_id)) 331 j_geom = ibs_geom(basis_id) 332 j2 = (Lj+1)*(Lj+2)/2 333 jhi = j_gen*j2 334c 335 do ksh = 1,ish 336 ucont = (sf_ibs_cn2ucn(ksh,basis_id)) 337 Lk = infbs_cont(CONT_TYPE ,ucont,basis_id) 338 k_prim = infbs_cont(CONT_NPRIM,ucont,basis_id) 339 k_gen = infbs_cont(CONT_NGEN ,ucont,basis_id) 340 k_iexp = infbs_cont(CONT_IEXP ,ucont,basis_id) 341 k_icfp = infbs_cont(CONT_ICFP ,ucont,basis_id) 342 k_cent = (sf_ibs_cn2ce(ish,basis_id)) 343 k_geom = ibs_geom(basis_id) 344 k2 = (Lk+1)*(Lk+2)/2 345 khi = k_gen*k2 346c 347 do lsh = 1,ksh 348 ucont = (sf_ibs_cn2ucn(lsh,basis_id)) 349 Ll = infbs_cont(CONT_TYPE ,ucont,basis_id) 350 l_prim = infbs_cont(CONT_NPRIM,ucont,basis_id) 351 l_gen = infbs_cont(CONT_NGEN ,ucont,basis_id) 352 l_iexp = infbs_cont(CONT_IEXP ,ucont,basis_id) 353 l_icfp = infbs_cont(CONT_ICFP ,ucont,basis_id) 354 l_cent = (sf_ibs_cn2ce(lsh,basis_id)) 355 l_geom = ibs_geom(basis_id) 356 l2 = (Ll+1)*(Ll+2)/2 357 lhi = l_gen*l2 358 359 n_cart = i2*j2*k2*l2 360 n_cont = i_gen*j_gen*k_gen*l_gen 361 Nints = ihi*jhi*khi*lhi 362 write (LuOut,*) ish,Li,i_prim,i_gen,i_cent 363 write (LuOut,*) jsh,Lj,j_prim,j_gen,j_cent 364 write (LuOut,*) ksh,Lk,k_prim,k_gen,k_cent 365 write (LuOut,*) lsh,Ll,l_prim,l_gen,l_cent 366C if (n_cont .gt. 1 .and. Li.ne.Lj .and. Lk.ne.Ll) then 367C write(luout,*)' ' 368C write(luout,*)' ============= (', 369C & ish,':',Li,',', 370C & jsh,':',Lj,'|', 371C & ksh,':',Lk,',', 372C & lsh,':',Ll,')', 373C & '====================' 374C write(luout,*)' ' 375 write (luout,*) 'Nints',Nints 376 377 call hf2( 378 & coords(1,i_cent,i_geom), 379 & dbl_mb(mb_exndcf(i_iexp,basis_id)), 380 & dbl_mb(mb_exndcf(i_icfp,basis_id)),i_prim,i_gen,Li, 381 & coords(1,j_cent,j_geom), 382 & dbl_mb(mb_exndcf(j_iexp,basis_id)), 383 & dbl_mb(mb_exndcf(j_icfp,basis_id)),j_prim,j_gen,Lj, 384 & coords(1,k_cent,k_geom), 385 & dbl_mb(mb_exndcf(k_iexp,basis_id)), 386 & dbl_mb(mb_exndcf(k_icfp,basis_id)),k_prim,k_gen,Lk, 387 & coords(1,l_cent,l_geom), 388 & dbl_mb(mb_exndcf(l_iexp,basis_id)), 389 & dbl_mb(mb_exndcf(l_icfp,basis_id)),l_prim,l_gen,Ll, 390 & dbl_mb(k_buf),Nints,.false.,.false.,.false., 391 & .false.,dbl_mb(k_scr),memscr) 392 393 i_eri = k_buf+Nints 394 i_seg = i_eri 395 i_c = i_icfp 396 do ic = 1,i_gen 397 j_c = j_icfp 398 do jc = 1,j_gen 399 k_c = k_icfp 400 do kc = 1,k_gen 401 l_c = l_icfp 402 do lc = 1,l_gen 403C write (luout,*) 'ic,jc,kc,lc',ic,jc,kc,lc 404 call hf2( 405 & coords(1,i_cent,i_geom), 406 & dbl_mb(mb_exndcf(i_iexp,basis_id)), 407 & dbl_mb(mb_exndcf(i_c,basis_id)), 408 & i_prim,1,Li, 409 & coords(1,j_cent,j_geom), 410 & dbl_mb(mb_exndcf(j_iexp,basis_id)), 411 & dbl_mb(mb_exndcf(j_c,basis_id)), 412 & j_prim,1,Lj, 413 & coords(1,k_cent,k_geom), 414 & dbl_mb(mb_exndcf(k_iexp,basis_id)), 415 & dbl_mb(mb_exndcf(k_c,basis_id)), 416 & k_prim,1,Lk, 417 & coords(1,l_cent,l_geom), 418 & dbl_mb(mb_exndcf(l_iexp,basis_id)), 419 & dbl_mb(mb_exndcf(l_c,basis_id)), 420 & l_prim,1,Ll, 421 & dbl_mb(i_eri),Nints,.false.,.false.,.false., 422 & .false.,dbl_mb(k_scr),memscr) 423 424 i_eri = i_eri+i2*j2*k2*l2 425 l_c = l_c+l_prim 426 end do 427 k_c = k_c+k_prim 428 end do 429 j_c = j_c+j_prim 430 end do 431 i_c = i_c+i_prim 432 end do 433 434 call reorder_2eints(dbl_mb(i_eri),dbl_mb(i_seg), 435 & l2,l_gen,k2,k_gen,j2,j_gen,i2,i_gen) 436 call dcopy (Nints,dbl_mb(i_eri),1,dbl_mb(i_seg),1) 437 call daxpy(Nints,-1.0d0,dbl_mb(k_buf),1,dbl_mb(i_eri),1) 438 ic = idamax(Nints,dbl_mb(i_eri),1)-1 439 write (luout,*) 'Maximum difference',dbl_mb(i_eri+ic) 440 difmax = max(difmax,abs(dbl_mb(i_eri+ic))) 441 if (abs(dbl_mb(i_eri+ic)) .gt. 1.0d-12) then 442 call ecp_matpr (dbl_mb(k_buf),1,khi*lhi,1,ihi*jhi, 443 & 1,khi*lhi,1,ihi*jhi,'General contraction', 444 & 'E',120,6) 445 call ecp_matpr (dbl_mb(i_seg),1,khi*lhi,1,ihi*jhi, 446 & 1,khi*lhi,1,ihi*jhi,'Segmented contraction', 447 & 'E',120,6) 448 end if 449C end if 450 end do 451 end do 452 end do 453 end do 454 write (luout,*) 'Maximum difference of all blocks',difmax 455 456 end 457************************************************************************ 458 subroutine reorder_2eints(bnew,bold,lc,lg,kc,kg,jc,jg,ic,ig) 459 double precision bnew(lc,lg,kc,kg,jc,jg,ic*ig) 460 double precision bold(lc,kc,jc,ic,lg,kg,jg*ig) 461 integer lc,lg,kc,kg,jc,jg,ic,ig 462 integer i,j,k,l,ii,jj,kk,ll,inew,iold 463 iold = 0 464 inew = 0 465 do i = 1,ig 466 do j = 1,jg 467 iold = iold+1 468 do k = 1,kg 469 do l = 1,lg 470 do ii = 1,ic 471 do jj = 1,jc 472 do kk = 1,kc 473 do ll = 1,lc 474 bnew(ll,l,kk,k,jj,j,inew+ii) = 475 & bold(ll,kk,jj,ii,l,k,iold) 476 end do 477 end do 478 end do 479 end do 480 end do 481 end do 482 end do 483 inew = inew+ic 484 end do 485 end 486 487************************************************************************ 488 subroutine kgdtest_soecp(rtdb) 489 implicit none 490#include "errquit.fh" 491#include "mafdecls.fh" 492#include "rtdb.fh" 493#include "context.fh" 494#include "geom.fh" 495#include "bas.fh" 496#include "nwc_const.fh" 497#include "basP.fh" 498#include "basdeclsP.fh" 499#include "geomP.fh" 500#include "geobasmapP.fh" 501#include "bas_exndcf_dec.fh" 502#include "bas_ibs_dec.fh" 503#include "ecp_nwc.fh" 504#include "stdio.fh" 505c 506 logical int_normalize 507 external int_normalize 508 integer idamax 509 external idamax 510c 511c test spin-orbit integrals 512c 513 integer rtdb 514 integer geom,basis, basis_id 515 integer nshell, memscr, membuf 516 integer h_scr, k_scr, h_buf, k_buf 517 integer ish, jsh, ucont 518 integer li, i_prim, i_gen, i_iexp, i_icfp, i_cent, i_geom 519 integer lj, j_prim, j_gen, j_iexp, j_icfp, j_cent, j_geom 520 integer ihi,jhi,i2,j2 521 integer ibug,n_blk 522 character*255 mo_basis, geom_name 523c 524#include "bas_exndcf_sfn.fh" 525#include "bas_ibs_sfn.fh" 526c 527 mo_basis = 'ao basis' 528 geom_name = 'geometry' 529c 530 if(.not.geom_create(geom,geom_name))call errquit 531 & ('kgdtest_gencon: geom create error',911, GEOM_ERR) 532 if(.not.bas_create(basis,mo_basis))call errquit 533 & ('kgdtest_gencon: basis create error',911, BASIS_ERR) 534c 535 if(.not.geom_rtdb_load(rtdb,geom,geom_name)) call errquit 536 & ('kgdtest_gencon: geom load ',911, RTDB_ERR) 537 if(.not.bas_rtdb_load(rtdb,geom,basis,mo_basis)) call errquit 538 & ('kgdtest_gencon: basis load ',911, RTDB_ERR) 539c 540 basis_id = basis + BASIS_HANDLE_OFFSET 541 nshell = ncont_tot_gb(basis_id) 542 if (.not.int_normalize(rtdb,basis)) call errquit 543 & ('kgdtest_gencon: error normalizing ',911, INT_ERR) 544c 545 call int_init(rtdb,1,basis) 546 memscr = 1 000 000 547 membuf = 10 000 548 if (.not.ma_push_get(mt_dbl,memscr,' scratch ', 549 & h_scr, k_scr)) call errquit 550 & (' ma error 1',911, MA_ERR) 551 if (.not.ma_push_get(mt_dbl,membuf,' buf ', 552 & h_buf, k_buf)) call errquit 553 & (' ma error 2',911, MA_ERR) 554c 555 if (.not. rtdb_get(rtdb,'ecp:ibug',mt_int,1,ibug)) 556 & ibug = 3 557 if (.not. rtdb_get(rtdb,'ecp:n_blk',mt_int,1,n_blk)) 558 & n_blk = 3 559c 560 do ish = 1,nshell 561 ucont = (sf_ibs_cn2ucn(ish,basis_id)) 562 Li = infbs_cont(CONT_TYPE ,ucont,basis_id) 563 i_prim = infbs_cont(CONT_NPRIM,ucont,basis_id) 564 i_gen = infbs_cont(CONT_NGEN ,ucont,basis_id) 565 i_iexp = infbs_cont(CONT_IEXP ,ucont,basis_id) 566 i_icfp = infbs_cont(CONT_ICFP ,ucont,basis_id) 567 i_cent = (sf_ibs_cn2ce(ish,basis_id)) 568 i_geom = ibs_geom(basis_id) 569 i2 = (Li+1)*(Li+2)/2 570 ihi = i_gen*i2 571c 572 do jsh = 1,ish 573 ucont = (sf_ibs_cn2ucn(jsh,basis_id)) 574 Lj = infbs_cont(CONT_TYPE ,ucont,basis_id) 575 j_prim = infbs_cont(CONT_NPRIM,ucont,basis_id) 576 j_gen = infbs_cont(CONT_NGEN ,ucont,basis_id) 577 j_iexp = infbs_cont(CONT_IEXP ,ucont,basis_id) 578 j_icfp = infbs_cont(CONT_ICFP ,ucont,basis_id) 579 j_cent = (sf_ibs_cn2ce(jsh,basis_id)) 580 j_geom = ibs_geom(basis_id) 581 j2 = (Lj+1)*(Lj+2)/2 582 jhi = j_gen*j2 583 584 write(luout,*)' ' 585 write(luout,*)' ============= ',ish,jsh, 586 & '====================' 587 write(luout,*)' ' 588 589 call ecp_integral ( 590 & coords(1,i_cent,i_geom), 591 & dbl_mb(mb_exndcf(i_iexp,basis_id)), 592 & dbl_mb(mb_exndcf(i_icfp,basis_id)),i_prim,i_gen,Li,i_cent, 593 & coords(1,j_cent,j_geom), 594 & dbl_mb(mb_exndcf(j_iexp,basis_id)), 595 & dbl_mb(mb_exndcf(j_icfp,basis_id)),j_prim,j_gen,Lj,j_cent, 596 & dbl_mb(k_xyzecp), 597 & dbl_mb(k_ecp_e),dbl_mb(k_ecp_c), 598 & int_mb(k_ecp_nprim_c), 599 & int_mb(k_ecp_ncoef_c), ! new name is n_colc_C 600 & int_mb(k_ecp_ind_z), 601 & int_mb(k_ecp_ind_c), 602 & n_zeta_c, 603 & n_zeta_c, 604 & int_mb(k_ecp_l_c), 605 & int_mb(k_ecp_lip), 606 & n_ecp,l_ecp, 607 & 0, 608 & dbl_mb(k_ecp_c2s),mem_c2s, 609 & dbl_mb(k_buf),ihi*jhi,n_blk, 610 & .false.,dbl_mb(k_scr),memscr, 611 & ibug) 612 613 end do 614 end do 615 call int_terminate() 616 end 617************************************************************************ 618 subroutine kgdtest_rel2e(rtdb) 619 implicit none 620#include "errquit.fh" 621#include "mafdecls.fh" 622#include "rtdb.fh" 623#include "context.fh" 624#include "geom.fh" 625#include "bas.fh" 626#include "nwc_const.fh" 627#include "basP.fh" 628#include "basdeclsP.fh" 629#include "geomP.fh" 630#include "geobasmapP.fh" 631#include "bas_exndcf_dec.fh" 632#include "bas_ibs_dec.fh" 633#include "apiP.fh" 634#include "rel_nwc.fh" 635#include "stdio.fh" 636c 637 logical int_normalize 638 external int_normalize 639 integer idamax 640 external idamax 641c 642c test relativistic 2e integrals for correct magnitude 643c 644 integer rtdb 645 integer geom, basis, basis_id 646 integer nshell, memscr, membuf 647 integer h_scr, k_scr, h_buf, k_buf 648 integer ish, jsh, ksh,lsh, ucont 649 integer li, ipr, igc, iex, icf, icfs, ice, igm 650 integer lj, jpr, jgc, jex, jcf, jcfs, jce, jgm 651 integer lk, kpr, kgc, kex, kcf, kcfs, kce, kgm 652 integer ll, lpr, lgc, lex, lcf, lcfs, lce, lgm 653 integer Nints 654 integer i_eri,i,abas,sbas 655 integer ihi,jhi,khi,lhi,i2,j2,k2,l2,n_cart,n_cont 656 character*255 mo_basis, geom_name 657 double precision difmax,errmax,dif,erimax 658c 659#include "bas_exndcf_sfn.fh" 660#include "bas_ibs_sfn.fh" 661c 662 mo_basis = 'ao basis' 663 geom_name = 'geometry' 664c 665 if(.not.geom_create(geom,geom_name))call errquit 666 & ('kgdtest_rel2e: geom create error',911, GEOM_ERR) 667 if(.not.bas_create(basis,mo_basis))call errquit 668 & ('kgdtest_rel2e: basis create error',911, BASIS_ERR) 669c 670 if(.not.geom_rtdb_load(rtdb,geom,geom_name)) call errquit 671 & ('kgdtest_rel2e: geom load ',911, GEOM_ERR) 672 if(.not.bas_rtdb_load(rtdb,geom,basis,mo_basis)) call errquit 673 & ('kgdtest_rel2e: basis load ',911, BASIS_ERR) 674c 675 basis_id = basis + BASIS_HANDLE_OFFSET 676 nshell = ncont_tot_gb(basis_id) 677 if (.not.int_normalize(rtdb,basis)) call errquit 678 & ('kgdtest_rel2e: error normalizing ',911, INT_ERR) 679c 680 call int_init(rtdb,1,basis) 681 memscr = 5 000 000 682 membuf = 100 000 683 if (.not.ma_push_get(mt_dbl,memscr,' scratch ', 684 & h_scr, k_scr)) call errquit 685 & (' ma error 1',911, MA_ERR) 686 if (.not.ma_push_get(mt_dbl,membuf,' buf ', 687 & h_buf, k_buf)) call errquit 688 & (' ma error 2',911, MA_ERR) 689c 690 sbas = sc_bsh + BASIS_HANDLE_OFFSET 691 abas = ao_bsh + BASIS_HANDLE_OFFSET 692c 693 difmax = 0.0d0 694 errmax = 0.0d0 695 rel_dbg = 3 696C do ish = 1,nshell 697C do ish = 2,12,10 698 ish = 12 699 ucont = (sf_ibs_cn2ucn(ish,abas)) 700 Li = infbs_cont(CONT_TYPE ,ucont,abas) 701 ipr = infbs_cont(CONT_NPRIM,ucont,abas) 702 igc = infbs_cont(CONT_NGEN ,ucont,abas) 703 iex = infbs_cont(CONT_IEXP ,ucont,abas) 704 icf = infbs_cont(CONT_ICFP ,ucont,abas) 705 ice = (sf_ibs_cn2ce(ish,abas)) 706 igm = ibs_geom(abas) 707 i2 = (Li+1)*(Li+2)/2 708 ihi = igc*i2 709 ucont = ao_to_ls(ucont) 710 icfs = infbs_cont(CONT_ICFP ,ucont,sbas) 711c 712C do jsh = 1,ish 713 jsh = 11 714C jsh = ish-1 715 ucont = (sf_ibs_cn2ucn(jsh,abas)) 716 Lj = infbs_cont(CONT_TYPE ,ucont,abas) 717 jpr = infbs_cont(CONT_NPRIM,ucont,abas) 718 jgc = infbs_cont(CONT_NGEN ,ucont,abas) 719 jex = infbs_cont(CONT_IEXP ,ucont,abas) 720 jcf = infbs_cont(CONT_ICFP ,ucont,abas) 721 jce = (sf_ibs_cn2ce(jsh,abas)) 722 jgm = ibs_geom(abas) 723 j2 = (Lj+1)*(Lj+2)/2 724 jhi = jgc*j2 725 ucont = ao_to_ls(ucont) 726 jcfs = infbs_cont(CONT_ICFP ,ucont,sbas) 727c 728C do ksh = 1,ish 729 ksh = 1 730 ucont = (sf_ibs_cn2ucn(ksh,abas)) 731 Lk = infbs_cont(CONT_TYPE ,ucont,abas) 732 kpr = infbs_cont(CONT_NPRIM,ucont,abas) 733 kgc = infbs_cont(CONT_NGEN ,ucont,abas) 734 kex = infbs_cont(CONT_IEXP ,ucont,abas) 735 kcf = infbs_cont(CONT_ICFP ,ucont,abas) 736 kce = (sf_ibs_cn2ce(ksh,abas)) 737 kgm = ibs_geom(abas) 738 k2 = (Lk+1)*(Lk+2)/2 739 khi = kgc*k2 740 ucont = ao_to_ls(ucont) 741 kcfs = infbs_cont(CONT_ICFP ,ucont,sbas) 742c 743C do lsh = 1,ksh 744 lsh = 1 745 ucont = (sf_ibs_cn2ucn(lsh,abas)) 746 Ll = infbs_cont(CONT_TYPE ,ucont,abas) 747 lpr = infbs_cont(CONT_NPRIM,ucont,abas) 748 lgc = infbs_cont(CONT_NGEN ,ucont,abas) 749 lex = infbs_cont(CONT_IEXP ,ucont,abas) 750 lcf = infbs_cont(CONT_ICFP ,ucont,abas) 751 lce = (sf_ibs_cn2ce(lsh,abas)) 752 lgm = ibs_geom(abas) 753 l2 = (Ll+1)*(Ll+2)/2 754 lhi = lgc*l2 755 ucont = ao_to_ls(ucont) 756 lcfs = infbs_cont(CONT_ICFP ,ucont,sbas) 757 758 n_cart = i2*j2*k2*l2 759 n_cont = igc*jgc*kgc*lgc 760 Nints = ihi*jhi*khi*lhi 761 i_eri = k_buf+Nints 762C if (n_cont .gt. 1 .and. Li.ne.Lj .and. Lk.ne.Ll) then 763 write(luout,*)' ' 764 write(luout,*)' ============= (', 765 & ish,':',Li,',', 766 & jsh,':',Lj,'|', 767 & ksh,':',Lk,',', 768 & lsh,':',Ll,')', 769 & '====================' 770 write(luout,*)' ' 771C write (luout,*) 'Nints',Nints 772 773 call hf2( 774 & coords(1,ice,igm), 775 & dbl_mb(mb_exndcf(iex,abas)), 776 & dbl_mb(mb_exndcf(icf,abas)),ipr,igc,Li, 777 & coords(1,jce,jgm), 778 & dbl_mb(mb_exndcf(jex,abas)), 779 & dbl_mb(mb_exndcf(jcf,abas)),jpr,jgc,Lj, 780 & coords(1,kce,kgm), 781 & dbl_mb(mb_exndcf(kex,abas)), 782 & dbl_mb(mb_exndcf(kcf,abas)),kpr,kgc,Lk, 783 & coords(1,lce,lgm), 784 & dbl_mb(mb_exndcf(lex,abas)), 785 & dbl_mb(mb_exndcf(lcf,abas)),lpr,lgc,Ll, 786 & dbl_mb(k_buf),Nints,.false.,.false.,.false., 787 & .false.,dbl_mb(k_scr),memscr) 788 call ecp_matpr (dbl_mb(k_buf),1,khi*lhi,1,ihi*jhi, 789 & 1,khi*lhi,1,ihi*jhi,'Large component only', 790 & 'E',120,6) 791C write (luout,*) 'hf2 completed' 792C call util_flush(6) 793 call rel_2e4c_sf ( 794 & coords(1,ice,igm),dbl_mb(mb_exndcf(iex,abas)), 795 & dbl_mb(mb_exndcf(icf,abas)), 796 & dbl_mb(mb_exndcf(icfs,sbas)),ipr,igc,Li,ice, 797 & coords(1,jce,jgm),dbl_mb(mb_exndcf(jex,abas)), 798 & dbl_mb(mb_exndcf(jcf,abas)), 799 & dbl_mb(mb_exndcf(jcfs,sbas)),jpr,jgc,Lj,jce, 800 & coords(1,kce,kgm),dbl_mb(mb_exndcf(kex,abas)), 801 & dbl_mb(mb_exndcf(kcf,abas)), 802 & dbl_mb(mb_exndcf(kcfs,sbas)),kpr,kgc,Lk,kce, 803 & coords(1,lce,lgm),dbl_mb(mb_exndcf(lex,abas)), 804 & dbl_mb(mb_exndcf(lcf,abas)), 805 & dbl_mb(mb_exndcf(lcfs,sbas)),lpr,lgc,Ll,lce, 806c...................... canAB canCD canPQ DryRun 807 & dbl_mb(i_eri),Nints,.false.,.false.,.false.,.false., 808 & dbl_mb(k_scr),memscr, 809 & .true.,.true.,ss_one_cent,do_ssss,rel_dbg) 810 811 call ecp_matpr (dbl_mb(i_eri),1,khi*lhi,1,ihi*jhi, 812 & 1,khi*lhi,1,ihi*jhi,'Relativistic integrals', 813 & 'E',120,6) 814C write (luout,*) 'rel_2e4c_sf completed' 815C call util_flush(6) 816 call daxpy(Nints,-1.0d0,dbl_mb(k_buf),1,dbl_mb(i_eri),1) 817 dif = 0.0d0 818 erimax = 0.0d0 819 do i = 0,Nints-1 820 erimax = max(erimax,abs(dbl_mb(k_buf+i))) 821 if (abs(dbl_mb(k_buf+i)) .gt. 1.0d-12) then 822 dif = max(dif,abs(dbl_mb(i_eri+i)/dbl_mb(k_buf+i))) 823 else if (abs(dbl_mb(i_eri+i)) .gt. 1.0d-12) then 824 if (abs(dbl_mb(k_buf+i)) .ne. 0.0d0) then 825 dif = max(dif,abs(dbl_mb(i_eri+i)/dbl_mb(k_buf+i))) 826 else 827 errmax = max(errmax,abs(dbl_mb(i_eri+i))) 828 end if 829 end if 830 end do 831 difmax = max(difmax,dif) 832 if (dif .gt. 0.1d0) then 833 write (LuOut,*) ish,Li,ipr,igc,ice 834 write (LuOut,*) jsh,Lj,jpr,jgc,jce 835 write (LuOut,*) ksh,Lk,kpr,kgc,kce 836 write (LuOut,*) lsh,Ll,lpr,lgc,lce 837 write (LuOut,*) dif,erimax 838 call util_flush(LuOut) 839 end if 840C end do 841C end do 842C end do 843C end do 844 write (luout,*) 'Maximum difference of all blocks',difmax 845 write (luout,*) 'Maximum difference from zero integrals',errmax 846 call int_terminate() 847 848 end 849************************************************************************ 850 subroutine kgdtest_ecpmem(rtdb) 851 implicit none 852#include "errquit.fh" 853#include "mafdecls.fh" 854#include "rtdb.fh" 855#include "context.fh" 856#include "geom.fh" 857#include "bas.fh" 858#include "nwc_const.fh" 859#include "basP.fh" 860#include "basdeclsP.fh" 861#include "geomP.fh" 862#include "geobasmapP.fh" 863#include "bas_exndcf_dec.fh" 864#include "bas_ibs_dec.fh" 865#include "ecp_nwc.fh" 866#include "stdio.fh" 867c 868 logical int_normalize 869 external int_normalize 870 integer idamax 871 external idamax 872c 873c test spin-orbit integrals 874c 875 integer rtdb 876 integer geom,basis, basis_id 877 integer nshell, memscr, membuf, maxscr 878 integer h_scr, k_scr, h_buf, k_buf 879 integer ish, jsh, ucont 880 integer li, i_prim, i_gen, i_iexp, i_icfp, i_cent, i_geom 881 integer lj, j_prim, j_gen, j_iexp, j_icfp, j_cent, j_geom 882 integer ihi,jhi,i2,j2 883 integer ibug,n_blk 884 character*255 mo_basis, geom_name 885c 886#include "bas_exndcf_sfn.fh" 887#include "bas_ibs_sfn.fh" 888c 889 mo_basis = 'ao basis' 890 geom_name = 'geometry' 891c 892 if(.not.geom_create(geom,geom_name))call errquit 893 & ('kgdtest_gencon: geom create error',911, GEOM_ERR) 894 if(.not.bas_create(basis,mo_basis))call errquit 895 & ('kgdtest_gencon: basis create error',911, BASIS_ERR) 896c 897 if(.not.geom_rtdb_load(rtdb,geom,geom_name)) call errquit 898 & ('kgdtest_gencon: geom load ',911, RTDB_ERR) 899 if(.not.bas_rtdb_load(rtdb,geom,basis,mo_basis)) call errquit 900 & ('kgdtest_gencon: basis load ',911, RTDB_ERR) 901c 902 basis_id = basis + BASIS_HANDLE_OFFSET 903 nshell = ncont_tot_gb(basis_id) 904 if (.not.int_normalize(rtdb,basis)) call errquit 905 & ('kgdtest_gencon: error normalizing ',911, INT_ERR) 906c 907 call int_init(rtdb,1,basis) 908 memscr = 1 000 000 909 membuf = 10 000 910 if (.not.ma_push_get(mt_dbl,memscr,' scratch ', 911 & h_scr, k_scr)) call errquit 912 & (' ma error 1',911, MA_ERR) 913 if (.not.ma_push_get(mt_dbl,membuf,' buf ', 914 & h_buf, k_buf)) call errquit 915 & (' ma error 2',911, MA_ERR) 916c 917 if (.not. rtdb_get(rtdb,'ecp:ibug',mt_int,1,ibug)) 918 & ibug = 1 919 if (.not. rtdb_get(rtdb,'ecp:n_blk',mt_int,1,n_blk)) 920 & n_blk = 1 921c 922 maxscr = 0 923 do ish = 1,nshell 924 ucont = (sf_ibs_cn2ucn(ish,basis_id)) 925 Li = infbs_cont(CONT_TYPE ,ucont,basis_id) 926 i_prim = infbs_cont(CONT_NPRIM,ucont,basis_id) 927 i_gen = infbs_cont(CONT_NGEN ,ucont,basis_id) 928 i_iexp = infbs_cont(CONT_IEXP ,ucont,basis_id) 929 i_icfp = infbs_cont(CONT_ICFP ,ucont,basis_id) 930 i_cent = (sf_ibs_cn2ce(ish,basis_id)) 931 i_geom = ibs_geom(basis_id) 932 i2 = (Li+1)*(Li+2)/2 933 ihi = i_gen*i2 934c 935 do jsh = 1,ish 936 ucont = (sf_ibs_cn2ucn(jsh,basis_id)) 937 Lj = infbs_cont(CONT_TYPE ,ucont,basis_id) 938 j_prim = infbs_cont(CONT_NPRIM,ucont,basis_id) 939 j_gen = infbs_cont(CONT_NGEN ,ucont,basis_id) 940 j_iexp = infbs_cont(CONT_IEXP ,ucont,basis_id) 941 j_icfp = infbs_cont(CONT_ICFP ,ucont,basis_id) 942 j_cent = (sf_ibs_cn2ce(jsh,basis_id)) 943 j_geom = ibs_geom(basis_id) 944 j2 = (Lj+1)*(Lj+2)/2 945 jhi = j_gen*j2 946 947 write(luout,*)' ' 948 write(luout,*)' ============= ',ish,jsh, 949 & '====================' 950 write(luout,*)' ' 951 952 call ecp_integral ( 953 & coords(1,i_cent,i_geom), 954 & dbl_mb(mb_exndcf(i_iexp,basis_id)), 955 & dbl_mb(mb_exndcf(i_icfp,basis_id)),i_prim,i_gen,Li,i_cent, 956 & coords(1,j_cent,j_geom), 957 & dbl_mb(mb_exndcf(j_iexp,basis_id)), 958 & dbl_mb(mb_exndcf(j_icfp,basis_id)),j_prim,j_gen,Lj,j_cent, 959 & dbl_mb(k_xyzecp), 960 & dbl_mb(k_ecp_e),dbl_mb(k_ecp_c), 961 & int_mb(k_ecp_nprim_c), 962 & int_mb(k_ecp_ncoef_c), ! new name is n_colc_C 963 & int_mb(k_ecp_ind_z), 964 & int_mb(k_ecp_ind_c), 965 & n_zeta_c, 966 & n_zeta_c, 967 & int_mb(k_ecp_l_c), 968 & int_mb(k_ecp_lip), 969 & n_ecp,l_ecp, 970 & 0, 971 & dbl_mb(k_ecp_c2s),mem_c2s, 972 & dbl_mb(k_buf),ihi*jhi,n_blk, 973 & .true.,dbl_mb(k_scr),memscr, 974 & ibug) 975 maxscr = max(memscr,maxscr) 976 977 end do 978 end do 979 write (luout,*) 'Maximum scratch needed:',maxscr 980 call int_terminate() 981 end 982************************************************************************ 983 subroutine kgdtest_relgen(rtdb) 984 implicit none 985#include "errquit.fh" 986#include "mafdecls.fh" 987#include "rtdb.fh" 988#include "context.fh" 989#include "geom.fh" 990#include "bas.fh" 991#include "nwc_const.fh" 992#include "basP.fh" 993#include "basdeclsP.fh" 994#include "geomP.fh" 995#include "geobasmapP.fh" 996#include "bas_exndcf_dec.fh" 997#include "bas_ibs_dec.fh" 998#include "apiP.fh" 999#include "rel_nwc.fh" 1000#include "stdio.fh" 1001c 1002 logical int_normalize 1003 external int_normalize 1004 integer idamax 1005 external idamax 1006c 1007c test general contraction in relativistic 2e integrals. 1008c 1009 integer rtdb 1010 integer geom, basis, basid 1011 integer nshell, memscr, membuf 1012 integer h_scr, k_scr, h_buf, k_buf 1013 integer ish, jsh, ksh,lsh, ucont 1014 integer li, ipr, igc, iex, icf, icfs, ice, ige 1015 integer lj, jpr, jgc, jex, jcf, jcfs, jce, jge 1016 integer lk, kpr, kgc, kex, kcf, kcfs, kce, kge 1017 integer ll, lpr, lgc, lex, lcf, lcfs, lce, lge 1018 integer Nints 1019 integer i_eri,i_c,j_c,k_c,l_c,ic,jc,kc,lc,i_seg 1020 integer ihi,jhi,khi,lhi,i2,j2,k2,l2,n_cart,n_cont 1021 integer abas,sbas,i_cs,j_cs,k_cs,l_cs 1022 character*255 basis_name, geom_name 1023 double precision difmax 1024c 1025#include "bas_exndcf_sfn.fh" 1026#include "bas_ibs_sfn.fh" 1027c 1028 write (luout,*) 'in kgdtest_relgen' 1029 basis_name = 'ao basis' 1030 geom_name = 'geometry' 1031c 1032 if(.not.geom_create(geom,geom_name))call errquit 1033 & ('kgdtest_relgen: geom create error',911, GEOM_ERR) 1034 if(.not.bas_create(basis,basis_name))call errquit 1035 & ('kgdtest_relgen: basis create error',911, BASIS_ERR) 1036c 1037 if(.not.geom_rtdb_load(rtdb,geom,geom_name)) call errquit 1038 & ('kgdtest_relgen: geom load ',911, GEOM_ERR) 1039 if(.not.bas_rtdb_load(rtdb,geom,basis,basis_name)) call errquit 1040 & ('kgdtest_relgen: basis load ',911, BASIS_ERR) 1041c 1042 basid = basis + BASIS_HANDLE_OFFSET 1043 nshell = ncont_tot_gb(basid) 1044 if (.not.int_normalize(rtdb,basis)) call errquit 1045 & ('kgdtest_relgen: error normalizing ',911, INT_ERR) 1046c 1047 call int_init(rtdb,1,basis) 1048 memscr = 5 000 000 1049 membuf = 100 000 1050 if (.not.ma_push_get(mt_dbl,memscr,' scratch ', 1051 & h_scr, k_scr)) call errquit 1052 & (' ma error 1',911, MA_ERR) 1053 if (.not.ma_push_get(mt_dbl,membuf,' buf ', 1054 & h_buf, k_buf)) call errquit 1055 & (' ma error 2',911, MA_ERR) 1056c 1057 sbas = sc_bsh + BASIS_HANDLE_OFFSET 1058 abas = ao_bsh + BASIS_HANDLE_OFFSET 1059 write (luout,*) ' starting shell loops ...' 1060c 1061 difmax = 0.0d0 1062 do ish = 1,nshell 1063 ucont = (sf_ibs_cn2ucn(ish,abas)) 1064 Li = infbs_cont(CONT_TYPE ,ucont,abas) 1065 ipr = infbs_cont(CONT_NPRIM,ucont,abas) 1066 igc = infbs_cont(CONT_NGEN ,ucont,abas) 1067 iex = infbs_cont(CONT_IEXP ,ucont,abas) 1068 icf = infbs_cont(CONT_ICFP ,ucont,abas) 1069 ice = (sf_ibs_cn2ce(ish,abas)) 1070 ige = ibs_geom(abas) 1071 i2 = (Li+1)*(Li+2)/2 1072 ihi = igc*i2 1073 ucont = ao_to_ls(ucont) 1074 icfs = infbs_cont(CONT_ICFP ,ucont,sbas) 1075c 1076 do jsh = 1,ish 1077 ucont = (sf_ibs_cn2ucn(jsh,abas)) 1078 Lj = infbs_cont(CONT_TYPE ,ucont,abas) 1079 jpr = infbs_cont(CONT_NPRIM,ucont,abas) 1080 jgc = infbs_cont(CONT_NGEN ,ucont,abas) 1081 jex = infbs_cont(CONT_IEXP ,ucont,abas) 1082 jcf = infbs_cont(CONT_ICFP ,ucont,abas) 1083 jce = (sf_ibs_cn2ce(jsh,abas)) 1084 jge = ibs_geom(abas) 1085 j2 = (Lj+1)*(Lj+2)/2 1086 jhi = jgc*j2 1087 ucont = ao_to_ls(ucont) 1088 jcfs = infbs_cont(CONT_ICFP ,ucont,sbas) 1089c 1090 do ksh = 1,ish 1091 ucont = (sf_ibs_cn2ucn(ksh,abas)) 1092 Lk = infbs_cont(CONT_TYPE ,ucont,abas) 1093 kpr = infbs_cont(CONT_NPRIM,ucont,abas) 1094 kgc = infbs_cont(CONT_NGEN ,ucont,abas) 1095 kex = infbs_cont(CONT_IEXP ,ucont,abas) 1096 kcf = infbs_cont(CONT_ICFP ,ucont,abas) 1097 kce = (sf_ibs_cn2ce(ish,abas)) 1098 kge = ibs_geom(abas) 1099 k2 = (Lk+1)*(Lk+2)/2 1100 khi = kgc*k2 1101 ucont = ao_to_ls(ucont) 1102 kcfs = infbs_cont(CONT_ICFP ,ucont,sbas) 1103c 1104 do lsh = 1,ksh 1105 ucont = (sf_ibs_cn2ucn(lsh,abas)) 1106 Ll = infbs_cont(CONT_TYPE ,ucont,abas) 1107 lpr = infbs_cont(CONT_NPRIM,ucont,abas) 1108 lgc = infbs_cont(CONT_NGEN ,ucont,abas) 1109 lex = infbs_cont(CONT_IEXP ,ucont,abas) 1110 lcf = infbs_cont(CONT_ICFP ,ucont,abas) 1111 lce = (sf_ibs_cn2ce(lsh,abas)) 1112 lge = ibs_geom(abas) 1113 l2 = (Ll+1)*(Ll+2)/2 1114 lhi = lgc*l2 1115 ucont = ao_to_ls(ucont) 1116 lcfs = infbs_cont(CONT_ICFP ,ucont,sbas) 1117 1118 n_cart = i2*j2*k2*l2 1119 n_cont = igc*jgc*kgc*lgc 1120 Nints = ihi*jhi*khi*lhi 1121C write (LuOut,*) ish,Li,ipr,igc,ice,ihi 1122C write (LuOut,*) jsh,Lj,jpr,jgc,jce,jhi 1123C write (LuOut,*) ksh,Lk,kpr,kgc,kce,khi 1124C write (LuOut,*) lsh,Ll,lpr,lgc,lce,lhi 1125C if (n_cont .gt. 1 .and. Li.ne.Lj .and. Lk.ne.Ll) then 1126C write(luout,*)' ' 1127C write(luout,*)' ============= (', 1128C & ish,':',Li,',', 1129C & jsh,':',Lj,'|', 1130C & ksh,':',Lk,',', 1131C & lsh,':',Ll,')', 1132C & '====================' 1133C write(luout,*)' ' 1134C write (luout,*) 'Nints',Nints 1135C write (luout,*) 'calling rel_2e4c_sf gen' 1136C call util_flush(6) 1137 1138 call rel_2e4c_sf ( 1139 & coords(1,ice,ige),dbl_mb(mb_exndcf(iex,abas)), 1140 & dbl_mb(mb_exndcf(icf,abas)), 1141 & dbl_mb(mb_exndcf(icfs,sbas)),ipr,igc,Li,ice, 1142 & coords(1,jce,jge),dbl_mb(mb_exndcf(jex,abas)), 1143 & dbl_mb(mb_exndcf(jcf,abas)), 1144 & dbl_mb(mb_exndcf(jcfs,sbas)),jpr,jgc,Lj,jce, 1145 & coords(1,kce,kge),dbl_mb(mb_exndcf(kex,abas)), 1146 & dbl_mb(mb_exndcf(kcf,abas)), 1147 & dbl_mb(mb_exndcf(kcfs,sbas)),kpr,kgc,Lk,kce, 1148 & coords(1,lce,lge),dbl_mb(mb_exndcf(lex,abas)), 1149 & dbl_mb(mb_exndcf(lcf,abas)), 1150 & dbl_mb(mb_exndcf(lcfs,sbas)),lpr,lgc,Ll,lce, 1151 & dbl_mb(k_buf),Nints,.false.,.false.,.false., 1152 & .false.,dbl_mb(k_scr),memscr, 1153 & .true.,.true.,.false.,.true.,0) 1154C call ecp_matpr (dbl_mb(k_buf),1,khi*lhi,1,ihi*jhi, 1155C & 1,khi*lhi,1,ihi*jhi,'General contraction', 1156C & 'E',120,6) 1157C write (luout,*) 'calling rel_2e4c_sf seg' 1158C call util_flush(6) 1159 1160 i_eri = k_buf+Nints 1161 i_seg = i_eri 1162 i_c = icf 1163 i_cs = icfs 1164 do ic = 1,igc 1165 j_c = jcf 1166 j_cs = jcfs 1167 do jc = 1,jgc 1168 k_c = kcf 1169 k_cs = kcfs 1170 do kc = 1,kgc 1171 l_c = lcf 1172 l_cs = lcfs 1173 do lc = 1,lgc 1174C write (luout,*) 'ic,jc,kc,lc',ic,jc,kc,lc 1175 1176 call rel_2e4c_sf ( 1177 & coords(1,ice,ige),dbl_mb(mb_exndcf(iex,abas)), 1178 & dbl_mb(mb_exndcf(i_c,abas)), 1179 & dbl_mb(mb_exndcf(i_cs,sbas)),ipr,1,Li,ice, 1180 & coords(1,jce,jge),dbl_mb(mb_exndcf(jex,abas)), 1181 & dbl_mb(mb_exndcf(j_c,abas)), 1182 & dbl_mb(mb_exndcf(j_cs,sbas)),jpr,1,Lj,jce, 1183 & coords(1,kce,kge),dbl_mb(mb_exndcf(kex,abas)), 1184 & dbl_mb(mb_exndcf(k_c,abas)), 1185 & dbl_mb(mb_exndcf(k_cs,sbas)),kpr,1,Lk,kce, 1186 & coords(1,lce,lge),dbl_mb(mb_exndcf(lex,abas)), 1187 & dbl_mb(mb_exndcf(l_c,abas)), 1188 & dbl_mb(mb_exndcf(l_cs,sbas)),lpr,1,Ll,lce, 1189 & dbl_mb(i_eri),Nints,.false.,.false.,.false., 1190 & .false.,dbl_mb(k_scr),memscr, 1191 & .true.,.true.,.false.,.true.,0) 1192C call ecp_matpr (dbl_mb(k_buf),1,k2*l2,1,i2*j2, 1193C & 1,k2*l2,1,i2*j2,'Segmented contraction', 1194C & 'E',120,6) 1195 1196 i_eri = i_eri+i2*j2*k2*l2 1197 l_c = l_c+lpr 1198 l_cs = l_cs+lpr 1199 end do 1200 k_c = k_c+kpr 1201 k_cs = k_cs+kpr 1202 end do 1203 j_c = j_c+jpr 1204 j_cs = j_cs+jpr 1205 end do 1206 i_c = i_c+ipr 1207 i_cs = i_cs+ipr 1208 end do 1209 1210 call reorder_2eints(dbl_mb(i_eri),dbl_mb(i_seg), 1211 & l2,lgc,k2,kgc,j2,jgc,i2,igc) 1212 call dcopy (Nints,dbl_mb(i_eri),1,dbl_mb(i_seg),1) 1213C call ecp_matpr (dbl_mb(i_seg),1,khi*lhi,1,ihi*jhi, 1214C & 1,khi*lhi,1,ihi*jhi,'Segmented contraction', 1215C & 'E',120,6) 1216 call daxpy(Nints,-1.0d0,dbl_mb(k_buf),1,dbl_mb(i_eri),1) 1217C call ecp_matpr (dbl_mb(i_eri),1,khi*lhi,1,ihi*jhi, 1218C & 1,khi*lhi,1,ihi*jhi,'Differences', 1219C & 'E',120,6) 1220 ic = idamax(Nints,dbl_mb(i_eri),1)-1 1221C write (luout,*) 'Maximum difference',dbl_mb(i_eri+ic) 1222 difmax = max(difmax,abs(dbl_mb(i_eri+ic))) 1223 if (abs(dbl_mb(i_eri+ic)) .gt. 1.0d-12) then 1224 write (luout,*) ish,jsh,ksh,lsh 1225 call ecp_matpr (dbl_mb(k_buf),1,khi*lhi,1,ihi*jhi, 1226 & 1,khi*lhi,1,ihi*jhi,'General contraction', 1227 & 'E',120,6) 1228 call ecp_matpr (dbl_mb(i_seg),1,khi*lhi,1,ihi*jhi, 1229 & 1,khi*lhi,1,ihi*jhi,'Segmented contraction', 1230 & 'E',120,6) 1231 end if 1232C end if 1233 end do 1234 end do 1235 end do 1236 end do 1237 write (luout,*) 'Maximum difference of all blocks',difmax 1238 1239 end 1240************************************************************************ 1241 subroutine kgdtest_gen1d(rtdb) 1242 implicit none 1243#include "errquit.fh" 1244#include "mafdecls.fh" 1245#include "rtdb.fh" 1246#include "context.fh" 1247#include "geom.fh" 1248#include "bas.fh" 1249#include "nwc_const.fh" 1250#include "basP.fh" 1251#include "basdeclsP.fh" 1252#include "geomP.fh" 1253#include "geobasmapP.fh" 1254#include "bas_exndcf_dec.fh" 1255#include "bas_ibs_dec.fh" 1256#include "apiP.fh" 1257#include "rel_nwc.fh" 1258#include "stdio.fh" 1259c 1260 logical int_normalize 1261 external int_normalize 1262 integer idamax 1263 external idamax 1264c 1265c test general contraction in derivative 1e integrals 1266c 1267 integer rtdb 1268 integer geom, basis, basid 1269 integer nshell, memscr, membuf 1270 integer h_scr, k_scr, h_buf, k_buf 1271 integer ish, jsh, ucont 1272 integer li, ipr, igc, iex, icf, ice, ige 1273 integer lj, jpr, jgc, jex, jcf, jce, jge 1274 integer Nints 1275 integer i_o2i,i_kei,i_nai 1276 integer ihi,jhi,i2,j2,n_cart,n_cont 1277 character*255 basis_name, geom_name 1278 double precision difmax 1279c 1280#include "bas_exndcf_sfn.fh" 1281#include "bas_ibs_sfn.fh" 1282c 1283 write (luout,*) 'in kgdtest_relgen' 1284 basis_name = 'ao basis' 1285 geom_name = 'geometry' 1286c 1287 if(.not.geom_create(geom,geom_name))call errquit 1288 & ('kgdtest_relgen: geom create error',911, GEOM_ERR) 1289 if(.not.bas_create(basis,basis_name))call errquit 1290 & ('kgdtest_relgen: basis create error',911, BASIS_ERR) 1291c 1292 if(.not.geom_rtdb_load(rtdb,geom,geom_name)) call errquit 1293 & ('kgdtest_relgen: geom load ',911, RTDB_ERR) 1294 if(.not.bas_rtdb_load(rtdb,geom,basis,basis_name)) call errquit 1295 & ('kgdtest_relgen: basis load ',911, RTDB_ERR) 1296c 1297 basid = basis + BASIS_HANDLE_OFFSET 1298 nshell = ncont_tot_gb(basid) 1299 if (.not.int_normalize(rtdb,basis)) call errquit 1300 & ('kgdtest_relgen: error normalizing ',911, INT_ERR) 1301c 1302 call int_init(rtdb,1,basis) 1303 memscr = 5 000 000 1304 membuf = 100 000 1305 if (.not.ma_push_get(mt_dbl,memscr,' scratch ', 1306 & h_scr, k_scr)) call errquit 1307 & (' ma error 1',911, MA_ERR) 1308 if (.not.ma_push_get(mt_dbl,membuf,' buf ', 1309 & h_buf, k_buf)) call errquit 1310 & (' ma error 2',911, MA_ERR) 1311c 1312 difmax = 0.0d0 1313 do ish = 1,nshell 1314 ucont = (sf_ibs_cn2ucn(ish,basid)) 1315 Li = infbs_cont(CONT_TYPE ,ucont,basid) 1316 ipr = infbs_cont(CONT_NPRIM,ucont,basid) 1317 igc = infbs_cont(CONT_NGEN ,ucont,basid) 1318 iex = infbs_cont(CONT_IEXP ,ucont,basid) 1319 icf = infbs_cont(CONT_ICFP ,ucont,basid) 1320 ice = (sf_ibs_cn2ce(ish,basid)) 1321 ige = ibs_geom(basid) 1322 i2 = (Li+1)*(Li+2)/2 1323 ihi = igc*i2 1324c 1325 do jsh = 1,ish 1326 ucont = (sf_ibs_cn2ucn(jsh,basid)) 1327 Lj = infbs_cont(CONT_TYPE ,ucont,basid) 1328 jpr = infbs_cont(CONT_NPRIM,ucont,basid) 1329 jgc = infbs_cont(CONT_NGEN ,ucont,basid) 1330 jex = infbs_cont(CONT_IEXP ,ucont,basid) 1331 jcf = infbs_cont(CONT_ICFP ,ucont,basid) 1332 jce = (sf_ibs_cn2ce(jsh,basid)) 1333 jge = ibs_geom(basid) 1334 j2 = (Lj+1)*(Lj+2)/2 1335 jhi = jgc*j2 1336 1337 n_cart = i2*j2 1338 n_cont = igc*jgc 1339 Nints = ihi*jhi 1340 write (LuOut,*) ish,Li,ipr,igc,ice,ihi 1341 write (LuOut,*) jsh,Lj,jpr,jgc,jce,jhi 1342 write (luout,*) 'Nints',Nints 1343 i_o2i = k_scr 1344 i_kei = i_o2i+Nints*6 1345 i_nai = i_kei+Nints*6 1346 call util_flush(6) 1347c 1348 call hf1d( 1349 & coords(1,ice,ige),dbl_mb(mb_exndcf(iex,basid)), 1350 & dbl_mb(mb_exndcf(icf,basid)),ipr,igc,Li,ice, 1351 & coords(1,jce,jge),dbl_mb(mb_exndcf(jex,basid)), 1352 & dbl_mb(mb_exndcf(jcf,basid)),jpr,jgc,Lj,jce, 1353 & coords(1,1,ige),charge(1,ige), 1354 & geom_invnucexp(1,ige),ncenter(ige), 1355 & dbl_mb(i_o2i),dbl_mb(i_kei),dbl_mb(i_nai),Nints, 1356 & .true.,.true.,.true.,.false.,.false., 1357 & dbl_mb(k_scr),memscr) 1358c 1359C call hf1d_new( 1360C & coords(1,ice,ige),dbl_mb(mb_exndcf(iex,basid)), 1361C & dbl_mb(mb_exndcf(icf,basid)),ipr,igc,Li,ice, 1362C & coords(1,jce,jge),dbl_mb(mb_exndcf(jex,basid)), 1363C & dbl_mb(mb_exndcf(jcf,basid)),jpr,jgc,Lj,jce, 1364C & coords(1,1,ige),charge(1,ige), 1365C & geom_invnucexp(1,ige),ncenter(ige), 1366C & dbl_mb(i_o2i),dbl_mb(i_kei),dbl_mb(i_nai),Nints, 1367C & .true.,.true.,.true.,.false.,.false., 1368C & dbl_mb(k_scr),memscr) 1369C 1370C write (6,*) k_buf,i_seg,i_eri 1371C call reorder_2eints(dbl_mb(i_eri),dbl_mb(i_seg), 1372C & l2,lgc,k2,kgc,j2,jgc,i2,igc) 1373C call dcopy (Nints,dbl_mb(i_eri),1,dbl_mb(i_seg),1) 1374C call ecp_matpr (dbl_mb(i_seg),1,khi*lhi,1,ihi*jhi, 1375C & 1,khi*lhi,1,ihi*jhi,'Segmented contraction', 1376C & 'E',120,6) 1377C call daxpy(Nints,-1.0d0,dbl_mb(k_buf),1,dbl_mb(i_eri),1) 1378C call ecp_matpr (dbl_mb(i_eri),1,khi*lhi,1,ihi*jhi, 1379C & 1,khi*lhi,1,ihi*jhi,'Differences', 1380C & 'E',120,6) 1381C ic = idamax(Nints,dbl_mb(i_eri),1)-1 1382C write (luout,*) 'Maximum difference',dbl_mb(i_eri+ic) 1383C difmax = max(difmax,abs(dbl_mb(i_eri+ic))) 1384C if (abs(dbl_mb(i_eri+ic)) .gt. 1.0d-12) then 1385C call ecp_matpr (dbl_mb(k_buf),1,khi*lhi,1,ihi*jhi, 1386C & 1,khi*lhi,1,ihi*jhi,'General contraction', 1387C & 'E',120,6) 1388C call ecp_matpr (dbl_mb(i_seg),1,khi*lhi,1,ihi*jhi, 1389C & 1,khi*lhi,1,ihi*jhi,'Segmented contraction', 1390C & 'E',120,6) 1391C end if 1392C end do 1393C end do 1394 end do 1395 end do 1396 write (luout,*) 'Maximum difference of all blocks',difmax 1397 1398 end 1399