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