1 logical function jvltest2(rtdb) 2* 3* $Id$ 4* 5 implicit none 6#include "errquit.fh" 7#include "global.fh" 8#include "mafdecls.fh" 9#include "bas.fh" 10#include "geom.fh" 11#include "rtdb.fh" 12#include "inp.fh" 13#include "util.fh" 14 integer rtdb 15c 16 integer basis, geom, nbf 17 character*255 movecs ! Name of movector file 18 character*80 title, name_of_basis, scftype 19 integer nbf_file, nsets, nmo_file(2) 20 logical cct_uhf 21 logical movecs_read, movecs_read_header 22 external movecs_read, movecs_read_header 23c 24 integer mult(8,8) 25c 26 integer nmo, nalpha, nbeta, g_tmp, l_occ, k_occ, 27 $ l_aeval, k_aeval, l_beval, k_beval, l_mos_t, k_mos_t, 28 $ l_eval, k_eval, l_tmp, k_tmp, nocca, noccb, 29 $ nvirta, nvirtb, 30 $ l_irs, k_irs 31c 32 integer nfrozen, nactive 33 logical status 34c 35 logical int_normalize 36 external int_normalize 37c 38 call ifill(8*8, 99999999, mult, 1) 39 mult(1,1) = 1 ! All we are using now! 40c 41 write(6,*) 42 write(6,*) 43 call util_print_centered(6,'Joop''s cool triples program', 44 $ 40, .true.) 45 write(6,*) 46 write(6,*) 47 title = ' ' 48 status = rtdb_cget(rtdb,'title',1,title) 49 if (title .ne. ' ') then 50 call util_print_centered(6,title,40, .false.) 51 write(6,*) 52 write(6,*) 53 endif 54c 55c load the geometry/basis set and get info 56c 57 if (.not. geom_create(geom, 'geometry')) 58 $ call errquit('scf_init: geom_create?', 0, GEOM_ERR) 59 if (.not. geom_rtdb_load(rtdb, geom, 'geometry')) 60 $ call errquit('scf_init: no geometry ', 0, RTDB_ERR) 61 if (.not. bas_create(basis, 'ao basis')) 62 $ call errquit('scf_init: bas_create?', 0, BASIS_ERR) 63 if (.not. bas_rtdb_load(rtdb, geom, basis, 'ao basis')) 64 $ call errquit('scf_init: no ao basis set', 0, RTDB_ERR) 65 if (.not. bas_numbf(basis, nbf)) call errquit 66 $ ('scf_init: basis info',0, BASIS_ERR) 67 if (.not.int_normalize(rtdb,basis)) 68 $ call errquit('scf:int_normalize failed', 0, INT_ERR) 69c 70c Read the MO vectors and evals from a UHF calculation 71c 72 call util_file_name('movecs',.false.,.false.,movecs) 73 if (.not. movecs_read_header(movecs, title, name_of_basis, 74 $ scftype, nbf_file, nsets, nmo_file, 2)) call errquit 75 $ ('jvltest2: failed to read movecs header',0, DISK_ERR) 76c 77 if (nsets.ne.2) call errquit('Joop says: UHF ONLY!',0, INPUT_ERR) 78c 79 write(6,*) ' Read movecs header from ', 80 $ movecs(1:inp_strlen(movecs)) 81 write(6,*) ' Job title : ', 82 $ title(1:inp_strlen(title)) 83 write(6,*) ' Basis name: ', 84 $ name_of_basis(1:inp_strlen(name_of_basis)) 85 nmo = nmo_file(1) 86 if (nmo_file(1).ne.nmo_file(2)) call errquit('Uh?',0, INPUT_ERR) 87 if (.not. rtdb_get(rtdb, 'scf:nalpha', mt_int, 1, nalpha)) 88 $ call errquit('nalpha?',0, RTDB_ERR) 89 if (.not. rtdb_get(rtdb, 'scf:nbeta', mt_int, 1, nbeta)) 90 $ call errquit('nbeta?',0, RTDB_ERR) 91c 92c Frozen core info 93c 94 if (rtdb_get(rtdb,'ccsd:frozen core:freeze by atoms',mt_log, 1, 95 $ status)) then 96 if (.not. geom_num_core(geom, nfrozen)) 97 $ call errquit('Joop says: geom_num_core?',0, GEOM_ERR) 98 else if (rtdb_get(rtdb, 'ccsd:frozen core', MT_INT, 1, 99 $ nfrozen)) then 100 else 101 nfrozen = 0 102 endif 103c 104 nactive = nmo - nfrozen 105 nocca = nalpha - nfrozen 106 noccb = nbeta - nfrozen 107 nvirta = nactive - nocca 108 nvirtb = nactive - noccb 109c 110 write(6,*) ' No. of alpha electrons: ', nalpha 111 write(6,*) ' No. of beta electrons: ', nbeta 112 write(6,*) ' No. of frozen core: ', nfrozen 113 write(6,*) ' No. of molecular orbitals: ', nmo 114 write(6,*) ' No. of active MOs: ', nactive 115 write(6,*) ' No. of basis functions: ', nbf 116c 117c Generate local copies of the TRANSPOSED, ACTIVE molecular orbitals. 118c 119c The active MOs and eigenvalues will be combined to form single 120c arrays with the first nactive entries being alpha and the next beta. 121c 122 if (.not. ma_push_get(mt_dbl, nbf*nactive*2,'mos_t', 123 $ l_mos_t, k_mos_t)) call errquit('Joop: mos', 2*nbf*nactive, 124 & MA_ERR) 125 if (.not. ma_push_get(mt_dbl, nactive*2,'evals', 126 $ l_eval, k_eval)) call errquit('Joop: evals', 2*nactive, 127 & MA_ERR) 128 if (.not. ma_push_get(mt_int, nactive*2, 'irs', 129 $ l_irs, k_irs)) call errquit('Joop: irs', 2*nactive, 130 & MA_ERR) 131 call ifill(2*nactive, 1, int_mb(k_irs), 1) 132c 133 if (.not. ma_push_get(mt_dbl, nbf,'occ',l_occ, k_occ)) 134 $ call errquit('ma occ', nbf, MA_ERR) 135 if (.not. ma_push_get(mt_dbl, nbf,'aeval',l_aeval, k_aeval)) 136 $ call errquit('ma eval', nbf, MA_ERR) 137 if (.not. ma_push_get(mt_dbl, nbf,'beval',l_beval, k_beval)) 138 $ call errquit('ma eval', nbf, MA_ERR) 139 if (.not. ma_push_get(mt_dbl, nbf*nactive,'tmp', 140 $ l_tmp, k_tmp)) call errquit('Joop: tmp', nbf*nactive, MA_ERR) 141c 142*ga:1:0 143 if (.not. ga_create(mt_dbl, nbf, nmo, 'tmp', 0, 0, g_tmp)) 144 & call errquit('scf_v_g: tmp', 0, GA_ERR) 145c 146 if (.not. movecs_read(movecs, 1, dbl_mb(k_occ), dbl_mb(k_aeval), 147 $ g_tmp)) call errquit('movecs_read of amos failed ',0, 148 & DISK_ERR) 149 call ga_get(g_tmp, 1, nbf, nfrozen+1, nmo, dbl_mb(k_tmp), nbf) 150 call util_transpose(dbl_mb(k_tmp),nbf,dbl_mb(k_mos_t),nactive, 151 $ nbf,nactive) 152 call dcopy(nactive, dbl_mb(k_aeval+nfrozen), 1, dbl_mb(k_eval), 1) 153c 154 if (.not. movecs_read(movecs, 2, dbl_mb(k_occ), dbl_mb(k_beval), 155 $ g_tmp)) call errquit('movecs_read of amos failed ',0, 156 & DISK_ERR) 157 call ga_get(g_tmp, 1, nbf, 1, nmo, dbl_mb(k_tmp), nbf) 158 call util_transpose(dbl_mb(k_tmp), nbf, 159 $ dbl_mb(k_mos_t+nactive*nbf), nactive, nbf, nactive) 160 call dcopy(nactive, dbl_mb(k_beval+nfrozen), 1, 161 $ dbl_mb(k_eval+nactive), 1) 162c 163 if (.not. ga_destroy(g_tmp)) call errquit(' ga bad?',0, GA_ERR) 164 if (.not. ma_chop_stack(l_occ)) call errquit(' ma chop?', 0, 165 & MA_ERR) 166c 167 write(6,*) ' Alpha active eigenvalues ' 168 call output(dbl_mb(k_eval),1,nactive,1,1,nactive,1,1) 169 write(6,*) ' Beta active eigenvalues ' 170 call output(dbl_mb(k_eval+nactive),1,nactive,1,1,nactive,1,1) 171 write(6,*) ' Alpha MOs transposed' 172 call output(dbl_mb(k_mos_t),1,nactive,1,nbf,nactive,nbf,1) 173 write(6,*) ' Beta active MOs transposed' 174 call output(dbl_mb(k_mos_t+nbf*nactive),1,nactive,1,nbf, 175 $ nactive,nbf,1) 176c 177c Now make the darn integral files 178c 179 call rjh_init_file_info() 180 call rjh_make_integral_files( 181 $ rtdb, basis, 182 $ nactive, nbf, nocca, noccb, 183 $ dbl_mb(k_mos_t), dbl_mb(k_eval)) 184c 185c Wake up Joop 186c 187 call cct_uhf_address(1, mult, int_mb(k_irs), 188 1 nocca, noccb, nactive) 189c 190c Trying reading something from a file ... all 3-c integrals 191c 192 if (.not. ga_create(mt_dbl, nvirta*nocca, nvirtb*nvirtb, 193 $ 'tmp', 0, 0, g_tmp)) 194 & call errquit('scf_v_g: tmp', 0, GA_ERR) 195c 196 call get_cct(g_tmp, '<ba||ek>(bkea)', 197 $ 1, nvirta*nocca, 1, 1, 198 $ nactive+noccb+1, nactive+nactive, 2, 199 $ nactive+noccb+1, nactive+nactive, 2) 200c 201 write(6,*) ' spins b k a e = ', 1, 1, 2, 2 202 call ga_print(g_tmp) 203c 204 call get_cct(g_tmp, '<ba||ek>(bkea)', 205 $ 1, nvirta*nocca, 2, 2, 206 $ noccb+1, nactive, 1, 207 $ noccb+1, nactive, 1) 208c 209 write(6,*) ' spins b k a e = ', 2, 2, 1, 1 210 call ga_print(g_tmp) 211c 212 if (.not. ga_destroy(g_tmp)) call errquit('ga?',0, GA_ERR) 213c 214 status = cct_uhf(rtdb) 215c 216c Tidy up 217c 218 if (.not. bas_destroy(basis)) call errquit(' bas ?',0, BASIS_ERR) 219 if (.not. geom_destroy(geom)) call errquit(' geom ?',0, GEOM_ERR) 220 if (.not. ma_chop_stack(l_mos_t)) call errquit(' ma ? ',0, MA_ERR) 221c 222 call rjh_tidy_file_info 223c 224 jvltest2 = status 225c 226 end 227c$$$ subroutine rjh_test_trans(rtdb, basis, amos, bmos, 228c$$$ $ nmo, nalpha, nbeta, nbf) 229c$$$ implicit none 230c$$$#include "errquit.fh" 231c$$$#include "mafdecls.fh" 232c$$$c 233c$$$ integer rtdb, basis 234c$$$ integer nmo, nalpha, nbeta, nbf 235c$$$ double precision amos(nbf,nmo), bmos(nbf,nmo), sum 236c$$$ 237c$$$c 238c$$$ integer lenfull, k_full, l_full, k_full2, l_full2 239c$$$ integer l_aoint, k_aoint, i 240c$$$c 241c$$$ if (nmo .ne. nbf) call errquit(' test hack ',0) 242c$$$c 243c$$$ lenfull = nbf**4 244c$$$ if (.not. ma_push_get(mt_dbl,lenfull,'full',l_full, k_full)) 245c$$$ $ call errquit('ma full', lenfull) 246c$$$ if (.not. ma_push_get(mt_dbl,lenfull,'full2',l_full2, k_full2)) 247c$$$ $ call errquit('ma full2', lenfull) 248c$$$ if (.not. ma_push_get(mt_dbl,lenfull,'aoint',l_aoint, k_aoint)) 249c$$$ $ call errquit('ma aoint', lenfull) 250c$$$c 251c$$$ call util_inplace_transpose(nbf, amos) ! TRANSPOSE 252c$$$ call util_inplace_transpose(nbf, bmos) 253c$$$c 254c$$$ call rjh_full_transform(rtdb, basis, nbf, nbf, nbf, nbf, 255c$$$ $ nmo, nmo, nmo, nmo, 256c$$$ $ amos, amos, amos, amos, dbl_mb(k_full), 'chargecloud') 257c$$$c 258c$$$ call rjh_make_integral_files(rtdb, basis, 259c$$$ $ nmo, nalpha, nbeta, amos, bmos) 260c$$$ call rjh_tidy_file_info() 261c$$$c 262c$$$ call util_inplace_transpose(nbf, amos) !UNDO TRANSPOSE 263c$$$ call util_inplace_transpose(nbf, bmos) 264c$$$c 265c$$$c$$$ call rjh_debug_print('full', 266c$$$c$$$ $ dbl_mb(k_full), nbf, nbf, nbf, nbf) 267c$$$c 268c$$$ call rjh_all_ao_integrals(rtdb,basis,nbf,'mulliken', 269c$$$ $ dbl_mb(k_aoint)) 270c$$$c 271c$$$ call rjh_uhf_trans_lo_mem(nbf, nbf, int_mb(k_full2), 272c$$$ $ amos, dbl_mb(k_aoint), dbl_mb(k_full2)) 273c$$$c$$$ call rjh_debug_print('full2', 274c$$$c$$$ $ dbl_mb(k_full2), nbf, nbf, nbf, nbf) 275c$$$c 276c$$$ sum = 0.0d0 277c$$$ do i = 0, lenfull-1 278c$$$ sum = sum + abs(dbl_mb(k_full2+i)-dbl_mb(k_full2+i)) 279c$$$ enddo 280c$$$ write(6,*) ' Difference between old and new ', sum 281c$$$c 282c$$$ if (.not. ma_chop_stack(l_full)) call errquit('ma99',0) 283c$$$c 284c$$$ end 285 subroutine rjh_debug_print(string,full, n1, n2, n3, n4) 286 implicit none 287 character*(*) string 288 integer n1, n2, n3, n4 289 double precision full(n1, n2, n3, n4) 290c 291 integer p, q, r, s 292 write(6,*) ' DEBUG FOR ', string, n1, n2, n3, n4 293 do s = 1, n4 294 do r = 1, n3 295 do q = 1, n2 296 do p = 1, n1 297 if (abs(full(p,q,r,s)).gt.1e-6) then 298 write(6,1) p,q,r,s,full(p,q,r,s) 299 1 format(4i5,2x,f12.6) 300 end if 301 end do 302 end do 303 end do 304 end do 305c 306 end 307 subroutine rjh_init_file_info() 308 implicit none 309#include "rjhfileinfo.fh" 310c 311 integer i 312c 313 do i = 1, maxfile 314 fds(i) = -1 315 enddo 316c 317 end 318 subroutine rjh_tidy_file_info() 319 implicit none 320#include "errquit.fh" 321#include "rjhfileinfo.fh" 322#include "eaf.fh" 323c 324 integer i 325c 326 do i = 1, maxfile 327 if (fds(i).ne.-1) then 328 if (eaf_close(fds(i)) .ne. 0) call errquit 329 $ ('rjh_tidy_file_info: eaf_close failed',i, UNKNOWN_ERR) 330 endif 331 enddo 332c 333 end 334 subroutine rjh_make_integral_files(rtdb, basis, 335 $ nmo, nbf, nalpha, nbeta, mos_t, evals) 336 implicit none 337#include "errquit.fh" 338#include "mafdecls.fh" 339#include "rjhfileinfo.fh" 340 integer k_all, norb 341 common /c_all/norb, k_all 342c 343 integer rtdb, basis 344 integer nmo, nalpha, nbeta, nbf 345 double precision mos_t(nmo,nbf,2), evals(nmo,2) 346c 347c <ba||ek> = <ba|ek> - <ba|ke> 348c Write the next 6 as file(b,k,e,a) 349c 350c 1. <ba||ek> b=a=e=k=alpha baek = 1111 351c 2. <ba||ek> b=a=e=k=beta baek = 2222 352c 3. <ba|ek> b=e=alpha, a=k=beta baek = 1212 353c 4. <ba|ek> b=e=beta, a=k=alpha baek = 2121 354c 5. -<ba|ke> = (bk|ae) b=k=alpha, a=e=beta baek = 1221 sign=-1 355c 6. -<ba|ke> = (bk|ae) b=k=beta, a=e=alpha baek = 2112 sign=-1 356c 357c <ij||ab> 358c Write the next set as (i,j,a,b) restricting i>=j, a>=b if possible. 359c 360c 7. <ij||ab> i=j=a=b=alpha ijab = 1111 361c 8. <ij||ab> i=j=a=b=beta ijab = 2222 362c 9. <ij|ab> i=a=alpha j=b=beta ijab = 1212 363c 364c T(ij,ab) written same as (i,j,a,b) 365c Write the next set as (i,j,a,b) restricting i>=j, a>=b if possible. 366c 367c 10. <ij||ab> i=j=a=b=alpha ijab = 1111 368c 11. <ij||ab> i=j=a=b=beta ijab = 2222 369c 12. <ij|ab> i=a=alpha j=b=beta ijab = 1212 370c 371c No use of spatial symmetry at present. Assume can hold 372c everything in memory. 373c 374 integer lenfull, l_full, k_full 375 integer occa, occb ! No. of occupied of each spin 376 integer virta, virtb ! No. of virtuals of each spin 377c 378 call rjh_init_file_info() 379 occa = nalpha 380 occb = nbeta 381 virta = nmo-occa 382 virtb = nmo-occb 383c 384* lenfull = max(occa,occb)*max(virta,virtb)**3 385 lenfull = nmo**4 386c 387c Slowly move to have the following using Joop-like spin cases. 388c 389 if (.not. ma_push_get(mt_dbl,lenfull,'full',l_full, k_full)) 390 $ call errquit('ma full', lenfull, MA_ERR) 391 392* call rjh_full_transform(rtdb, basis, 393* $ nmo, nmo, nmo, nmo, 394* $ nmo, nmo, nmo, nmo, 395* $ mos_t(1,1,1), mos_t(1,1,1), 396* $ mos_t(1,1,1), mos_t(1,1,1), 397* $ dbl_mb(k_full), 'Chargecloud') 398* k_all = k_full 399* norb = nmo 400 401 if (.not. ma_push_get(mt_dbl,lenfull,'full',l_full, k_full)) 402 $ call errquit('ma full', lenfull, MA_ERR) 403 404c 405c 1. <ba||ek> b=a=e=k=alpha baek = 1111 406c 407 call rjh_full_transform(rtdb, basis, 408 $ virta, virta, virta, occa, 409 $ nmo, nmo, nmo, nmo, 410 $ mos_t(occa+1,1,1), mos_t(occa+1,1,1), 411 $ mos_t(occa+1,1,1), mos_t(1,1,1), 412 $ dbl_mb(k_full), 'LeftAsymDirac') 413 call rjh_integ_write(dbl_mb(k_full), virta, virta, virta, occa, 414 $ '1432', '<ba||ek>', 1111, 1.0d0) 415c 416c 2. <ba||ek> b=a=e=k=beta baek = 2222 417c 418 call rjh_full_transform(rtdb, basis, 419 $ virtb, virtb, virtb, occb, 420 $ nmo, nmo, nmo, nmo, 421 $ mos_t(occb+1,1,2), mos_t(occb+1,1,2), 422 $ mos_t(occb+1,1,2), mos_t(1,1,2), 423 $ dbl_mb(k_full), 'LeftAsymDirac') 424 call rjh_integ_write(dbl_mb(k_full), virtb, virtb, virtb, occb, 425 $ '1432', '<ba||ek>', 2222, 1.0d0) 426c 427c 3. <ba|ek> b=e=alpha, a=k=beta baek = 1212 428c 429 call rjh_full_transform(rtdb, basis, 430 $ virta, virtb, virta, occb, 431 $ nmo, nmo, nmo, nmo, 432 $ mos_t(occa+1,1,1), mos_t(occb+1,1,2), 433 $ mos_t(occa+1,1,1), mos_t(1,1,2), 434 $ dbl_mb(k_full), 'Dirac') 435 call rjh_integ_write(dbl_mb(k_full), virta, virtb, virta, occb, 436 $ '1432', '<ba||ek>', 1212, 1.0d0) 437c 438c 4. <ba|ek> b=e=beta, a=k=alpha baek = 2121 439c 440 call rjh_full_transform(rtdb, basis, 441 $ virtb, virta, virtb, occa, 442 $ nmo, nmo, nmo, nmo, 443 $ mos_t(occb+1,1,2), mos_t(occa+1,1,1), 444 $ mos_t(occb+1,1,2), mos_t(1,1,1), 445 $ dbl_mb(k_full), 'Dirac') 446 call rjh_integ_write(dbl_mb(k_full), virtb, virta, virtb, occa, 447 $ '1432', '<ba||ek>', 2121, 1.0d0) 448c 449c 5. <ba|ke> = (bk|ae) b=k=alpha, a=e=beta baek = 1221 sign=-1 450c 451 call rjh_full_transform(rtdb, basis, 452 $ virta, occa, virtb, virtb, 453 $ nmo, nmo, nmo, nmo, 454 $ mos_t(occa+1,1,1), mos_t(1,1,1), 455 $ mos_t(occb+1,1,2), mos_t(occb+1,1,2), 456 $ dbl_mb(k_full), 'Chargecloud') 457 call rjh_integ_write(dbl_mb(k_full), virta, occa, virtb, virtb, 458 $ '1234', '<ba||ek>', 1221, -1.0d0) 459c 460c 6. <ba|ke> = (bk|ae) b=k=beta, a=e=alpha baek = 2112 sign=-1 461c 462 call rjh_full_transform(rtdb, basis, 463 $ virtb, occb, virta, virta, 464 $ nmo, nmo, nmo, nmo, 465 $ mos_t(occb+1,1,2), mos_t(1,1,2), 466 $ mos_t(occa+1,1,1), mos_t(occa+1,1,1), 467 $ dbl_mb(k_full), 'Chargecloud') 468 call rjh_integ_write(dbl_mb(k_full), virtb, occb, virta, virta, 469 $ '1234', '<ba||ek>', 2112, -1.0d0) 470c 471c 7. <ij||ab> i=j=a=b=alpha ijab = 1111 472c i>=j, a>=b 473c 474 call rjh_full_transform(rtdb, basis, 475 $ occa, occa, virta, virta, 476 $ nmo, nmo, nmo, nmo, 477 $ mos_t(1,1,1), mos_t(1,1,1), 478 $ mos_t(occa+1,1,1), mos_t(occa+1,1,1), 479 $ dbl_mb(k_full), 'LeftAsymDirac') 480 call rjh_integ_write(dbl_mb(k_full), occa, occa, virta, virta, 481 $ '1<=234', '<ij||ab>', 1111, 1.0d0) 482c 483c 8. <ij||ab> i=j=a=b=beta ijab = 2222 484c 485 call rjh_full_transform(rtdb, basis, 486 $ occb, occb, virtb, virtb, 487 $ nmo, nmo, nmo, nmo, 488 $ mos_t(1,1,2), mos_t(1,1,2), 489 $ mos_t(occb+1,1,2), mos_t(occb+1,1,2), 490 $ dbl_mb(k_full), 'LeftAsymDirac') 491 call rjh_integ_write(dbl_mb(k_full), occb, occb, virtb, virtb, 492 $ '1<=234', '<ij||ab>', 2222, 1.0d0) 493c 494c 9. <ij|ab> i=a=alpha j=b=beta ijab = 1212 495c 496 call rjh_full_transform(rtdb, basis, 497 $ occa, occb, virta, virtb, 498 $ nmo, nmo, nmo, nmo, 499 $ mos_t(1,1,1), mos_t(1,1,2), 500 $ mos_t(occa+1,1,1), mos_t(occb+1,1,2), 501 $ dbl_mb(k_full), 'Dirac') 502 call rjh_integ_write(dbl_mb(k_full), occa, occb, virta, virtb, 503 $ '1234', '<ij||ab>', 1212, 1.0d0) 504c 505c 10. T(ij,ab) i=j=a=b=alpha ijab = 1111 506c i>=j, a>=b 507c 508 call rjh_full_transform(rtdb, basis, 509 $ occa, occa, virta, virta, 510 $ nmo, nmo, nmo, nmo, 511 $ mos_t(1,1,1), mos_t(1,1,1), 512 $ mos_t(occa+1,1,1), mos_t(occa+1,1,1), 513 $ dbl_mb(k_full), 'LeftAsymDirac') 514 call rjh_denoms(dbl_mb(k_full), 515 $ occa, occa, virta, virta, 516 $ evals(1,1), evals(1,1), evals(occa+1,1), evals(occa+1,1)) 517 call rjh_integ_write(dbl_mb(k_full), occa, occa, virta, virta, 518 $ '1<=234', 'T(ij,ec)', 1111, 1.0d0) 519c 520c 11. T(ij,ab) i=j=a=b=beta ijab = 2222 521c 522 call rjh_full_transform(rtdb, basis, 523 $ occb, occb, virtb, virtb, 524 $ nmo, nmo, nmo, nmo, 525 $ mos_t(1,1,2), mos_t(1,1,2), 526 $ mos_t(occb+1,1,2), mos_t(occb+1,1,2), 527 $ dbl_mb(k_full), 'LeftAsymDirac') 528 call rjh_denoms(dbl_mb(k_full), 529 $ occb, occb, virtb, virtb, 530 $ evals(1,2), evals(1,2), evals(occb+1,2), evals(occb+1,2)) 531 call rjh_integ_write(dbl_mb(k_full), occb, occb, virtb, virtb, 532 $ '1<=234', 'T(ij,ec)', 2222, 1.0d0) 533c 534c 12. T(ij,ab) i=a=alpha j=b=beta ijab = 1212 535c 536 call rjh_full_transform(rtdb, basis, 537 $ occa, occb, virta, virtb, 538 $ nmo, nmo, nmo, nmo, 539 $ mos_t(1,1,1), mos_t(1,1,2), 540 $ mos_t(occa+1,1,1), mos_t(occb+1,1,2), 541 $ dbl_mb(k_full), 'Dirac') 542 call rjh_denoms(dbl_mb(k_full), 543 $ occa, occb, virta, virtb, 544 $ evals(1,1), evals(1,2), evals(occa+1,1), evals(occb+1,2)) 545 call rjh_integ_write(dbl_mb(k_full), occa, occb, virta, virtb, 546 $ '1234', 'T(ij,ec)', 1212, 1.0d0) 547c 548c 13. T(ij,ab) i=b=alpha j=a=beta ijab = 1221 549c 550 call rjh_full_transform(rtdb, basis, 551 $ occa, occb, virtb, virta, 552 $ nmo, nmo, nmo, nmo, 553 $ mos_t(1,1,1), mos_t(1,1,2), 554 $ mos_t(occb+1,1,2), mos_t(occa+1,1,1), 555 $ dbl_mb(k_full), 'Dirac') 556 call rjh_denoms(dbl_mb(k_full), 557 $ occa, occb, virta, virtb, 558 $ evals(1,1), evals(1,2), evals(occb+1,2), evals(occa+1,1)) 559 call rjh_integ_write(dbl_mb(k_full), occa, occb, virtb, virta, 560 $ '1234', 'T(ij,ec)', 1221, 1.0d0) 561c 562 if (.not. ma_chop_stack(l_full)) call errquit 563 $ ('rjh_write_integ: ma chop?',0, MA_ERR) 564c 565 end 566 subroutine rjh_denoms(full, dim1, dim2, dim3, dim4, 567 $ e1, e2, e3, e4) 568 implicit none 569 integer dim1, dim2, dim3, dim4 570 double precision full(dim1, dim2, dim3, dim4) 571 double precision e1(*), e2(*), e3(*), e4(*) 572c 573 integer i, j, a, b 574 double precision energy 575c 576 energy = 0.0d0 577 do b = 1, dim4 578 do a = 1, dim3 579 do j = 1, dim2 580 do i = 1, dim1 581 full(i,j,a,b) = full(i,j,a,b) / 582 $ (e1(i) + e2(j) - e3(a) - e4(b)) 583 energy = energy + full(i,j,a,b)**2* 584 $ (e1(i) + e2(j) - e3(a) - e4(b)) 585 enddo 586 enddo 587 enddo 588 enddo 589c 590 write(6,*) ' ENERGY from T2(1) before writing ', energy 591c 592 end 593 subroutine rjh_integ_write(full, dim1, dim2, dim3, dim4, 594 $ order, partstub, spin, sign) 595 implicit none 596#include "errquit.fh" 597#include "eaf.fh" 598#include "inp.fh" 599#include "rjhfileinfo.fh" 600 integer dim1, dim2, dim3, dim4, spin 601 double precision full(dim1, dim2, dim3, dim4), sign 602 character*(*) order, partstub 603c 604 integer fd, ierr, i2, i3, i4, filenum 605 double precision offset 606 character*255 filename, errmsg 607 character*16 stub 608c 609 stub = ' ' 610 write(stub,'(a,''-'',i4)') partstub, spin 611 call util_file_name(stub, .true., .false., filename) 612c 613 call dscal(dim1*dim2*dim3*dim4, sign, full, 1) 614 if (eaf_open(filename, eaf_rw, fd) .ne. 0) 615 $ call errquit('rjh_integ_write: opening file?',0, DISK_ERR) 616 offset = 0.0d0 617c 618** call rjh_debug_print(stub,full,dim1,dim2,dim3,dim4) 619c 620 if (order .eq. '1234') then ! full(i1,i2,i3,i4) -> file(i1,i2,i3,i4) 621 do i4 = 1, dim4 622 do i3 = 1, dim3 623 do i2 = 1, dim2 624 ierr = eaf_write(fd, offset, full(1,i2,i3,i4), 8*dim1) 625 if (ierr .ne. 0) goto 100 626 offset = offset + dim1*8 627 enddo 628 enddo 629 enddo 630 else if (order .eq. '1<=234') then 631 if (dim1.ne.dim2) call errquit('rjh_integ_write: order?',0, 632 & INT_ERR) 633 do i4 = 1, dim4 634 do i3 = 1, dim3 635 do i2 = 1, dim2 636 ierr = eaf_write(fd, offset, full(1,i2,i3,i4), 8*i2) 637 if (ierr .ne. 0) goto 100 638 offset = offset + 8*i2 639 enddo 640 enddo 641 enddo 642 else if (order .eq. '1432') then ! full(i1,i2,i3,i4) -> file(i1,i4,i3,i2) 643 do i2 = 1, dim2 644 do i3 = 1, dim3 645 do i4 = 1, dim4 646 ierr = eaf_write(fd, offset, full(1,i2,i3,i4), 8*dim1) 647 if (ierr .ne. 0) goto 100 648 offset = offset + 8*dim1 649 enddo 650 enddo 651 enddo 652 else 653 call errquit('rjh_integ_write: unknown order',0, INT_ERR) 654 endif 655c 656c Save all of the info in common 657c 658 do filenum = 1, maxfile 659 if (fds(filenum).eq.-1) goto 10 660 enddo 661 call errquit('rjh_integ_write: too many open files?', 0, INT_ERR) 662 10 fds(filenum) = fd 663 stubs(filenum) = stub 664 filenames(filenum) = filename 665 write(6,11) filenum, stub(1:inp_strlen(stub)), 666 $ filename(1:inp_strlen(filename)), order, dim1, dim2, 667 $ dim3, dim4, offset 668 11 format(' wrote file=',i2,2x,a,2x,a,' ',a,2x,4i5,f10.1) 669c 670 return 671c 672 100 call eaf_errmsg(ierr, errmsg) 673 write(6,*) ' filename ', filename(1:inp_strlen(filename)) 674 write(6,*) ' IO offset ', offset 675 write(6,*) ' IO error message ',errmsg(1:inp_strlen(errmsg)) 676 call errquit('rij_integ_write: write failed',0, DISK_ERR) 677c 678 end 679 subroutine get_cct(g_a, type, 680 $ i1i2_lo, i1i2_hi, i1_spin, i2_spin, 681 $ i3_lo, i3_hi, i3_spin, 682 $ i4_lo, i4_hi, i4_spin) 683 implicit none 684#include "errquit.fh" 685#include "eaf.fh" 686#include "rjhfileinfo.fh" 687#include "cct_table_uhf.fh" 688#include "mafdecls.fh" 689#include "inp.fh" 690 691 integer k_all, norb 692 common /c_all/norb, k_all 693 694c 695 integer g_a, i1i2_lo, i1i2_hi, i1_spin, i2_spin, 696 $ i3_lo, i3_hi, i3_spin, i4_lo, i4_hi, i4_spin 697 character*(*) type 698c 699c Read into the global array which is dimensioned 700c . ga(i1i2_lo:i1i2_hi, (i3_lo:i3_hi)*(i4_lo:i4_hi)) 701c from the list described by type which is logically stored 702c as an array (i1, i2, i3, i4) as follows for given spatial 703c symmetries of i3 and i4. 704c 705c do i4 = i4_lo, i4_hi ! All the same symmetry 706c . do i3 = i3_lo, i3_hi ! All the same symmetry 707c . do = sym of i2 708c . -> sym of i1 709c . do i2 in symmetry block 710c . do i1 in symmetry block 711c . ... 712c 713c Spatial symmetry is not yet used. 714c 715c Indices 1 and 2 will be summed over and must be weighted. 716c Index 3 will always be a non-weighted index (e/f/m/n) 717c 718c The input ranges are over SPIN orbitals 719c 720 integer filenum, fd 721 character*8 order 722 character*80 errmsg 723 integer dim1, dim2, dim3, dim12, i4_lo_spatial, i4_hi_spatial, 724 $ i3_lo_spatial, i3_hi_spatial, i1_lo, i2_lo, i1_hi, i2_hi, 725 $ reclen, i1, i2, i3, i4, i3i4, ierr 726 integer i1_occ, i2_occ, i3_occ, i4_occ, ind, l_buff, k_buff 727 integer i1_top, i3_base, i4_base 728 double precision offset, xxx, test 729c 730c 0) determine which physical file 731c 1) determine order of indices in file on disk 732c 2) determine dimensions of each index 733c 3) read the data 734c 4) apply denominators if reading first order T2 735c 5) apply weights to summation indices 736c 737 if (type .eq. '<ba||ek>(bkea)') then 738c 739c i1=b, i2=k, i3=e, i4=a 740c 741 call rjh_which_file('<ba||ek>', 742 $ i1_spin, i4_spin, i3_spin, i2_spin, filenum) 743 order = '1234' 744 i1_occ = 2 ! b 745 i2_occ = 1 ! k 746 i3_occ = 2 ! e 747 i4_occ = 2 ! a 748 else if (type .eq. 'T(ij,ce)(ijec)') then ! ? ORDERING OF INDICES 749c 750c i1=i, i2=j, i3=e, i4=c 751c 752 call rjh_which_file('T(ij,ec)', 753 $ i1_spin, i2_spin, i3_spin, i4_spin, filenum) 754 if (i1_spin .eq. i2_spin) then 755 order = '1<=234' 756 else 757 order = '1234' 758 endif 759 i1_occ = 1 ! i 760 i2_occ = 1 ! j 761 i3_occ = 2 ! c/e 762 i4_occ = 2 ! c/e 763 else 764 call errquit('get_cct: unknown type',0, INPUT_ERR) 765 endif 766c 767 fd = fds(filenum) 768c 769 dim12 = lenxy(1,i1_spin,i2_spin,i1_occ+i2_occ) 770 dim1 = n_occ_virt(1,i1_spin,i1_occ) 771 dim2 = n_occ_virt(1,i2_spin,i2_occ) 772 dim3 = n_occ_virt(1,i3_spin,i3_occ) 773c 774 i4_lo_spatial = i4_lo 775 if (i4_spin .eq. 2) i4_lo_spatial = i4_lo - orb 776 i4_hi_spatial = i4_lo_spatial + (i4_hi - i4_lo) 777 i3_lo_spatial = i3_lo 778 if (i3_spin .eq. 2) i3_lo_spatial = i3_lo - orb 779 i3_hi_spatial = i3_lo_spatial + (i3_hi - i3_lo) 780c 781 i1_lo = s_occ_virt(1,i1_spin,i1_occ) 782 i1_hi = s_occ_virt(2,i1_spin,i1_occ)-1 783 i2_lo = s_occ_virt(1,i2_spin,i2_occ) 784 i2_hi = s_occ_virt(2,i2_spin,i2_occ)-1 785 i3_base = s_occ_virt(1,i3_spin,i3_occ) 786 if (i3_spin .eq. 2) i3_base = i3_base - orb 787 i4_base = s_occ_virt(1,i4_spin,i4_occ) 788 if (i4_spin .eq. 2) i4_base = i4_base - orb 789 790c 791 reclen = i1i2_hi - i1i2_lo + 1 792 if (.not. ma_push_get(mt_dbl,reclen,'buffer',l_buff, k_buff)) 793 $ call errquit('ma buff', reclen, MA_ERR) 794c 795* write(6,*) ' GET_CCT: ', type, i1i2_lo, i1i2_hi, i1_spin, i2_spin, 796* $ i3_lo, i3_hi, i3_spin, 797* $ i4_lo, i4_hi, i4_spin, ' ',order 798c 799* write(6,*) ' FAC' 800* call output(fac,1,orb,1,2,orb,2,1) 801c 802 if (order .eq. '1234' .or. order .eq. '1<=234') then 803 do i4 = i4_lo_spatial, i4_hi_spatial 804 do i3 = i3_lo_spatial, i3_hi_spatial 805 offset = (i1i2_lo - 1 + dim12*(i3-i3_base + 806 $ dim3*(i4-i4_base)))*8.0d0 807 ierr = eaf_read(fd, offset, dbl_mb(k_buff), 8*reclen) 808 if (ierr .ne. 0) goto 100 809* write(6,*) ' THIS IS WHAT I READ ', offset, reclen 810* call output(dbl_mb(k_buff), 1, reclen, 1, 1, reclen,1,1) 811c 812 if (i1i2_lo .ne. 1) stop 99999 813 ind = 0 814 do i2 = i2_lo, i2_hi 815 i1_top = i1_hi 816 if (order .eq. '1<=234') i1_top = i2 817 do i1 = i1_lo, i1_top 818 dbl_mb(k_buff+ind) = dbl_mb(k_buff+ind) * 819 $ fac(i1) * fac(i2) 820 ind = ind + 1 821 enddo 822 enddo 823c 824 if (ind .ne. reclen) call errquit('ind?',ind, 825 & UNKNOWN_ERR) 826c 827 i3i4 = 1 + i3-i3_lo_spatial + 828 $ (i3_hi-i3_lo+1)*(i4-i4_lo_spatial) 829c 830 call ga_put(g_a, 1, reclen, i3i4, i3i4, 831 $ dbl_mb(k_buff), 1) 832c 833* write(6,*) ' ga_put ', 1, reclen, i3i4, i3i4 834* call output(dbl_mb(k_buff), 1, reclen, 1, 1, reclen,1,1) 835c 836 enddo 837 enddo 838 endif 839c 840 if (i1_spin .eq. 2) then 841 i1_lo = i1_lo - norb 842 i1_hi = i1_hi - norb 843 endif 844 if (i2_spin .eq. 2) then 845 i2_lo = i2_lo - norb 846 i2_hi = i2_hi - norb 847 endif 848c$$$ if (type .eq. '<ba||ek>(bkea)' ) then 849c$$$ write(6,*) ' SPINS ', i1_spin, i2_spin, i3_spin, i4_spin 850c$$$ write(6,*) ' LOS ', i1_lo, i2_lo, i3_lo_spatial, 851c$$$ $ i4_lo_spatial 852c$$$ write(6,*) ' HIS ', i1_hi, i2_hi, i3_hi_spatial, i4_hi_spatial 853c$$$ do i4 = i4_lo_spatial, i4_hi_spatial 854c$$$ do i3 = i3_lo_spatial, i3_hi_spatial 855c$$$ i3i4 = 1 + i3-i3_lo_spatial + 856c$$$ $ (i3_hi-i3_lo+1)*(i4-i4_lo_spatial) 857c$$$ ind = 1 858c$$$ do i2 = i2_lo, i2_hi 859c$$$ i1_top = i1_hi 860c$$$ if (order .eq. '1<=234') i1_top = i2 861c$$$ do i1 = i1_lo, i1_top 862c$$$ call ga_get(g_a,ind,ind,i3i4,i3i4,xxx,1) 863c$$$ ind = ind + 1 864c$$$c 865c$$$c XXX should be <i1 i4||i3 i2> = (i1 i3|i4 i2) - (i1 i2|i4 i3) 866c$$$c . b a e k = b e a k b k a e 867c$$$c 868c$$$ test = 0.0d0 869c$$$ if (i1_spin .eq. i3_spin) test = test + 870c$$$ $ dbl_mb(k_all + 871c$$$ $ i1-1+norb*(i3-1+norb*(i2-1+norb*(i4-1)))) 872c$$$ if (i1_spin .eq. i2_spin) test = test - 873c$$$ $ dbl_mb(k_all + 874c$$$ $ i1-1+norb*(i2-1+norb*(i4-1+norb*(i3-1)))) 875c$$$ if (abs(test-xxx).gt.1d-8) then 876c$$$ write(6,1) i1,i2,i3,i4, test, xxx 877c$$$ 1 format(4i5,2f12.8) 878c$$$ else if (abs(test).gt.1d-8) then 879c$$$ write(6,2) i1,i2,i3,i4 880c$$$ 2 format(4i5,' ok') 881c$$$ endif 882c$$$ enddo 883c$$$ enddo 884c$$$ enddo 885c$$$ enddo 886c$$$ endif 887c 888 if (.not. ma_pop_stack(l_buff)) call errquit 889 $ ('get_cct: ma corrupted', 0, MA_ERR) 890c 891 return 892c 893 100 call eaf_errmsg(ierr, errmsg) 894 write(6,*) ' filename ', 895 $ filenames(filenum)(1:inp_strlen(filenames(filenum))) 896 write(6,*) ' IO offset ', offset 897 write(6,*) ' IO error message ',errmsg(1:inp_strlen(errmsg)) 898 call errquit('get_cct: read failed',0, DISK_ERR) 899c 900 end 901 subroutine rjh_which_file(partstub, 902 $ i1_spin, i2_spin, i3_spin, i4_spin, 903 $ filenum) 904 implicit none 905#include "errquit.fh" 906#include "rjhfileinfo.fh" 907#include "inp.fh" 908c 909 integer i1_spin, i2_spin, i3_spin, i4_spin, filenum 910 character*(*) partstub 911c 912c Determine which physical file needs to be read by examining 913c the stubs and spin labels 914c 915 character*16 stub 916 integer ls, spin, i 917c 918 spin = i1_spin*1000 + i2_spin*100 + i3_spin*10 + i4_spin 919 write(stub,'(a,''-'',i4)') partstub, spin 920 ls = inp_strlen(stub) 921c 922 do i = 1, maxfile 923 if (fds(i) .ne. -1) then 924 if (stub(1:ls).eq.stubs(i)(1:ls)) then 925 filenum = i 926 write(6,*) stub, ' -> ', i, filenames(i) 927 return 928 endif 929 endif 930 enddo 931c 932 write(6,*) ' STUB ', stub 933 call errquit(' rjh_which_file: failed to find ', 0, DISK_ERR) 934c 935 end 936 subroutine rjh_full_transform( 937 $ rtdb, basis, 938 $ n1, n2, n3, n4, 939 $ ld1, ld2, ld3, ld4, 940 $ c1t, c2t, c3t, c4t, 941 $ full, order) 942 implicit none 943#include "errquit.fh" 944#include "schwarz.fh" 945#include "bas.fh" 946#include "mafdecls.fh" 947#include "inp.fh" 948c 949 integer rtdb 950 integer basis ! AO basis handle 951 integer n1, n2, n3, n4 ! Dimension of each MO set 952 integer ld1, ld2, ld3, ld4 953 double precision c1t(ld1,*), c2t(ld2,*), ! Transposed MO coeffs 954 $ c3t(ld3,*), c4t(ld4,*) 955 double precision full(n1,n2,n3,n4) 956 character*(*) order 957c 958c Generate the specified block of MO integrals with 959c no assumptions of equivalence between the sets of coefficients. 960c 961c Order can be either 962c . ChargeCloud -> full(p,q,r,s) = (pq|rs) 963c or 964c . Dirac -> full(p,q,r,s) = <pq|rs> 965c or 966c . LeftAsymDirac -> full(p,q,r,s) = <pq|rs>-<qp|rs> 967c . (must have c1t=c2t, n1=n2) 968c or 969c . RightAsymDirac -> full(p,q,r,s) = <pq|rs>-<pq|sr> 970c . (must have c3t=c4t, n3=n4) 971c 972c Presently the antisymmetrization is done at the top level 973c and the storage of full is not reduced to use the symmetry. 974c 975c Memory requirements are 976c . n1*n2*n3*n4 + S*n2*n3*n4 + S*S*n3*n4 + S*S*S*n4 + 977c . maxs2 + maxg2*(1 + 4*integer) 978c 979c Index 4 is the first transformed so there is advantage in 980c making it the smallest range. 981c 982 double precision tol2e 983 parameter (tol2e = 1d-12) 984C 985 integer nsh, maxbfsh, lenhalf, lenthird, geom 986 integer l_half, k_half, l_third, k_third 987 integer lsh, ksh, llo, lhi, klo, khi 988 integer p, q, r, s 989 logical ochargecloud, oasym 990 character*8 side 991c 992 if (inp_compare(.false.,order,'chargecloud')) then 993 ochargecloud = .true. 994 oasym = .false. 995 side = ' ' 996 else if (inp_compare(.false.,order,'dirac')) then 997 ochargecloud = .false. 998 oasym = .false. 999 side = ' ' 1000 else if (inp_compare(.false.,order,'leftasymdirac')) then 1001 ochargecloud = .false. 1002 oasym = .true. 1003 side = 'left' 1004 else if (inp_compare(.false.,order,'rightasymdirac')) then 1005 ochargecloud = .false. 1006 oasym = .true. 1007 side = 'right' 1008 else 1009 call errquit('rjh_full_trans: unkown integral option',0, 1010 & INPUT_ERR) 1011 endif 1012c 1013c Initialize integrals and Schwarz screening 1014c 1015 if (.not. bas_geom(basis, geom)) 1016 $ call errquit('rjh_transform: basis ', basis, BASIS_ERR) 1017 call int_init(rtdb, 1, basis) 1018 call schwarz_init(geom, basis) 1019c 1020 if (.not. bas_numcont(basis, nsh)) call errquit( 1021 $ 'rjh_transform: bas_numcont', basis, BASIS_ERR) 1022 if (.not. bas_nbf_cn_max(basis,maxbfsh)) call errquit( 1023 $ 'rjh_transform: bas_nbf_cn_max', basis, BASIS_ERR) 1024c 1025 lenhalf = n3*n4*maxbfsh**2 1026 lenthird= n2*n3*n4*maxbfsh 1027c 1028 if (.not. ma_push_get(mt_dbl,lenhalf,'half',l_half, k_half)) 1029 $ call errquit('ma half', lenhalf, MA_ERR) 1030 if (.not. ma_push_get(mt_dbl,lenthird,'third',l_third, k_third)) 1031 $ call errquit('ma third', lenthird, MA_ERR) 1032c 1033 call dfill(n1*n2*n3*n4, 0.0d0, full, 1) 1034 do ksh = 1, nsh 1035 if (.not. bas_cn2bfr(basis, ksh, klo, khi)) 1036 $ call errquit('rjh_transform: bas_cn2bfr',basis, BASIS_ERR) 1037 call dfill(n2*n3*n4*(khi-klo+1), 0.0d0, dbl_mb(k_third), 1) 1038 do lsh = 1, nsh 1039 if (.not. bas_cn2bfr(basis, lsh, llo, lhi)) 1040 $ call errquit('rjh_transform: bas_cn2bfr',basis, 1041 & BASIS_ERR) 1042 if (schwarz_shell(ksh,lsh)*schwarz_max() 1043 $ .gt. tol2e) then 1044c 1045c Make (rs|kl) all rs (indices 3 and 4) given shells k and l 1046c 1047 call rjh_half_transform(basis, ksh, lsh, n3, n4, 1048 $ c3t, c4t, ld3, ld4, 1049 $ dbl_mb(k_half), ochargecloud, tol2e) 1050 1051* write(6,*) ' ksh, lsh ', ksh, lsh 1052* call rjh_debug_print('half', 1053* $ dbl_mb(k_half), n3, n4, khi-klo+1, lhi-llo+1) 1054 1055c 1056 call rjh_third_transform(llo, lhi, klo, khi, 1057 $ n2, n3, n4, dbl_mb(k_half), dbl_mb(k_third), 1058 $ c2t, ld2, tol2e) 1059 end if 1060 end do 1061* write(6,*) ' lsh ', lsh 1062* call rjh_debug_print('third', 1063* $ dbl_mb(k_third), n2, n3, n4, lhi-llo+1) 1064 call rjh_final_transform(klo, khi, n1, n2, n3, n4, 1065 $ dbl_mb(k_third), full, c1t, ld1, tol2e) 1066 end do 1067c 1068 if (oasym) call rjh_asym_trans(full,n1,n2,n3,n4,side) 1069c 1070 do s = 1, n4 1071 do r = 1, n3 1072 do q = 1, n2 1073 do p = 1, n1 1074 if (abs(full(p,q,r,s)).lt.1d-10) 1075 $ full(p,q,r,s) = 0.0d0 1076 end do 1077 end do 1078 end do 1079 end do 1080c 1081 if (.not. ma_pop_stack(l_third)) call errquit('ma third',0, 1082 & MA_ERR) 1083 if (.not. ma_pop_stack(l_half)) call errquit('ma half',0, MA_ERR) 1084c 1085 call schwarz_tidy() 1086 call int_terminate 1087c 1088 end 1089 subroutine rjh_asym_trans(full,n1,n2,n3,n4,side) 1090 implicit none 1091#include "errquit.fh" 1092c 1093 integer n1, n2, n3, n4 1094 double precision full(n1,n2,n3,n4) 1095 character*(*) side 1096c 1097 integer p, q, r, s 1098 double precision tmp 1099c 1100 if (side .eq. 'left') then 1101 if (n1 .ne. n2) call errquit('rjh_asym_trans: left', n1, 1102 & UNKNOWN_ERR) 1103 do s = 1, n4 1104 do r = 1, n3 1105 do q = 1, n2 1106 do p = 1, q 1107 tmp = full(p,q,r,s) - full(q,p,r,s) 1108 full(p,q,r,s) = tmp 1109 full(q,p,r,s) =-tmp 1110 enddo 1111 enddo 1112 enddo 1113 enddo 1114 else 1115 if (n3 .ne. n4) call errquit('rjh_asym_trans: right', n3, 1116 & UNKNOWN_ERR) 1117 do s = 1, n4 1118 do r = 1, s 1119 do q = 1, n2 1120 do p = 1, n1 1121 tmp = full(p,q,r,s) - full(q,p,s,r) 1122 full(p,q,r,s) = tmp 1123 full(p,q,s,r) =-tmp 1124 enddo 1125 enddo 1126 enddo 1127 enddo 1128 endif 1129c 1130 end 1131 subroutine rjh_final_transform(klo, khi, 1132 $ n1, n2, n3, n4, third, full, c1t, ld1, tol2e) 1133 implicit none 1134c 1135 integer klo, khi, n1, n2, n3, n4, ld1 1136 double precision third(n2,n3,n4,klo:khi) 1137 double precision full(n1,n2,n3,n4) 1138 double precision c1t(ld1,*) 1139 double precision tol2e 1140c 1141 integer k, s, r, q, p 1142 double precision g 1143c 1144 do k = klo, khi 1145 do s = 1, n4 1146 do r = 1, n3 1147 do q = 1, n2 1148 g = third(q,r,s,k) 1149 if (abs(g) .gt. tol2e) then 1150 do p = 1, n1 1151 full(p,q,r,s) = full(p,q,r,s) + g*c1t(p,k) 1152 end do 1153 end if 1154 end do 1155 end do 1156 end do 1157 end do 1158c 1159 end 1160 subroutine rjh_third_transform(llo, lhi, klo, khi, 1161 $ n2, n3, n4, half, third, c2t, ld2, tol2e) 1162 implicit none 1163c 1164 integer llo, lhi, klo, khi, n2, n3, n4, ld2 1165 double precision half(n3,n4,klo:khi,llo:lhi) 1166 double precision third(n2,n3,n4,klo:khi) 1167 double precision c2t(ld2,*) 1168 double precision tol2e 1169c 1170 integer k, l, s, r, q 1171 double precision g 1172c 1173 do l = llo, lhi 1174 do k = klo, khi 1175 do s = 1, n4 1176 do r = 1, n3 1177 g = half(r,s,k,l) 1178 if (abs(g) .gt. tol2e) then 1179 do q = 1, n2 1180 third(q,r,s,k) = third(q,r,s,k) + g*c2t(q,l) 1181 end do 1182 end if 1183 end do 1184 end do 1185 end do 1186 end do 1187c 1188 end 1189 subroutine rjh_half_transform(basis, ksh, lsh, n1, n2, 1190 $ c1t, c2t, ld1, ld2, half, ochargecloud, tol2e) 1191 implicit none 1192#include "errquit.fh" 1193#include "bas.fh" 1194#include "mafdecls.fh" 1195 integer basis, ksh, lsh, n1, n2, ld1, ld2 1196 double precision c1t(*), c2t(*) ! Transposed MO coeffs 1197 double precision half(*) ! n1*n2*kdim*ldim (pq|kl) 1198 logical ochargecloud 1199 double precision tol2e 1200c 1201c For a pair of shells k and l fill 1202c 1203c . half(p,q,k,l) = (pq|kl) 1204c 1205c for all p, q, and k, l within their respective shells 1206c where p and q are transformed into the new bases and 1207c k and l are AO indices. 1208c 1209c Eventually exploiting sparsity and abelian symmetry 1210c ... should also eventually use the texas integrals 1211c 1212c Assumes that integrals and schwarz have been initialized. 1213c 1214 integer nbf, nsh, k_i, k_j, k_k, k_l, l_i, l_j, l_k, l_l, 1215 $ maxg2, maxs2, k_buf, l_buf, k_scr, l_scr, maxbfsh 1216 integer llo, lhi, klo, khi, lenaobuf, l_aobuf, k_aobuf 1217c 1218c Get dimensions and required scratch space info 1219c 1220 if (.not. bas_numcont(basis, nsh)) call errquit( 1221 $ 'rjh_transform: bas_numcont', basis, BASIS_ERR) 1222 if (.not. bas_nbf_cn_max(basis,maxbfsh)) call errquit( 1223 $ 'rjh_transform: bas_nbf_cn_max', basis, BASIS_ERR) 1224 if (.not. bas_numbf(basis, nbf)) 1225 $ call errquit('rjh_transform: nbf',basis, BASIS_ERR) 1226 if (.not. bas_cn2bfr(basis, ksh, klo, khi)) 1227 $ call errquit('rjh_transform: bas_cn2bfr',basis, BASIS_ERR) 1228 if (.not. bas_cn2bfr(basis, lsh, llo, lhi)) 1229 $ call errquit('rjh_transform: bas_cn2bfr',basis, BASIS_ERR) 1230 call int_mem_2e4c(maxg2,maxs2) 1231 lenaobuf = (khi-klo+1)*(lhi-llo+1)*maxbfsh*n2 ! (iq|kl) 1232c 1233c Allocate scratch space for integrals and buffers for 1234c transformation 1235c 1236 if (.not. ma_push_get(mt_dbl,maxs2,'scr',l_scr, k_scr)) 1237 $ call errquit('ma scr',maxs2, MA_ERR) 1238 if (.not. ma_push_get(mt_dbl,maxg2,'buf',l_buf, k_buf)) 1239 $ call errquit('ma buf',maxg2, MA_ERR) 1240 if (.not. ma_push_get(mt_int,maxg2,'i',l_i, k_i)) 1241 $ call errquit('ma ibuf',maxg2, MA_ERR) 1242 if (.not. ma_push_get(mt_int,maxg2,'j',l_j, k_j)) 1243 $ call errquit('ma jbuf',maxg2, MA_ERR) 1244 if (.not. ma_push_get(mt_int,maxg2,'k',l_k, k_k)) 1245 $ call errquit('ma kbuf',maxg2, MA_ERR) 1246 if (.not. ma_push_get(mt_int,maxg2,'l',l_l, k_l)) 1247 $ call errquit('ma lbuf',maxg2, MA_ERR) 1248 if (.not. ma_push_get(mt_dbl,lenaobuf,'aobuf',l_aobuf, k_aobuf)) 1249 $ call errquit('ma aobuf',lenaobuf, MA_ERR) 1250c 1251 call rjh_do_half_transform( 1252 $ basis, 1253 $ dbl_mb(k_buf), dbl_mb(k_scr), 1254 $ int_mb(k_i), int_mb(k_j), int_mb(k_k), int_mb(k_l), 1255 $ maxg2, maxs2, 1256 $ nbf, nsh, maxbfsh, n1, n2, tol2e, 1257 $ ksh, lsh, klo, khi, llo, lhi, 1258 $ half, c1t, c2t, ld1, ld2, dbl_mb(k_aobuf), ochargecloud) 1259c 1260 if (.not. ma_chop_stack(l_scr)) call errquit 1261 $ ('rjh_transform: chopping stack', 0, MA_ERR) 1262c 1263 end 1264 subroutine rjh_do_half_transform( 1265 $ basis, 1266 $ buf, scr, ilab, jlab, klab, llab, maxg2, maxs2, 1267 $ nbf, nsh, maxbfsh, n1, n2, tol2e, 1268 $ ksh, lsh, klo, khi, llo, lhi, 1269 $ half, c1t, c2t, ld1, ld2, aobuf, ochargecloud) 1270 implicit none 1271#include "errquit.fh" 1272#include "schwarz.fh" 1273#include "bas.fh" 1274 integer basis 1275 integer maxg2, maxs2 1276 integer nbf, nsh, n1, n2, maxbfsh, ld1, ld2 1277 double precision buf(maxg2), scr(maxs2) 1278 double precision c1t(ld1,nbf), c2t(ld2,nbf) 1279 integer ilab(maxg2), jlab(maxg2), klab(maxg2), llab(maxg2) 1280 integer ksh, lsh, klo, khi, llo, lhi 1281 double precision aobuf(n2,maxbfsh,klo:khi,llo:lhi) 1282 double precision half(n1,n2,klo:khi,llo:lhi) 1283 double precision tol2e 1284 logical ochargecloud 1285c 1286c For a pair of shells k and l, fill aobuf with integrals (ij|kl) 1287c for all i>=j ... eventually exploiting sparsity and abelian symmetry 1288c ... should also use the texas integrals 1289c 1290 double precision skl, g 1291 integer ish, jsh, i, j, k, l, p, q, ijkl, nint, ilo, ihi, 1292 $ idim, kdim, ldim 1293c 1294 kdim = khi - klo + 1 1295 ldim = lhi - llo + 1 1296 skl = schwarz_shell(ksh,lsh) 1297 call dfill(n1*n2*kdim*ldim,0.0d0,half,1) 1298c 1299 do ish = 1, nsh 1300 if (.not. bas_cn2bfr(basis,ish,ilo,ihi)) 1301 $ call errquit('rjh_do_half_transform',ish, BASIS_ERR) 1302 idim = ihi-ilo+1 1303 do l = llo,lhi 1304 do k = klo, khi 1305 do i = 1, idim 1306 do q = 1, n2 1307 aobuf(q,i,k,l) = 0.0d0 1308 end do 1309 end do 1310 end do 1311 end do 1312c 1313 do jsh = 1, nsh 1314 if (ochargecloud) then ! (ij|kl) 1315 if (schwarz_shell(ish,jsh)*skl .gt. tol2e) then 1316 call int_l2e4c(basis, ish, jsh, basis, ksh, lsh, 1317 & tol2e, .false., maxg2, buf, nint, 1318 $ ilab, jlab, klab, llab, maxs2, scr) 1319 do ijkl = 1, nint 1320 i = ilab(ijkl)-ilo+1 1321 j = jlab(ijkl) 1322 k = klab(ijkl) 1323 l = llab(ijkl) 1324 g = buf(ijkl) 1325 if (abs(g) .gt. tol2e) then 1326 do q = 1, n2 1327 aobuf(q,i,k,l) = aobuf(q,i,k,l) + g*c2t(q,j) 1328 end do 1329 end if 1330 end do 1331 end if 1332 else ! <ij|kl> = (ik|jl) 1333 if (schwarz_shell(ish,ksh)*schwarz_shell(jsh,lsh) 1334 $ .gt. tol2e) then 1335 call int_l2e4c(basis, ish, ksh, basis, jsh, lsh, 1336 & tol2e, .false., maxg2, buf, nint, 1337 $ ilab, klab, jlab, llab, maxs2, scr) 1338 do ijkl = 1, nint 1339 i = ilab(ijkl)-ilo+1 1340 j = jlab(ijkl) 1341 k = klab(ijkl) 1342 l = llab(ijkl) 1343 g = buf(ijkl) 1344 if (abs(g) .gt. tol2e) then 1345 do q = 1, n2 1346 aobuf(q,i,k,l) = aobuf(q,i,k,l) + g*c2t(q,j) 1347 end do 1348 end if 1349 end do 1350 end if 1351 endif 1352 end do 1353 do l = llo, lhi 1354 do k = klo, khi 1355 do i = ilo, ihi 1356 do q = 1, n2 1357 g = aobuf(q,i-ilo+1,k,l) 1358 if (abs(g) .gt. tol2e) then 1359 do p = 1, n1 1360 half(p,q,k,l) = half(p,q,k,l) + g*c1t(p,i) 1361 end do 1362 end if 1363 end do 1364 end do 1365 end do 1366 end do 1367 end do 1368c 1369 end 1370