1 function i2addr(iorb,jorb,korb,lorb,ijklof,noccsym,no12sym) 2* 3* obtain address of integral (iorb jorb ! korb lorb) in molcas order 4* iorb jorb korb lorb corresponds to symmetry ordered indeces !! 5* integrals assumed in core 6* 7* ITRA_ROUTE switch added May 2011 8* 9* 10 implicit real*8(a-h,o-z) 11* 12 include 'mxpdim.inc' 13 include 'orbinp.inc' 14 include 'lucinp.inc' 15 include 'cintfo.inc' 16 include 'multd2h.inc' 17 INCLUDE 'crun.inc' 18#include "errquit.fh" 19#include "mafdecls.fh" 20#include "global.fh" 21 22* 23 dimension ijklof(nsmob,nsmob,nsmob) 24 logical isymj,ksyml,isymk,jsyml,ijsymkl,iksymjl 25*. 26 ntest = 00 27 noccsym_l = noccsym 28 if (i12s.eq.0.and.i34s.eq.0) noccsym_l = 1 29* 30 iabs = iorb 31 ism = ismfto(ireost(iorb)) 32 ioff = ibso(ism) 33* 34 jabs = jorb 35 jsm = ismfto(ireost(jorb)) 36 joff = ibso(jsm) 37* 38 kabs = korb 39 ksm = ismfto(ireost(korb)) 40 koff = ibso(ksm) 41* 42 labs = lorb 43 lsm = ismfto(ireost(lorb)) 44 loff = ibso(lsm) 45* 46 if( ntest.ge. 100) then 47 write(6,*) ' gmijkl at your service ' 48 write(6,*) ' iorb iabs ism ioff ',iorb,iabs,ism,ioff 49 write(6,*) ' jorb jabs jsm joff ',jorb,jabs,jsm,joff 50 write(6,*) ' korb kabs ksm koff ',korb,kabs,ksm,koff 51 write(6,*) ' lorb labs lsm loff ',lorb,labs,lsm,loff 52 end if 53* 54c test 55 ijsm = multd2h(ism,jsm) 56 klsm = multd2h(ksm,lsm) 57 ijklsm = multd2h(ijsm,klsm) 58 if (ijklsm.ne.1) then 59 i2addr = -1 60 return 61 end if 62c test end 63* 64 if (noccsym_l.eq.0.and. 65 & (jsm.gt.ism .or. ( ism.eq.jsm .and. jabs.gt.iabs))) then 66 isym=jsm 67 jsym=ism 68 i = jabs - joff + 1 69 j = iabs - ioff + 1 70 else 71 isym=ism 72 jsym=jsm 73 i = iabs - ioff + 1 74 j = jabs - joff + 1 75 end if 76 if (noccsym_l.eq.0) then 77 ijblk=jsym+isym*(isym-1)/2 78 else 79 ijblk = (isym-1)*nsmob + jsym 80 end if 81 if ( noccsym_l.eq.0 .and. 82 & (lsm.gt.ksm .or. ( ksm.eq.lsm .and. labs.gt.kabs)) ) then 83 ksym=lsm 84 lsym=ksm 85 k = labs -loff + 1 86 l = kabs - koff + 1 87 else 88 ksym=ksm 89 lsym=lsm 90 k = kabs - koff + 1 91 l = labs -loff + 1 92 end if 93 if ( noccsym_l.eq.0) then 94 klblk=lsym+ksym*(ksym-1)/2 95 else 96 klblk = (ksym-1)*nsmob + lsym 97 end if 98* 99 if ( klblk.gt.ijblk .and. no12sym.eq.0 ) then 100 itemp=isym 101 isym=ksym 102 ksym=itemp 103 itemp=jsym 104 jsym=lsym 105 lsym=itemp 106 itemp=ijblk 107 ijblk=klblk 108 klblk=itemp 109* 110 itemp = i 111 i = k 112 k = itemp 113 itemp = j 114 j = l 115 l = itemp 116 end if 117 if(ntest .ge. 100 ) then 118 write(6,*) ' i j k l ',i,j,k,l 119 write(6,*) ' isym,jsym,ksym,lsym',isym,jsym,ksym,lsym 120 end if 121* 122* define offset for given symmetry block 123 ibloff = ijklof(isym,jsym,ksym) 124 if(ntest .ge. 100 ) 125 &write(6,*) ' ibloff isym jsym ksym ', ibloff,isym,jsym,ksym 126 if (noccsym_l.eq.0) then 127 isymj=isym.eq.jsym 128 ksyml=ksym.eq.lsym 129 isymk=(no12sym.eq.0).and.isym.eq.ksym 130 jsyml=(no12sym.eq.0).and.jsym.eq.lsym 131 ijsymkl=isymj.and.ksyml 132 iksymjl=isymk.and.jsyml 133 else 134 isymj=.false. 135 ksyml=.false. 136 isymk=.false. 137 jsyml=.false. 138 ijsymkl=.false. 139 iksymjl=(no12sym.eq.0).and.(isym.eq.ksym).and.(jsym.eq.lsym) 140 end if 141* 142 itorb=ntoobs(isym) 143 jtorb=ntoobs(jsym) 144 ktorb=ntoobs(ksym) 145 ltorb=ntoobs(lsym) 146c? print *,' itorb,jtorb,ktorb,ltorb',itorb,jtorb,ktorb,ltorb 147 if ( isymj ) then 148 ijpairs=itorb*(itorb+1)/2 149 ij=j+i*(i-1)/2 150 else 151 ijpairs=itorb*jtorb 152 IF(ITRA_ROUTE.EQ.1) THEN 153 ij=j + (i-1)*jtorb 154 ELSE 155 ij=i + (j-1)*itorb 156 END IF 157 end if 158* 159 if(ksyml ) then 160 klpairs=ktorb*(ktorb+1)/2 161 kl=l+k*(k-1)/2 162 else 163 klpairs=ktorb*ltorb 164 IF(ITRA_ROUTE.EQ.1) THEN 165 kl=l+(k-1)*ltorb 166 ELSE 167 KL = K + (L-1)*KTORB 168 END IF 169 end if 170c? print *,' ijpairs,klpairs',ijpairs,klpairs 171* 172 if ( iksymjl ) then 173 if ( ij.gt.kl ) then 174 kloff=kl+(kl-1)*(kl-2)/2-1 175 ijkl=ij+(kl-1)*ijpairs-kloff 176 else 177 ijoff=ij+(ij-1)*(ij-2)/2-1 178 ijkl=kl+(ij-1)*klpairs-ijoff 179 end if 180 else 181 ijkl=ij+(kl-1)*ijpairs 182 end if 183 if( ntest .ge. 100 ) 184 & write(6,*) ' ijkl ', ijkl 185* 186 i2addr = ibloff-1+ijkl 187 if ( ntest .ge. 100 ) then 188 write(6,*) 'i j k l ', i,j,k,l 189 write(6,*) ' ibloff ijkl ',ibloff,ijkl 190 write(6,*) ' i2addr = ', i2addr 191 end if 192* 193 return 194 end 195* 196 function i2addr2(iorb,jorb,korb,lorb,ijklof,i12,i34,i1234) 197* 198* obtain address of 4index quantity (iorb,jorb;korb,lorb) 199* iorb jorb korb lorb corresponds to symmetry ordered indeces !! 200* 201* i12 (0,-1,1) symmetry wrt permutation of iorb and jorb 202* i34 (0,-1,1) symmetry wrt permutation of korb and lorb 203* i1234 (0,-1,1) symmetry wrt permutation of pairs iorb,jorb and korb,lorb 204* 205* ITRA_ROUTE switch added May 2011 206* 207 implicit real*8(a-h,o-z) 208* 209 include 'mxpdim.inc' 210 include 'orbinp.inc' 211 include 'lucinp.inc' 212 include 'cintfo.inc' 213 include 'multd2h.inc' 214 include 'frorbs.inc' 215 INCLUDE 'crun.inc' 216 217* 218 dimension ijklof(nsmob,nsmob,nsmob) 219 logical isymj,ksyml,isymk,jsyml,ilsymjk,iksymjl,swapij,swapkl 220*. 221 ntest = 00 222* 223 iabs = iorb 224 isym = ismfto(ireost(iorb)) 225 ioff = ibso(isym) 226 i = iabs - ioff + 1 - nfrobs(isym) 227* 228 jabs = jorb 229 jsym = ismfto(ireost(jorb)) 230 joff = ibso(jsym) 231 j = jabs - joff + 1 - nfrobs(jsym) 232* 233 kabs = korb 234 ksym = ismfto(ireost(korb)) 235 koff = ibso(ksym) 236 k = kabs - koff + 1 - nfrobs(ksym) 237* 238 labs = lorb 239 lsym = ismfto(ireost(lorb)) 240 loff = ibso(lsym) 241 l = labs -loff + 1 - nfrobs(lsym) 242* 243 if (i.le.0.or.j.le.0.or.k.le.0.or.l.le.0) then 244 i2addr2 = -1 245 return 246 end if 247 248* 249 if( ntest.ge. 100) then 250 write(6,*) ' gmijkl at your service ' 251 write(6,*) ' iorb iabs isym ioff ',iorb,iabs,isym,ioff 252 write(6,*) ' jorb jabs jsym joff ',jorb,jabs,jsym,joff 253 write(6,*) ' korb kabs ksym koff ',korb,kabs,ksym,koff 254 write(6,*) ' lorb labs lsym loff ',lorb,labs,lsym,loff 255 end if 256* 257c test 258 ijsym = multd2h(isym,jsym) 259 klsym = multd2h(ksym,lsym) 260 ijklsym = multd2h(ijsym,klsym) 261 if (ijklsym.ne.1) then 262 if (ntest.ge.100) 263 & write(6,*) 'WARNING: symmetry error in i2addr2' 264 i2addr2 = -2 265 return 266 end if 267c test end 268 269 if(ntest .ge. 100 ) then 270 write(6,*) ' as entered:' 271 write(6,*) ' i j k l ',i,j,k,l 272 write(6,*) ' isym,jsym,ksym,lsym',isym,jsym,ksym,lsym 273 end if 274* 275* resort index quadruple to standard sequence 276* 277 ijblk = (min(isym,jsym)-1)*nsmob + max(isym,jsym) 278 klblk = (min(ksym,lsym)-1)*nsmob + max(ksym,lsym) 279 280 ijdx = (min(iabs,jabs)-1)*ntoob + max(iabs,jabs) 281 kldx = (min(kabs,labs)-1)*ntoob + max(kabs,labs) 282 if ( i1234.ne.0 .and. 283 & (ijblk.gt.klblk .or. 284 & (ijblk.eq.klblk .and. ijdx.gt.kldx) ) ) then 285 itemp = isym 286 isym = ksym 287 ksym = itemp 288 itemp = jsym 289 jsym = lsym 290 lsym = itemp 291 itemp = i 292 i = k 293 k = itemp 294 itemp = j 295 j = l 296 l = itemp 297 end if 298 299 swapij = ( i12.ne.0 .and. (isym.gt.jsym .or. 300 & (isym.eq.jsym .and. i.gt.j) ) ) 301 swapkl = ( i34.ne.0 .and. (ksym.gt.lsym .or. 302 & (ksym.eq.lsym .and. k.gt.l) ) ) 303 swapij = swapij.or.( i34.ne.0 .and. i12.eq.0 .and. swapkl ) 304 swapkl = swapkl.or.( i12.ne.0 .and. i34.eq.0 .and. swapij ) 305 306 swapij = swapij.or.( i34.ne.0.and.i12.eq.0 .and. 307 & ksym.eq.lsym .and. k.eq.l 308 & .and.(isym.gt.jsym .or.(isym.eq.jsym .and. i.gt.j))) 309 310 swapkl = swapkl.or.( i12.ne.0.and.i34.eq.0 .and. 311 & isym.eq.jsym .and. i.eq.j 312 & .and.(ksym.gt.lsym .or.(ksym.eq.lsym .and. k.gt.l))) 313 314 if (swapij) then 315 itemp = isym 316 isym = jsym 317 jsym = itemp 318 itemp = i 319 i = j 320 j = itemp 321 end if 322 323 if (swapkl) then 324 itemp = ksym 325 ksym = lsym 326 lsym = itemp 327 itemp = k 328 k = l 329 l = itemp 330 end if 331* 332 if(ntest .ge. 100 ) then 333 write(6,*) ' resorted to:' 334 write(6,*) ' i j k l ',i,j,k,l 335 write(6,*) ' isym,jsym,ksym,lsym',isym,jsym,ksym,lsym 336 end if 337* 338* offset for given symmetry block 339 ibloff = ijklof(isym,jsym,ksym) 340 341 if(ntest .ge. 100 ) 342 & write(6,*) ' ibloff isym jsym ksym ', ibloff,isym,jsym,ksym 343 344 isymj=(i12.ne.0).and.isym.eq.jsym 345 ksyml=(i34.ne.0).and.ksym.eq.lsym 346 isymk=(i1234.ne.0).and.isym.eq.ksym 347 jsyml=(i1234.ne.0).and.jsym.eq.lsym 348 ilsymjk=(i12.ne.0.and.i34.eq.0.and.i1234.ne.0).and. 349 & isym.eq.lsym.and.jsym.eq.ksym 350 iksymjl=isymk.and.jsyml 351* 352 itorb=ntaobs(isym) 353 jtorb=ntaobs(jsym) 354 ktorb=ntaobs(ksym) 355 ltorb=ntaobs(lsym) 356c? print *,' itorb,jtorb,ktorb,ltorb',itorb,jtorb,ktorb,ltorb 357 if ( isymj.and.i12.eq.1 ) then 358 ijpairs=itorb*(itorb+1)/2 359 ij=i+j*(j-1)/2 360 ijdia = 1 361 else if ( isymj.and.i12.eq.-1 ) then 362 if (i.eq.j) then 363 i2addr2 = -1 364 return 365 end if 366 ijpairs=itorb*(itorb-1)/2 367 ij=i+(j-1)*(j-2)/2 368 ijdia = -1 369 else 370 ijpairs=itorb*jtorb 371 IF(ITRA_ROUTE.EQ.1) THEN 372 ij=i + (j-1)*itorb 373 ELSE 374 ij=j + (i-1)*jtorb 375 END IF 376 ijdia = 0 377 end if 378* 379 if(ksyml.and.i34.eq.1) then 380 klpairs=ktorb*(ktorb+1)/2 381 kl=k+l*(l-1)/2 382 kldia = 1 383 else if(ksyml.and.i34.eq.-1) then 384 if (k.eq.l) then 385 i2addr2 = -1 386 return 387 end if 388 klpairs=ktorb*(ktorb-1)/2 389 kl=k+(l-1)*(l-2)/2 390 kldia = -1 391 else 392 klpairs=ktorb*ltorb 393 IF(ITRA_ROUTE.EQ.1) THEN 394 kl=k+(l-1)*ktorb 395 ELSE 396 kl=l+(k-1)*ltorb 397 END IF 398 kldia = 0 399 end if 400c? print *,' ijpairs,klpairs',ijpairs,klpairs 401* 402 if ( iksymjl .and. i1234.eq.1 ) then 403 if (ijdia.eq.kldia) then 404 if ( ij.le.kl ) then 405 ijkl=(kl-1)*kl/2 + ij 406 else 407 ijkl=(ij-1)*ij/2 + kl 408 end if 409 else 410 stop 'incomplete i2addr2' 411 end if 412 else if ( iksymjl .and. i1234.eq.-1 ) then 413 if (ijdia.eq.kldia) then 414 if ( ij.lt.kl ) then 415 ijkl=(kl-2)*(kl-1)/2 + ij 416 else if (ij.gt.kl) then 417 ijkl=(ij-2)*(ij-1)/2 + kl 418 else 419 ! not allowed 420 i2addr2 = -1 421 return 422 end if 423 else if (ijdia.eq.1.and.kldia.eq.0) then 424 ! ij is upper-triangular packed 425 ! kl is full matrix 426 ! case 1: 427 ! i.le.j and k.le.l 428 if (k.le.l) then 429 kl_tri = k+l*(l-1)/2 430 if (ij.gt.kl_tri) then 431 ijkl = (ij-1)*(ij-2)/2 + kl_tri 432 else if (ij.lt.kl_tri) then 433 ijkl = (kl_tri-1)*(kl_tri-2)/2 + ij 434 else 435 ! not allowed 436 i2addr2 = -1 437 return 438 end if 439 ! case 2: 440 ! i.lt.j and k.gt.l 441 else if (i.lt.j.and.k.gt.l) then 442 ijkl_off = ijpairs*(ijpairs-1)/2 443 ij_tri = (j-1)*(j-2)/2+i 444 kl_tri = (k-1)*(k-2)/2+l 445 if (ij_tri.le.kl_tri) then 446 ijkl = ijkl_off+kl_tri*(kl_tri-1)/2+ij_tri 447 else 448 ijkl = ijkl_off+ij_tri*(ij_tri-1)/2+kl_tri 449 end if 450 ! forbidden case: 451 else 452 i2addr2 = -1 453 return 454 end if 455 456 else if (ijdia.eq.0.and.kldia.eq.1) then 457 nnmx = itorb 458 stop 'adapt' 459 ! case 1: 460 ! i.le.j and k.le.l 461 if (k.le.l) then 462 kl_tri = l+k*(k-1)/2 463 if (ij.gt.kl_tri) then 464 ijkl = ij*(ij-1)/2 + kl_tri 465 else if (ij.lt.kl_tri) then 466 ijkl = kl_tri*(kl_tri-1)/2 + ij 467 else 468 ! not allowed 469 i2addr2 = -1 470 return 471 end if 472 ! case 2: 473 ! i.lt.j and k.gt.l 474 else if (i.lt.j.and.k.gt.l) then 475 ijkl_off = ijpairs*(ijpairs-1)/2 476 ij_tri = (j-1)*(j-2)/2+i 477 kl_tri = (k-1)*(k-2)/2+l 478 if (ij_tri.le.kl_tri) then 479 ijkl = ijkl_off+kl_tri*(kl_tri-1)/2+ij_tri 480 else 481 ijkl = ijkl_off+ij_tri*(ij_tri-1)/2+kl_tri 482 end if 483 ! forbidden case: 484 else 485 i2addr2 = -1 486 return 487 end if 488 489 else 490 stop 'incomplete i2addr2' 491 end if 492 else if (ksym.eq.lsym.and.ijdia.eq.1.and.kldia.eq.0) then 493 if (i.lt.j) then 494 ij_tri = (j-1)*(j-2)/2+i 495 ij_tripairs = itorb*(itorb-1)/2 496 ijkl = ij_tri + (kl-1)*ij_tripairs 497 else if (i.eq.j) then 498 ijkl_off = itorb*(itorb-1)/2*klpairs 499 kl_tri = l*(l-1)/2+k 500 ijkl = ijkl_off + i + (kl_tri-1)*itorb 501 else 502 stop 'unexpected' 503 end if 504 else if (isym.eq.jsym.and.ijdia.eq.0.and.kldia.eq.1) then 505 if (k.lt.l) then 506 kl_tri = (l-1)*(l-2)/2+k 507 kl_tripairs = ktorb*(ktorb-1)/2 508 ijkl = kl_tri + (ij-1)*kl_tripairs 509 else if (k.eq.l) then 510 ijkl_off = ktorb*(ktorb-1)/2*ijpairs 511 ij_tri = j*(j-1)/2+i 512 ijkl = ijkl_off + k + (ij_tri-1)*ktorb 513 else 514 stop 'unexpected' 515 end if 516 517 else if (ilsymjk) then 518 ! i<j, k>l 519 IF(ITRA_ROUTE.EQ.1) THEN 520 ij = (j-1)*itorb + i 521 kl = (k-1)*ltorb + l 522 ELSE 523 ij = (i-1)*jtorb + j 524 kl = (l-1)*ktorb + k 525 END IF 526 527 if (ij.le.kl) then 528 ijkl = kl*(kl-1)/2+ij 529 else 530 ijkl = ij*(ij-1)/2+kl 531 end if 532 else 533 ijkl=ij+(kl-1)*ijpairs 534 end if 535 if( ntest .ge. 100 ) 536 & write(6,*) ' ijkl ', ijkl 537* 538 i2addr2 = ibloff-1+ijkl 539 if( ntest .ge. 100 ) then 540 write(6,*) 'i j k l ', i,j,k,l 541 write(6,*) ' ibloff ijkl ',ibloff,ijkl 542 write(6,*) ' i2addr2 = ', i2addr2 543 end if 544* 545 return 546 end 547 548 subroutine ijkl2iadr(ijkl,iadr,nadr,ntoob,ireost, 549 & ijklof,i12,i34,i1234) 550c 551c return the orbital quadruples for a list of addresses listed on iadr() 552c 553 implicit none 554 555 integer, intent(in) :: 556 & nadr, iadr(nadr), ntoob, 557 & i12, i34, i1234, ijklof(*), ireost(*) 558 559 integer, intent(out) :: 560 & ijkl(4,nadr) 561 562 integer :: 563 & idx, jdx, kdx, ldx, ii, iad 564 565 integer, external :: 566 & i2addr2 567 568 ijkl(1:4,1:nadr) = 0 569 570 do idx = 1, ntoob 571 do jdx = 1, ntoob 572 do kdx = 1, ntoob 573 do ldx = 1, ntoob 574 iad = i2addr2(idx,jdx,kdx,ldx,ijklof,i12,i34,i1234) 575 if (iad.lt.0) cycle 576 do ii = 1, nadr 577 if (iad.eq.iadr(ii).and.ijkl(1,ii).eq.0) then 578 ijkl(1,ii) = ireost(idx) 579 ijkl(2,ii) = ireost(jdx) 580 ijkl(3,ii) = ireost(kdx) 581 ijkl(4,ii) = ireost(ldx) 582 end if 583 end do 584 end do 585 end do 586 end do 587 end do 588 589 return 590 591 end 592 FUNCTION GTIJKL(I,J,K,L) 593* 594* Obtain integral (I J ! K L ) 595* 596* reads oper and will get similarity transformed integrals if 597* I_USE_SIMTRH==1 598* if i_unrorb.eq.1 it will obey ispcas (both on oper.inc) 599* 600* I,J,K L refers to orbitals in Type ordering 601* ============== 602* 603* from GTIJKL_SM_AB, replacing old GTIJKL 604* 605c IMPLICIT REAL*8(A-H,O-Z) 606c INCLUDE 'mxpdim.inc' 607 INCLUDE 'wrkspc.inc' 608#include "errquit.fh" 609#include "mafdecls.fh" 610#include "global.fh" 611 INCLUDE 'glbbas.inc' 612 INCLUDE 'lucinp.inc' 613 INCLUDE 'orbinp.inc' 614 INCLUDE 'crun.inc' 615 INCLUDE 'oper.inc' 616* 617 XIJKL = 0.0D0 618* 619 ! set up offset and pointer array 620 IF (I_USE_SIMTRH.EQ.0) THEN 621 NOCCSYM=0 622 IF (I_UNRORB.EQ.0.OR.ISPCAS.EQ.1) THEN 623 NO12SYM=0 624 K2ADR = KINT2 625 KP2ADR= KPINT2 626 ELSE IF (I_UNRORB.EQ.1.AND.ISPCAS.EQ.2) THEN 627 NO12SYM=0 628 K2ADR = KINT2BB 629 KP2ADR= KPINT2 630 ELSE IF (I_UNRORB.EQ.1.AND.(ISPCAS.EQ.3.OR.ISPCAS.EQ.4)) THEN 631 NO12SYM=1 632 K2ADR = KINT2AB 633 KP2ADR= KPINT2AB 634 ELSE 635 WRITE(6,*) 'unknown case: i_unrorb, ispcas: ',i_unrorb, ispcas 636 STOP 'GTIJKL' 637 END IF 638 ELSE 639 NOCCSYM=1 640 IF (I_UNRORB.EQ.0) THEN 641 NO12SYM=0 642 K2ADR = KINT2_SIMTRH 643 KP2ADR= KPINT2_SIMTRH 644 ELSE IF (I_UNRORB.EQ.1.AND.ISPCAS.EQ.1) THEN 645 NO12SYM=0 646 K2ADR = KINT2_SIMTRH_AA 647 KP2ADR= KPINT2_SIMTRH 648 ELSE IF (I_UNRORB.EQ.1.AND.ISPCAS.EQ.2) THEN 649 NO12SYM=0 650 K2ADR = KINT2_SIMTRH_BB 651 KP2ADR= KPINT2_SIMTRH 652 ELSE IF (I_UNRORB.EQ.1.AND.(ISPCAS.EQ.3.OR.ISPCAS.EQ.4)) THEN 653 NO12SYM=1 654 K2ADR = KINT2_SIMTRH_AB 655 KP2ADR= KPINT2_SIMTRH_AB 656 END IF 657 658 END IF 659 IF (.NOT.(I_UNRORB.EQ.1.AND.ISPCAS.EQ.4)) THEN 660 IADR = I2ADDR(IREOTS(I),IREOTS(J), 661 & IREOTS(K),IREOTS(L), 662 & int_mb(KP2ADR),NOCCSYM,NO12SYM) 663 ELSE 664 IADR = I2ADDR(IREOTS(K),IREOTS(L), 665 & IREOTS(I),IREOTS(J), 666 & int_mb(KP2ADR),NOCCSYM,NO12SYM) 667 END IF 668 669 IF (IADR.GT.0) THEN 670 XIJKL = dbl_mb(K2ADR-1+IADR) 671 ELSE 672 XIJKL = 0D0 673 END IF 674* 675 GTIJKL = XIJKL 676* 677 NTEST = 00 678 IF(NTEST.GE.100) THEN 679 WRITE(6,*) ' 2e integral for I,J,K,L = ', I,J,K,L 680 WRITE(6,*) ' is ', XIJKL 681 END IF 682 RETURN 683 END 684 SUBROUTINE PTIJKL(I,J,K,L,XINT,XLIST) 685* 686* Put integral (I J ! K L ) on XINT to its correct place in XLIST 687* 688* a quick hack: never tested, never used :-)) 689* 690* reads oper and will behave correspondingly 691* 692* I,J,K L refers to active orbitals in Type ordering 693* ============== 694* 695c IMPLICIT REAL*8(A-H,O-Z) 696c INCLUDE 'mxpdim.inc' 697 INCLUDE 'wrkspc.inc' 698#include "errquit.fh" 699#include "mafdecls.fh" 700#include "global.fh" 701 INCLUDE 'glbbas.inc' 702 INCLUDE 'lucinp.inc' 703 INCLUDE 'orbinp.inc' 704 INCLUDE 'crun.inc' 705 INCLUDE 'oper.inc' 706 707 DIMENSION XLIST(*) 708* 709 ! set up offset and pointer array 710 IF (I_USE_SIMTRH.EQ.0) THEN 711 NOCCSYM=0 712 IF (I_UNRORB.EQ.0.OR.ISPCAS.EQ.1) THEN 713 NO12SYM=0 714 KP2ADR= KPINT2 715 ELSE IF (I_UNRORB.EQ.1.AND.ISPCAS.EQ.2) THEN 716 NO12SYM=0 717 KP2ADR= KPINT2 718 ELSE IF (I_UNRORB.EQ.1.AND.(ISPCAS.EQ.3.OR.ISPCAS.EQ.4)) THEN 719 NO12SYM=1 720 KP2ADR= KPINT2AB 721 ELSE 722 WRITE(6,*) 'unknown case: i_unrorb, ispcas: ',i_unrorb, ispcas 723 STOP 'GTIJKL' 724 END IF 725 ELSE 726 NOCCSYM=1 727 IF (I_UNRORB.EQ.0) THEN 728 NO12SYM=0 729 KP2ADR= KPINT2_SIMTRH 730 ELSE IF (I_UNRORB.EQ.1.AND.ISPCAS.EQ.1) THEN 731 NO12SYM=0 732 KP2ADR= KPINT2_SIMTRH 733 ELSE IF (I_UNRORB.EQ.1.AND.ISPCAS.EQ.2) THEN 734 NO12SYM=0 735 KP2ADR= KPINT2_SIMTRH 736 ELSE IF (I_UNRORB.EQ.1.AND.(ISPCAS.EQ.3.OR.ISPCAS.EQ.4)) THEN 737 NO12SYM=1 738 KP2ADR= KPINT2_SIMTRH_AB 739 END IF 740 741 END IF 742 IF (.NOT.(I_UNRORB.EQ.1.AND.ISPCAS.EQ.4)) THEN 743 IADR = I2ADDR(IREOTS(I),IREOTS(J), 744 & IREOTS(K),IREOTS(L), 745 & dbl_mb(KP2ADR),NOCCSYM,NO12SYM) 746 ELSE 747 IADR = I2ADDR(IREOTS(K),IREOTS(L), 748 & IREOTS(I),IREOTS(J), 749 & dbl_mb(KP2ADR),NOCCSYM,NO12SYM) 750 END IF 751 752 IF (IADR.GT.0) THEN 753 XLIST(IADR) = XINT 754 ELSE 755 IF (XINT.GT.1D3*EPSILON(1D0)) THEN 756 WRITE(6,*) 'WARNING: INTEGRAL ',I,J,K,L, 757 & ' IS NOT TOTALLY SYMMETRIC, THUS IGNORED!' 758 WRITE(6,*) 'VALUE = ',XINT 759 END IF 760 END IF 761* 762 RETURN 763 END 764 FUNCTION GTIJKL_GN(IORB,JORB,KORB,LORB) 765* 766* Obtain integral (IORB JORB ! KORB LORB) from current 767* active array of integrals. It is assumed that the 768* current list is a complete integral list 769* 770* IST allows switching between symmetry(1)- and type(2)-ordered orbitals 771* 772* Jeppe Olsen, July 2011 773* 774 INCLUDE 'implicit.inc' 775 INCLUDE 'mxpdim.inc' 776 INCLUDE 'cintfo.inc' 777 INCLUDE 'glbbas.inc' 778 INCLUDE 'wrkspc-static.inc' 779#include "errquit.fh" 780#include "mafdecls.fh" 781#include "global.fh" 782 783* 784 NTEST = 000 785 IF(NTEST.GE.200) THEN 786 WRITE(6,*) ' Output from GT_IJKL_GN ' 787 WRITE(6,'(A,4I4)') ' IORB, JORB, KORB, LORB = ', 788 & IORB, JORB, KORB, LORB 789 END IF 790* 791* 792 KINT2_LA = KINT2_A(IE2ARRAY_A) 793 KPINT2_LA = KPINT2_A(IE2ARRAY_A) 794* 795 IST = 2 796 I2ADDR = I2ADDR_GN(IORB,JORB,KORB,LORB, 797 & int_mb(KPINT2_LA),IST) 798* 799 IF(I2ADDR.GT.0) THEN 800 XINT = dbl_mb(KINT2_LA-1+I2ADDR) 801 ELSE 802 WRITE(6,*) ' Negative integral address ' 803 WRITE(6,*) ' IORB, JORB, KORB, LORB, I2ADDR = ', 804 & IORB, JORB, KORB, LORB, I2ADDR 805 WRITE(6,*) ' IE2ARRAY_A, KPINT2_LA = ', 806 & IE2ARRAY_A, KPINT2_LA 807 STOP ' Negative integral address ' 808 END IF 809* 810 GTIJKL_GN = XINT 811* 812 IF(NTEST.GE.100) THEN 813 WRITE(6,*) ' Output from GT_IJKL_GN ' 814 WRITE(6,'(A,4I4)') ' IORB, JORB, KORB, LORB = ', 815 & IORB, JORB, KORB, LORB 816 WRITE(6,'(A,I8,E15.8)') 817 & 'Address and value of integral ', I2ADDR, GTIJKL_GN 818 END IF 819* 820 RETURN 821 END 822 823 824* 825* Obtain 826 FUNCTION I2ADDR_GN(IORB,JORB,KORB,LORB,IJKLOF,IST) 827* 828* obtain address of integral (iorb jorb ! korb lorb) using 829* new integral order.. 830* IST = 1 => Symmetry ordered orbital indices 831* IST = 2 => Type ordered orbital indices 832* 833* Permutational symmetry is defined by I12S_A,I34S_A,I1234S_A coming through 834* CINTFO 835* 836* 837 INCLUDE 'implicit.inc' 838* 839 include 'mxpdim.inc' 840 include 'orbinp.inc' 841 include 'lucinp.inc' 842 include 'cintfo.inc' 843 include 'multd2h.inc' 844 INCLUDE 'crun.inc' 845 846* 847 dimension ijklof(nsmob,nsmob,nsmob) 848 logical isymj,ksyml,ijsymkl 849*. 850 ntest = 00 851*. Reform to symmetry ordered indices if required 852 IF(IST.EQ.2) THEN 853 IABS = IREOTS(IORB) 854 JABS = IREOTS(JORB) 855 KABS = IREOTS(KORB) 856 LABS = IREOTS(LORB) 857 ELSE 858 IABS = IORB 859 JABS = JORB 860 KABS = KORB 861 LABS = LORB 862 END IF 863* 864 IF(IST.EQ.1) THEN 865 ism = ismfto(ireost(iorb)) 866 ELSE 867 ISM = ISMFTO(IORB) 868 END IF 869 ioff = ibso(ism) 870* 871 IF(IST.EQ.1) THEN 872 jsm = ismfto(ireost(jorb)) 873 ELSE 874 JSM = ISMFTO(JORB) 875 END IF 876 joff = ibso(jsm) 877* 878 IF(IST.EQ.1) THEN 879 ksm = ismfto(ireost(korb)) 880 ELSE 881 KSM = ISMFTO(KORB) 882 END IF 883 koff = ibso(ksm) 884* 885 IF(IST.EQ.1) THEN 886 lsm = ismfto(ireost(lorb)) 887 ELSE 888 LSM = ISMFTO(LORB) 889 END IF 890 loff = ibso(lsm) 891* 892 if( ntest.ge. 100) then 893 write(6,*) ' gmijkl at your service ' 894 write(6,*) ' iorb iabs ism ioff ',iorb,iabs,ism,ioff 895 write(6,*) ' jorb jabs jsm joff ',jorb,jabs,jsm,joff 896 write(6,*) ' korb kabs ksm koff ',korb,kabs,ksm,koff 897 write(6,*) ' lorb labs lsm loff ',lorb,labs,lsm,loff 898* 899 WRITE(6,*) ' I12S_A,I34S_A,I1234S_A = ', 900 & I12S_A,I34S_A,I1234S_A 901 end if 902*. Test symmetry 903 ijsm = multd2h(ism,jsm) 904 klsm = multd2h(ksm,lsm) 905 ijklsm = multd2h(ijsm,klsm) 906 if (ijklsm.ne.1) then 907 i2addr = -1 908 return 909 end if 910* 911 if (I12S_A.EQ.1.and. 912 & (jsm.gt.ism .or. ( ism.eq.jsm .and. jabs.gt.iabs))) then 913 isym=jsm 914 jsym=ism 915 i = jabs - joff + 1 916 j = iabs - ioff + 1 917 else 918 isym=ism 919 jsym=jsm 920 i = iabs - ioff + 1 921 j = jabs - joff + 1 922 end if 923 if (I12S_A.eq.1) then 924 ijblk=jsym+isym*(isym-1)/2 925 else 926 ijblk = (isym-1)*nsmob + jsym 927 end if 928 if ( I34S_A.EQ.1 .and. 929 & (lsm.gt.ksm .or. ( ksm.eq.lsm .and. labs.gt.kabs)) ) then 930 ksym=lsm 931 lsym=ksm 932 k = labs -loff + 1 933 l = kabs - koff + 1 934 else 935 ksym=ksm 936 lsym=lsm 937 k = kabs - koff + 1 938 l = labs -loff + 1 939 end if 940 if (I34S_A .eq.1) then 941 klblk=lsym+ksym*(ksym-1)/2 942 else 943 klblk = (ksym-1)*nsmob + lsym 944 end if 945* 946 if ( klblk.gt.ijblk .and. I1234S_A.EQ.1 ) then 947 itemp=isym 948 isym=ksym 949 ksym=itemp 950 itemp=jsym 951 jsym=lsym 952 lsym=itemp 953 itemp=ijblk 954 ijblk=klblk 955 klblk=itemp 956* 957 itemp = i 958 i = k 959 k = itemp 960 itemp = j 961 j = l 962 l = itemp 963 end if 964 if(ntest .ge. 100 ) then 965 write(6,*) ' i j k l ',i,j,k,l 966 write(6,*) ' isym,jsym,ksym,lsym',isym,jsym,ksym,lsym 967 end if 968* 969* define offset for given symmetry block 970 ibloff = ijklof(isym,jsym,ksym) 971 if(ntest .ge. 100 ) 972 &write(6,*) ' ibloff isym jsym ksym ', ibloff,isym,jsym,ksym 973 isymj=.false. 974 ksyml=.false. 975 ijsymkl=.false. 976 if (I12S_A.eq.1) isymj=isym.eq.jsym 977 IF(I34S_A.EQ.1) KSYML = KSYM.EQ.LSYM 978 IF(I1234S_A.EQ.1) IJSYMKL = (ISYM.EQ.KSYM).AND.(JSYM.EQ.LSYM) 979* 980 itorb=ntoobs(isym) 981 jtorb=ntoobs(jsym) 982 ktorb=ntoobs(ksym) 983 ltorb=ntoobs(lsym) 984c? print *,' itorb,jtorb,ktorb,ltorb',itorb,jtorb,ktorb,ltorb 985 if ( isymj ) then 986 ijpairs=itorb*(itorb+1)/2 987 ij=j+i*(i-1)/2 988 else 989 ijpairs=itorb*jtorb 990 ij=i + (j-1)*itorb 991 end if 992* 993 if(ksyml ) then 994 klpairs=ktorb*(ktorb+1)/2 995 kl=l+k*(k-1)/2 996 else 997 klpairs=ktorb*ltorb 998 KL = K + (L-1)*KTORB 999 end if 1000c? print *,' ijpairs,klpairs',ijpairs,klpairs 1001* 1002 if ( IJSYMKL ) then 1003 if ( ij.gt.kl ) then 1004 kloff=kl+(kl-1)*(kl-2)/2-1 1005 ijkl=ij+(kl-1)*ijpairs-kloff 1006 else 1007 ijoff=ij+(ij-1)*(ij-2)/2-1 1008 ijkl=kl+(ij-1)*klpairs-ijoff 1009 end if 1010 else 1011 ijkl=ij+(kl-1)*ijpairs 1012 end if 1013 if( ntest .ge. 100 ) 1014 & write(6,*) ' ijkl ', ijkl 1015* 1016 I2ADDR_GN = ibloff-1+ijkl 1017 if( ntest .ge. 100 ) then 1018 write(6,*) 'i j k l ', i,j,k,l 1019 write(6,*) ' ibloff ijkl ',ibloff,ijkl 1020 write(6,*) ' I2ADDR_GN = ', I2ADDR_GN 1021 end if 1022* 1023 return 1024 end 1025c $Id$ 1026