1 logical function jvltest(rtdb) 2* 3* $Id$ 4* calling routine for CCSD(T) trials 5* 6 implicit none 7 integer rtdb 8c... rtdb is handle to database 9c 10 call triple_so(rtdb) 11c 12 jvltest = .true. 13c 14 end 15 16 subroutine triple_so(rtdb) 17c 18c... first stab at a triple code for CCSD(T) 19c... using spin-orbital integrals sitting in core 20c... for now a 4-dimensional array of <ij||ab> 21c 22 implicit none 23#include "errquit.fh" 24#include "global.fh" 25#include "mafdecls.fh" 26#include "rtdb.fh" 27c 28 integer nn_almlof,nn_laguer 29 parameter (nn_almlof=8, nn_laguer=200) 30 integer rtdb,n_almlof,np_almlof(nn_almlof), 31 1 n_laguer,np_laguer(nn_laguer) 32c 33 integer occ,virt,nbasis 34 logical oexist 35 integer l_rint,k_rint,l_t1,k_t1,l_t2,k_t2,l_orb,k_orb 36 integer l_orbn,k_orbn 37 double precision delta,energy 38 integer i 39 character*26 date 40c 41c... check existence of required files 42c 43 inquire(file='T2',exist=oexist) 44 if (.not. oexist) 45 $ call errquit('jvltest: no T2 ', 0, UNKNOWN_ERR) 46 inquire(file='ASOINTS',exist=oexist) 47 if (.not. oexist) 48 $ call errquit('jvltest: no ASOINTS ', 0, UNKNOWN_ERR) 49 inquire(file='EVALS',exist=oexist) 50 if (.not. oexist) 51 $ call errquit('jvltest: no EVALS ', 0, UNKNOWN_ERR) 52c 53c... get basis size from T2 54c 55 open(9,file='T2',form='formatted') 56 read(9,*) occ,virt 57 nbasis = occ + virt 58 close(9) 59c 60c.. integrals 61c double precision rint(nbasis,nbasis,nbasis,nbasis) 62 if (.not. ma_push_get(mt_dbl, nbasis**4,'rint',l_rint, k_rint)) 63 $ call errquit('rint',0, MA_ERR) 64c.. amplitudes 65c double precision t1(occ,virt), t2(occ,occ,virt,virt) 66 if (.not. ma_push_get(mt_dbl, occ**2*virt**2,'t2',l_t2, k_t2)) 67 $ call errquit('t2',0, MA_ERR) 68 if (.not. ma_push_get(mt_dbl, occ*virt,'t1',l_t1, k_t1)) 69 $ call errquit('t1',0, MA_ERR) 70c.. the orb (orbital energies) 71c double precision orb(nbasis) 72 if (.not. ma_push_get(mt_dbl, nbasis,'orb',l_orb, k_orb)) 73 $ call errquit('orb',0, MA_ERR) 74c 75 call get_INTS(dbl_mb(k_rint),nbasis) 76 call get_T1(dbl_mb(k_t1),occ,virt) 77 call get_T2(dbl_mb(k_t2),occ,virt) 78 call get_EVALS(dbl_mb(k_orb),nbasis,occ,virt) 79c 80* call calc_robert(dbl_mb(k_rint),dbl_mb(k_t1),dbl_mb(k_t2), 81* 1 dbl_mb(k_orb),nbasis,occ,virt) 82c 83* call calc_pople(dbl_mb(k_rint),dbl_mb(k_t1),dbl_mb(k_t2), 84* 1 dbl_mb(k_orb),nbasis,occ,virt) 85c 86c... now make modified (homo-lumo) orbital energies and 87c... call all + cullen 88c 89* if (.not. ma_push_get(mt_dbl, nbasis,'orbn',l_orbn, k_orbn)) 90* $ call errquit('orbn',0) 91* call mod_EVALS(dbl_mb(k_orb),dbl_mb(k_orbn),nbasis,occ,virt) 92c 93* print *,' ORBITAL ENERGIES TO HOMO-LUMO ' 94* delta = dbl_mb(k_orb+occ-1) - dbl_mb(k_orb+occ) 95* print *,' delta e used ',delta 96* delta = delta * 3.0d0 97c 98* call calc_robert(dbl_mb(k_rint),dbl_mb(k_t1),dbl_mb(k_t2), 99* 1 dbl_mb(k_orbn),nbasis,occ,virt) 100c 101* call calc_pople(dbl_mb(k_rint),dbl_mb(k_t1),dbl_mb(k_t2), 102* 1 dbl_mb(k_orbn),nbasis,occ,virt) 103c 104c... set orbital energies to 1 to get almlof to do cullen 105c... delta contains de delta e (*3) 106c 107 call dfill(nbasis,1.0d0,dbl_mb(k_orbn),1) 108 delta = 1.0d0 109 call calc_almlof(dbl_mb(k_rint),dbl_mb(k_t1),dbl_mb(k_t2), 110 1 delta,dbl_mb(k_orbn), 111 2 nbasis,occ,virt,energy,.true.) 112c 113c try integrating 114c 115c call mod_EVALS(dbl_mb(k_orb),dbl_mb(k_orbn),nbasis,occ,virt) 116c call dcopy(nbasis,dbl_mb(k_orbn),1,dbl_mb(k_orb),1) 117c print *,' cw keep ORBITAL ENERGIES as HOMO-LUMO ' 118c 119* if (.not.rtdb_get_info(rtdb,'cct_almlof',i,n_almlof,date)) 120* 1 then 121* n_almlof = 1 122* np_almlof(1) = 5 123* else 124* if (.not.rtdb_get(rtdb,'cct_almlof',MT_INT,nn_almlof,np_almlof)) 125* 1 call errquit(' cct_almlof ',0) 126* end if 127c 128* do i=1,n_almlof 129* call int_almlof(dbl_mb(k_rint),dbl_mb(k_t1),dbl_mb(k_t2), 130* 1 dbl_mb(k_orb),dbl_mb(k_orbn), 131* 2 nbasis,occ,virt,energy,np_almlof(i)) 132* print *,np_almlof(i),' point almlof integrated energy ',energy 133* end do 134c 135* if (.not.rtdb_get_info(rtdb,'cct_laguer',i,n_laguer,date)) 136* 1 then 137* n_laguer = 1 138* np_laguer(1) = 5 139* else 140* if (.not.rtdb_get(rtdb,'cct_laguer',MT_INT,nn_laguer,np_laguer)) 141* 1 call errquit(' cct_laguer ',0) 142* end if 143c 144* do i=1,n_laguer 145* call int_laguer(dbl_mb(k_rint),dbl_mb(k_t1),dbl_mb(k_t2), 146* 1 dbl_mb(k_orb),dbl_mb(k_orbn), 147* 2 nbasis,occ,virt,energy,np_laguer(i),rtdb) 148* print *,np_laguer(i),' point laguer integrated energy ',energy 149* end do 150c 151c... free all core beyond rint 152c 153 if (.not. ma_chop_stack(l_rint)) call errquit(' ma chop?', 0, 154 & MA_ERR) 155c 156c could also use ma_pop_stack(l_occ) which is like gmem_free 157c (but using handle instead of pointer) 158c 159 end 160 subroutine int_almlof(rint,t1,t2,orb,fac,nbasis,occ,virt, 161 1 energy,np) 162c 163c... control routine for Almloef integration 164c... energy returns the finally computed energy 165c... np is thew # points to be used 166c 167 implicit none 168#include "errquit.fh" 169c 170 integer nbasis,occ,virt,np 171 double precision rint(nbasis,nbasis,nbasis,nbasis) 172 double precision t2(occ,occ,virt,virt),t1(occ,virt) 173 double precision orb(nbasis),fac(nbasis),energy 174c 175 integer ip,itw 176 double precision enerp 177 common/flop/ term1,term2,term3,term4 178 double precision term1,term2,term3,term4,tt1,tt2,tt3,tt4 179c 180 integer npoint,nd 181 parameter (npoint=8, nd=npoint*(npoint+1)/2) 182 double precision t_alpha(nd), weight(nd) 183c 184 data t_alpha/0.2543066d0, 185 2 0.1057374d0, 0.8263985d0, 186 3 0.0664618d0, 0.4380785d0, 1.6101925d0, 187 4 0.0135967d0, 0.1646374d0, 0.6569702d0, 2.0357815d0, 188 5 0.0064785d0, 0.0669154d0, 0.3089156d0, 0.9478808d0, 189 5 2.5786393d0, 190 6 0.0060050d0, 0.0553016d0, 0.2346500d0, 0.6587781d0, 191 6 1.5760606d0, 3.6502575d0, 192 7 0.0035243d0, 0.0251939d0, 0.1193583d0, 0.3653936d0, 193 7 0.8850526d0, 1.9454617d0, 4.2401468d0, 194 8 0.0034250d0, 0.0238627d0, 0.1023430d0, 0.2947359d0, 195 8 0.6736046d0, 1.3821837d0, 2.7133497d0, 5.4051090d0/ 196 data weight /0.7443778d0, 197 2 0.2998529d0, 1.3750404d0, 198 3 0.1793919d0, 0.6329602d0, 1.9996848d0, 199 4 0.0527171d0, 0.2757647d0, 0.7858024d0, 2.2786679d0, 200 5 0.0177806d0, 0.1278500d0, 0.3878375d0, 0.9741098d0, 201 5 2.6125732d0, 202 6 0.0164421d0, 0.0995694d0, 0.2759458d0, 0.6106806d0, 203 6 1.3178664d0, 3.1857008d0, 204 7 0.0095904d0, 0.0437062d0, 0.1566528d0, 0.3542937d0, 205 7 0.7264528d0, 1.4928836d0, 3.4634348d0, 206 8 0.0093035d0, 0.0390022d0, 0.1267465d0, 0.2692945d0, 207 8 0.5112993d0, 0.9512673d0, 1.8163835d0, 3.9522146d0/ 208c 209 tt1 = 0.0d0 210 tt2 = 0.0d0 211 tt3 = 0.0d0 212 tt4 = 0.0d0 213c 214 if (np.gt.npoint) call errquit('almlof overflow',0, UNKNOWN_ERR) 215 energy=0.0d0 216 do ip=1,np 217 energy=energy+dexp(-5.0d0*t_alpha(np*(np-1)/2 + ip)) 218 1 *weight(np*(np-1)/2 + ip) 219 end do 220 print *,np,' points 0.2 is ', energy 221c 222 energy = 0.0d0 223 do ip=1,np 224 itw = np*(np-1)/2 + ip 225 call cct_fac(orb,fac,nbasis,occ,virt,t_alpha(itw)) 226c 227 call calc_almlof(rint,t1,t2,-1.0d0,fac, 228 1 nbasis,occ,virt,enerp,.false.) 229c 230 energy = energy + enerp * weight(itw) 231c print *,' point ',ip,' of ',np,' contr ',enerp,' energy ',energy 232 tt1 = tt1 + term1*weight(itw) 233 tt2 = tt2 + term2*weight(itw) 234 tt3 = tt3 + term3*weight(itw) 235 tt4 = tt4 + term4*weight(itw) 236 end do 237c print *,' t ',tt1,tt2,tt3,tt4 238c 239 return 240 end 241 subroutine int_laguer(rint,t1,t2,orb,fac,nbasis,occ,virt, 242 1 energy,np,rtdb) 243c 244c... control routine for Almloef integration 245c... using straight gauss-laguerre 246c... energy returns the finally computed energy 247c... np is thew # points to be used 248c 249 implicit none 250#include "errquit.fh" 251#include "global.fh" 252#include "mafdecls.fh" 253#include "rtdb.fh" 254c 255 integer nbasis,occ,virt,np,rtdb 256 double precision rint(nbasis,nbasis,nbasis,nbasis) 257 double precision t2(occ,occ,virt,virt),t1(occ,virt) 258 double precision orb(nbasis),fac(nbasis),energy 259c 260 integer ip,i,n 261 double precision enerp,sum,avocc,avirt,delta 262c 263 integer npoint 264 parameter (npoint= 200) 265 double precision t_alpha(npoint), weight(npoint) 266 double precision scra(npoint), dum(2),elimit(2) 267 character*26 date 268c 269 if (np.gt.npoint) call errquit(' enlarge int_laguer ',0, 270 & UNKNOWN_ERR) 271 call gaussq(6,np,0.0d0,0.0d0,0,dum,scra,t_alpha,weight) 272c 273c... dermine integration scale required 274c default : average delta e (*3) 275c might exclude extremes or all 276c 277 if (.not.rtdb_get_info(rtdb,'cct_laguer_scale',i,n,date)) 278 1 then 279c.. default (-5 to +5 for now) 280 elimit(1) = -5.0d0 281 elimit(2) = +5.0d0 282 else 283 if (n.eq.0) then 284c.. next default (all) 285 elimit(1) = -9.0d99 286 elimit(2) = +9.0d99 287 else 288 if (.not.rtdb_get(rtdb,'cct_laguer_scale',MT_DBL,2,elimit)) 289 1 call errquit(' cct_laguer_scale ',0, RTDB_ERR) 290 if (n.eq.1) elimit(2) = elimit(1) 291 elimit(1) = -elimit(1) 292 end if 293 end if 294c 295c... compute scale 296c 297 sum = 0.0d0 298 n = 0 299 do i=1,occ 300 if (orb(i).gt.elimit(1)) then 301 n = n + 1 302 sum = sum + orb(i) 303 end if 304 end do 305 avocc = sum/max0(n,1) 306 ip = n 307 sum = 0.0d0 308 n = 0 309 do i=occ+1,nbasis 310 if (orb(i).lt.elimit(2)) then 311 n = n + 1 312 sum = sum + orb(i) 313 end if 314 end do 315 avirt = sum/max0(n,1) 316c... scale 317 delta = 3.0d0*(avirt-avocc) 318 if (elimit(1).eq.0.0d0.and.elimit(2).ne.0.0d0) delta = elimit(2) 319 if (delta .lt.1.0d-2) delta = 1.0d0 320 print *,' laguerre scaled by ',delta,' kept ',ip,'/',n 321c 322c... add the ex(xi) to weights 323c 324 energy=0.0d0 325 do ip=1,np 326 energy=energy+dexp(-5.0d0*t_alpha(ip)/delta)*weight(ip) 327 1 *dexp(t_alpha(ip))/delta 328 end do 329 print *,np,' points 0.2 is ', energy 330c 331 energy = 0.0d0 332 do ip=1,np 333 call cct_fac(orb,fac,nbasis,occ,virt,t_alpha(ip)/delta) 334c 335 call calc_almlof(rint,t1,t2,-1.0d0,fac, 336 1 nbasis,occ,virt,enerp,.false.) 337c 338 energy = energy + enerp * weight(ip)*dexp(t_alpha(ip))/delta 339 end do 340c 341 return 342 end 343 344 subroutine get_EVALS(orb,nbasis,occ,virt) 345c 346c... read orbital energies from EVALS 347c 348 implicit none 349#include "errquit.fh" 350c 351 integer nbasis,occ,virt 352 double precision orb(nbasis) 353c 354 integer nbin,i 355 logical oexist 356c 357 inquire(file='EVALS',exist=oexist) 358 if (.not.oexist) call errquit(' no EVAL ',0, UNKNOWN_ERR) 359 open(9,file='EVALS',form='formatted') 360 read(9,*) nbin 361 if (nbin.ne.nbasis) call errquit('nbasis wrong',1, BASIS_ERR) 362 do i=1,nbasis 363 read(9,'(5x,f20.10)') orb(i) 364 end do 365 close(9) 366c 367 print *,' orbital energies read ' 368c print *, orb 369 return 370 end 371 subroutine mod_EVALS(orb,orbn,nbasis,occ,virt) 372c 373c... make modified orbital energies (homo/lumo) 374c 375 implicit none 376c 377 integer nbasis,occ,virt 378 double precision orb(nbasis),orbn(nbasis) 379 integer i 380c 381 do i=1,occ 382 orbn(i) = orb(occ) 383 end do 384 do i=occ+1,nbasis 385 orbn(i) = orb(occ+1) 386 end do 387c 388 return 389 end 390 subroutine get_T2(t2,occ,virt) 391c 392c... real t2 amplitudes from T2 393c 394 implicit none 395#include "errquit.fh" 396c 397 integer occ,virt 398 double precision t2(occ,occ,virt,virt) 399c 400 integer ocin,vin,i,j,a,b 401 logical oexist 402c 403 inquire(file='T2',exist=oexist) 404 if (.not.oexist) call errquit(' no T2 ',0, UNKNOWN_ERR) 405 open(9,file='T2',form='formatted') 406 read(9,*) ocin,vin 407 if (ocin.ne.occ.or.vin.ne.virt) call errquit(' occ,virt wrong',0, 408 & UNKNOWN_ERR) 409c 410 call dfill(occ*occ*virt*virt,0.0d0,t2,1) 411c 41210 read(9,'(4i5,f20.10)',end=20) i,j,a,b,t2(i,j,a,b) 413 go to 10 414c 41520 close(9) 416 print *,' T2 read ' 417c print *,((((t2(i,j,a,b),i=1,3),j=1,3),a=1,3),b=1,3) 418c 419 return 420 end 421 subroutine get_T1(t1,occ,virt) 422c 423c... real t1 amplitudes from T1 424c 425 implicit none 426#include "errquit.fh" 427c 428 integer occ,virt 429 double precision t1(occ,virt) 430c 431 integer ocin,vin,i,a 432 logical oexist 433c 434 inquire(file='T1',exist=oexist) 435 if (.not.oexist) then 436 print *,' no T1 ' 437 call dfill(occ*virt,0.0d0,t1,1) 438 return 439 end if 440 open(9,file='T1',form='formatted') 441 read(9,*) ocin,vin 442 if (ocin.ne.occ.or.vin.ne.virt) call errquit(' occ,virt wrong',1, 443 & UNKNOWN_ERR) 444c 445 call dfill(occ*virt,0.0d0,t1,1) 446c 44710 read(9,'(2i5,f20.10)',end=20) i,a,t1(i,a) 448 go to 10 449c 45020 close(9) 451 print *,' T1 read ' 452c 453 return 454 end 455 subroutine get_INTS(rint,nbasis) 456c 457c... real integrals from ASOINTS 458c 459 implicit none 460#include "errquit.fh" 461c 462 integer nbasis 463 double precision rint(nbasis,nbasis,nbasis,nbasis) 464c 465 integer in,i,j,k,l 466 logical oexist 467c 468 inquire(file='ASOINTS',exist=oexist) 469 if (.not.oexist) call errquit(' no INTS ',0, DISK_ERR) 470 open(9,file='ASOINTS',form='formatted') 471 read(9,*) in 472 if (in.ne.nbasis) call errquit(' nbasis on INTS wrong',0, 473 & BASIS_ERR) 474c 475 call dfill(nbasis**4,0.0d0,rint,1) 476c 47710 read(9,'(4i5,f20.10)',end=20) i,j,k,l,rint(i,j,k,l) 478 go to 10 479c 48020 close(9) 481 print *,' integrals read ' 482c print *,((((rint(i,j,k,l),i=1,3),j=1,3),k=1,3),l=1,3) 483c 484 return 485 end 486 subroutine cct_fac(orb,fac,nbasis,occ,virt,factor) 487c 488c... generatewd factores for laplace integration 489c... calculated as exp(-orb)*factor for occ 490c... exp(orb)*factor for virt 491c 492 implicit none 493#include "errquit.fh" 494 integer nbasis,occ,virt 495 double precision orb(nbasis),fac(nbasis),factor 496c 497 integer i 498c 499 do i=1,occ 500 fac(i) = dexp(orb(i)*factor) 501 end do 502 do i=occ+1,nbasis 503 fac(i) = dexp(-orb(i)*factor) 504 end do 505c 506 if (orb(occ).gt.0.0d0.or.orb(occ+1).lt.0.0d0) 507 1 call errquit(' OEPS orbital energies ',0, UNKNOWN_ERR) 508c print *,' factor ... ',orb(occ),orb(occ+1) 509c print *,fac(1),fac(nbasis) 510c print *,(fac(1)*fac(nbasis))**3 511c 512 return 513 end 514 subroutine calc_robert(rint,t1,t2,orb,nbasis,occ,virt) 515c 516c... use Robert's formulae 517c 518 implicit none 519c 520 integer nbasis,occ,virt 521 double precision rint(nbasis,nbasis,nbasis,nbasis) 522 double precision t1(occ,virt), t2(occ,occ,virt,virt) 523 double precision orb(nbasis) 524c.. intermediates 525 double precision w1,w2,w3,w4,v,w1w1,w1w2,w1w3,w1w4 526 double precision vw1,vw2,vw3,vw4 527c 528 double precision e1,e2,energy,delta 529 integer i,j,k,a,b,c,e,m 530c 531 w1w1 = 0.0d0 532 w1w2 = 0.0d0 533 w1w3 = 0.0d0 534 w1w4 = 0.0d0 535 vw1 = 0.0d0 536 vw2 = 0.0d0 537 vw3 = 0.0d0 538 vw4 = 0.0d0 539c 540 do i=1,occ 541 do j=1,occ 542 do k=1,occ 543 do a=1,virt 544 do b=1,virt 545 do c=1,virt 546 w1 = 0.0d0 547 w2 = 0.0d0 548 w3 = 0.0d0 549 w4 = 0.0d0 550 do e=1,virt 551 w1 = w1 + t2(i,j,c,e)*rint(a+occ,b+occ,e+occ,k) 552 w2 = w2 + t2(k,j,c,e)*rint(a+occ,b+occ,e+occ,i) 553 w3 = w3 + t2(i,j,a,e)*rint(c+occ,b+occ,e+occ,k) 554 w4 = w4 + t2(k,j,a,e)*rint(c+occ,b+occ,e+occ,i) 555 end do 556 do m=1,occ 557 w1 = w1 + t2(k,m,a,b)*rint(c+occ,m,i,j) 558 w2 = w2 + t2(i,m,a,b)*rint(c+occ,m,k,j) 559 w3 = w3 + t2(k,m,c,b)*rint(a+occ,m,i,j) 560 w4 = w4 + t2(i,m,c,b)*rint(a+occ,m,k,j) 561 end do 562 v = t1(k,c)*rint(i,j,a,b) 563 delta = (orb(i)+orb(j)+orb(k) 564 1 - orb(a+occ)-orb(b+occ)-orb(c+occ)) 565 w1w1 = w1w1 + w1*w1 / delta 566 w1w2 = w1w2 + w1*w2 / delta 567 w1w3 = w1w3 + w1*w3 / delta 568 w1w4 = w1w4 + w1*w4 / delta 569 vw1 = vw1 + v*w1 / delta 570 vw2 = vw2 + v*w2 / delta 571 vw3 = vw3 + v*w3 / delta 572 vw4 = vw4 + v*w4 / delta 573 end do 574 end do 575 end do 576 end do 577 end do 578 end do 579c 580c... calculate energies 581c 582 e1 = (vw1 - 2.0d0*vw2 -2.0d0*vw3 + 4.0d0*vw4) /4.0d0 583 e2 = (w1w1 - 2.0d0*w1w2 -2.0d0*w1w3 + 4.0d0*w1w4) /4.0d0 584 print *,' w1w1 ',w1w1 585 print *,' w1w2 ',w1w2 586 print *,' w1w3 ',w1w3 587 print *,' w1w4 ',w1w4 588 energy = e1 + e2 589c 590 print *,' Robert : e1= ',e1,' e2= ',e2,' e= ',energy 591c 592 end 593 subroutine calc_pople(rint,t1,t2,orb,nbasis,occ,virt) 594c 595c... use Pople's formulae (who can argue with a Nobel-prize) 596c... BUT ... typed in from Laplace triple document (jvl) 597c... ( + formulae(QCI ...)) 598c 599 implicit none 600c 601 integer nbasis,occ,virt 602 double precision rint(nbasis,nbasis,nbasis,nbasis) 603 double precision t1(occ,virt), t2(occ,occ,virt,virt) 604 double precision orb(nbasis) 605c.. intermediates 606 double precision u1,u2,u3,u4,u5,u6,u7,u8,u9,v1 607c 608 double precision e1,e2,energy,delta 609 integer i,j,k,a,b,c,e,m 610c 611 e1 = 0.0d0 612 e2 = 0.0d0 613c 614 do i=1,occ 615 do j=1,occ 616 do k=1,occ 617 do a=1,virt 618 do b=1,virt 619 do c=1,virt 620 u1 = 0.0d0 621 u2 = 0.0d0 622 u3 = 0.0d0 623 u4 = 0.0d0 624 u5 = 0.0d0 625 u6 = 0.0d0 626 u7 = 0.0d0 627 u8 = 0.0d0 628 u9 = 0.0d0 629 do e=1,virt 630 u1 = u1 + t2(i,j,a,e)*rint(b+occ,c+occ,e+occ,k) 631 u2 = u2 + t2(i,j,b,e)*rint(c+occ,a+occ,e+occ,k) 632 u3 = u3 + t2(i,j,c,e)*rint(a+occ,b+occ,e+occ,k) 633 u4 = u4 + t2(k,i,a,e)*rint(b+occ,c+occ,e+occ,j) 634 u5 = u5 + t2(k,i,b,e)*rint(c+occ,a+occ,e+occ,j) 635 u6 = u6 + t2(k,i,c,e)*rint(a+occ,b+occ,e+occ,j) 636 u7 = u7 + t2(j,k,a,e)*rint(b+occ,c+occ,e+occ,i) 637 u8 = u8 + t2(j,k,b,e)*rint(c+occ,a+occ,e+occ,i) 638 u9 = u9 + t2(j,k,c,e)*rint(a+occ,b+occ,e+occ,i) 639 end do 640 do m=1,occ 641 u1 = u1 + t2(i,m,a,b)*rint(c+occ,m,j,k) 642 u2 = u2 + t2(i,m,b,c)*rint(a+occ,m,j,k) 643 u3 = u3 + t2(i,m,c,a)*rint(b+occ,m,j,k) 644 u4 = u4 + t2(j,m,a,b)*rint(c+occ,m,k,i) 645 u5 = u5 + t2(j,m,b,c)*rint(a+occ,m,k,i) 646 u6 = u6 + t2(j,m,c,a)*rint(b+occ,m,k,i) 647 u7 = u7 + t2(k,m,a,b)*rint(c+occ,m,i,j) 648 u8 = u8 + t2(k,m,b,c)*rint(a+occ,m,i,j) 649 u9 = u9 + t2(k,m,c,a)*rint(b+occ,m,i,j) 650 end do 651 v1 = t1(i,a)*rint(j,k,b+occ,c+occ) 652 1 + t1(i,b)*rint(j,k,c+occ,a+occ) 653 2 + t1(i,c)*rint(j,k,a+occ,b+occ) 654 3 + t1(j,a)*rint(k,i,b+occ,c+occ) 655 4 + t1(j,b)*rint(k,i,c+occ,a+occ) 656 5 + t1(j,c)*rint(k,i,a+occ,b+occ) 657 6 + t1(k,a)*rint(i,j,b+occ,c+occ) 658 7 + t1(k,b)*rint(i,j,c+occ,a+occ) 659 8 + t1(k,c)*rint(i,j,a+occ,b+occ) 660 delta = (orb(i)+orb(j)+orb(k) 661 1 - orb(a+occ)-orb(b+occ)-orb(c+occ)) 662 e1 = e1+ v1*(u1+u2+u3+u4+u5+u6+u7+u8+u9) / delta 663 e2 = e2+ (u1+u2+u3+u4+u5+u6+u7+u8+u9)**2 / delta 664 end do 665 end do 666 end do 667 end do 668 end do 669 end do 670c 671c... calculate energies 672c 673 e1 = e1 / 36.0d0 674 e2 = e2 / 36.0d0 675 energy = e1 + e2 676c 677 print *,' Pople : e1= ',e1,' e2= ',e2,' e= ',energy 678c 679 end 680 subroutine calc_almlof(rint,t1,t2,delta,fac,nbasis,occ,virt, 681 1 energy,opr) 682c 683c... this routine can do Cullen+Zerner's approach by 684c... setting all fac'a ,to 1 and delta to delta e 685c... for almlof delta is set to 1 686c 687c... starting from "Robert's" formula 688c notation v(vv) : running index (virtual), o(oo)= running occupied) 689c 690c 691 implicit none 692#include "errquit.fh" 693c 694 integer nbasis,occ,virt 695 double precision rint(nbasis,nbasis,nbasis,nbasis) 696 double precision t1(occ,virt), t2(occ,occ,virt,virt) 697 double precision delta,fac(nbasis),energy 698 logical opr 699#include "global.fh" 700#include "mafdecls.fh" 701c 702c... functions 703c 704 double precision cct_so_a4a,cct_so_a4b,cct_so_b4b 705 double precision cct_so_a3a,cct_so_a3b,cct_so_b3b 706 double precision cct_so_a2a,cct_so_a2b,cct_so_b2b 707 double precision cct_so_a1a,cct_so_a1b,cct_so_b1b 708c 709 double precision term4,term_a4a,term_a4b,term_b4b 710 double precision term3,term_a3a,term_a3b,term_b3b 711 double precision term2,term_a2a,term_a2b,term_b2b 712 double precision term1,term_a1a,term_a1b,term_b1b 713 double precision e1,e2 714 common/flop/ term1,term2,term3,term4 715c 716 integer lenabcd,leniabc,lenijab,lenijka,lenijkl,lenab,lenia,lenij 717 integer l_w,k_w 718c 719 lenabcd = virt*virt*virt*virt 720 leniabc = occ*virt*virt*virt 721 lenijab = occ*occ*virt*virt 722 lenijka = occ*occ*occ*virt 723 lenijkl = occ*occ*occ*occ 724 lenab = virt*virt 725 lenia = virt*occ 726 lenij = occ*occ 727 if (opr) print *,' triple delta e used ',delta 728c 729c term 1 730c 731 if (.not. ma_push_get(mt_dbl,lenab*2,'a2a', l_w, k_w)) 732 1 call errquit('push a1a',0, MA_ERR) 733 term_a1a = cct_so_a1a(rint,nbasis,t2,occ,virt,fac, 734 1 dbl_mb(k_w),dbl_mb(k_w+lenab)) 735 if (.not.ma_pop_stack(l_w)) call errquit('pop a2a',0, MA_ERR) 736 if (.not. ma_push_get(mt_dbl,lenia*2,'a1b', l_w, k_w)) 737 1 call errquit('push a1b',0, MA_ERR) 738 term_a1b = cct_so_a1b(rint,nbasis,t2,occ,virt,fac, 739 1 dbl_mb(k_w),dbl_mb(k_w+lenia)) 740 if (.not.ma_pop_stack(l_w)) call errquit('pop a1b',0, MA_ERR) 741 if (.not. ma_push_get(mt_dbl,lenij*2,'b1b', l_w, k_w)) 742 1 call errquit('push b1b',0, MA_ERR) 743 term_b1b = cct_so_b1b(rint,nbasis,t2,occ,virt,fac, 744 1 dbl_mb(k_w),dbl_mb(k_w+lenij)) 745 if (.not.ma_pop_stack(l_w)) call errquit('pop b1b',0, MA_ERR) 746c 747 term1 = term_a1a + 2.0d0*term_a1b + term_b1b 748c 749 if (opr) then 750c print *,' term a1a ',term_a1a/ delta 751c print *,' term a1b=b1a ',term_a1b/ delta 752c print *,' term b1b ',term_b1b/ delta 753 print *,' term1 ',term1/ delta 754 end if 755 756c 757c term 2 758c 759 if (.not. ma_push_get(mt_dbl,lenijab+lenijab,'a2a', l_w, k_w)) 760 1 call errquit('push a2a',0, MA_ERR) 761 term_a2a = cct_so_a2a(rint,nbasis,t2,occ,virt,fac, 762 1 dbl_mb(k_w),dbl_mb(k_w+lenijab)) 763 if (.not.ma_pop_stack(l_w)) call errquit('pop a2a',0, MA_ERR) 764 if (.not. ma_push_get(mt_dbl,lenijka*2,'a2a', l_w, k_w)) 765 1 call errquit('push a2b',0, MA_ERR) 766 term_a2b = cct_so_a2b(rint,nbasis,t2,occ,virt,fac, 767 1 dbl_mb(k_w),dbl_mb(k_w+lenijka)) 768 if (.not.ma_pop_stack(l_w)) call errquit('pop a2b',0, MA_ERR) 769 if (.not. ma_push_get(mt_dbl,lenijka+lenijkl,'b2b', l_w, k_w)) 770 1 call errquit('push b2b',0, MA_ERR) 771 term_b2b = cct_so_b2b(rint,nbasis,t2,occ,virt,fac, 772 1 dbl_mb(k_w),dbl_mb(k_w+lenijka)) 773 if (.not.ma_pop_stack(l_w)) call errquit('pop b2b',0, MA_ERR) 774c 775 term2 = term_a2a + 2.0d0*term_a2b + term_b2b 776c 777 if (opr) then 778c print *,' term a2a ',term_a2a/ delta 779c print *,' term a2b=b2a ',term_a2b/ delta 780c print *,' term b2b ',term_b2b/ delta 781 print *,' term2 ',term2/ delta 782 end if 783c 784c term 3 785c 786 if (.not. ma_push_get(mt_dbl,lenabcd+leniabc,'a3a', l_w, k_w)) 787 1 call errquit('push a3a',0, MA_ERR) 788 term_a3a = cct_so_a3a(rint,nbasis,t2,occ,virt,fac, 789 1 dbl_mb(k_w),dbl_mb(k_w+lenabcd)) 790 if (.not.ma_pop_stack(l_w)) call errquit('pop a3a',0, MA_ERR) 791 if (.not. ma_push_get(mt_dbl,leniabc+lenijka,'a3b', l_w, k_w)) 792 1 call errquit('push a3b',0, MA_ERR) 793 term_a3b = cct_so_a3b(rint,nbasis,t2,occ,virt,fac, 794 1 dbl_mb(k_w),dbl_mb(k_w+leniabc)) 795 if (.not.ma_pop_stack(l_w)) call errquit('pop a3b',0, MA_ERR) 796 if (.not. ma_push_get(mt_dbl,lenijab+lenijka,'b3b', l_w, k_w)) 797 1 call errquit('push b3b',0, MA_ERR) 798 term_b3b = cct_so_b3b(rint,nbasis,t2,occ,virt,fac, 799 1 dbl_mb(k_w),dbl_mb(k_w+lenijab)) 800 if (.not.ma_pop_stack(l_w)) call errquit('pop b3b',0, MA_ERR) 801c 802 print *,' term_a3a *.5 ',term_a3a*0.5d0 803 print *,' term_a3b *1, ',term_a3b 804 print *,' term_b3b *.5 ',term_b3b*0.5d0 805 term3 = term_a3a + 2.0d0*term_a3b + term_b3b 806c 807 if (opr) then 808c print *,' term a3a ',term_a3a/ delta 809c print *,' term a3b=b3a ',term_a3b/ delta 810c print *,' term b3b ',term_b3b/ delta 811 print *,' term3 ',term3/ delta 812 end if 813 814c 815c term 4 816c 817 if (.not. ma_push_get(mt_dbl,leniabc+lenijab,'a4a', l_w, k_w)) 818 1 call errquit('push a4a',0, MA_ERR) 819 term_a4a = cct_so_a4a(rint,nbasis,t2,occ,virt,fac, 820 1 dbl_mb(k_w),dbl_mb(k_w+leniabc)) 821 if (.not.ma_pop_stack(l_w)) call errquit('pop a4a',0, MA_ERR) 822c 823 if (.not. ma_push_get(mt_dbl,lenijab+lenijka,'a4b', l_w, k_w)) 824 1 call errquit('push a4b',0, MA_ERR) 825 term_a4b = cct_so_a4b(rint,nbasis,t2,occ,virt,fac, 826 1 dbl_mb(k_w),dbl_mb(k_w+lenijab)) 827 if (.not.ma_pop_stack(l_w)) call errquit('pop a4b',0, MA_ERR) 828c 829 if (.not. ma_push_get(mt_dbl,lenijka*2,'b4b', l_w, k_w)) 830 1 call errquit('push b4b',0, MA_ERR) 831 term_b4b = cct_so_b4b(rint,nbasis,t2,occ,virt,fac, 832 1 dbl_mb(k_w),dbl_mb(k_w+lenijka)) 833 if (.not.ma_pop_stack(l_w)) call errquit('pop b4b',0, MA_ERR) 834c 835 term4 = term_a4a + 2.0d0*term_a4b + term_b4b 836c 837 if (opr) then 838c print *,' term a4a ',term_a4a/ delta 839c print *,' term a4b=b4a ',term_a4b/ delta 840c print *,' term b4b ',term_b4b/ delta 841 print *,' term4 ',term4/ delta 842 end if 843c 844 e1 = 0.0d0 845 e2 = (0.25d0*term1 - 0.5d0*term2 - 0.5d0*term3 + term4)/ delta 846 energy = e1 + e2 847 if (opr) print *,' Cullen : e1= ',e1,' e2= ',e2,' e= ',energy 848c 849 end 850 double precision function cct_so_a4a(rint,nbasis,t2,occ,virt,fac, 851 1 wiabc,wijab) 852c 853c term A4A 854c t(ij,ce) (ab||ek) t(kj,af) (cb||fi) 855c = t(ij,ce) (cb||fi) (ab||ek) t(kj,af) 856c weighting 1. 857c 858c... fac has to be added to i,j,k,a,b,c (there is choice) 859c... early is done / late might be better ? 860c 861 implicit none 862c 863 integer nbasis,occ,virt 864 double precision rint(nbasis,nbasis,nbasis,nbasis) 865 double precision t2(occ,occ,virt,virt),fac(nbasis) 866c 867 integer j,e,b,f,i,c,k,a 868 double precision wiabc(virt,virt,virt,occ), 869 1 wijab(virt,virt,occ,occ),sum 870c 871c... t2 * integral (ic) t(ij,ce) (cb||fi) => (j,e,b,f) 872c 873 do j=1,occ 874 do e=1,virt 875 do b=1,virt 876 do f=1,virt 877 wiabc(f,b,e,j) = 0.0d0 878 do i=1,occ 879 do c=1,virt 880 wiabc(f,b,e,j) = wiabc(f,b,e,j) + t2(i,j,c,e) 881 1 * rint(c+occ,b+occ,f+occ,i) 882 2 * fac(j)*fac(b+occ)*fac(i)*fac(c+occ) 883 end do 884 end do 885 end do 886 end do 887 end do 888 end do 889c 890c... (eb) (j,e,b,f) (ab||ek) => (j,k,f,a) (choice) 891c 892 do j=1,occ 893 do k=1,occ 894 do f=1,virt 895 do a=1,virt 896 wijab(a,f,k,j) = 0.0d0 897 do e=1,virt 898 do b = 1,virt 899 wijab(a,f,k,j) = wijab(a,f,k,j) + wiabc(f,b,e,j) 900 1 * rint(a+occ,b+occ,e+occ,k) 901 2 * fac(k)*fac(a+occ) 902 end do 903 end do 904 end do 905 end do 906 end do 907 end do 908c 909c... (jkfa) (j,k,f,a) t(kj,af) 910c 911 sum = 0.0d0 912 do j=1,occ 913 do k=1,occ 914 do f=1,virt 915 do a=1,virt 916 sum= sum + wijab(a,f,k,j) * t2(k,j,a,f) 917 end do 918 end do 919 end do 920 end do 921c 922 cct_so_a4a = sum 923c 924 return 925 end 926 double precision function cct_so_a4b(rint,nbasis,t2,occ,virt,fac, 927 1 wijab,wijka) 928c 929c term A4B = B4A 930c A4B = t(ij,ce) (ab||ek) t(in,cb) (an||kj) 931c = t(ij,ce) t(in,cb) (ab||ek) (an||kj) 932c 933c B4A = t(km,ab) (cm||ij) t(kj,af) (cb||fi) 934c relabel in cb an kj ij ce ab ek 935c weighting 1. 936c 937 implicit none 938c 939 integer nbasis,occ,virt 940 double precision rint(nbasis,nbasis,nbasis,nbasis) 941 double precision t2(occ,occ,virt,virt),fac(nbasis) 942c 943 integer j,n,e,b,i,c,k,a 944 double precision wijab(virt,virt,occ,occ), 945 1 wijka(virt,occ,occ,occ),sum 946c 947c.. contract t's (ic) t(ij,ce) t(in,cb) => (j,n,e,b) choice 948c 949 do j=1,occ 950 do n=1,occ 951 do e=1,virt 952 do b=1,virt 953 wijab(b,e,n,j) = 0.0d0 954 do i=1,occ 955 do c=1,virt 956 wijab(b,e,n,j) = wijab(b,e,n,j) + t2(i,j,c,e) 957 1 * t2(i,n,c,b) 958 2 * fac(j)*fac(b+occ)*fac(i)*fac(c+occ) 959 end do 960 end do 961 end do 962 end do 963 end do 964 end do 965c 966c.. (be) (j,n,e,b) (ab||ek) => (j,n,k,a) (choice) 967c 968 do j=1,occ 969 do n=1,occ 970 do k=1,occ 971 do a=1,virt 972 wijka(a,k,n,j) = 0.0d0 973 do b=1,virt 974 do e=1,virt 975 wijka(a,k,n,j) = wijka(a,k,n,j) + wijab(b,e,n,j) 976 1 * rint(a+occ,b+occ,e+occ,k) 977 2 * fac(k)*fac(a+occ) 978 end do 979 end do 980 end do 981 end do 982 end do 983 end do 984c 985c (jnka) (j,n,k,a) ((an||kj) => sum 986c 987 sum = 0.0d0 988 do j=1,occ 989 do n=1,occ 990 do k=1,occ 991 do a=1,virt 992 sum = sum + wijka(a,k,n,j) * rint(a+occ,n,k,j) 993 end do 994 end do 995 end do 996 end do 997c 998 cct_so_a4b = sum 999c 1000 return 1001 end 1002 double precision function cct_so_b4b(rint,nbasis,t2,occ,virt,fac, 1003 1 wijka,wijka2) 1004c 1005c term B4B 1006c t(km,ab)(cm||ij) t(in,cb) (an||kj) 1007c = t(km,ab)(an||kj) (cm||ij) t(in,cb) 1008c 1009c weighting 1. 1010c 1011 implicit none 1012c 1013 integer nbasis,occ,virt 1014 double precision rint(nbasis,nbasis,nbasis,nbasis) 1015 double precision t2(occ ,occ,virt,virt) 1016 double precision fac(nbasis) 1017c 1018 integer m,n,j,b,k,a,i,c 1019 double precision wijka(virt,occ,occ,occ), 1020 1 wijka2(virt,occ,occ,occ),sum 1021c 1022c... (ka) .. t*rint (left) t(km,ab)(an||kj => (m,n,j,b) 1023c 1024 do m=1,occ 1025 do n=1,occ 1026 do j=1,occ 1027 do b=1,virt 1028 wijka(b,j,n,m) = 0.0d0 1029 do k=1,occ 1030 do a=1,virt 1031 wijka(b,j,n,m) = wijka(b,j,n,m) + t2(k,m,a,b) 1032 1 * rint(a+occ,n,k,j) 1033 2 * fac(j)*fac(b+occ)*fac(k)*fac(a+occ) 1034 end do 1035 end do 1036 end do 1037 end do 1038 end do 1039 end do 1040c 1041c... (nb) .. (m,n,j,b) t(in,cb) => (m,j,i,c) 1042c 1043 do m=1,occ 1044 do j=1,occ 1045 do i=1,occ 1046 do c=1,virt 1047 wijka2(c,i,j,m) = 0.0d0 1048 do n=1,occ 1049 do b=1,virt 1050 wijka2(c,i,j,m) = wijka2(c,i,j,m) + wijka(b,j,n,m) 1051 1 * t2(i,n,c,b) 1052 2 * fac(i)*fac(c+occ) 1053 end do 1054 end do 1055 end do 1056 end do 1057 end do 1058 end do 1059c 1060c... (mjic) .. (m,j,i,c) ((cm||ij) => sum 1061c 1062 sum = 0.0d0 1063 do m=1,occ 1064 do j=1,occ 1065 do i=1,occ 1066 do c=1,virt 1067 sum = sum + wijka2(c,i,j,m) * rint(c+occ,m,i,j) 1068 end do 1069 end do 1070 end do 1071 end do 1072c 1073 cct_so_b4b = sum 1074c 1075 return 1076 end 1077 double precision function cct_so_a3a(rint,nbasis,t2,occ,virt,fac, 1078 1 wabcd,wiabc) 1079c 1080c term A3A 1081c 1082c t(ij,ce) (ab||ek) t(ij,af) (cb||fk) 1083c = t(ij,ce) t(ij,af) (ab||ek) (cb||fk) 1084c 1085c weighting 0.5 1086c 1087 implicit none 1088c 1089 integer nbasis,occ,virt 1090 double precision rint(nbasis,nbasis,nbasis,nbasis) 1091 double precision t2(occ,occ,virt,virt) 1092 double precision fac(nbasis) 1093c 1094 integer c,e,a,f,i,j,k,b 1095 double precision wabcd(virt,virt,virt,virt), 1096 1 wiabc(virt,virt,virt,occ),sum 1097c 1098c... (ij) t(ij,ce) t(ij,af) => (c,e,a,f) 1099c 1100 do c=1,virt 1101 do e=1,virt 1102 do a=1,virt 1103 do f=1,virt 1104 wabcd(f,a,e,c) = 0.0d0 1105 do i=1,occ 1106 do j=1,occ 1107 wabcd(f,a,e,c) = wabcd(f,a,e,c) + t2(i,j,c,e)*t2(i,j,a,f) 1108 2 * fac(i)*fac(j)*fac(a+occ)*fac(c+occ) 1109 end do 1110 end do 1111 end do 1112 end do 1113 end do 1114 end do 1115c 1116c... (ae) (c,e,a,f) (ab||ek) => (k,c,f,b) 1117c 1118 do k=1,occ 1119 do c=1,virt 1120 do f=1,virt 1121 do b=1,virt 1122 wiabc(b,f,c,k) = 0.0d0 1123 do a=1,virt 1124 do e=1,virt 1125 wiabc(b,f,c,k) = wiabc(b,f,c,k) + wabcd(f,a,e,c) 1126 1 * rint(a+occ,b+occ,e+occ,k) 1127 2 * fac(k)*fac(b+occ) 1128 end do 1129 end do 1130 end do 1131 end do 1132 end do 1133 end do 1134c 1135c.. (kcfb) (k,c,f,b) (cb||fk) => sum 1136c 1137 sum = 0.0d0 1138 do k=1,occ 1139 do c=1,virt 1140 do f=1,virt 1141 do b=1,virt 1142 sum = sum + wiabc(b,f,c,k)*rint(c+occ,b+occ,f+occ,k) 1143 end do 1144 end do 1145 end do 1146 end do 1147c 1148 cct_so_a3a = sum 1149c 1150 return 1151 end 1152 double precision function cct_so_a3b(rint,nbasis,t2,occ,virt,fac, 1153 1 wiabc,wijka) 1154c 1155c term A3B = B3A 1156c 1157c t(ij,ce) (ab||ek) t(kn,cb) (an||ij) 1158c = (ab||ek) t(kn,cb) t(ij,ce) (an||ij) 1159c 1160c B3A = t(km,ab) (cm||ij) t(ij,af) (cb||fk) 1161c relabel kn cb an ij ij ce ab ek 1162c weighting 0.5 1163c 1164 implicit none 1165c 1166 integer nbasis,occ,virt 1167 double precision rint(nbasis,nbasis,nbasis,nbasis) 1168 double precision t2(occ,occ,virt,virt) 1169 double precision fac(nbasis) 1170c 1171 integer n,c,a,e,k,b,i,j 1172 double precision wiabc(virt,virt,virt,occ), 1173 1 wijka(virt,occ,occ,occ),sum 1174c 1175c... (kb) t(kn,cb) (ab||ek) => (n,c,a,e) 1176c 1177 do n=1,occ 1178 do c=1,virt 1179 do a=1,virt 1180 do e=1,virt 1181 wiabc(e,a,c,n) = 0.0d0 1182 do k=1,occ 1183 do b=1,virt 1184 wiabc(e,a,c,n) = wiabc(e,a,c,n) + t2(k,n,c,b) 1185 1 * rint(a+occ,b+occ,e+occ,k) 1186 2 * fac(c+occ)*fac(a+occ)*fac(k)*fac(b+occ) 1187 end do 1188 end do 1189 end do 1190 end do 1191 end do 1192 end do 1193c 1194c... (ce) (n,c,a,e) t(i,j,c,e) => (n,i,j,a) 1195c 1196 do n=1,occ 1197 do i=1,occ 1198 do j=1,occ 1199 do a=1,virt 1200 wijka(a,j,i,n) = 0.0d0 1201 do c=1,virt 1202 do e=1,virt 1203 wijka(a,j,i,n) = wijka(a,j,i,n) + wiabc(e,a,c,n)*t2(i,j,c,e) 1204 2 * fac(i)*fac(j) 1205 end do 1206 end do 1207 end do 1208 end do 1209 end do 1210 end do 1211c 1212c... (nija) (n,i,j,a) (an||ij) => sum 1213c 1214 sum = 0.0d0 1215 do n=1,occ 1216 do i=1,occ 1217 do j=1,occ 1218 do a=1,virt 1219 sum = sum + wijka(a,j,i,n) * rint(a+occ,n,i,j) 1220 end do 1221 end do 1222 end do 1223 end do 1224c 1225 cct_so_a3b = sum 1226c 1227 return 1228 end 1229 double precision function cct_so_b3b(rint,nbasis,t2,occ,virt,fac, 1230 1 wijab,wijka) 1231c 1232c term B3B 1233c 1234c t(km,ab) (cm||ij) t(kn,cb) (an||ij) 1235c = t(km,ab) t(kn,cb) (cm||ij) (an||ij) 1236c 1237c weighting 0.5 1238c 1239 implicit none 1240c 1241 integer nbasis,occ,virt 1242 double precision rint(nbasis,nbasis,nbasis,nbasis) 1243 double precision t2(occ,occ,virt,virt) 1244 double precision fac(nbasis) 1245c 1246 integer m,n,j,b,k,a,i,c 1247 double precision wijab(virt,virt,occ,occ), 1248 1 wijka(virt,occ,occ,occ),sum 1249c 1250c... (kb) t(km,ab) t(kn,cb) => (m,n,a,c) 1251c 1252 do m=1,occ 1253 do n=1,occ 1254 do a=1,virt 1255 do c=1,virt 1256 wijab(c,a,n,m) = 0.0d0 1257 do k=1,occ 1258 do b=1,virt 1259 wijab(c,a,n,m) = wijab(c,a,n,m) + t2(k,m,a,b) * t2(k,n,c,b) 1260 2 * fac(a+occ)*fac(c+occ)*fac(k)*fac(b+occ) 1261 end do 1262 end do 1263 end do 1264 end do 1265 end do 1266 end do 1267c 1268c... (mc) (m,n,a,c) (cm||ij) => (n,i,j,a) 1269c 1270 do n=1,occ 1271 do i=1,occ 1272 do j=1,occ 1273 do a=1,virt 1274 wijka(a,j,i,n) = 0.0d0 1275 do m=1,occ 1276 do c=1,virt 1277 wijka(a,j,i,n) = wijka(a,j,i,n) + wijab(c,a,n,m) 1278 1 * rint(c+occ,m,i,j) 1279 2 * fac(i)*fac(j) 1280 end do 1281 end do 1282 end do 1283 end do 1284 end do 1285 end do 1286c 1287c.. (nija) (n,i,j,a) (an||ij) => sum 1288c 1289 sum = 0.0d0 1290 do n=1,occ 1291 do i=1,occ 1292 do j=1,occ 1293 do a=1,virt 1294 sum = sum + wijka(a,j,i,n)*rint(a+occ,n,i,j) 1295 end do 1296 end do 1297 end do 1298 end do 1299c 1300 cct_so_b3b = sum 1301c 1302 return 1303 end 1304 double precision function cct_so_a2a(rint,nbasis,t2,occ,virt,fac, 1305 1 wijab,wijab2) 1306c 1307c term A2A 1308c 1309c t(ij,ce) (ab||ek) t(kj,cf) (ab||fi) 1310c = (ab||ek) (ab||fi) t(ij,ce) t(kj,cf) 1311c 1312c weighting 0.5 1313c 1314c 1315 implicit none 1316#include "errquit.fh" 1317#include "mafdecls.fh" 1318c 1319 integer nbasis,occ,virt 1320 double precision rint(nbasis,nbasis,nbasis,nbasis) 1321 double precision t2(occ,occ,virt,virt) 1322 double precision fac(nbasis) 1323c 1324 integer k,i,e,f,a,b,j,c 1325 double precision wijab(virt,virt,occ,occ), 1326 1 wijab2(virt,virt,occ,occ),sum 1327c 1328c... (ab) (ab||ek) (ab||fi) => (k,i,e,f) 1329c 1330 do k=1,occ 1331 do i=1,occ 1332 do e=1,virt 1333 do f=1,virt 1334 wijab(f,e,i,k) = 0.0d0 1335 do a=1,virt 1336 do b=1,virt 1337 wijab(f,e,i,k) = wijab(f,e,i,k)+rint(a+occ,b+occ,e+occ,k) 1338 1 *rint(a+occ,b+occ,f+occ,i) 1339 2 * fac(k)*fac(i)*fac(a+occ)*fac(b+occ) 1340 if (.not. ma_verify_allocator_stuff()) call errquit('a',0, 1341 & MA_ERR) 1342 end do 1343 end do 1344 end do 1345 end do 1346 end do 1347 end do 1348c 1349c... (ie) (k,i,e,f) t(ij,ce) => (k,j,c,f) 1350c 1351 do k=1,occ 1352 do j=1,occ 1353 do c=1,virt 1354 do f=1,virt 1355 wijab2(f,c,j,k) = 0.0d0 1356 do i=1,occ 1357 do e=1,virt 1358 wijab2(f,c,j,k) = wijab2(f,c,j,k) + wijab(f,e,i,k) 1359 1 * t2(i,j,c,e) 1360 2 * fac(j)*fac(c+occ) 1361 if (.not. ma_verify_allocator_stuff()) call errquit('b',0, 1362 & MA_ERR) 1363 end do 1364 end do 1365 end do 1366 end do 1367 end do 1368 end do 1369c 1370c.. (kjcf) (k,j,c,f) t(kj,cf) => sum 1371c 1372 sum = 0.0d0 1373 do k=1,occ 1374 do j=1,occ 1375 do c=1,virt 1376 do f=1,virt 1377 sum = sum + wijab2(f,c,j,k) * t2(k,j,c,f) 1378 end do 1379 end do 1380 end do 1381 end do 1382c 1383 cct_so_a2a = sum 1384c 1385 return 1386 end 1387 double precision function cct_so_a2b(rint,nbasis,t2,occ,virt,fac, 1388 1 wijka,wijka2) 1389c 1390c term A2B=B2A 1391c 1392c t(ij,ce) (ab||ek) t(in,ab) (cn||kj) 1393c = (ab||ek) t(in,ab) t(ij,ce) (cn||kj) 1394c B2A = t(km,ab) (cm||ij) t(kj,cf) (ab||fi) 1395c rename in ab cn kj ij ce ab ek 1396c 1397c weighting 0.5 1398c 1399 implicit none 1400c 1401 integer nbasis,occ,virt 1402 double precision rint(nbasis,nbasis,nbasis,nbasis) 1403 double precision t2(occ,occ,virt,virt) 1404 double precision fac(nbasis) 1405c 1406 integer i,n,k,e,a,b,j,c 1407 double precision wijka(virt,occ,occ,occ), 1408 1 wijka2(virt,occ,occ,occ),sum 1409c 1410c... (ab) t(in,ab) (ab||ek) => (i,n,k,e) 1411c 1412 do i=1,occ 1413 do n=1,occ 1414 do k=1,occ 1415 do e=1,virt 1416 wijka(e,k,n,i) = 0.0d0 1417 do a=1,virt 1418 do b=1,virt 1419 wijka(e,k,n,i) = wijka(e,k,n,i) + t2(i,n,a,b) 1420 1 * rint(a+occ,b+occ,e+occ,k) 1421 2 * fac(i)*fac(k)*fac(a+occ)*fac(b+occ) 1422 end do 1423 end do 1424 end do 1425 end do 1426 end do 1427 end do 1428c 1429c... (ie) (i,n,k,e) t(ij,ce) => (n,k,j,c) 1430c 1431 do n=1,occ 1432 do k=1,occ 1433 do j=1,occ 1434 do c=1,virt 1435 wijka2(c,j,k,n) = 0.0d0 1436 do i=1,occ 1437 do e=1,virt 1438 wijka2(c,j,k,n) = wijka2(c,j,k,n) + wijka(e,k,n,i) 1439 1 * t2(i,j,c,e) 1440 2 * fac(j)*fac(c+occ) 1441 end do 1442 end do 1443 end do 1444 end do 1445 end do 1446 end do 1447c 1448c.. (nkjc) (n,k,j,c) (cn||kj) => sum 1449c 1450 sum = 0.0d0 1451 do n=1,occ 1452 do k=1,occ 1453 do j=1,occ 1454 do c=1,virt 1455 sum = sum + wijka2(c,j,k,n)*rint(c+occ,n,k,j) 1456 end do 1457 end do 1458 end do 1459 end do 1460c 1461 cct_so_a2b = sum 1462c 1463 return 1464 end 1465 double precision function cct_so_b2b(rint,nbasis,t2,occ,virt,fac, 1466 1 wijka,wijkl) 1467c 1468c term B2B 1469c 1470c t(km,ab) (cm||ij) t(in,ab) (cn||kj) 1471c = t(km,ab) t(in,ab) (cm||ij) (cn||kj) 1472c 1473c weighting 0.5 1474c 1475 implicit none 1476c 1477 integer nbasis,occ,virt 1478 double precision rint(nbasis,nbasis,nbasis,nbasis) 1479 double precision t2(occ,occ,virt,virt) 1480 double precision fac(nbasis) 1481c 1482 integer k,m,i,n,a,b,j,c 1483 double precision wijka(virt,occ,occ,occ), 1484 1 wijkl(occ,occ,occ,occ),sum 1485c 1486c... (ab) t(km,ab) t(in,ab) => (k,m,i,n) 1487c 1488 do k=1,occ 1489 do m=1,occ 1490 do i=1,occ 1491 do n=1,occ 1492 wijkl(n,i,m,k) = 0.0d0 1493 do a=1,virt 1494 do b=1,virt 1495 wijkl(n,i,m,k) = wijkl(n,i,m,k) + t2(k,m,a,b)*t2(i,n,a,b) 1496 2 * fac(k)*fac(i)*fac(a+occ)*fac(b+occ) 1497 end do 1498 end do 1499 end do 1500 end do 1501 end do 1502 end do 1503c 1504c... (mi) (k,m,i,n) (cm||ij) => (k,n,j,c) 1505c 1506 do k=1,occ 1507 do n=1,occ 1508 do j=1,occ 1509 do c=1,virt 1510 wijka(c,j,n,k) = 0.0d0 1511 do m=1,occ 1512 do i=1,occ 1513 wijka(c,j,n,k) = wijka(c,j,n,k) + wijkl(n,i,m,k) 1514 1 * rint(c+occ,m,i,j) 1515 2 * fac(j)*fac(c+occ) 1516 end do 1517 end do 1518 end do 1519 end do 1520 end do 1521 end do 1522c 1523c.. (knjc) (k,n,j,c) (cn||kj) => sum 1524c 1525 sum = 0.0d0 1526 do k=1,occ 1527 do n=1,occ 1528 do j=1,occ 1529 do c=1,virt 1530 sum = sum + wijka(c,j,n,k)*rint(c+occ,n,k,j) 1531 end do 1532 end do 1533 end do 1534 end do 1535c 1536 cct_so_b2b = sum 1537c 1538 return 1539 end 1540 double precision function cct_so_a1a(rint,nbasis,t2,occ,virt,fac, 1541 1 wab,wab2) 1542c 1543c term A1A 1544c *seems might be combined with A2A* 1545c 1546c t(ij,ce) (ab||ek) t(ij,cf) (ab||fk) 1547c = (ab||ek) (ab||fk) t(ij,ce) t(ij,cf) 1548c 1549c weighting 0.25 1550c 1551c 1552 implicit none 1553c 1554 integer nbasis,occ,virt 1555 double precision rint(nbasis,nbasis,nbasis,nbasis) 1556 double precision t2(occ,occ,virt,virt) 1557 double precision fac(nbasis) 1558c 1559 integer e,f,k,a,b,i,j,c 1560 double precision wab(virt,virt), 1561 1 wab2(virt,virt),sum,ddot 1562c 1563c... (abk) (ab||ek) (ab||fk) => (e,f) 1564c 1565 do e=1,virt 1566 do f=1,virt 1567 wab(f,e) = 0.0d0 1568 do k=1,occ 1569 do a=1,virt 1570 do b=1,virt 1571 wab(f,e) = wab(f,e) + rint(a+occ,b+occ,e+occ,k) 1572 1 * rint(a+occ,b+occ,f+occ,k) 1573 2 * fac(k)*fac(a+occ)*fac(b+occ) 1574 end do 1575 end do 1576 end do 1577 end do 1578 end do 1579c 1580c... (ijc) t(ij,ce) t(ijcf) => (e,f) 1581c 1582 do e=1,virt 1583 do f=1,virt 1584 wab2(f,e) = 0.0d0 1585 do i=1,occ 1586 do j=1,occ 1587 do c=1,virt 1588 wab2(f,e) = wab2(f,e) + t2(i,j,c,e) * t2(i,j,c,f) 1589 2 * fac(i)*fac(j)*fac(c+occ) 1590 end do 1591 end do 1592 end do 1593 end do 1594 end do 1595c 1596 sum = ddot(virt*virt,wab,1,wab2,1) 1597c 1598 cct_so_a1a = sum 1599c 1600 return 1601 end 1602 double precision function cct_so_a1b(rint,nbasis,t2,occ,virt,fac, 1603 1 wia,wia2) 1604c 1605c term A1B=B1A 1606c 1607c t(ij,ce) (ab||ek) t(kn,ab) (cn||ij) 1608c = (ab||ek) t(kn,ab) t(ij,ce) (cn||ij) 1609c 1610c weighting 0.25 1611c 1612 implicit none 1613c 1614 integer nbasis,occ,virt 1615 double precision rint(nbasis,nbasis,nbasis,nbasis) 1616 double precision t2(occ,occ,virt,virt) 1617 double precision fac(nbasis) 1618c 1619 integer n,e,k,a,b,i,j,c 1620 double precision wia(virt,occ), 1621 1 wia2(virt,occ),sum,ddot 1622c 1623c... (abk) t(kn,ab) (ab||ek) => (n,e) 1624c 1625 do n=1,occ 1626 do e=1,virt 1627 wia(e,n) = 0.0d0 1628 do k=1,occ 1629 do a=1,virt 1630 do b=1,virt 1631 wia(e,n) = wia(e,n) + t2(k,n,a,b) * rint(a+occ,b+occ,e+occ,k) 1632 2 * fac(k)*fac(a+occ)*fac(b+occ) 1633 end do 1634 end do 1635 end do 1636 end do 1637 end do 1638c 1639c... (ijc) (cn||ij) t(ij,ce) => (n,e) 1640c 1641 do n=1,occ 1642 do e=1,virt 1643 wia2(e,n) = 0.0d0 1644 do i=1,occ 1645 do j=1,occ 1646 do c=1,virt 1647 wia2(e,n) = wia2(e,n) + t2(i,j,c,e) * rint(c+occ,n,i,j) 1648 2 * fac(i)*fac(j)*fac(c+occ) 1649 end do 1650 end do 1651 end do 1652 end do 1653 end do 1654c 1655 sum = ddot(virt*occ,wia,1,wia2,1) 1656c 1657 cct_so_a1b = sum 1658c 1659 return 1660 end 1661 double precision function cct_so_b1b(rint,nbasis,t2,occ,virt,fac, 1662 1 wij,wij2) 1663c 1664c term B1B 1665c 1666c t(km,ab) (cm||ij) t(kn,ab) (cn||ij) 1667c = t(km,ab) t(kn,ab) (cm||ij) (cn||ij) 1668c 1669c weighting 0.25 1670c 1671 implicit none 1672c 1673 integer nbasis,occ,virt 1674 double precision rint(nbasis,nbasis,nbasis,nbasis) 1675 double precision t2(occ,occ,virt,virt) 1676 double precision fac(nbasis) 1677c 1678 integer m,n,k,a,b,j,i,c 1679 double precision wij(occ,occ), 1680 1 wij2(occ,occ),sum,ddot 1681c 1682c... (abk) t(km,ab) t(kn,ab) => (m,n) 1683c 1684 do m=1,occ 1685 do n=1,occ 1686 wij(n,m) = 0.0d0 1687 do k=1,occ 1688 do a=1,virt 1689 do b=1,virt 1690 wij(n,m) = wij(n,m) + t2(k,m,a,b) * t2(k,n,a,b) 1691 2 * fac(k)*fac(a+occ)*fac(b+occ) 1692 end do 1693 end do 1694 end do 1695 end do 1696 end do 1697c 1698c... (cij) (cn||ij) (cm||ij) => (m,n) 1699c 1700 do m=1,occ 1701 do n=1,occ 1702 wij2(n,m) = 0.0d0 1703 do j=1,occ 1704 do i=1,occ 1705 do c=1,virt 1706 wij2(n,m) = wij2(n,m) + rint(c+occ,n,i,j) * rint(c+occ,m,i,j) 1707 2 * fac(i)*fac(j)*fac(c+occ) 1708 end do 1709 end do 1710 end do 1711 end do 1712 end do 1713c 1714 sum = ddot(occ*occ,wij,1,wij2,1) 1715c 1716 cct_so_b1b = sum 1717c 1718 return 1719 end 1720#if(0) 1721 subroutine w_ii(rint,nbasis,occ,virt, 1722 1 w,con) 1723c 1724c... calculated the integral only contracted intermediate 1725c... (ab||ei)*(ab||fj) => (ij,ef) ('ab') 1726c 1727 implicit none 1728#include "errquit.fh" 1729c 1730 integer nbasis,occ,virt 1731 double precision rint(nbasis,nbasis,nbasis,nbasis) 1732 character*(*) con 1733c 1734 double precision w(occ,occ,virt,virt) 1735c 1736 if (con.eq.'ab') then 1737c 1738 do e=1,virt 1739 do f=1,virt 1740 do i=1,occ 1741 do j=1,occ 1742c 1743c.. wab(occ,occ,virt,virt)=sum(ab)(ab||ei)*(ab||fj) 1744c 1745 w(i,j,e,f) = 0.0d0 1746 do a=1,virt 1747 do b=1,virt 1748 w(i,j,e,f) = rint(a+occ,b+occ,e+occ,i) 1749 1 * rint(a+occ,b+occ,f+occ,j) 1750 end do 1751 end do 1752c 1753 end do 1754 end do 1755 end do 1756 end do 1757c 1758 else if (con.eq.'ij') then 1759c 1760 do e=1,virt 1761 do f=1,virt 1762 do i=1,occ 1763 do j=1,occ 1764c 1765c.. w(occ,occ,virt,virt)=sum(ab)(ab||ei)*(ab||fj) 1766c 1767 w(i,j,e,f) = 0.0d0 1768 do a=1,occ 1769 do b=1,occ 1770 w(i,j,e,f) = rint(a,b,e+occ,i) 1771 1 * rint(a,b,f+occ,j) 1772 end do 1773 end do 1774c 1775 end do 1776 end do 1777 end do 1778 end do 1779c 1780 else 1781 call errquit(' what ?',0, UNKNOWN_ERR) 1782 end if 1783c 1784 end 1785 subroutine cct_tr_ij_vv_kc(rint,nbasis,t2,occ,virt,w) 1786c 1787c... calculated t2(ij,vv) (vv//kc) contracted intermediate 1788c 1789 implicit none 1790c 1791 integer nbasis,occ,virt 1792 double precision rint(nbasis,nbasis,nbasis,nbasis) 1793 double precision t2(occ,occ,virt,virt) 1794 double precision w(i,j,k,c) 1795c 1796c 1797 do i=1,occ 1798 do j=1,occ 1799 do k=1,occ 1800 do c=1,virt 1801c 1802 w(i,j,k,c) = 0.0d0 1803 do a=1,virt 1804 do b=1,virt 1805 w(i,j,k,c) = rint(a+occ,b+occ,k,c+occ) 1806 1 * t2(i,j,a,b) 1807 end do 1808 end do 1809c 1810 end do 1811 end do 1812 end do 1813 end do 1814c 1815 end 1816 subroutine cct_rt_avbo_oicv(rint,nbasis,t2,occ,virt,wiabc) 1817c 1818c... calculated intermediate 1819c 1820 implicit none 1821c 1822 integer nbasis,occ,virt 1823 double precision rint(nbasis,nbasis,nbasis,nbasis) 1824 double precision t2(occ,occ,virt,virt) 1825 double precision w(i,j,k,c) 1826c 1827 integer i,a,b,c,v,o 1828c 1829 do i=1,occ 1830 do a=1,virt 1831 do b=1,virt 1832 do c=1,virt 1833c 1834 w(i,a,b,c) = 0.0d0 1835 do o=1,occ 1836 do v=1,virt 1837 w(i,a,b,c) = rint(a+occ,v+occ,b+occ,o) 1838 1 * t2(o,i,c,v) 1839 end do 1840 end do 1841c 1842 end do 1843 end do 1844 end do 1845 end do 1846c 1847 end 1848#endif 1849