1 logical function jantest(rtdb) 2* 3* $$ 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 integer rtdb 14c 15 integer basis, geom, nbf 16 integer l_aoint, k_aoint ! nbf**4 array of AO integrals 17 integer l_moint, k_moint ! nbf**4 array of AO integrals 18 character*255 movecs ! Name of movector file 19 character*80 title, name_of_basis, scftype 20 integer nbf_file, nsets, nmo_file(2) 21 logical movecs_read, movecs_read_header 22 external movecs_read, movecs_read_header 23c 24 integer nmo, g_tmp, l_occ, k_occ, 25 $ l_eval, k_eval, l_mos, k_mos, l_most, k_most 26 integer nocc, nvirt, nopen, nclosed, nso, noso 27 integer l_t, k_t, l_t2, k_t2, l_f, k_f 28 integer l_t1, k_t1, l_r1, k_r1 29 integer l_tg, k_tg, l_tf, k_tf 30 integer l_w, k_w, l_v, k_v 31 logical int_normalize 32 external int_normalize 33c 34c load the geometry/basis set and get info 35c 36 if (.not. geom_create(geom, 'geometry')) 37 $ call errquit('scf_init: geom_create?', 0, GEOM_ERR) 38 if (.not. geom_rtdb_load(rtdb, geom, 'geometry')) 39 $ call errquit('scf_init: no geometry ', 0, RTDB_ERR) 40 if (.not. bas_create(basis, 'ao basis')) 41 $ call errquit('scf_init: bas_create?', 0, BASIS_ERR) 42 if (.not. bas_rtdb_load(rtdb, geom, basis, 'ao basis')) 43 $ call errquit('scf_init: no ao basis set', 0, RTDB_ERR) 44 if (.not.int_normalize(rtdb,basis)) 45 $ call errquit('scf:int_normalize failed', 0, INT_ERR) 46 if (.not. bas_numbf(basis, nbf)) call errquit 47 $ ('scf_init: basis info',0, BASIS_ERR) 48c 49c Read the MO vectors and evals from a RHF calculation 50c 51 call util_file_name('movecs',.false.,.false.,movecs) 52 if (.not. movecs_read_header(movecs, title, name_of_basis, 53 $ scftype, nbf_file, nsets, nmo_file, 2)) call errquit 54 $ ('jantest: failed to read movecs header',911, DISK_ERR) 55 write(6,*) ' Read movecs header from ', movecs 56 write(6,*) ' Job title : ', 57 $ title(1:inp_strlen(title)) 58 write(6,*) ' Basis name: ', 59 $ name_of_basis(1:inp_strlen(name_of_basis)) 60 nmo = nmo_file(1) 61 if (.not. rtdb_get(rtdb, 'scf:nclosed', mt_int, 1, nocc)) 62 $ call errquit('nocc?',0, RTDB_ERR) 63 if (.not. rtdb_get(rtdb, 'scf:nopen', mt_int, 1, nopen)) 64 $ call errquit('nopen?',0, RTDB_ERR) 65 if (nopen .ne. 0) call errquit('asjdlfkadjsl',0, UNKNOWN_ERR) 66 nvirt= nmo - nocc 67 write(6,*) ' No. of closed shells ', nocc 68 write(6,*) ' No. of molecular orbitals: ', nmo 69 write(6,*) ' No. of basis functions: ', nbf 70c 71*ga:1:0 72 if (.not. ga_create(mt_dbl, nbf, nmo, 'tmp', 0, 0, g_tmp)) 73 & call errquit('scf_v_g: tmp', 0, GA_ERR) 74 if (.not. ma_push_get(mt_dbl, nbf,'occ',l_occ, k_occ)) 75 $ call errquit('ma occ', nbf, MA_ERR) 76 if (.not. ma_push_get(mt_dbl, nbf,'eval',l_eval, k_eval)) 77 $ call errquit('ma eval', nbf, MA_ERR) 78 if (.not. ma_push_get(mt_dbl, nbf*nbf,'mos', l_mos, k_mos)) 79 $ call errquit('ma mos', nbf*nbf, MA_ERR) 80 if (.not. ma_push_get(mt_dbl, nbf*nbf,'mos', l_most, k_most)) 81 $ call errquit('ma mos', nbf*nbf, MA_ERR) 82c 83 if (.not. movecs_read(movecs, 1, dbl_mb(k_occ), dbl_mb(k_eval), 84 $ g_tmp)) call errquit('movecs_read of amos failed ',0, 85 & DISK_ERR) 86 call ga_get(g_tmp, 1, nbf, 1, nmo, dbl_mb(k_mos), nbf) 87 call util_transpose(dbl_mb(k_mos),nbf,dbl_mb(k_most),nmo, 88 $ nbf,nmo) 89c 90 write(6,*) ' Orbital eigenvalues ' 91 call output(dbl_mb(k_eval),1,nmo,1,1,nmo,1,1) 92 write(6,*) ' MOs' 93 call output(dbl_mb(k_mos),1,nbf,1,nmo,nbf,nmo,1) 94 write(6,*) ' MOs T' 95 call output(dbl_mb(k_most),1,nmo,1,nbf,nmo,nbf,1) 96c 97 if (.not. ga_destroy(g_tmp)) call errquit(' ga bad?',0, GA_ERR) 98c 99c Make all AO integrals 100c 101 if (.not. ma_push_get(mt_dbl,nbf**4,'aoint',l_aoint,k_aoint)) 102 $ call errquit('allocation of AO integrals failed',nbf**4, 103 & MA_ERR) 104 call jan_all_ao_integrals(rtdb,basis,nbf,'dirac',dbl_mb(k_aoint)) 105c call jan_debug_print('AOINTS',dbl_mb(k_aoint), nbf, nbf, nbf, 106c $ nbf) 107c 108c Make all MO integrals in Dirac order 109c 110 if (.not. ma_push_get(mt_dbl,nmo**4,'moint',l_moint,k_moint)) 111 $ call errquit('allocation of MO integrals failed',nbf**4, 112 & MA_ERR) 113 call jan_full_transform( 114 $ rtdb, basis, 115 $ nmo, nmo, nmo, nmo, 116 $ nmo, nmo, nmo, nmo, 117 $ dbl_mb(k_most),dbl_mb(k_most),dbl_mb(k_most),dbl_mb(k_most), 118 $ dbl_mb(k_moint), 'Dirac') 119c call jan_debug_print('MOINTS',dbl_mb(k_moint), nmo, nmo, nmo, 120c $ nmo) 121c 122c do some incore cc 123c set nso=2*nmo, noso=2*nocc 124c 125 nso=2*nmo 126 noso=2*nocc 127 if (.not. ma_push_get(mt_dbl,nso**4,'t amps',l_t,k_t)) 128 $ call errquit('allocation of t amplitudes failed',nso**4, 129 & MA_ERR) 130 if (.not. ma_push_get(mt_dbl,nso**4,'t2 amps',l_t2,k_t2)) 131 $ call errquit('allocation of t2 amplitudes failed',nso**4, 132 & MA_ERR) 133 if (.not. ma_push_get(mt_dbl,nmo**2,'Fock Matrix',l_f,k_f)) 134 $ call errquit('allocation of t2 amplitudes failed',nmo**2, 135 & MA_ERR) 136 if (.not. ma_push_get(mt_dbl,nso**2,'t1 amps',l_t1,k_t1)) 137 $ call errquit('allocation of t1 amplitudes failed',nso**2, 138 & MA_ERR) 139 if (.not. ma_push_get(mt_dbl,nso**2,'t1 resi',l_r1,k_r1)) 140 $ call errquit('allocation of t1 residual failed',nso**2, 141 & MA_ERR) 142c 143 call ccsd_incore(rtdb, basis, dbl_mb(k_moint), 144 & dbl_mb(k_eval), dbl_mb(k_t), dbl_mb(k_t2), 145 & dbl_mb(k_f), dbl_mb(k_t1), dbl_mb(k_r1), 146 & nbf, nmo, nocc, nso, noso) 147c 148 if (.not. ma_push_get(mt_dbl,noso**3*nso**3,'t3 tg',l_tg,k_tg)) 149 $ call errquit('allocation of t3 tg failed',nso**6, 150 & MA_ERR) 151 if (.not. ma_push_get(mt_dbl,noso**3*nso**3,'t3 tf',l_tg,k_tf)) 152 $ call errquit('allocation of t3 tf failed',nso**6, MA_ERR) 153 if (.not. ma_push_get(mt_dbl,noso**3*nso**3,'t3 w',l_tg,k_w)) 154 $ call errquit('allocation of t3 w failed',nso**6, MA_ERR) 155 if (.not. ma_push_get(mt_dbl,noso**3*nso**3,'t3 v',l_tg,k_v)) 156 $ call errquit('allocation of t3 v failed',nso**6, MA_ERR) 157c 158 if (.not. ga_create(mt_dbl, nbf, nmo, 'tmp', 0, 0, g_tmp)) 159 & call errquit('scf_v_g: tmp', 0, GA_ERR) 160 if (.not. movecs_read(movecs, 1, dbl_mb(k_occ), dbl_mb(k_eval), 161 $ g_tmp)) call errquit('movecs_read of amos failed ',0, 162 & DISK_ERR) 163 call ga_get(g_tmp, 1, nbf, 1, nmo, dbl_mb(k_mos), nbf) 164 call util_transpose(dbl_mb(k_mos),nbf,dbl_mb(k_most),nmo, 165 $ nbf,nmo) 166 if (.not. ga_destroy(g_tmp)) call errquit(' ga bad?',0, GA_ERR) 167c 168 call jan_full_transform( 169 $ rtdb, basis, 170 $ nmo, nmo, nmo, nmo, 171 $ nmo, nmo, nmo, nmo, 172 $ dbl_mb(k_most),dbl_mb(k_most),dbl_mb(k_most),dbl_mb(k_most), 173 $ dbl_mb(k_moint), 'Dirac') 174c 175 call triples_incore(rtdb, basis, 176 & dbl_mb(k_moint), 177 & dbl_mb(k_eval), dbl_mb(k_t2), 178 & dbl_mb(k_t1), 179 & dbl_mb(k_tg), dbl_mb(k_tf), 180 & dbl_mb(k_w), dbl_mb(k_v), 181 & nbf, nmo, nocc, nso, noso) 182c 183c Tidy up 184c 185 if (.not. ma_chop_stack(l_occ)) call errquit(' ma chop?', 0, 186 & MA_ERR) 187 if (.not. bas_destroy(basis)) call errquit(' bas ?',0, BASIS_ERR) 188 if (.not. geom_destroy(geom)) call errquit(' geom ?',0, GEOM_ERR) 189 jantest = .true. 190c 191 end 192 subroutine jan_all_ao_integrals(rtdb, basis, nbf, order, ao) 193 implicit none 194#include "errquit.fh" 195#include "mafdecls.fh" 196#include "bas.fh" 197 integer rtdb, basis, nbf 198 double precision ao(nbf,nbf,nbf,nbf) 199 character*(*) order 200c 201 integer nsh, k_i, k_j, k_k, k_l, l_i, l_j, l_k, l_l, 202 $ maxg2, maxs2, k_buf, l_buf, k_scr, l_scr 203c 204 call int_init(rtdb, 1, basis) 205 if ( .not. bas_numcont(basis, nsh) ) call errquit( 206 $ 'ao_fock_2e: problem with call to bas_numcont', basis, 207 & BASIS_ERR) 208 call int_mem_2e4c(maxg2,maxs2) 209 if (.not. ma_push_get(mt_dbl,maxs2,'scr',l_scr, k_scr)) 210 $ call errquit('ma scr',maxg2, MA_ERR) 211 if (.not. ma_push_get(mt_dbl,maxg2,'buf',l_buf, k_buf)) 212 $ call errquit('ma buf',maxg2, MA_ERR) 213 if (.not. ma_push_get(mt_int,maxg2,'i',l_i, k_i)) 214 $ call errquit('ma ibuf',maxg2, MA_ERR) 215 if (.not. ma_push_get(mt_int,maxg2,'j',l_j, k_j)) 216 $ call errquit('ma jbuf',maxg2, MA_ERR) 217 if (.not. ma_push_get(mt_int,maxg2,'k',l_k, k_k)) 218 $ call errquit('ma kbuf',maxg2, MA_ERR) 219 if (.not. ma_push_get(mt_int,maxg2,'l',l_l, k_l)) 220 $ call errquit('ma lbuf',maxg2, MA_ERR) 221c 222 call jan_do_all_ao_integrals(basis, dbl_mb(k_buf), dbl_mb(k_scr), 223 $ int_mb(k_i), int_mb(k_j), int_mb(k_k), int_mb(k_l), 224 $ maxg2, maxs2, nbf, nsh, order, ao) 225c 226 if (.not. ma_chop_stack(l_scr)) call errquit('janallao: ma?',0, 227 & MA_ERR) 228c 229 call int_terminate 230c 231 end 232 subroutine jan_do_all_ao_integrals( 233 $ basis, buf, scr, ilab, jlab, klab, 234 $ llab, maxg2, maxs2, nbf, nsh, order, ao) 235 implicit none 236#include "errquit.fh" 237c 238 integer basis, nbf, nsh, maxg2, maxs2 239 double precision buf(maxg2), scr(maxs2) 240 integer ilab(maxg2), jlab(maxs2), klab(maxs2), llab(maxs2) 241 integer i, j, k ,l, ish, jsh, ksh, lsh, ijkl, nint 242 character*(*) order 243 double precision ao(nbf,nbf,nbf,nbf) 244 double precision zerotol 245 logical omulliken 246c 247 omulliken = .false. ! avoids compiler warning 248 if (order .eq. 'mulliken') then 249 omulliken = .true. 250 else if (order .eq. 'dirac') then 251 omulliken = .false. 252 else 253 call errquit(' unknown order',0, UNKNOWN_ERR) 254 end if 255c 256 call dfill(nbf**4, 0.0d0, ao, 1) 257 zerotol = 1d-12 258c 259 do ish = 1, nsh 260 do jsh = 1, nsh 261 do ksh = 1, nsh 262 do lsh = 1, nsh 263 call int_l2e4c(basis, ish, jsh, basis, ksh, lsh, 264 & zerotol, .false., maxg2, buf, nint, 265 $ ilab, jlab, klab, llab, maxs2, scr) 266 do ijkl = 1, nint 267 i = ilab(ijkl) 268 j = jlab(ijkl) 269 k = klab(ijkl) 270 l = llab(ijkl) 271 if (omulliken) then 272 ao(i,j,k,l) = buf(ijkl) 273 else 274 ao(i,k,j,l) = buf(ijkl) 275 end if 276 end do 277 end do 278 end do 279 end do 280 end do 281c 282c$$$ write(6,*) 283c$$$ write(6,*) ' AO integrals ' 284c$$$ write(6,*) 285c$$$ do i = 1, nbf 286c$$$ do j = 1, nbf 287c$$$ do k = 1, nbf 288c$$$ do l = 1, nbf 289c$$$ if ( abs(ao(i,j,k,l)) .gt. 1e-6 ) 290c$$$ $ write(6,7) i,j,k,l,ao(i,j,k,l) 291c$$$ 7 format(1x,4i5,2x,f12.6) 292c$$$ end do 293c$$$ end do 294c$$$ end do 295c$$$ end do 296c 297 end 298 subroutine jan_full_transform( 299 $ rtdb, basis, 300 $ n1, n2, n3, n4, 301 $ ld1, ld2, ld3, ld4, 302 $ c1t, c2t, c3t, c4t, 303 $ full, order) 304 implicit none 305#include "errquit.fh" 306#include "schwarz.fh" 307#include "bas.fh" 308#include "mafdecls.fh" 309#include "inp.fh" 310c 311 integer rtdb 312 integer basis ! AO basis handle 313 integer n1, n2, n3, n4 ! Dimension of each MO set 314 integer ld1, ld2, ld3, ld4 315 double precision c1t(ld1,*), c2t(ld2,*), ! Transposed MO coeffs 316 $ c3t(ld3,*), c4t(ld4,*) 317 double precision full(n1,n2,n3,n4) 318 character*(*) order 319c 320c Generate the specified block of MO integrals with 321c no assumptions of equivalence between the sets of coefficients. 322c 323c Order can be either 324c . ChargeCloud -> full(p,q,r,s) = (pq|rs) 325c or 326c . Dirac -> full(p,q,r,s) = <pq|rs> 327c or 328c . LeftAsymDirac -> full(p,q,r,s) = <pq|rs>-<qp|rs> 329c . (must have c1t=c2t, n1=n2) 330c or 331c . RightAsymDirac -> full(p,q,r,s) = <pq|rs>-<pq|sr> 332c . (must have c3t=c4t, n3=n4) 333c 334c Presently the antisymmetrization is done at the top level 335c and the storage of full is not reduced to use the symmetry. 336c 337c Memory requirements are 338c . n1*n2*n3*n4 + S*n2*n3*n4 + S*S*n3*n4 + S*S*S*n4 + 339c . maxs2 + maxg2*(1 + 4*integer) 340c 341c Index 4 is the first transformed so there is advantage in 342c making it the smallest range. 343c 344 double precision tol2e 345 parameter (tol2e = 1d-12) 346C 347 integer nsh, maxbfsh, lenhalf, lenthird, geom 348 integer l_half, k_half, l_third, k_third 349 integer lsh, ksh, llo, lhi, klo, khi 350 integer p, q, r, s 351 logical ochargecloud, oasym 352 character*8 side 353c 354 ochargecloud = .false. ! avoids compiler warning 355 oasym = .false. ! avoids compiler warning 356 if (inp_compare(.false.,order,'chargecloud')) then 357 ochargecloud = .true. 358 oasym = .false. 359 side = ' ' 360 else if (inp_compare(.false.,order,'dirac')) then 361 ochargecloud = .false. 362 oasym = .false. 363 side = ' ' 364 else if (inp_compare(.false.,order,'leftasymdirac')) then 365 ochargecloud = .false. 366 oasym = .true. 367 side = 'left' 368 else if (inp_compare(.false.,order,'rightasymdirac')) then 369 ochargecloud = .false. 370 oasym = .true. 371 side = 'right' 372 else 373 call errquit('jan_full_trans: unkown integral option',0, 374 & INT_ERR) 375 endif 376c 377c Initialize integrals and Schwarz screening 378c 379 if (.not. bas_geom(basis, geom)) 380 $ call errquit('jan_transform: basis ', basis, BASIS_ERR) 381 call int_init(rtdb, 1, basis) 382 call schwarz_init(geom, basis) 383c 384 if (.not. bas_numcont(basis, nsh)) call errquit( 385 $ 'jan_transform: bas_numcont', basis, BASIS_ERR) 386 if (.not. bas_nbf_cn_max(basis,maxbfsh)) call errquit( 387 $ 'jan_transform: bas_nbf_cn_max', basis, BASIS_ERR) 388c 389 lenhalf = n3*n4*maxbfsh**2 390 lenthird= n2*n3*n4*maxbfsh 391c 392 if (.not. ma_push_get(mt_dbl,lenhalf,'half',l_half, k_half)) 393 $ call errquit('ma half', lenhalf, MA_ERR) 394 if (.not. ma_push_get(mt_dbl,lenthird,'third',l_third, k_third)) 395 $ call errquit('ma third', lenthird, MA_ERR) 396c 397 call dfill(n1*n2*n3*n4, 0.0d0, full, 1) 398 do ksh = 1, nsh 399 if (.not. bas_cn2bfr(basis, ksh, klo, khi)) 400 $ call errquit('jan_transform: bas_cn2bfr',basis, BASIS_ERR) 401 call dfill(n2*n3*n4*(khi-klo+1), 0.0d0, dbl_mb(k_third), 1) 402 do lsh = 1, nsh 403 if (.not. bas_cn2bfr(basis, lsh, llo, lhi)) 404 $ call errquit('jan_transform: bas_cn2bfr',basis, 405 & BASIS_ERR) 406 if (schwarz_shell(ksh,lsh)*schwarz_max() 407 $ .gt. tol2e) then 408c 409c Make (rs|kl) all rs (indices 3 and 4) given shells k and l 410c 411 call jan_half_transform(basis, ksh, lsh, n3, n4, 412 $ c3t, c4t, ld3, ld4, 413 $ dbl_mb(k_half), ochargecloud, tol2e) 414 415* write(6,*) ' ksh, lsh ', ksh, lsh 416* call jan_debug_print('half', 417* $ dbl_mb(k_half), n3, n4, khi-klo+1, lhi-llo+1) 418 419c 420 call jan_third_transform(llo, lhi, klo, khi, 421 $ n2, n3, n4, dbl_mb(k_half), dbl_mb(k_third), 422 $ c2t, ld2, tol2e) 423 end if 424 end do 425* write(6,*) ' lsh ', lsh 426* call jan_debug_print('third', 427* $ dbl_mb(k_third), n2, n3, n4, lhi-llo+1) 428 call jan_final_transform(klo, khi, n1, n2, n3, n4, 429 $ dbl_mb(k_third), full, c1t, ld1, tol2e) 430 end do 431c 432 if (oasym) call jan_asym_trans(full,n1,n2,n3,n4,side) 433c 434 do s = 1, n4 435 do r = 1, n3 436 do q = 1, n2 437 do p = 1, n1 438 if (abs(full(p,q,r,s)).lt.1d-10) 439 $ full(p,q,r,s) = 0.0d0 440 end do 441 end do 442 end do 443 end do 444c 445 if (.not. ma_pop_stack(l_third)) call errquit('ma third',0, 446 & MA_ERR) 447 if (.not. ma_pop_stack(l_half)) call errquit('ma half',0, 448 & MA_ERR) 449c 450 call schwarz_tidy() 451 call int_terminate 452c 453 end 454 subroutine jan_full_transform_noinit( 455 $ rtdb, basis, 456 $ n1, n2, n3, n4, 457 $ ld1, ld2, ld3, ld4, 458 $ c1t, c2t, c3t, c4t, 459 $ full, order) 460 implicit none 461#include "errquit.fh" 462#include "schwarz.fh" 463#include "bas.fh" 464#include "mafdecls.fh" 465#include "inp.fh" 466c 467 integer rtdb 468 integer basis ! AO basis handle 469 integer n1, n2, n3, n4 ! Dimension of each MO set 470 integer ld1, ld2, ld3, ld4 471 double precision c1t(ld1,*), c2t(ld2,*), ! Transposed MO coeffs 472 $ c3t(ld3,*), c4t(ld4,*) 473 double precision full(n1,n2,n3,n4) 474 character*(*) order 475c 476c Generate the specified block of MO integrals with 477c no assumptions of equivalence between the sets of coefficients. 478c 479c Order can be either 480c . ChargeCloud -> full(p,q,r,s) = (pq|rs) 481c or 482c . Dirac -> full(p,q,r,s) = <pq|rs> 483c or 484c . LeftAsymDirac -> full(p,q,r,s) = <pq|rs>-<qp|rs> 485c . (must have c1t=c2t, n1=n2) 486c or 487c . RightAsymDirac -> full(p,q,r,s) = <pq|rs>-<pq|sr> 488c . (must have c3t=c4t, n3=n4) 489c 490c Presently the antisymmetrization is done at the top level 491c and the storage of full is not reduced to use the symmetry. 492c 493c Memory requirements are 494c . n1*n2*n3*n4 + S*n2*n3*n4 + S*S*n3*n4 + S*S*S*n4 + 495c . maxs2 + maxg2*(1 + 4*integer) 496c 497c Index 4 is the first transformed so there is advantage in 498c making it the smallest range. 499c 500 double precision tol2e 501 parameter (tol2e = 1d-12) 502C 503 integer nsh, maxbfsh, lenhalf, lenthird, geom 504 integer l_half, k_half, l_third, k_third 505 integer lsh, ksh, llo, lhi, klo, khi 506 integer p, q, r, s 507 logical ochargecloud, oasym 508 character*8 side 509c 510 ochargecloud = .false. ! avoids compiler warning 511 oasym = .false. ! avoids compiler warning 512 if (inp_compare(.false.,order,'chargecloud')) then 513 ochargecloud = .true. 514 oasym = .false. 515 side = ' ' 516 else if (inp_compare(.false.,order,'dirac')) then 517 ochargecloud = .false. 518 oasym = .false. 519 side = ' ' 520 else if (inp_compare(.false.,order,'leftasymdirac')) then 521 ochargecloud = .false. 522 oasym = .true. 523 side = 'left' 524 else if (inp_compare(.false.,order,'rightasymdirac')) then 525 ochargecloud = .false. 526 oasym = .true. 527 side = 'right' 528 else 529 call errquit('jan_full_trans: unkown integral option',0, 530 & INT_ERR) 531 endif 532c 533c Initialize integrals and Schwarz screening 534c 535 if (.not. bas_geom(basis, geom)) 536 $ call errquit('jan_transform: basis ', basis, BASIS_ERR) 537* call int_init(rtdb, 1, basis) 538* call schwarz_init(geom, basis) 539c 540 if (.not. bas_numcont(basis, nsh)) call errquit( 541 $ 'jan_transform: bas_numcont', basis, BASIS_ERR) 542 if (.not. bas_nbf_cn_max(basis,maxbfsh)) call errquit( 543 $ 'jan_transform: bas_nbf_cn_max', basis, BASIS_ERR) 544c 545 lenhalf = n3*n4*maxbfsh**2 546 lenthird= n2*n3*n4*maxbfsh 547c 548 if (.not. ma_push_get(mt_dbl,lenhalf,'half',l_half, k_half)) 549 $ call errquit('ma half', lenhalf, MA_ERR) 550 if (.not. ma_push_get(mt_dbl,lenthird,'third',l_third, k_third)) 551 $ call errquit('ma third', lenthird, MA_ERR) 552c 553 call dfill(n1*n2*n3*n4, 0.0d0, full, 1) 554 do ksh = 1, nsh 555 if (.not. bas_cn2bfr(basis, ksh, klo, khi)) 556 $ call errquit('jan_transform: bas_cn2bfr',basis, BASIS_ERR) 557 call dfill(n2*n3*n4*(khi-klo+1), 0.0d0, dbl_mb(k_third), 1) 558 do lsh = 1, nsh 559 if (.not. bas_cn2bfr(basis, lsh, llo, lhi)) 560 $ call errquit('jan_transform: bas_cn2bfr',basis, 561 & BASIS_ERR) 562 if (schwarz_shell(ksh,lsh)*schwarz_max() 563 $ .gt. tol2e) then 564c 565c Make (rs|kl) all rs (indices 3 and 4) given shells k and l 566c 567 call jan_half_transform(basis, ksh, lsh, n3, n4, 568 $ c3t, c4t, ld3, ld4, 569 $ dbl_mb(k_half), ochargecloud, tol2e) 570 571* write(6,*) ' ksh, lsh ', ksh, lsh 572* call jan_debug_print('half', 573* $ dbl_mb(k_half), n3, n4, khi-klo+1, lhi-llo+1) 574 575c 576 call jan_third_transform(llo, lhi, klo, khi, 577 $ n2, n3, n4, dbl_mb(k_half), dbl_mb(k_third), 578 $ c2t, ld2, tol2e) 579 end if 580 end do 581* write(6,*) ' lsh ', lsh 582* call jan_debug_print('third', 583* $ dbl_mb(k_third), n2, n3, n4, lhi-llo+1) 584 call jan_final_transform(klo, khi, n1, n2, n3, n4, 585 $ dbl_mb(k_third), full, c1t, ld1, tol2e) 586 end do 587c 588 if (oasym) call jan_asym_trans(full,n1,n2,n3,n4,side) 589c 590 do s = 1, n4 591 do r = 1, n3 592 do q = 1, n2 593 do p = 1, n1 594 if (abs(full(p,q,r,s)).lt.1d-10) 595 $ full(p,q,r,s) = 0.0d0 596 end do 597 end do 598 end do 599 end do 600c 601 if (.not. ma_pop_stack(l_third)) call errquit('ma third',0, 602 & MA_ERR) 603 if (.not. ma_pop_stack(l_half)) call errquit('ma half',0, 604 & MA_ERR) 605c 606* call schwarz_tidy() 607* call int_terminate 608c 609 end 610 subroutine jan_asym_trans(full,n1,n2,n3,n4,side) 611 implicit none 612#include "errquit.fh" 613c 614 integer n1, n2, n3, n4 615 double precision full(n1,n2,n3,n4) 616 character*(*) side 617c 618 integer p, q, r, s 619 double precision tmp 620c 621 if (side .eq. 'left') then 622 if (n1 .ne. n2) call errquit('jan_asym_trans: left', n1, 623 & UNKNOWN_ERR) 624 do s = 1, n4 625 do r = 1, n3 626 do q = 1, n2 627 do p = 1, q 628 tmp = full(p,q,r,s) - full(q,p,r,s) 629 full(p,q,r,s) = tmp 630 full(q,p,r,s) =-tmp 631 enddo 632 enddo 633 enddo 634 enddo 635 else 636 if (n3 .ne. n4) call errquit('jan_asym_trans: right', n3, 637 & UNKNOWN_ERR) 638 do s = 1, n4 639 do r = 1, s 640 do q = 1, n2 641 do p = 1, n1 642 full(p,q,r,s) = full(p,q,r,s) - full(p,q,s,r) 643 enddo 644 enddo 645 if (r .ne. s) then 646 do q = 1, n2 647 do p = 1, n1 648 full(p,q,s,r) = -full(p,q,r,s) 649 enddo 650 enddo 651 else 652 do q = 1, n2 653 do p = 1, n1 654 full(p,q,s,r) = 0.0d0 655 enddo 656 enddo 657 endif 658 enddo 659 enddo 660 endif 661c 662 end 663 subroutine jan_final_transform(klo, khi, 664 $ n1, n2, n3, n4, third, full, c1t, ld1, tol2e) 665 implicit none 666c 667 integer klo, khi, n1, n2, n3, n4, ld1 668 double precision third(n2,n3,n4,klo:khi) 669 double precision full(n1,n2,n3,n4) 670 double precision c1t(ld1,*) 671 double precision tol2e 672c 673 integer k, s, r, q, p 674 double precision g 675c 676 do k = klo, khi 677 do s = 1, n4 678 do r = 1, n3 679 do q = 1, n2 680 g = third(q,r,s,k) 681 if (abs(g) .gt. tol2e) then 682 do p = 1, n1 683 full(p,q,r,s) = full(p,q,r,s) + g*c1t(p,k) 684 end do 685 end if 686 end do 687 end do 688 end do 689 end do 690c 691 end 692 subroutine jan_third_transform(llo, lhi, klo, khi, 693 $ n2, n3, n4, half, third, c2t, ld2, tol2e) 694 implicit none 695c 696 integer llo, lhi, klo, khi, n2, n3, n4, ld2 697 double precision half(n3,n4,klo:khi,llo:lhi) 698 double precision third(n2,n3,n4,klo:khi) 699 double precision c2t(ld2,*) 700 double precision tol2e 701c 702 integer k, l, s, r, q 703 double precision g 704c 705 do l = llo, lhi 706 do k = klo, khi 707 do s = 1, n4 708 do r = 1, n3 709 g = half(r,s,k,l) 710 if (abs(g) .gt. tol2e) then 711 do q = 1, n2 712 third(q,r,s,k) = third(q,r,s,k) + g*c2t(q,l) 713 end do 714 end if 715 end do 716 end do 717 end do 718 end do 719c 720 end 721 subroutine jan_half_transform(basis, ksh, lsh, n1, n2, 722 $ c1t, c2t, ld1, ld2, half, ochargecloud, tol2e) 723 implicit none 724#include "errquit.fh" 725#include "bas.fh" 726#include "mafdecls.fh" 727 integer basis, ksh, lsh, n1, n2, ld1, ld2 728 double precision c1t(*), c2t(*) ! Transposed MO coeffs 729 double precision half(*) ! n1*n2*kdim*ldim (pq|kl) 730 logical ochargecloud 731 double precision tol2e 732c 733c For a pair of shells k and l fill 734c 735c . half(p,q,k,l) = (pq|kl) 736c 737c for all p, q, and k, l within their respective shells 738c where p and q are transformed into the new bases and 739c k and l are AO indices. 740c 741c Eventually exploiting sparsity and abelian symmetry 742c ... should also eventually use the texas integrals 743c 744c Assumes that integrals and schwarz have been initialized. 745c 746 integer nbf, nsh, k_i, k_j, k_k, k_l, l_i, l_j, l_k, l_l, 747 $ maxg2, maxs2, k_buf, l_buf, k_scr, l_scr, maxbfsh 748 integer llo, lhi, klo, khi, lenaobuf, l_aobuf, k_aobuf 749c 750c Get dimensions and required scratch space info 751c 752 if (.not. bas_numcont(basis, nsh)) call errquit( 753 $ 'jan_transform: bas_numcont', basis, BASIS_ERR) 754 if (.not. bas_nbf_cn_max(basis,maxbfsh)) call errquit( 755 $ 'jan_transform: bas_nbf_cn_max', basis, BASIS_ERR) 756 if (.not. bas_numbf(basis, nbf)) 757 $ call errquit('jan_transform: nbf',basis, BASIS_ERR) 758 if (.not. bas_cn2bfr(basis, ksh, klo, khi)) 759 $ call errquit('jan_transform: bas_cn2bfr',basis, BASIS_ERR) 760 if (.not. bas_cn2bfr(basis, lsh, llo, lhi)) 761 $ call errquit('jan_transform: bas_cn2bfr',basis, BASIS_ERR) 762 call int_mem_2e4c(maxg2,maxs2) 763 lenaobuf = (khi-klo+1)*(lhi-llo+1)*maxbfsh*n2 ! (iq|kl) 764c 765c Allocate scratch space for integrals and buffers for 766c transformation 767c 768 if (.not. ma_push_get(mt_dbl,maxs2,'scr',l_scr, k_scr)) 769 $ call errquit('ma scr',maxs2, MA_ERR) 770 if (.not. ma_push_get(mt_dbl,maxg2,'buf',l_buf, k_buf)) 771 $ call errquit('ma buf',maxg2, MA_ERR) 772 if (.not. ma_push_get(mt_int,maxg2,'i',l_i, k_i)) 773 $ call errquit('ma ibuf',maxg2, MA_ERR) 774 if (.not. ma_push_get(mt_int,maxg2,'j',l_j, k_j)) 775 $ call errquit('ma jbuf',maxg2, MA_ERR) 776 if (.not. ma_push_get(mt_int,maxg2,'k',l_k, k_k)) 777 $ call errquit('ma kbuf',maxg2, MA_ERR) 778 if (.not. ma_push_get(mt_int,maxg2,'l',l_l, k_l)) 779 $ call errquit('ma lbuf',maxg2, MA_ERR) 780 if (.not. ma_push_get(mt_dbl,lenaobuf,'aobuf',l_aobuf, k_aobuf)) 781 $ call errquit('ma aobuf',lenaobuf, MA_ERR) 782c 783 call jan_do_half_transform( 784 $ basis, 785 $ dbl_mb(k_buf), dbl_mb(k_scr), 786 $ int_mb(k_i), int_mb(k_j), int_mb(k_k), int_mb(k_l), 787 $ maxg2, maxs2, 788 $ nbf, nsh, maxbfsh, n1, n2, tol2e, 789 $ ksh, lsh, klo, khi, llo, lhi, 790 $ half, c1t, c2t, ld1, ld2, dbl_mb(k_aobuf), ochargecloud) 791c 792 if (.not. ma_chop_stack(l_scr)) call errquit 793 $ ('jan_transform: chopping stack', 0, MA_ERR) 794c 795 end 796 subroutine jan_do_half_transform( 797 $ basis, 798 $ buf, scr, ilab, jlab, klab, llab, maxg2, maxs2, 799 $ nbf, nsh, maxbfsh, n1, n2, tol2e, 800 $ ksh, lsh, klo, khi, llo, lhi, 801 $ half, c1t, c2t, ld1, ld2, aobuf, ochargecloud) 802 implicit none 803#include "errquit.fh" 804#include "schwarz.fh" 805#include "bas.fh" 806 integer basis 807 integer maxg2, maxs2 808 integer nbf, nsh, n1, n2, maxbfsh, ld1, ld2 809 double precision buf(maxg2), scr(maxs2) 810 double precision c1t(ld1,nbf), c2t(ld2,nbf) 811 integer ilab(maxg2), jlab(maxg2), klab(maxg2), llab(maxg2) 812 integer ksh, lsh, klo, khi, llo, lhi 813 double precision aobuf(n2,maxbfsh,klo:khi,llo:lhi) 814 double precision half(n1,n2,klo:khi,llo:lhi) 815 double precision tol2e 816 logical ochargecloud 817c 818c For a pair of shells k and l, fill aobuf with integrals (ij|kl) 819c for all i>=j ... eventually exploiting sparsity and abelian symmetry 820c ... should also use the texas integrals 821c 822 double precision skl, g 823 integer ish, jsh, i, j, k, l, p, q, ijkl, nint, ilo, ihi, 824 $ idim, kdim, ldim 825c 826 kdim = khi - klo + 1 827 ldim = lhi - llo + 1 828 skl = schwarz_shell(ksh,lsh) 829 call dfill(n1*n2*kdim*ldim,0.0d0,half,1) 830c 831 do ish = 1, nsh 832 if (.not. bas_cn2bfr(basis,ish,ilo,ihi)) 833 $ call errquit('jan_do_half_transform',ish, BASIS_ERR) 834 idim = ihi-ilo+1 835 do l = llo,lhi 836 do k = klo, khi 837 do i = 1, idim 838 do q = 1, n2 839 aobuf(q,i,k,l) = 0.0d0 840 end do 841 end do 842 end do 843 end do 844c 845 do jsh = 1, nsh 846 if (ochargecloud) then ! (ij|kl) 847 if (schwarz_shell(ish,jsh)*skl .gt. tol2e) then 848 call int_l2e4c(basis, ish, jsh, basis, ksh, lsh, 849 & tol2e, .false., maxg2, buf, nint, 850 $ ilab, jlab, klab, llab, maxs2, scr) 851 do ijkl = 1, nint 852 i = ilab(ijkl)-ilo+1 853 j = jlab(ijkl) 854 k = klab(ijkl) 855 l = llab(ijkl) 856 g = buf(ijkl) 857 if (abs(g) .gt. tol2e) then 858 do q = 1, n2 859 aobuf(q,i,k,l) = aobuf(q,i,k,l) + g*c2t(q,j) 860 end do 861 end if 862 end do 863 end if 864 else ! <ij|kl> = (ik|jl) 865 if (schwarz_shell(ish,ksh)*schwarz_shell(jsh,lsh) 866 $ .gt. tol2e) then 867 call int_l2e4c(basis, ish, ksh, basis, jsh, lsh, 868 & tol2e, .false., maxg2, buf, nint, 869 $ ilab, klab, jlab, llab, maxs2, scr) 870 do ijkl = 1, nint 871 i = ilab(ijkl)-ilo+1 872 j = jlab(ijkl) 873 k = klab(ijkl) 874 l = llab(ijkl) 875 g = buf(ijkl) 876 if (abs(g) .gt. tol2e) then 877 do q = 1, n2 878 aobuf(q,i,k,l) = aobuf(q,i,k,l) + g*c2t(q,j) 879 end do 880 end if 881 end do 882 end if 883 endif 884 end do 885 do l = llo, lhi 886 do k = klo, khi 887 do i = ilo, ihi 888 do q = 1, n2 889 g = aobuf(q,i-ilo+1,k,l) 890 if (abs(g) .gt. tol2e) then 891 do p = 1, n1 892 half(p,q,k,l) = half(p,q,k,l) + g*c1t(p,i) 893 end do 894 end if 895 end do 896 end do 897 end do 898 end do 899 end do 900c 901 end 902 subroutine jan_debug_print(string,full, n1, n2, n3, n4) 903 implicit none 904 character*(*) string 905 integer n1, n2, n3, n4 906 double precision full(n1, n2, n3, n4) 907c 908 integer p, q, r, s 909 write(6,*) ' DEBUG FOR ', string, n1, n2, n3, n4 910 do s = 1, n4 911 do r = 1, n3 912 do q = 1, n2 913 do p = 1, n1 914 if (abs(full(p,q,r,s)).gt.1e-6) then 915 write(6,1) p,q,r,s,full(p,q,r,s) 916 1 format(4i5,2x,f14.8) 917 end if 918 end do 919 end do 920 end do 921 end do 922c 923 end 924 subroutine ccsd_incore(rtdb, basis, g_mo, e, r2, t2, fock, 925 & t1, r1, nbf, nmo, nocc, nso, noso) 926 implicit none 927 integer rtdb, basis, nbf, nmo, nocc, nso, noso 928c 929c nbf = number of basis functions 930c nmo = number of molecular orbitals 931c nocc = number of occupied orbitals 932c nso = number of spin orbitals (for now set to 2*nmo) 933c noso = number of occupied spin orbitals (for now set to 2*nocc) 934c 935c toy coupled cluster program based on the green book 936c 937 double precision g_mo(nmo,nmo,nmo,nmo) 938 double precision e(nmo), fock(nmo,nmo) 939 double precision r2(nso,nso,nso,nso) 940 double precision t2(nso,nso,nso,nso) 941 double precision t1(nso,nso) 942 double precision r1(nso,nso) 943 integer n 944c 945c cc stuff (greek indices denote occupied spin-orbitals) 946c cc stuff (r, s, t, u, ... denote unoccupied spin-orbitals) 947c 948 write(6,*)' nbf, nmo, nocc, nso, noso: ', 949 & nbf, nmo, nocc, nso, noso 950c 951c Initially fill the Fock matrix with diagonal elements only 952c 953 call dfill(nmo**2, 0.0d0, fock, 1) 954 do n = 1, nmo 955 fock(n,n) = e(n) 956 enddo 957c write(6,*) ' Fock Matrix' 958c call output(fock,1,nmo,1,nmo,nmo,nmo,1) 959c 960c First generate an initial guess for t2 961c 962 call dfill(nso**4, 0.0d0, r2, 1) 963 call dfill(nso**4, 0.0d0, t2, 1) 964 call dfill(nso**2, 0.0d0, r1, 1) 965 call dfill(nso**2, 0.0d0, t1, 1) 966c 967c call t2_init(g_mo, e, t, nmo, nso, noso) 968c write(6,*)' T2 initial guess. ' 969c call writet2(t,nso) 970c call correlation(g_mo,t,nmo,nso,noso) 971c 972c 973c next put this guess into t2 expression keeping all terms linear in t2 974c 975c call t2_l_gb(g_mo, e, t, t2, nmo, nso, noso) 976c do n = 1, 3 977c call dfill(nso**4, 0.0d0, t2, 1) 978c call t2_l_rjh(g_mo, t, t2, fock, nmo, nso, noso) 979c call t2_update_rjh(g_mo, t, t2, fock, nmo, nso, noso) 980c write(6,*)' T2 with linear terms. ' 981c call writet2(t,nso) 982c call correlation(g_mo,t,nmo,nso,noso) 983c enddo 984c 985c next put this into t2 expression keeping all terms 986c 987 do n = 1, 50 988 call dfill(nso**4, 0.0d0, r2, 1) 989 call dfill(nso**2, 0.0d0, r1, 1) 990c call t2_l_q_gb(g_mo, e, t, t2, nmo, nso, noso) 991c call t2_l_rjh(g_mo, t, t2, fock, nmo, nso, noso) 992c 993c Generate delta T1 994c 995 call t1_rjh(g_mo, r1, t2, fock, nmo, nso, noso) 996c 997c Generate delta T2 998c 999 call t2_l_q_rjh(g_mo, t2, r2, fock, nmo, nso, noso) 1000c 1001c Update T1 1002c 1003 call t1_update_rjh(g_mo, t1, r1, fock, nmo, nso, noso) 1004c write(6,*)' T1: ' 1005c call writet1(t1,nso) 1006c 1007c Update T2 1008c 1009 call t2_update_rjh(g_mo, t2, r2, fock, nmo, nso, noso) 1010c write(6,*)' T2 with quadratic terms. ' 1011c call writet2(t2,nso) 1012c 1013c Compute Energy expression 1014c 1015 call correlation(g_mo,t1,t2,nmo,nso,noso) 1016c 1017c Use T1 to transform 1e- and 2e- ints 1018c 1019 call get_new_f_g(rtdb, basis, fock, g_mo, t1, nso, nmo) 1020 enddo 1021c 1022 return 1023 end 1024 subroutine triples_incore(rtdb, basis, g_mo, ener, t2, 1025 & t1, tg, tf, w, v, 1026 & nbf, nmo, nocc, nso, noso) 1027 implicit none 1028 integer rtdb, basis, nbf, nmo, nocc, nso, noso 1029c 1030c nbf = number of basis functions 1031c nmo = number of molecular orbitals 1032c nocc = number of occupied orbitals 1033c nso = number of spin orbitals (for now set to 2*nmo) 1034c noso = number of occupied spin orbitals (for now set to 2*nocc) 1035c 1036 double precision g_mo(nmo,nmo,nmo,nmo) 1037 double precision ener(nmo) 1038 double precision t2(nso,nso,nso,nso) 1039 double precision t1(nso,nso) 1040 double precision tg(nso,nso,nso,noso,noso,noso) 1041 double precision tf(nso,nso,nso,noso,noso,noso) 1042 double precision w(nso,nso,nso,noso,noso,noso) 1043 double precision v(nso,nso,nso,noso,noso,noso) 1044 double precision g, et 1045c occupied orbitals 1046 integer i, j, k, l, m, n 1047 integer i_orb, j_orb, k_orb, l_orb, m_orb, n_orb 1048 integer i_spin, j_spin, k_spin, l_spin, m_spin, n_spin 1049c unoccupied orbitals 1050 integer a, b, c, d, e, f 1051 integer a_orb, b_orb, c_orb, d_orb, e_orb, f_orb 1052 integer a_spin, b_spin, c_spin, d_spin, e_spin, f_spin 1053c 1054 write(6,*)' nbf, nmo, nocc, nso, noso: ', 1055 & nbf, nmo, nocc, nso, noso 1056c 1057 call dfill(noso**3*nso**3, 0.0d0, tg, 1) 1058 call dfill(noso**3*nso**3, 0.0d0, tf, 1) 1059 call dfill(noso**3*nso**3, 0.0d0, w, 1) 1060 call dfill(noso**3*nso**3, 0.0d0, v, 1) 1061c 1062 do a = noso+1, nso 1063 a_orb = (a-1)/2 + 1 1064 a_spin = 0 1065 if ((a-2*a_orb).ne.0)a_spin = 1 1066 do b = noso+1, nso 1067 b_orb = (b-1)/2 + 1 1068 b_spin = 0 1069 if ((b-2*b_orb).ne.0)b_spin = 1 1070 do c = noso+1, nso 1071 c_orb = (c-1)/2 + 1 1072 c_spin = 0 1073 if ((c-2*c_orb).ne.0)c_spin = 1 1074c 1075 do i = 1, noso 1076 i_orb = (i-1)/2 + 1 1077 i_spin = 0 1078 if ((i-2*i_orb).ne.0)i_spin = 1 1079 do j = 1, noso 1080 j_orb = (j-1)/2 + 1 1081 j_spin = 0 1082 if ((j-2*j_orb).ne.0)j_spin = 1 1083 do k = 1, noso 1084 k_orb = (k-1)/2 + 1 1085 k_spin = 0 1086 if ((k-2*k_orb).ne.0)k_spin = 1 1087c 1088 do e = noso+1, nso 1089 e_orb = (e-1)/2 + 1 1090 e_spin = 0 1091 if ((e-2*e_orb).ne.0)e_spin = 1 1092c 1093 g=0.0d0 1094 if (a_spin.eq.e_spin.and. 1095 & b_spin.eq.k_spin)then 1096 g = g_mo(a_orb,b_orb,e_orb,k_orb) 1097 endif 1098 if (a_spin.eq.k_spin.and. 1099 & b_spin.eq.e_spin)then 1100 g = g - g_mo(a_orb,b_orb,k_orb,e_orb) 1101 endif 1102c write(6,*)' g:',g 1103 tg(a,b,c,i,j,k) = tg(a,b,c,i,j,k) + 1104 & g*t2(c,e,i,j) 1105 enddo 1106c 1107 do m = 1, noso 1108 m_orb = (m-1)/2 + 1 1109 m_spin = 0 1110 if ((m-2*m_orb).ne.0)m_spin = 1 1111c 1112 g=0.0d0 1113 if (c_spin.eq.i_spin.and. 1114 & m_spin.eq.j_spin)then 1115 g = g_mo(c_orb,m_orb,i_orb,j_orb) 1116 endif 1117 if (c_spin.eq.j_spin.and. 1118 & m_spin.eq.i_spin)then 1119 g = g - g_mo(c_orb,m_orb,j_orb,i_orb) 1120 endif 1121c write(6,*)' g:',g 1122 tg(a,b,c,i,j,k) = tg(a,b,c,i,j,k) - 1123 & g*t2(a,b,m,k) 1124 enddo 1125c 1126 g=0.0d0 1127 if (i_spin.eq.a_spin.and.j_spin.eq.b_spin)then 1128 g = g_mo(i_orb,j_orb,a_orb,b_orb) 1129 endif 1130 if (i_spin.eq.b_spin.and.j_spin.eq.a_spin)then 1131 g = g - g_mo(i_orb,j_orb,b_orb,a_orb) 1132 endif 1133c write(6,*)' g:',g 1134 tf(a,b,c,i,j,k) = tg(a,b,c,i,j,k) + 1135 & g*t1(c,k) 1136 enddo 1137 enddo 1138 enddo 1139 enddo 1140 enddo 1141 enddo 1142c 1143 do a = noso+1, nso 1144 a_orb = (a-1)/2 + 1 1145 a_spin = 0 1146 if ((a-2*a_orb).ne.0)a_spin = 1 1147 do b = noso+1, nso 1148 b_orb = (b-1)/2 + 1 1149 b_spin = 0 1150 if ((b-2*b_orb).ne.0)b_spin = 1 1151 do c = noso+1, nso 1152 c_orb = (c-1)/2 + 1 1153 c_spin = 0 1154 if ((c-2*c_orb).ne.0)c_spin = 1 1155c 1156 do i = 1, noso 1157 i_orb = (i-1)/2 + 1 1158 i_spin = 0 1159 if ((i-2*i_orb).ne.0)i_spin = 1 1160 do j = 1, noso 1161 j_orb = (j-1)/2 + 1 1162 j_spin = 0 1163 if ((j-2*j_orb).ne.0)j_spin = 1 1164 do k = 1, noso 1165 k_orb = (k-1)/2 + 1 1166 k_spin = 0 1167 if ((k-2*k_orb).ne.0)k_spin = 1 1168c 1169 w(a,b,c,i,j,k) = w(a,b,c,i,j,k) + 1170 & tg(a,b,c,i,j,k) - 1171 & tg(a,b,c,k,j,i) - 1172 & tg(a,b,c,i,k,j) - 1173 & tg(c,b,a,i,j,k) - 1174 & tg(a,c,b,i,j,k) + 1175 & tg(c,b,a,k,j,i) + 1176 & tg(a,c,b,k,j,i) + 1177 & tg(c,b,a,i,k,j) + 1178 & tg(a,c,b,i,k,j) 1179c 1180 v(a,b,c,i,j,k) = v(a,b,c,i,j,k) + 1181 & tf(a,b,c,i,j,k) - 1182 & tf(a,b,c,k,j,i) - 1183 & tf(a,b,c,i,k,j) - 1184 & tf(c,b,a,i,j,k) - 1185 & tf(a,c,b,i,j,k) + 1186 & tf(c,b,a,k,j,i) + 1187 & tf(a,c,b,k,j,i) + 1188 & tf(c,b,a,i,k,j) + 1189 & tf(a,c,b,i,k,j) 1190c 1191 enddo 1192 enddo 1193 enddo 1194 enddo 1195 enddo 1196 enddo 1197c 1198 Et = 0.0d0 1199c 1200 do a = noso+1, nso 1201 a_orb = (a-1)/2 + 1 1202 a_spin = 0 1203 if ((a-2*a_orb).ne.0)a_spin = 1 1204 do b = noso+1, nso 1205 b_orb = (b-1)/2 + 1 1206 b_spin = 0 1207 if ((b-2*b_orb).ne.0)b_spin = 1 1208 do c = noso+1, nso 1209 c_orb = (c-1)/2 + 1 1210 c_spin = 0 1211 if ((c-2*c_orb).ne.0)c_spin = 1 1212c 1213 do i = 1, noso 1214 i_orb = (i-1)/2 + 1 1215 i_spin = 0 1216 if ((i-2*i_orb).ne.0)i_spin = 1 1217 do j = 1, noso 1218 j_orb = (j-1)/2 + 1 1219 j_spin = 0 1220 if ((j-2*j_orb).ne.0)j_spin = 1 1221 do k = 1, noso 1222 k_orb = (k-1)/2 + 1 1223 k_spin = 0 1224 if ((k-2*k_orb).ne.0)k_spin = 1 1225c 1226 Et = Et - 1.0d0/36.0d0*w(a,b,c,i,j,k)* 1227 & v(a,b,c,i,j,k)/ 1228 & ( ener(a_orb) + 1229 & ener(b_orb) + 1230 & ener(c_orb) - 1231 & ener(i_orb) - 1232 & ener(j_orb) - 1233 & ener(k_orb) ) 1234 1235 1236 enddo 1237 enddo 1238 enddo 1239 enddo 1240 enddo 1241 enddo 1242c 1243 write(6,*)' (t) energy = ',et 1244c 1245c do a = noso+1, nso 1246c a_orb = (a-1)/2 + 1 1247c a_spin = 0 1248c if ((a-2*a_orb).ne.0)a_spin = 1 1249c do b = noso+1, nso 1250c b_orb = (b-1)/2 + 1 1251c b_spin = 0 1252c if ((b-2*b_orb).ne.0)b_spin = 1 1253c do c = noso+1, nso 1254c c_orb = (c-1)/2 + 1 1255c c_spin = 0 1256c if ((c-2*c_orb).ne.0)c_spin = 1 1257cc 1258c do i = 1, noso 1259c i_orb = (i-1)/2 + 1 1260c i_spin = 0 1261c if ((i-2*i_orb).ne.0)i_spin = 1 1262c do j = 1, noso 1263c j_orb = (j-1)/2 + 1 1264c j_spin = 0 1265c if ((j-2*j_orb).ne.0)j_spin = 1 1266c do k = 1, noso 1267c k_orb = (k-1)/2 + 1 1268c k_spin = 0 1269c if ((k-2*k_orb).ne.0)k_spin = 1 1270cc 1271c if(dabs(w(a,b,c,i,j,k)).gt.1d-10)then 1272c write(6,*)' i, j, k, a, b, c, w, v: ', 1273c & i, j, k, a, b, c, 1274c & w(a,b,c,i,j,k), 1275c & v(a,b,c,i,j,k) 1276c endif 1277cc 1278c enddo 1279c enddo 1280c enddo 1281c enddo 1282c enddo 1283c enddo 1284c 1285 return 1286 end 1287 subroutine t2_init(g_mo, e, t, nmo, nso, noso) 1288 implicit none 1289 integer nmo, nso, noso 1290 double precision g_mo(nmo,nmo,nmo,nmo) 1291 double precision e(nmo) 1292 double precision t(nso,nso,nso,nso) 1293 double precision g, denom 1294c 1295c cc stuff (greek indices denote occupied spin-orbitals) 1296c cc stuff (r, s, t, u, ... denote unoccupied spin-orbitals) 1297c 1298c occupied orbitals 1299 integer alpha, beta 1300 integer alpha_orb, beta_orb 1301 integer alpha_spin, beta_spin 1302c unoccupied orbitals 1303 integer m, n 1304 integer m_orb, n_orb 1305 integer m_spin, n_spin 1306c 1307c First generate an initial guess for t2 1308c 1309 call dfill(nso**4, 0.0d0, t, 1) 1310 do m = noso+1, nso 1311c get spatial orbital and whether alpha or beta spin 1312 m_orb = (m-1)/2 + 1 1313 m_spin = 0 1314 if ((m-2*m_orb).ne.0)m_spin = 1 1315c write(6,*)' m, m_orb, m_spin: ', 1316c & m, m_orb, m_spin 1317 do n = noso+1, nso 1318 n_orb = (n-1)/2 + 1 1319 n_spin = 0 1320 if ((n-2*n_orb).ne.0)n_spin = 1 1321c write(6,*)' n, n_orb, n_spin: ', 1322c & n, n_orb, n_spin 1323 do alpha = 1, noso 1324 alpha_orb = (alpha-1)/2 + 1 1325 alpha_spin = 0 1326 if ((alpha-2*alpha_orb).ne.0)alpha_spin = 1 1327c write(6,*)' alpha, alpha_orb, alpha_spin: ', 1328c & alpha, alpha_orb, alpha_spin 1329 do beta = 1, noso 1330 beta_orb = (beta-1)/2 + 1 1331 beta_spin = 0 1332 if ((beta-2*beta_orb).ne.0)beta_spin = 1 1333c write(6,*)' beta, beta_orb, beta_spin: ', 1334c & beta, beta_orb, beta_spin 1335c 1336 denom = -e(m_orb)-e(n_orb)+e(alpha_orb)+e(beta_orb) 1337c write(6,*)' denom: ', denom 1338 g=0.0d0 1339 if (m_spin.eq.alpha_spin.and.n_spin.eq.beta_spin)then 1340 g = g_mo(m_orb,n_orb,alpha_orb,beta_orb) 1341 endif 1342 if (m_spin.eq.beta_spin.and.n_spin.eq.alpha_spin)then 1343 g = g - g_mo(m_orb,n_orb,beta_orb,alpha_orb) 1344 endif 1345c write(6,*)' g:',g 1346 t(m,n,alpha,beta) = g/denom 1347c write(6,*)' m,n,alpha,beta,t(m,n,alpha,beta):', 1348c & m,n,alpha,beta,t(m,n,alpha,beta) 1349 enddo 1350 enddo 1351 enddo 1352 enddo 1353 return 1354 end 1355 subroutine t2_l_gb(g_mo, e, t, t2, nmo, nso, noso) 1356 implicit none 1357 integer nmo, nso, noso 1358 double precision g_mo(nmo,nmo,nmo,nmo) 1359 double precision e(nmo) 1360 double precision t(nso,nso,nso,nso) 1361 double precision t2(nso,nso,nso,nso) 1362 double precision g, denom 1363c 1364c cc stuff (greek indices denote occupied spin-orbitals) 1365c cc stuff (r, s, t, u, ... denote unoccupied spin-orbitals) 1366c 1367c occupied orbitals 1368 integer alpha, beta, delta, gamma 1369 integer alpha_orb, beta_orb, delta_orb, gamma_orb 1370 integer alpha_spin, beta_spin, delta_spin, gamma_spin 1371c unoccupied orbitals 1372 integer m, n, p, q 1373 integer m_orb, n_orb, p_orb, q_orb 1374 integer m_spin, n_spin, p_spin, q_spin 1375c 1376c put t into t2 expression keeping all terms linear in t2 1377c 1378 do m = noso+1, nso 1379c get spatial orbital and whether alpha or beta spin 1380 m_orb = (m-1)/2 + 1 1381 m_spin = 0 1382 if ((m-2*m_orb).ne.0)m_spin = 1 1383 do n = noso+1, nso 1384 n_orb = (n-1)/2 + 1 1385 n_spin = 0 1386 if ((n-2*n_orb).ne.0)n_spin = 1 1387 do alpha = 1, noso 1388 alpha_orb = (alpha-1)/2 + 1 1389 alpha_spin = 0 1390 if ((alpha-2*alpha_orb).ne.0)alpha_spin = 1 1391 do beta = 1, noso 1392 beta_orb = (beta-1)/2 + 1 1393 beta_spin = 0 1394 if ((beta-2*beta_orb).ne.0)beta_spin = 1 1395 denom = e(m_orb)+e(n_orb)-e(alpha_orb)-e(beta_orb) 1396c write(6,*)' denom: ', denom 1397 g=0.0d0 1398 if (m_spin.eq.alpha_spin.and.n_spin.eq.beta_spin)then 1399 g = g_mo(m_orb,n_orb,alpha_orb,beta_orb) 1400 endif 1401 if (m_spin.eq.beta_spin.and.n_spin.eq.alpha_spin)then 1402 g = g - g_mo(m_orb,n_orb,beta_orb,alpha_orb) 1403 endif 1404c write(6,*)' g:',g 1405 t2(m,n,alpha,beta) = g 1406c 1407 do p = noso+1, nso 1408 p_orb = (p-1)/2 + 1 1409 p_spin = 0 1410 if ((p-2*p_orb).ne.0)p_spin = 1 1411 do q = noso+1, p-1 1412 q_orb = (q-1)/2 + 1 1413 q_spin = 0 1414 if ((q-2*q_orb).ne.0)q_spin = 1 1415c 1416 g=0.0d0 1417 if (m_spin.eq.p_spin.and.n_spin.eq.q_spin)then 1418 g = g_mo(m_orb,n_orb,p_orb,q_orb) 1419 endif 1420 if (m_spin.eq.q_spin.and.n_spin.eq.p_spin)then 1421 g = g - g_mo(m_orb,n_orb,q_orb,p_orb) 1422 endif 1423c write(6,*)' g:',g 1424 t2(m,n,alpha,beta) = t2(m,n,alpha,beta) - 1425 & g*t(p,q,alpha,beta) 1426c 1427 enddo 1428 enddo 1429c 1430 do gamma = 1, noso 1431 gamma_orb = (gamma-1)/2 + 1 1432 gamma_spin = 0 1433 if ((gamma-2*gamma_orb).ne.0)gamma_spin = 1 1434 do delta = 1, gamma-1 1435 delta_orb = (delta-1)/2 + 1 1436 delta_spin = 0 1437 if ((delta-2*delta_orb).ne.0)delta_spin = 1 1438c 1439 g=0.0d0 1440 if (gamma_spin.eq.alpha_spin.and. 1441 & delta_spin.eq.beta_spin)then 1442 g = 1443 & g_mo(gamma_orb,delta_orb,alpha_orb,beta_orb) 1444 endif 1445 if (gamma_spin.eq.beta_spin.and. 1446 & delta_spin.eq.alpha_spin)then 1447 g = g - 1448 & g_mo(gamma_orb,delta_orb,beta_orb,alpha_orb) 1449 endif 1450c write(6,*)' g:',g 1451 t2(m,n,alpha,beta) = t2(m,n,alpha,beta) - 1452 & g*t(m,n,gamma,delta) 1453c 1454 enddo 1455 enddo 1456c 1457 do p = noso+1, nso 1458 p_orb = (p-1)/2 + 1 1459 p_spin = 0 1460 if ((p-2*p_orb).ne.0)p_spin = 1 1461 do gamma = 1, noso 1462 gamma_orb = (gamma-1)/2 + 1 1463 gamma_spin = 0 1464 if ((gamma-2*gamma_orb).ne.0)gamma_spin = 1 1465c 1466 g=0.0d0 1467 if (gamma_spin.eq.beta_spin.and. 1468 & n_spin.eq.p_spin)then 1469 g = g_mo(gamma_orb,n_orb,beta_orb,p_orb) 1470 endif 1471 if (gamma_spin.eq.p_spin.and. 1472 & n_spin.eq.beta_spin)then 1473 g = g - g_mo(gamma_orb,n_orb,p_orb,beta_orb) 1474 endif 1475c write(6,*)' g:',g 1476 t2(m,n,alpha,beta) = t2(m,n,alpha,beta) + 1477 & g*t(m,p,alpha,gamma) 1478c 1479 g=0.0d0 1480 if (gamma_spin.eq.beta_spin.and. 1481 & m_spin.eq.p_spin)then 1482 g = g_mo(gamma_orb,m_orb,beta_orb,p_orb) 1483 endif 1484 if (gamma_spin.eq.p_spin.and. 1485 & m_spin.eq.beta_spin)then 1486 g = g - g_mo(gamma_orb,m_orb,p_orb,beta_orb) 1487 endif 1488c write(6,*)' g:',g 1489 t2(m,n,alpha,beta) = t2(m,n,alpha,beta) - 1490 & g*t(n,p,alpha,gamma) 1491c 1492 g=0.0d0 1493 if (gamma_spin.eq.alpha_spin.and. 1494 & n_spin.eq.p_spin)then 1495 g = g_mo(gamma_orb,n_orb,alpha_orb,p_orb) 1496 endif 1497 if (gamma_spin.eq.p_spin.and. 1498 & n_spin.eq.alpha_spin)then 1499 g = g - g_mo(gamma_orb,n_orb,p_orb,alpha_orb) 1500 endif 1501c write(6,*)' g:',g 1502 t2(m,n,alpha,beta) = t2(m,n,alpha,beta) - 1503 & g*t(m,p,beta,gamma) 1504c 1505 g=0.0d0 1506 if (gamma_spin.eq.alpha_spin.and. 1507 & m_spin.eq.p_spin)then 1508 g = g_mo(gamma_orb,m_orb,alpha_orb,p_orb) 1509 endif 1510 if (gamma_spin.eq.p_spin.and. 1511 & m_spin.eq.alpha_spin)then 1512 g = g - g_mo(gamma_orb,m_orb,p_orb,alpha_orb) 1513 endif 1514c write(6,*)' g:',g 1515 t2(m,n,alpha,beta) = t2(m,n,alpha,beta) + 1516 & g*t(n,p,beta,gamma) 1517 enddo 1518 enddo 1519 t2(m,n,alpha,beta) = t2(m,n,alpha,beta)/denom 1520c write(6,*)' t2(m,n,alpha,beta):', 1521c & t2(m,n,alpha,beta) 1522 enddo 1523 enddo 1524 enddo 1525 enddo 1526c write(6,*)' T2 with linear terms. ' 1527c call writet2(t2,nso) 1528c call correlation(g_mo,t2,nmo,nso,noso) 1529 return 1530 end 1531 subroutine t2_l_q_gb(g_mo, e, t, t2, nmo, nso, noso) 1532 implicit none 1533 integer nmo, nso, noso 1534 double precision g_mo(nmo,nmo,nmo,nmo) 1535 double precision e(nmo) 1536 double precision t(nso,nso,nso,nso) 1537 double precision t2(nso,nso,nso,nso) 1538 double precision g, denom 1539c 1540c cc stuff (greek indices denote occupied spin-orbitals) 1541c cc stuff (r, s, t, u, ... denote unoccupied spin-orbitals) 1542c 1543c occupied orbitals 1544 integer alpha, beta, delta, gamma 1545 integer alpha_orb, beta_orb, delta_orb, gamma_orb 1546 integer alpha_spin, beta_spin, delta_spin, gamma_spin 1547c unoccupied orbitals 1548 integer m, n, p, q 1549 integer m_orb, n_orb, p_orb, q_orb 1550 integer m_spin, n_spin, p_spin, q_spin 1551c 1552c put t into t2 expression keeping all terms (linear and quadratic) 1553c 1554 do m = noso+1, nso 1555c get spatial orbital and whether alpha or beta spin 1556 m_orb = (m-1)/2 + 1 1557 m_spin = 0 1558 if ((m-2*m_orb).ne.0)m_spin = 1 1559 do n = noso+1, nso 1560 n_orb = (n-1)/2 + 1 1561 n_spin = 0 1562 if ((n-2*n_orb).ne.0)n_spin = 1 1563 do alpha = 1, noso 1564 alpha_orb = (alpha-1)/2 + 1 1565 alpha_spin = 0 1566 if ((alpha-2*alpha_orb).ne.0)alpha_spin = 1 1567 do beta = 1, noso 1568 beta_orb = (beta-1)/2 + 1 1569 beta_spin = 0 1570 if ((beta-2*beta_orb).ne.0)beta_spin = 1 1571 denom = e(m_orb)+e(n_orb)-e(alpha_orb)-e(beta_orb) 1572c write(6,*)' denom: ', denom 1573 g=0.0d0 1574 if (m_spin.eq.alpha_spin.and.n_spin.eq.beta_spin)then 1575 g = g_mo(m_orb,n_orb,alpha_orb,beta_orb) 1576 endif 1577 if (m_spin.eq.beta_spin.and.n_spin.eq.alpha_spin)then 1578 g = g - g_mo(m_orb,n_orb,beta_orb,alpha_orb) 1579 endif 1580c write(6,*)' g:',g 1581 t2(m,n,alpha,beta) = g 1582c 1583 do p = noso+1, nso 1584 p_orb = (p-1)/2 + 1 1585 p_spin = 0 1586 if ((p-2*p_orb).ne.0)p_spin = 1 1587 do q = noso+1, p-1 1588 q_orb = (q-1)/2 + 1 1589 q_spin = 0 1590 if ((q-2*q_orb).ne.0)q_spin = 1 1591c 1592 g=0.0d0 1593 if (m_spin.eq.p_spin.and.n_spin.eq.q_spin)then 1594 g = g_mo(m_orb,n_orb,p_orb,q_orb) 1595 endif 1596 if (m_spin.eq.q_spin.and.n_spin.eq.p_spin)then 1597 g = g - g_mo(m_orb,n_orb,q_orb,p_orb) 1598 endif 1599c write(6,*)' g:',g 1600 t2(m,n,alpha,beta) = t2(m,n,alpha,beta) - 1601 & g*t(p,q,alpha,beta) 1602c 1603 enddo 1604 enddo 1605c 1606 do gamma = 1, noso 1607 gamma_orb = (gamma-1)/2 + 1 1608 gamma_spin = 0 1609 if ((gamma-2*gamma_orb).ne.0)gamma_spin = 1 1610 do delta = 1, gamma-1 1611 delta_orb = (delta-1)/2 + 1 1612 delta_spin = 0 1613 if ((delta-2*delta_orb).ne.0)delta_spin = 1 1614c 1615 g=0.0d0 1616 if (gamma_spin.eq.alpha_spin.and. 1617 & delta_spin.eq.beta_spin)then 1618 g = 1619 & g_mo(gamma_orb,delta_orb,alpha_orb,beta_orb) 1620 endif 1621 if (gamma_spin.eq.beta_spin.and. 1622 & delta_spin.eq.alpha_spin)then 1623 g = g - 1624 & g_mo(gamma_orb,delta_orb,beta_orb,alpha_orb) 1625 endif 1626c write(6,*)' g:',g 1627 t2(m,n,alpha,beta) = t2(m,n,alpha,beta) - 1628 & g*t(m,n,gamma,delta) 1629c 1630 enddo 1631 enddo 1632c 1633 do p = noso+1, nso 1634 p_orb = (p-1)/2 + 1 1635 p_spin = 0 1636 if ((p-2*p_orb).ne.0)p_spin = 1 1637 do gamma = 1, noso 1638 gamma_orb = (gamma-1)/2 + 1 1639 gamma_spin = 0 1640 if ((gamma-2*gamma_orb).ne.0)gamma_spin = 1 1641c 1642 g=0.0d0 1643 if (gamma_spin.eq.beta_spin.and. 1644 & n_spin.eq.p_spin)then 1645 g = g_mo(gamma_orb,n_orb,beta_orb,p_orb) 1646 endif 1647 if (gamma_spin.eq.p_spin.and. 1648 & n_spin.eq.beta_spin)then 1649 g = g - g_mo(gamma_orb,n_orb,p_orb,beta_orb) 1650 endif 1651c write(6,*)' g:',g 1652 t2(m,n,alpha,beta) = t2(m,n,alpha,beta) + 1653 & g*t(m,p,alpha,gamma) 1654c 1655 g=0.0d0 1656 if (gamma_spin.eq.beta_spin.and. 1657 & m_spin.eq.p_spin)then 1658 g = g_mo(gamma_orb,m_orb,beta_orb,p_orb) 1659 endif 1660 if (gamma_spin.eq.p_spin.and. 1661 & m_spin.eq.beta_spin)then 1662 g = g - g_mo(gamma_orb,m_orb,p_orb,beta_orb) 1663 endif 1664c write(6,*)' g:',g 1665 t2(m,n,alpha,beta) = t2(m,n,alpha,beta) - 1666 & g*t(n,p,alpha,gamma) 1667c 1668 g=0.0d0 1669 if (gamma_spin.eq.alpha_spin.and. 1670 & n_spin.eq.p_spin)then 1671 g = g_mo(gamma_orb,n_orb,alpha_orb,p_orb) 1672 endif 1673 if (gamma_spin.eq.p_spin.and. 1674 & n_spin.eq.alpha_spin)then 1675 g = g - g_mo(gamma_orb,n_orb,p_orb,alpha_orb) 1676 endif 1677c write(6,*)' g:',g 1678 t2(m,n,alpha,beta) = t2(m,n,alpha,beta) - 1679 & g*t(m,p,beta,gamma) 1680c 1681 g=0.0d0 1682 if (gamma_spin.eq.alpha_spin.and. 1683 & m_spin.eq.p_spin)then 1684 g = g_mo(gamma_orb,m_orb,alpha_orb,p_orb) 1685 endif 1686 if (gamma_spin.eq.p_spin.and. 1687 & m_spin.eq.alpha_spin)then 1688 g = g - g_mo(gamma_orb,m_orb,p_orb,alpha_orb) 1689 endif 1690c write(6,*)' g:',g 1691 t2(m,n,alpha,beta) = t2(m,n,alpha,beta) + 1692 & g*t(n,p,beta,gamma) 1693 enddo 1694 enddo 1695c 1696 do gamma = 1, noso 1697 gamma_orb = (gamma-1)/2 + 1 1698 gamma_spin = 0 1699 if ((gamma-2*gamma_orb).ne.0)gamma_spin = 1 1700 do delta = 1, gamma-1 1701 delta_orb = (delta-1)/2 + 1 1702 delta_spin = 0 1703 if ((delta-2*delta_orb).ne.0)delta_spin = 1 1704 do p = noso+1, nso 1705 p_orb = (p-1)/2 + 1 1706 p_spin = 0 1707 if ((p-2*p_orb).ne.0)p_spin = 1 1708 do q = noso+1, p-1 1709 q_orb = (q-1)/2 + 1 1710 q_spin = 0 1711 if ((q-2*q_orb).ne.0)q_spin = 1 1712c 1713 g=0.0d0 1714 if (gamma_spin.eq.p_spin.and. 1715 & delta_spin.eq.q_spin)then 1716 g = g_mo(gamma_orb,delta_orb,p_orb,q_orb) 1717 endif 1718 if (gamma_spin.eq.q_spin.and. 1719 & delta_spin.eq.p_spin)then 1720 g = g - g_mo(gamma_orb,delta_orb,q_orb,p_orb) 1721 endif 1722c write(6,*)' g:',g 1723c 1724 t2(m,n,alpha,beta) = t2(m,n,alpha,beta) + 1725 & g*(t(p,q,alpha,beta)*t(m,n,gamma,delta) - 1726 & 2.0d0*(t(m,p,alpha,beta)*t(n,q,gamma,delta) + 1727 & t(n,q,alpha,beta)*t(m,p,gamma,delta))- 1728 & 2.0d0*(t(m,n,alpha,gamma)*t(p,q,beta,delta) + 1729 & t(p,q,alpha,gamma)*t(m,n,beta,delta))+ 1730 & 4.0d0*(t(m,p,alpha,gamma)*t(n,q,beta,delta) + 1731 & t(n,q,alpha,gamma)*t(m,p,beta,delta))) 1732 enddo 1733 enddo 1734 enddo 1735 enddo 1736 t2(m,n,alpha,beta) = t2(m,n,alpha,beta)/denom 1737c write(6,*)' t2(m,n,alpha,beta):', 1738c & t2(m,n,alpha,beta) 1739 enddo 1740 enddo 1741 enddo 1742 enddo 1743c write(6,*)' T2 with all terms. ' 1744c call writet2(t2,nso) 1745c call correlation(g_mo,t2,nmo,nso,noso) 1746 return 1747 end 1748 subroutine writet2(t,len) 1749 implicit none 1750 integer len, i, j, k, l 1751 double precision t(len,len,len,len) 1752 do i = 1, len 1753 do j = 1, len 1754 do k = 1, len 1755 do l = 1, len 1756 if (abs(t(i,j,k,l)).gt.1d-6) 1757 & write(6,*)i,j,k,l,t(i,j,k,l) 1758 enddo 1759 enddo 1760 enddo 1761 enddo 1762 return 1763 end 1764 subroutine writet1(t,len) 1765 implicit none 1766 integer len, i, j 1767 double precision t(len,len) 1768 do i = 1, len 1769 do j = 1, len 1770 if (abs(t(i,j)).gt.1d-6) 1771 & write(6,*)i,j,t(i,j) 1772 enddo 1773 enddo 1774 return 1775 end 1776 subroutine correlation(g_mo,t1,t2,nmo,nso,noso) 1777 implicit none 1778 integer nmo, nso, noso 1779 double precision g_mo(nmo,nmo,nmo,nmo) 1780 double precision t2(nso,nso,nso,nso) 1781 double precision t1(nso,nso) 1782 double precision g, e2 1783c occupied orbitals 1784 integer alpha, beta 1785 integer alpha_orb, beta_orb 1786 integer alpha_spin, beta_spin 1787c unoccupied orbitals 1788 integer m, n 1789 integer m_orb, n_orb 1790 integer m_spin, n_spin 1791c 1792 e2 = 0.0d0 1793 do m = noso+1, nso 1794c get spatial orbital and whether alpha or beta spin 1795 m_orb = (m-1)/2 + 1 1796 m_spin = 0 1797 if ((m-2*m_orb).ne.0)m_spin = 1 1798 do n = noso+1, nso 1799 n_orb = (n-1)/2 + 1 1800 n_spin = 0 1801 if ((n-2*n_orb).ne.0)n_spin = 1 1802 do alpha = 1, noso 1803 alpha_orb = (alpha-1)/2 + 1 1804 alpha_spin = 0 1805 if ((alpha-2*alpha_orb).ne.0)alpha_spin = 1 1806 do beta = 1, noso 1807 beta_orb = (beta-1)/2 + 1 1808 beta_spin = 0 1809 if ((beta-2*beta_orb).ne.0)beta_spin = 1 1810c 1811 g=0.0d0 1812 if (m_spin.eq.beta_spin.and.n_spin.eq.alpha_spin)then 1813 g = g_mo(beta_orb,alpha_orb,m_orb,n_orb) 1814 endif 1815 if (m_spin.eq.alpha_spin.and.n_spin.eq.beta_spin)then 1816 g = g - g_mo(beta_orb,alpha_orb,n_orb,m_orb) 1817 endif 1818c write(6,*)' g:',g 1819 e2 = e2 - 0.25d0*g*(t2(m,n,alpha,beta)+ 1820 & t1(m,alpha)*t1(n,beta) - 1821 & t1(n,alpha)*t1(m,beta)) 1822c 1823 enddo 1824 enddo 1825 enddo 1826 enddo 1827 write(6,*)' e2 = ',e2 1828c 1829 return 1830 end 1831 subroutine t2_l_rjh(g_mo, t, t2, fock, nmo, nso, noso) 1832 implicit none 1833 integer nmo, nso, noso 1834 double precision g_mo(nmo,nmo,nmo,nmo) 1835 double precision fock(nmo,nmo) 1836 double precision t(nso,nso,nso,nso) 1837 double precision t2(nso,nso,nso,nso) 1838 double precision g, f_t, denom 1839c occupied orbitals 1840 integer i, j, k, l, m, n 1841 integer i_orb, j_orb, k_orb, l_orb, m_orb, n_orb 1842 integer i_spin, j_spin, k_spin, l_spin, m_spin, n_spin 1843c unoccupied orbitals 1844 integer a, b, c, d, e, f 1845 integer a_orb, b_orb, c_orb, d_orb, e_orb, f_orb 1846 integer a_spin, b_spin, c_spin, d_spin, e_spin, f_spin 1847c 1848c put t into t2 expression keeping all terms linear in t2 1849c 1850 do e = noso+1, nso 1851c get spatial orbital and whether alpha or beta spin 1852 e_orb = (e-1)/2 + 1 1853 e_spin = 0 1854 if ((e-2*e_orb).ne.0)e_spin = 1 1855 do f = noso+1, nso 1856 f_orb = (f-1)/2 + 1 1857 f_spin = 0 1858 if ((f-2*f_orb).ne.0)f_spin = 1 1859 do m = 1, noso 1860 m_orb = (m-1)/2 + 1 1861 m_spin = 0 1862 if ((m-2*m_orb).ne.0)m_spin = 1 1863 do n = 1, noso 1864 n_orb = (n-1)/2 + 1 1865 n_spin = 0 1866 if ((n-2*n_orb).ne.0)n_spin = 1 1867c 1868 g=0.0d0 1869 if (e_spin.eq.m_spin.and.f_spin.eq.n_spin)then 1870 g = g_mo(e_orb,f_orb,m_orb,n_orb) 1871 endif 1872 if (e_spin.eq.n_spin.and.f_spin.eq.m_spin)then 1873 g = g - g_mo(e_orb,f_orb,n_orb,m_orb) 1874 endif 1875c write(6,*)' g:',g 1876 t2(e,f,m,n) = g 1877c 1878 do a = noso+1, nso 1879 a_orb = (a-1)/2 + 1 1880 a_spin = 0 1881 if ((a-2*a_orb).ne.0)a_spin = 1 1882 f_t = 0.0d0 1883 if (f_spin.eq.a_spin)then 1884 f_t = fock(f_orb,a_orb) 1885 endif 1886 t2(e,f,m,n) = t2(e,f,m,n) + 1887 & f_t*t(e,a,m,n) 1888 f_t = 0.0d0 1889 if (e_spin.eq.a_spin)then 1890 f_t = fock(e_orb,a_orb) 1891 endif 1892 t2(e,f,m,n) = t2(e,f,m,n) - 1893 & f_t*t(f,a,m,n) 1894 enddo 1895c 1896 do i = 1, noso 1897 i_orb = (i-1)/2 + 1 1898 i_spin = 0 1899 if ((i-2*i_orb).ne.0)i_spin = 1 1900 f_t = 0.0d0 1901 if (i_spin.eq.n_spin)then 1902 f_t = fock(i_orb,n_orb) 1903 endif 1904 t2(e,f,m,n) = t2(e,f,m,n) - 1905 & f_t*t(e,f,m,i) 1906 f_t = 0.0d0 1907 if (i_spin.eq.m_spin)then 1908 f_t = fock(i_orb,m_orb) 1909 endif 1910 t2(e,f,m,n) = t2(e,f,m,n) + 1911 & f_t*t(e,f,n,i) 1912 enddo 1913c 1914 do a = noso+1, nso 1915 a_orb = (a-1)/2 + 1 1916 a_spin = 0 1917 if ((a-2*a_orb).ne.0)a_spin = 1 1918 do i = 1, noso 1919 i_orb = (i-1)/2 + 1 1920 i_spin = 0 1921 if ((i-2*i_orb).ne.0)i_spin = 1 1922 g=0.0d0 1923 if (i_spin.eq.a_spin.and.f_spin.eq.n_spin)then 1924 g = g_mo(i_orb,f_orb,a_orb,n_orb) 1925 endif 1926 if (i_spin.eq.n_spin.and.a_spin.eq.f_spin)then 1927 g = g - g_mo(i_orb,f_orb,n_orb,a_orb) 1928 endif 1929c write(6,*)' g:',g 1930 t2(e,f,m,n) = t2(e,f,m,n) + 1931 & g*t(e,a,m,i) 1932c ef flip 1933 g=0.0d0 1934 if (i_spin.eq.a_spin.and.e_spin.eq.n_spin)then 1935 g = g_mo(i_orb,e_orb,a_orb,n_orb) 1936 endif 1937 if (i_spin.eq.n_spin.and.a_spin.eq.e_spin)then 1938 g = g - g_mo(i_orb,e_orb,n_orb,a_orb) 1939 endif 1940c write(6,*)' g:',g 1941 t2(e,f,m,n) = t2(e,f,m,n) - 1942 & g*t(f,a,m,i) 1943c mn flip 1944 g=0.0d0 1945 if (i_spin.eq.a_spin.and.f_spin.eq.m_spin)then 1946 g = g_mo(i_orb,f_orb,a_orb,m_orb) 1947 endif 1948 if (i_spin.eq.m_spin.and.a_spin.eq.f_spin)then 1949 g = g - g_mo(i_orb,f_orb,m_orb,a_orb) 1950 endif 1951c write(6,*)' g:',g 1952 t2(e,f,m,n) = t2(e,f,m,n) - 1953 & g*t(e,a,n,i) 1954c ef+mn flip 1955 g=0.0d0 1956 if (i_spin.eq.a_spin.and.e_spin.eq.m_spin)then 1957 g = g_mo(i_orb,e_orb,a_orb,m_orb) 1958 endif 1959 if (i_spin.eq.m_spin.and.a_spin.eq.e_spin)then 1960 g = g - g_mo(i_orb,e_orb,m_orb,a_orb) 1961 endif 1962c write(6,*)' g:',g 1963 t2(e,f,m,n) = t2(e,f,m,n) + 1964 & g*t(f,a,n,i) 1965 enddo 1966 enddo 1967c 1968 do i = 1, noso 1969 i_orb = (i-1)/2 + 1 1970 i_spin = 0 1971 if ((i-2*i_orb).ne.0)i_spin = 1 1972 do j = 1, noso 1973 j_orb = (j-1)/2 + 1 1974 j_spin = 0 1975 if ((j-2*j_orb).ne.0)j_spin = 1 1976 g=0.0d0 1977 if (i_spin.eq.m_spin.and.j_spin.eq.n_spin)then 1978 g = g_mo(i_orb,j_orb,m_orb,n_orb) 1979 endif 1980 if (i_spin.eq.n_spin.and.j_spin.eq.m_spin)then 1981 g = g - g_mo(i_orb,j_orb,n_orb,m_orb) 1982 endif 1983c write(6,*)' g:',g 1984 t2(e,f,m,n) = t2(e,f,m,n) + 1985 & 0.5d0*g*t(e,f,i,j) 1986 enddo 1987 enddo 1988 do a = noso+1, nso 1989 a_orb = (a-1)/2 + 1 1990 a_spin = 0 1991 if ((a-2*a_orb).ne.0)a_spin = 1 1992 do b = noso+1, nso 1993 b_orb = (b-1)/2 + 1 1994 b_spin = 0 1995 if ((b-2*b_orb).ne.0)b_spin = 1 1996 g=0.0d0 1997 if (e_spin.eq.a_spin.and.f_spin.eq.b_spin)then 1998 g = g_mo(e_orb,f_orb,a_orb,b_orb) 1999 endif 2000 if (e_spin.eq.b_spin.and.f_spin.eq.a_spin)then 2001 g = g - g_mo(e_orb,f_orb,b_orb,a_orb) 2002 endif 2003c write(6,*)' g:',g 2004 t2(e,f,m,n) = t2(e,f,m,n) + 2005 & 0.5d0*g*t(a,b,m,n) 2006 enddo 2007 enddo 2008 enddo 2009 enddo 2010 enddo 2011 enddo 2012c write(6,*)' t2(e,f,m,n) L residual:' 2013c call writet2(t2,nso) 2014 return 2015 end 2016 subroutine t2_l_q_rjh(g_mo, t2, r2, fock, nmo, nso, noso) 2017 implicit none 2018 integer nmo, nso, noso 2019 double precision g_mo(nmo,nmo,nmo,nmo) 2020 double precision fock(nmo,nmo) 2021 double precision r2(nso,nso,nso,nso) 2022 double precision t2(nso,nso,nso,nso) 2023 double precision g, f_t, denom 2024c occupied orbitals 2025 integer i, j, k, l, m, n 2026 integer i_orb, j_orb, k_orb, l_orb, m_orb, n_orb 2027 integer i_spin, j_spin, k_spin, l_spin, m_spin, n_spin 2028c unoccupied orbitals 2029 integer a, b, c, d, e, f 2030 integer a_orb, b_orb, c_orb, d_orb, e_orb, f_orb 2031 integer a_spin, b_spin, c_spin, d_spin, e_spin, f_spin 2032c 2033c put t2 into r2 expression keeping all terms quadratic in t2 2034c 2035 do e = noso+1, nso 2036c get spatial orbital and whether alpha or beta spin 2037 e_orb = (e-1)/2 + 1 2038 e_spin = 0 2039 if ((e-2*e_orb).ne.0)e_spin = 1 2040 do f = noso+1, nso 2041 f_orb = (f-1)/2 + 1 2042 f_spin = 0 2043 if ((f-2*f_orb).ne.0)f_spin = 1 2044 do m = 1, noso 2045 m_orb = (m-1)/2 + 1 2046 m_spin = 0 2047 if ((m-2*m_orb).ne.0)m_spin = 1 2048 do n = 1, noso 2049 n_orb = (n-1)/2 + 1 2050 n_spin = 0 2051 if ((n-2*n_orb).ne.0)n_spin = 1 2052c 2053 g=0.0d0 2054 if (e_spin.eq.m_spin.and.f_spin.eq.n_spin)then 2055 g = g_mo(e_orb,f_orb,m_orb,n_orb) 2056 endif 2057 if (e_spin.eq.n_spin.and.f_spin.eq.m_spin)then 2058 g = g - g_mo(e_orb,f_orb,n_orb,m_orb) 2059 endif 2060c write(6,*)' g:',g 2061 r2(e,f,m,n) = g 2062c 2063 do a = noso+1, nso 2064 a_orb = (a-1)/2 + 1 2065 a_spin = 0 2066 if ((a-2*a_orb).ne.0)a_spin = 1 2067 f_t = 0.0d0 2068 if (f_spin.eq.a_spin)then 2069 f_t = fock(f_orb,a_orb) 2070 endif 2071 r2(e,f,m,n) = r2(e,f,m,n) + 2072 & f_t*t2(e,a,m,n) 2073 f_t = 0.0d0 2074 if (e_spin.eq.a_spin)then 2075 f_t = fock(e_orb,a_orb) 2076 endif 2077 r2(e,f,m,n) = r2(e,f,m,n) - 2078 & f_t*t2(f,a,m,n) 2079 enddo 2080c 2081 do i = 1, noso 2082 i_orb = (i-1)/2 + 1 2083 i_spin = 0 2084 if ((i-2*i_orb).ne.0)i_spin = 1 2085 f_t = 0.0d0 2086 if (i_spin.eq.n_spin)then 2087 f_t = fock(i_orb,n_orb) 2088 endif 2089 r2(e,f,m,n) = r2(e,f,m,n) - 2090 & f_t*t2(e,f,m,i) 2091 f_t = 0.0d0 2092 if (i_spin.eq.m_spin)then 2093 f_t = fock(i_orb,m_orb) 2094 endif 2095 r2(e,f,m,n) = r2(e,f,m,n) + 2096 & f_t*t2(e,f,n,i) 2097 enddo 2098c 2099 do a = noso+1, nso 2100 a_orb = (a-1)/2 + 1 2101 a_spin = 0 2102 if ((a-2*a_orb).ne.0)a_spin = 1 2103 do i = 1, noso 2104 i_orb = (i-1)/2 + 1 2105 i_spin = 0 2106 if ((i-2*i_orb).ne.0)i_spin = 1 2107 g=0.0d0 2108 if (i_spin.eq.a_spin.and.f_spin.eq.n_spin)then 2109 g = g_mo(i_orb,f_orb,a_orb,n_orb) 2110 endif 2111 if (i_spin.eq.n_spin.and.a_spin.eq.f_spin)then 2112 g = g - g_mo(i_orb,f_orb,n_orb,a_orb) 2113 endif 2114c write(6,*)' g:',g 2115 r2(e,f,m,n) = r2(e,f,m,n) + 2116 & g*t2(e,a,m,i) 2117c ef flip 2118 g=0.0d0 2119 if (i_spin.eq.a_spin.and.e_spin.eq.n_spin)then 2120 g = g_mo(i_orb,e_orb,a_orb,n_orb) 2121 endif 2122 if (i_spin.eq.n_spin.and.a_spin.eq.e_spin)then 2123 g = g - g_mo(i_orb,e_orb,n_orb,a_orb) 2124 endif 2125c write(6,*)' g:',g 2126 r2(e,f,m,n) = r2(e,f,m,n) - 2127 & g*t2(f,a,m,i) 2128c mn flip 2129 g=0.0d0 2130 if (i_spin.eq.a_spin.and.f_spin.eq.m_spin)then 2131 g = g_mo(i_orb,f_orb,a_orb,m_orb) 2132 endif 2133 if (i_spin.eq.m_spin.and.a_spin.eq.f_spin)then 2134 g = g - g_mo(i_orb,f_orb,m_orb,a_orb) 2135 endif 2136c write(6,*)' g:',g 2137 r2(e,f,m,n) = r2(e,f,m,n) - 2138 & g*t2(e,a,n,i) 2139c ef+mn flip 2140 g=0.0d0 2141 if (i_spin.eq.a_spin.and.e_spin.eq.m_spin)then 2142 g = g_mo(i_orb,e_orb,a_orb,m_orb) 2143 endif 2144 if (i_spin.eq.m_spin.and.a_spin.eq.e_spin)then 2145 g = g - g_mo(i_orb,e_orb,m_orb,a_orb) 2146 endif 2147c write(6,*)' g:',g 2148 r2(e,f,m,n) = r2(e,f,m,n) + 2149 & g*t2(f,a,n,i) 2150 enddo 2151 enddo 2152c 2153 do i = 1, noso 2154 i_orb = (i-1)/2 + 1 2155 i_spin = 0 2156 if ((i-2*i_orb).ne.0)i_spin = 1 2157 do j = 1, noso 2158 j_orb = (j-1)/2 + 1 2159 j_spin = 0 2160 if ((j-2*j_orb).ne.0)j_spin = 1 2161 g=0.0d0 2162 if (i_spin.eq.m_spin.and.j_spin.eq.n_spin)then 2163 g = g_mo(i_orb,j_orb,m_orb,n_orb) 2164 endif 2165 if (i_spin.eq.n_spin.and.j_spin.eq.m_spin)then 2166 g = g - g_mo(i_orb,j_orb,n_orb,m_orb) 2167 endif 2168c write(6,*)' g:',g 2169 r2(e,f,m,n) = r2(e,f,m,n) + 2170 & 0.5d0*g*t2(e,f,i,j) 2171 enddo 2172 enddo 2173 do a = noso+1, nso 2174 a_orb = (a-1)/2 + 1 2175 a_spin = 0 2176 if ((a-2*a_orb).ne.0)a_spin = 1 2177 do b = noso+1, nso 2178 b_orb = (b-1)/2 + 1 2179 b_spin = 0 2180 if ((b-2*b_orb).ne.0)b_spin = 1 2181 g=0.0d0 2182 if (e_spin.eq.a_spin.and.f_spin.eq.b_spin)then 2183 g = g_mo(e_orb,f_orb,a_orb,b_orb) 2184 endif 2185 if (e_spin.eq.b_spin.and.f_spin.eq.a_spin)then 2186 g = g - g_mo(e_orb,f_orb,b_orb,a_orb) 2187 endif 2188c write(6,*)' g:',g 2189 r2(e,f,m,n) = r2(e,f,m,n) + 2190 & 0.5d0*g*t2(a,b,m,n) 2191 enddo 2192 enddo 2193c 2194 do i = 1, noso 2195 i_orb = (i-1)/2 + 1 2196 i_spin = 0 2197 if ((i-2*i_orb).ne.0)i_spin = 1 2198 do j = 1, noso 2199 j_orb = (j-1)/2 + 1 2200 j_spin = 0 2201 if ((j-2*j_orb).ne.0)j_spin = 1 2202 do a = noso+1, nso 2203 a_orb = (a-1)/2 + 1 2204 a_spin = 0 2205 if ((a-2*a_orb).ne.0)a_spin = 1 2206 do b = noso+1, nso 2207 b_orb = (b-1)/2 + 1 2208 b_spin = 0 2209 if ((b-2*b_orb).ne.0)b_spin = 1 2210 g=0.0d0 2211 if (i_spin.eq.a_spin.and. 2212 & j_spin.eq.b_spin)then 2213 g = g_mo(i_orb,j_orb,a_orb,b_orb) 2214 endif 2215 if (i_spin.eq.b_spin.and. 2216 & j_spin.eq.a_spin)then 2217 g = g - g_mo(i_orb,j_orb,b_orb,a_orb) 2218 endif 2219c write(6,*)' g:',g 2220 r2(e,f,m,n) = r2(e,f,m,n) + 2221 & 0.5d0*g* ( t2(f,b,j,n)*t2(e,a,i,m) - 2222 & t2(e,b,j,n)*t2(f,a,i,m) - 2223 & t2(f,b,j,m)*t2(e,a,i,n) + 2224 & t2(e,b,j,m)*t2(f,a,i,n) - 2225 & t2(f,b,m,n)*t2(e,a,i,j) + 2226 & t2(e,b,m,n)*t2(f,a,i,j) - 2227 & t2(e,f,m,i)*t2(a,b,n,j) + 2228 & t2(e,f,n,i)*t2(a,b,m,j) + 2229 & 0.5d0*t2(e,f,i,j)*t2(a,b,m,n) ) 2230 2231 enddo 2232 enddo 2233 enddo 2234 enddo 2235 enddo 2236 enddo 2237 enddo 2238 enddo 2239c write(6,*)' t2(e,f,m,n) L+Q residual:' 2240c call writet2(r2,nso) 2241 return 2242 end 2243 subroutine t2_update_rjh(g_mo, t2, r2, fock, nmo, nso, noso) 2244 implicit none 2245 integer nmo, nso, noso 2246 double precision g_mo(nmo,nmo,nmo,nmo) 2247 double precision fock(nmo,nmo) 2248 double precision t2(nso,nso,nso,nso) 2249 double precision r2(nso,nso,nso,nso) 2250 double precision g, denom, r2norm 2251c occupied orbitals 2252 integer i, j, k, l, m, n 2253 integer i_orb, j_orb, k_orb, l_orb, m_orb, n_orb 2254 integer i_spin, j_spin, k_spin, l_spin, m_spin, n_spin 2255c unoccupied orbitals 2256 integer a, b, c, d, e, f 2257 integer a_orb, b_orb, c_orb, d_orb, e_orb, f_orb 2258 integer a_spin, b_spin, c_spin, d_spin, e_spin, f_spin 2259 double precision ddot 2260 external ddot 2261c 2262c update t with delta t (in t2) 2263c 2264 do e = noso+1, nso 2265c get spatial orbital and whether alpha or beta spin 2266 e_orb = (e-1)/2 + 1 2267 e_spin = 0 2268 if ((e-2*e_orb).ne.0)e_spin = 1 2269 do f = noso+1, nso 2270 f_orb = (f-1)/2 + 1 2271 f_spin = 0 2272 if ((f-2*f_orb).ne.0)f_spin = 1 2273 do m = 1, noso 2274 m_orb = (m-1)/2 + 1 2275 m_spin = 0 2276 if ((m-2*m_orb).ne.0)m_spin = 1 2277 do n = 1, noso 2278 n_orb = (n-1)/2 + 1 2279 n_spin = 0 2280 if ((n-2*n_orb).ne.0)n_spin = 1 2281c 2282 denom = fock(e_orb,e_orb) + fock(f_orb,f_orb) - 2283 & fock(m_orb,m_orb) - fock(n_orb,n_orb) 2284c write(6,*)' denom: ', denom 2285c 2286 t2(e,f,m,n) = t2(e,f,m,n) - r2(e,f,m,n)/denom 2287 enddo 2288 enddo 2289 enddo 2290 enddo 2291 r2norm = dsqrt(ddot(nso**4,r2,1,r2,1)) 2292 write(6,*)' r2norm = ', r2norm 2293 return 2294 end 2295 subroutine t1_rjh(g_mo, r1, t2, fock, nmo, nso, noso) 2296 implicit none 2297 integer nmo, nso, noso 2298 double precision g_mo(nmo,nmo,nmo,nmo) 2299 double precision fock(nmo,nmo) 2300 double precision r1(nso,nso) 2301 double precision t2(nso,nso,nso,nso) 2302 double precision g, f_t, denom 2303c occupied orbitals 2304 integer i, j, k, l, m, n 2305 integer i_orb, j_orb, k_orb, l_orb, m_orb, n_orb 2306 integer i_spin, j_spin, k_spin, l_spin, m_spin, n_spin 2307c unoccupied orbitals 2308 integer a, b, c, d, e, f 2309 integer a_orb, b_orb, c_orb, d_orb, e_orb, f_orb 2310 integer a_spin, b_spin, c_spin, d_spin, e_spin, f_spin 2311c 2312c put t into t2 expression keeping all terms linear in t2 2313c 2314c write(6,*)' debug inside t1_rjh ' 2315c call jan_debug_print('MOINTS', g_mo, nmo, nmo, 2316c $ nmo, nmo) 2317c write(6,*)' T2 with quadratic terms. ' 2318c call writet2(t2,nso) 2319c write(6,*) ' Fock Matrix' 2320c call output(fock,1,nmo,1,nmo,nmo,nmo,1) 2321c 2322 do e = noso+1, nso 2323c get spatial orbital and whether alpha or beta spin 2324 e_orb = (e-1)/2 + 1 2325 e_spin = 0 2326 if ((e-2*e_orb).ne.0)e_spin = 1 2327 do m = 1, noso 2328 m_orb = (m-1)/2 + 1 2329 m_spin = 0 2330 if ((m-2*m_orb).ne.0)m_spin = 1 2331c 2332 r1(e,m) = 0.0d0 2333 if (e_spin.eq.m_spin)then 2334 r1(e,m) = fock(e_orb,m_orb) 2335 endif 2336c 2337 do j = 1, noso 2338 j_orb = (j-1)/2 + 1 2339 j_spin = 0 2340 if ((j-2*j_orb).ne.0)j_spin = 1 2341 do b = noso+1, nso 2342 b_orb = (b-1)/2 + 1 2343 b_spin = 0 2344 if ((b-2*b_orb).ne.0)b_spin = 1 2345c 2346 f_t = 0.0d0 2347 if (j_spin.eq.b_spin)then 2348 f_t = fock(j_orb,b_orb) 2349 endif 2350c 2351 r1(e,m) = r1(e,m) + f_t*t2(e,b,m,j) 2352c 2353 do i = 1, noso 2354 i_orb = (i-1)/2 + 1 2355 i_spin = 0 2356 if ((i-2*i_orb).ne.0)i_spin = 1 2357c 2358 g=0.0d0 2359 if (i_spin.eq.m_spin.and.j_spin.eq.b_spin)then 2360 g = g_mo(i_orb,j_orb,m_orb,b_orb) 2361 endif 2362 if (i_spin.eq.b_spin.and.j_spin.eq.m_spin)then 2363 g = g - g_mo(i_orb,j_orb,b_orb,m_orb) 2364 endif 2365 2366c write(6,*)' e,m,j,b,i,g:',e,m,j,b,i,g 2367c 2368 r1(e,m) = r1(e,m) - 0.5d0*g*t2(e,b,i,j) 2369c 2370 enddo 2371 enddo 2372 enddo 2373 do b = noso+1, nso 2374 b_orb = (b-1)/2 + 1 2375 b_spin = 0 2376 if ((b-2*b_orb).ne.0)b_spin = 1 2377c 2378 do i = 1, noso 2379 i_orb = (i-1)/2 + 1 2380 i_spin = 0 2381 if ((i-2*i_orb).ne.0)i_spin = 1 2382c 2383 do a = noso+1, nso 2384 a_orb = (a-1)/2 + 1 2385 a_spin = 0 2386 if ((a-2*a_orb).ne.0)a_spin = 1 2387 g=0.0d0 2388 if (i_spin.eq.a_spin.and.e_spin.eq.b_spin)then 2389 g = g_mo(i_orb,e_orb,a_orb,b_orb) 2390 endif 2391 if (i_spin.eq.b_spin.and.e_spin.eq.a_spin)then 2392 g = g - g_mo(i_orb,e_orb,b_orb,a_orb) 2393 endif 2394c 2395c write(6,*)' g:',g 2396c 2397 r1(e,m) = r1(e,m) + 0.5d0*g*t2(a,b,i,m) 2398c 2399 enddo 2400 enddo 2401 enddo 2402 enddo 2403 enddo 2404c write(6,*)' t1(e,m) residual:' 2405c call writet1(r1,nso) 2406 return 2407 end 2408 subroutine t1_update_rjh(g_mo, t1, r1, fock, nmo, nso, noso) 2409 implicit none 2410 integer nmo, nso, noso 2411 double precision g_mo(nmo,nmo,nmo,nmo) 2412 double precision fock(nmo,nmo) 2413 double precision r1(nso,nso) 2414 double precision t1(nso,nso) 2415 double precision g, denom, r1norm 2416c occupied orbitals 2417 integer i, j, k, l, m, n 2418 integer i_orb, j_orb, k_orb, l_orb, m_orb, n_orb 2419 integer i_spin, j_spin, k_spin, l_spin, m_spin, n_spin 2420c unoccupied orbitals 2421 integer a, b, c, d, e, f 2422 integer a_orb, b_orb, c_orb, d_orb, e_orb, f_orb 2423 integer a_spin, b_spin, c_spin, d_spin, e_spin, f_spin 2424 double precision ddot 2425 external ddot 2426c 2427c update t with delta t (in t2) 2428c 2429 do e = noso+1, nso 2430c get spatial orbital and whether alpha or beta spin 2431 e_orb = (e-1)/2 + 1 2432 e_spin = 0 2433 if ((e-2*e_orb).ne.0)e_spin = 1 2434 do m = 1, noso 2435 m_orb = (m-1)/2 + 1 2436 m_spin = 0 2437 if ((m-2*m_orb).ne.0)m_spin = 1 2438c 2439 denom = fock(e_orb,e_orb) - fock(m_orb,m_orb) 2440c write(6,*)' denom: ', denom 2441c 2442 t1(e,m) = t1(e,m) - r1(e,m)/denom 2443 enddo 2444 enddo 2445 r1norm = dsqrt(ddot(nso**2,r1,1,r1,1)) 2446 write(6,*)' r1norm = ', r1norm 2447 return 2448 end 2449 subroutine jan_h( 2450 $ rtdb, basis, 2451 $ n1, n2, 2452 $ lda1, lda2, 2453 $ c1t, c2t, 2454 $ h) 2455 implicit none 2456#include "errquit.fh" 2457#include "bas.fh" 2458#include "global.fh" 2459#include "mafdecls.fh" 2460 integer rtdb, basis 2461 integer lda1, lda2, n1, n2 2462 double precision c1t(lda1, n1), c2t(lda2, n2) 2463 double precision h(n1, n2) 2464c 2465c Return hij = sum(kl) c1t(i,k) hAO(k,l) c2t(j,l) 2466c 2467c Note that the transposed MO coeffs are passed in 2468c to be compatible with jan_full_transform 2469c 2470 integer geom, nbf, g_tmp 2471 integer l_tmp1, k_tmp1, l_tmp2, k_tmp2 2472c 2473 call int_init(rtdb, 1, basis) 2474 if (.not. bas_geom(basis, geom)) 2475 $ call errquit('jan_transform: basis ', basis, BASIS_ERR) 2476 call schwarz_init(geom, basis) 2477 if (.not. bas_numbf(basis, nbf)) 2478 $ call errquit('jan_h: nbf',basis, BASIS_ERR) 2479c 2480 if (.not. ma_push_get(mt_dbl, nbf*nbf,'tmp1', l_tmp1, k_tmp1)) 2481 $ call errquit('tmp1', nbf*nbf, MA_ERR) 2482 if (.not. ma_push_get(mt_dbl, nbf*nbf,'tmp2', l_tmp2, k_tmp2)) 2483 $ call errquit('tmp2', nbf*nbf, MA_ERR) 2484c 2485 if (.not. ga_create(mt_dbl, nbf, nbf, 'tmp', 1, 1, g_tmp)) 2486 & call errquit('scf_v_g: tmp', 0, GA_ERR) 2487 call ga_zero(g_tmp) 2488 call int_1e_ga(basis, basis, g_tmp, 'kinetic', .false.) 2489 call int_1e_ga(basis, basis, g_tmp, 'potential', .false.) 2490 call ga_get(g_tmp, 1, nbf, 1, nbf, dbl_mb(k_tmp1), nbf) 2491 if (.not. ga_destroy(g_tmp)) call errquit('ga?',0, GA_ERR) 2492c 2493 call dgemm('n', 'n', n1, nbf, nbf, 2494 $ 1.0d0, c1t, lda1, dbl_mb(k_tmp1), nbf, 2495 $ 0.0d0, dbl_mb(k_tmp2), n1) 2496 call dgemm('n', 't', n1, n2, nbf, 2497 $ 1.0d0, dbl_mb(k_tmp2), n1, c2t, lda2, 2498 $ 0.0d0, h, n1) 2499c 2500c write(6,*) ' Transformed H' 2501c call output(h, 1, n1, 1, n2, n1, n2, 1) 2502c 2503 call schwarz_tidy() 2504 call int_terminate 2505c 2506 if (.not. ma_chop_stack(l_tmp1)) call errquit('ma',0, MA_ERR) 2507c 2508 return 2509 end 2510 subroutine uccsdtest_lambda(nbf, nmo, no, nv, 2511 $ t1, ct, pt, ht, work) 2512 implicit none 2513c 2514 integer nbf, nmo, no, nv 2515 double precision t1(nv, no) 2516 double precision ct(nmo, nbf) 2517 double precision pt(nmo, nbf) 2518 double precision ht(nmo, nbf) 2519 double precision work(nmo, nmo) 2520c 2521 integer i, a 2522c 2523c Particle transformation 2524c 2525 call dfill(nmo*nmo, 0.0d0, work, 1) 2526 call dfill(nmo, 1.0d0, work, nmo+1) 2527 do i = 1, no 2528 do a = 1, nv 2529 work(i,a+no) = work(i,a+no) - t1(a,i) 2530 end do 2531 end do 2532 call dgemm('t','n',nmo,nbf,nmo, 2533 $ 1.0d0, work, nmo, ct, nmo, 2534 $ 0.0d0, pt, nmo) 2535c 2536c Hole transformation 2537c 2538 call dfill(nmo*nmo, 0.0d0, work, 1) 2539 call dfill(nmo, 1.0d0, work, nmo+1) 2540 do i = 1, no 2541 do a = 1, nv 2542 work(a+no,i) = work(a+no,i) + t1(a,i) 2543 end do 2544 end do 2545 call dgemm('t','n',nmo,nbf,nmo, 2546 $ 1.0d0, work, nmo, ct, nmo, 2547 $ 0.0d0, ht, nmo) 2548c 2549 end 2550 subroutine get_new_f_g(rtdb, basis, fock, g_mo, t1, nso, nmo) 2551 implicit none 2552#include "errquit.fh" 2553#include "global.fh" 2554#include "mafdecls.fh" 2555#include "bas.fh" 2556#include "geom.fh" 2557#include "rtdb.fh" 2558#include "inp.fh" 2559 integer rtdb 2560c 2561 integer nmo, g_tmp, l_occ, k_occ, 2562 $ l_eval, k_eval, l_mos, k_mos, l_most, k_most 2563 integer nocc, nvirt, nopen, nclosed, nso, noso 2564 integer l_t, k_t, l_t2, k_t2, l_f, k_f 2565 integer l_t1, k_t1, l_r1, k_r1 2566 integer l_pmost, k_pmost, l_hmost, k_hmost 2567 integer l_work, k_work, l_h, k_h 2568 integer l_t1_vomos, k_t1_vomos 2569c 2570 double precision g_mo(nmo,nmo,nmo,nmo) 2571 double precision fock(nmo,nmo) 2572 double precision t1(nso,nso) 2573 integer basis, nbf 2574 character*255 movecs ! Name of movector file 2575 character*80 title, name_of_basis, scftype 2576 integer nbf_file, nsets, nmo_file(2) 2577 logical movecs_read, movecs_read_header 2578 external movecs_read, movecs_read_header 2579c 2580 logical int_normalize 2581 external int_normalize 2582c 2583c Read the MO vectors and evals from a RHF calculation 2584c 2585 call util_file_name('movecs',.false.,.false.,movecs) 2586 if (.not. movecs_read_header(movecs, title, name_of_basis, 2587 $ scftype, nbf_file, nsets, nmo_file, 2)) call errquit 2588 $ ('jantest: failed to read movecs header',911, 2589 & DISK_ERR) 2590c write(6,*) ' Read movecs header from ', movecs 2591c write(6,*) ' Job title : ', 2592c $ title(1:inp_strlen(title)) 2593c write(6,*) ' Basis name: ', 2594c $ name_of_basis(1:inp_strlen(name_of_basis)) 2595 nmo = nmo_file(1) 2596 if (.not. rtdb_get(rtdb, 'scf:nclosed', mt_int, 1, nocc)) 2597 $ call errquit('nocc?',0, RTDB_ERR) 2598 if (.not. rtdb_get(rtdb, 'scf:nopen', mt_int, 1, nopen)) 2599 $ call errquit('nopen?',0, RTDB_ERR) 2600 if (nopen .ne. 0) call errquit('asjdlfkadjsl',0, UNKNOWN_ERR) 2601 if (.not. bas_numbf(basis, nbf)) call errquit 2602 $ ('scf_init: basis info',0, BASIS_ERR) 2603 nvirt= nmo - nocc 2604c write(6,*) ' No. of closed shells ', nocc 2605c write(6,*) ' No. of molecular orbitals: ', nmo 2606c write(6,*) ' No. of basis functions: ', nbf 2607c 2608 if (.not. ga_create(mt_dbl, nbf, nmo, 'tmp', 0, 0, g_tmp)) 2609 & call errquit('scf_v_g: tmp', 0, GA_ERR) 2610 if (.not. ma_push_get(mt_dbl, nbf,'occ',l_occ, k_occ)) 2611 $ call errquit('ma occ', nbf, MA_ERR) 2612 if (.not. ma_push_get(mt_dbl, nbf,'eval',l_eval, k_eval)) 2613 $ call errquit('ma eval', nbf, MA_ERR) 2614 if (.not. ma_push_get(mt_dbl, nbf*nbf,'mos', l_mos, k_mos)) 2615 $ call errquit('ma mos', nbf*nbf, MA_ERR) 2616 if (.not. ma_push_get(mt_dbl, nbf*nbf,'mos', l_most, k_most)) 2617 $ call errquit('ma mos', nbf*nbf, MA_ERR) 2618c 2619 if (.not. movecs_read(movecs, 1, dbl_mb(k_occ), dbl_mb(k_eval), 2620 $ g_tmp)) call errquit('movecs_read of amos failed ',0, 2621 & DISK_ERR) 2622 call ga_get(g_tmp, 1, nbf, 1, nmo, dbl_mb(k_mos), nbf) 2623 call util_transpose(dbl_mb(k_mos),nbf,dbl_mb(k_most),nmo, 2624 $ nbf,nmo) 2625c 2626c write(6,*) ' Orbital eigenvalues ' 2627c call output(dbl_mb(k_eval),1,nmo,1,1,nmo,1,1) 2628c write(6,*) ' MOs' 2629c call output(dbl_mb(k_mos),1,nbf,1,nmo,nbf,nmo,1) 2630c write(6,*) ' MOs T' 2631c call output(dbl_mb(k_most),1,nmo,1,nbf,nmo,nbf,1) 2632c 2633 if (.not. ga_destroy(g_tmp)) call errquit(' ga bad?',0, GA_ERR) 2634c 2635c Transform MO integrals in Dirac order to new T1 bases 2636c 2637 if (.not. ma_push_get(mt_dbl, nbf*nbf,'pmost', l_pmost, k_pmost)) 2638 $ call errquit('ma pmost', nbf*nbf, MA_ERR) 2639 if (.not. ma_push_get(mt_dbl, nbf*nbf,'hmost', l_hmost, k_hmost)) 2640 $ call errquit('ma hmost', nbf*nbf, MA_ERR) 2641 if (.not. ma_push_get(mt_dbl, nbf*nbf,'work', l_work, k_work)) 2642 $ call errquit('ma work', nbf*nbf, MA_ERR) 2643 if (.not. ma_push_get(mt_dbl, nbf*nbf,'1e-', l_h, k_h)) 2644 $ call errquit('ma h', nbf*nbf, MA_ERR) 2645 if (.not. ma_push_get(mt_dbl, nvirt*nocc,'t1_vomos', 2646 $ l_t1_vomos, k_t1_vomos)) 2647 $ call errquit('ma t1_vomos', nvirt*nocc, MA_ERR) 2648 call pack_t1(t1, dbl_mb(k_t1_vomos), nso, nmo, nocc, nvirt) 2649 call uccsdtest_lambda(nbf, nmo, nocc, nvirt, 2650 $ dbl_mb(k_t1_vomos), dbl_mb(k_most), 2651 $ dbl_mb(k_pmost), 2652 $ dbl_mb(k_hmost), dbl_mb(k_work)) 2653c 2654c write(6,*) ' P mos' 2655c call output(dbl_mb(k_pmost), 1, nmo, 1, nmo, nmo, nmo, 1) 2656c 2657c write(6,*) ' H mos' 2658c call output(dbl_mb(k_hmost), 1, nmo, 1, nmo, nmo, nmo, 1) 2659c 2660 call jan_h(rtdb, basis, nmo, nmo, nmo, nmo, 2661 $ dbl_mb(k_pmost), dbl_mb(k_hmost), dbl_mb(k_h)) 2662 call jan_full_transform( 2663 $ rtdb, basis, 2664 $ nmo, nmo, nmo, nmo, 2665 $ nmo, nmo, nmo, nmo, 2666 $ dbl_mb(k_pmost), dbl_mb(k_pmost), dbl_mb(k_hmost), 2667 $ dbl_mb(k_hmost), g_mo, 'Dirac') 2668c call jan_debug_print('MOINTS', g_mo, nmo, nmo, 2669c $ nmo, nmo) 2670 call build_f(fock, g_mo, dbl_mb(k_h), nmo, nocc) 2671c 2672c Tidy up 2673c 2674 if (.not. ma_chop_stack(l_occ)) call errquit(' ma chop?', 0, 2675 & MA_ERR) 2676c 2677 return 2678c 2679 end 2680 subroutine pack_t1(t1, t1_vomos, nso, nmo, nocc, nvirt) 2681 implicit none 2682c 2683 integer nso, nmo, nocc, nvirt 2684 integer i, j, i_orb, j_orb, i_so, j_so 2685 double precision t1_vomos(nvirt, nocc) 2686 double precision t1(nso,nso) 2687c 2688c write(6,*)' t1 from pack_t1 ' 2689c do i = 1, nso 2690c do j = 1, nso 2691c if (abs(t1(i,j)).gt.1d-6) 2692c & write(6,*)i,j,t1(i,j) 2693c enddo 2694c enddo 2695c 2696 call dfill(nvirt*nocc, 0.0d0, t1_vomos, 1) 2697c 2698 do i = 1, nvirt 2699 do j = 1, nocc 2700 i_orb = nocc + i 2701 j_orb = j 2702 i_so = 2*(i_orb-1)+1 2703 j_so = 2*(j_orb-1)+1 2704 t1_vomos(i,j) = t1(i_so,j_so) 2705 enddo 2706 enddo 2707c write(6,*)' t1_vomos ' 2708c do i = 1, nvirt 2709c do j = 1, nocc 2710c if (abs(t1_vomos(i,j)).gt.1d-6) 2711c & write(6,*)i,j,t1_vomos(i,j) 2712c enddo 2713c enddo 2714 return 2715 end 2716 subroutine build_f(fock, g_mo, h_mo, nmo, nocc) 2717 implicit none 2718c 2719 integer nmo, nocc 2720 double precision g_mo(nmo,nmo,nmo,nmo) 2721 double precision fock(nmo,nmo), h_mo(nmo,nmo) 2722c 2723c fij = hij + 2*<ik|jk> - <ik|kj> 2724c 2725 integer i, j, k 2726c 2727 call dfill(nmo*nmo, 0.0d0, fock, 1) 2728c 2729 do j = 1, nmo 2730 do i = 1, nmo 2731 fock(i,j) = h_mo(i,j) 2732 do k = 1, nocc 2733 fock(i,j) = fock(i,j) + 2.0d0*g_mo(i,k,j,k) 2734 & - g_mo(i,k,k,j) 2735 end do 2736 end do 2737 end do 2738c 2739c write(6,*) ' FOCK' 2740c call output(fock, 1, nmo, 1, nmo, nmo, nmo, 1) 2741c 2742 return 2743 end 2744 2745c $Id$ 2746