1      subroutine ccsd_ht2pm(basis,nsh,ncor,nocc,nvir,nact,nbf,g_sht2,
2     &                      g_nht2,cmo,scra,scrb,offsh)
3C     $Id$
4      implicit none
5#include "errquit.fh"
6      integer basis,nsh,ncor,nocc,nvir,nact,nbf,g_sht2,g_nht2,
7     &        offsh(nsh,nsh)
8      double precision cmo(nbf,nbf),scra(nbf*nbf),scrb(nbf*nbf)
9
10#include "mafdecls.fh"
11#include "global.fh"
12#include "bas.fh"
13#include "rtdb.fh"
14#include "ccsd_debug.fh"
15#include "ccsdps.fh"
16c
17      integer g_sht2t,g_jlo,g_jhi,g_ilo,g_ihi,
18     &        ish,ilo,ihi,jsh,jlo,jhi,ksh,klo,khi,lsh,llo,lhi,
19     &        x,y,xy,ymax,nfi,nfj,ipp,imm,ii,jj,iijj,
20     &        jjii,i,j,k,l,ad1,ad2,ad3,iam,lnoo,a,b
21      integer nxtask
22      external nxtask
23c
24      if (occsdps) then
25         call pstat_on(ps_ht2pm)
26      else
27         call qenter('ht2pm',0)
28      endif
29c
30      iam=ga_nodeid()
31      lnoo=nocc*nocc
32c
33      call ga_distribution(g_sht2,iam,g_jlo,g_jhi,g_ilo,g_ihi)
34      do ish=1,nsh
35        if (.not. bas_cn2bfr(basis,ish,ilo,ihi))
36     $      call errquit('vvvv: bas_cn2bfr',ish, BASIS_ERR)
37        nfi=ihi-ilo+1
38        do jsh=1,ish
39          if (.not. bas_cn2bfr(basis,jsh,jlo,jhi))
40     $        call errquit('vvvv: bas_cn2bfr',jsh, BASIS_ERR)
41          nfj=jhi-jlo+1
42          do x=1,nfi
43            ymax=nfj
44            if (ish.eq.jsh)ymax=x
45            do y=1,ymax
46              xy=offsh(ish,jsh)+(x-1)*nfj+y
47              if (xy.ge.g_ilo.and.xy.le.g_ihi)then
48                call ga_get(g_sht2,1,lnoo,xy,xy,scra,lnoo)
49                ipp=0
50                imm=nocc*(nocc+1)/2
51                do ii=1,nocc
52                  do jj=1,ii-1
53                    ipp=ipp+1
54                    imm=imm+1
55                    iijj=(ii-1)*nocc+jj
56                    jjii=(jj-1)*nocc+ii
57                    scrb(iijj)=scra(ipp)+scra(imm)
58                    scrb(jjii)=scra(ipp)-scra(imm)
59                  enddo
60                  iijj=(ii-1)*nocc+ii
61                  ipp=ipp+1
62                  scrb(iijj)=scra(ipp)
63                enddo
64                call ga_put(g_sht2,1,lnoo,xy,xy,scrb,lnoo)
65                if (ish.ne.jsh.or.x.ne.y)then
66                  xy=offsh(jsh,ish)+(y-1)*nfi+x
67                  do ii=1,nocc
68                    do jj=1,nocc
69                     iijj=(ii-1)*nocc+jj
70                     jjii=(jj-1)*nocc+ii
71                     scra(iijj)=scrb(jjii)
72                    enddo
73                  enddo
74                  call ga_put(g_sht2,1,lnoo,xy,xy,scra,lnoo)
75                endif
76              endif
77            enddo
78          enddo
79        enddo
80      enddo
81      call ga_sync()
82c
83c ------------------------------------------------------------
84c - transform ao indices of ht2 array into the mo basis
85c ------------------------------------------------------------
86
87*ga:1:0
88      if (.not.ga_create(MT_DBL,nbf*nbf,lnoo,'sht2t',
89     &                   nbf*nbf,0,g_sht2t))
90     &     call errquit('ga_create g_sht2t failed',0, GA_ERR)
91
92      call ga_transpose(g_sht2,g_sht2t)
93c - redefine g_sht2
94      if (.not.ga_destroy(g_sht2))
95     &    call errquit('ga_dest g_sht2 fail',0, GA_ERR)
96*ga:1:0
97      if (.not.ga_create(MT_DBL,nact*nact,lnoo,'sht2',
98     &                   nact*nact,0,g_sht2))
99     &     call errquit('ga_create g_sht2 failed',0, GA_ERR)
100
101      call ga_distribution(g_nht2,iam,g_jlo,g_jhi,g_ilo,g_ihi)
102      do i=1,nocc
103        ad1=(i-1)*nvir
104        if (ad1+1.ge.g_ilo.and.ad1+1.le.g_ihi)then
105          do j=1,nocc
106            ad2=(j-1)*nvir
107            if (ad2+1.ge.g_jlo.and.ad2+1.le.g_jhi)then
108            ad3=(i-1)*nocc+j
109              call ga_get(g_sht2t,1,nbf*nbf,ad3,ad3,
110     &                    scrb,nbf*nbf)
111              ad3=0
112              do ksh=1,nsh
113                if (.not. bas_cn2bfr(basis,ksh,klo,khi))
114     &              call errquit('vvvv: bas_cn2bfr',ksh, BASIS_ERR)
115                do lsh=1,nsh
116                  if (.not. bas_cn2bfr(basis,lsh,llo,lhi))
117     &                call errquit('vvvv: bas_cn2bfr',lsh, BASIS_ERR)
118                  do k=klo,khi
119                    do l=llo,lhi
120                      ad3=ad3+1
121                      scra((k-1)*nbf+l)=scrb(ad3)
122                    enddo
123                  enddo
124                enddo
125              enddo
126              call dgemm('n','n',nbf,nact,nbf,1.0d00,scra,nbf,
127     &                   cmo(1,ncor+1),nbf,0.0d00,scrb,nbf)
128              call dgemm('t','n',nact,nact,nbf,1.0d00,cmo(1,ncor+1),nbf,
129     &                   scrb,nbf,0.0d00,scra,nact)
130              ad3=(i-1)*nocc+j
131              call ga_put(g_sht2,1,nact*nact,ad3,ad3,
132     &                    scra,nact*nact)
133c
134              if (dob(1).eq.2)then
135                ad3=0
136                do a=nocc+1,nact
137                  do b=nocc+1,nact
138                    ad3=ad3+1
139                    scrb(ad3)=scra((a-1)*nact+b)
140                  enddo
141                enddo
142                call ga_acc(g_nht2,ad2+1,ad2+nvir,ad1+1,ad1+nvir,
143     &                      scrb,nvir,1.0d00)
144              endif
145c
146            endif
147          enddo
148        endif
149      enddo
150      if (.not.ga_destroy(g_sht2t))
151     &    call errquit('ga_dest g_sht2t fail',0, GA_ERR)
152c
153      if (occsdps) then
154         call pstat_off(ps_ht2pm)
155      else
156         call qexit('ht2pm',0)
157      endif
158c
159      return
160      end
161      subroutine ccsd_t2pm(basis,nsh,ncor,nocc,nvir,nbf,g_st2,g_nt2,cmo,
162     &                     t1,scra,scrb,nbfdim)
163      implicit none
164#include "errquit.fh"
165      integer basis,nsh,ncor,nocc,nvir,nbf,g_st2,g_nt2,nbfdim
166      double precision cmo(nbf,nbf),scra(nbf*nbf),scrb(nbf*nbf),
167     &                 t1(nocc,nvir)
168#include "mafdecls.fh"
169#include "global.fh"
170#include "tcgmsg.fh"
171#include "bas.fh"
172#include "rtdb.fh"
173#include "ccsdps.fh"
174c
175      integer g_st2t,lnoo,iam
176      integer i,j,k,l,a,b,ad1,ad2,ad3,g_jlo,g_jhi,g_ilo,g_ihi,
177     &        ish,ilo,ihi,jsh,jlo,jhi,ksh,klo,khi,lsh,llo,lhi,
178     &        xy,x,y,ipp,imm,ii,jj,iijj,jjii,ad4,nodes
179c
180      if (occsdps) then
181         call pstat_on(ps_t2pm)
182      else
183         call qenter('t2pm',0)
184      endif
185c
186      iam=ga_nodeid()
187      nodes=ga_nnodes()
188      lnoo=nocc*nocc
189c
190*ga:1:0
191      do ad4=iam+1,lnoo,nodes
192         i=(ad4-1)/nocc+1
193         j=ad4-(i-1)*nocc
194              ad1=(i-1)*nvir
195              ad2=(j-1)*nvir
196              call ga_get(g_nt2,ad2+1,ad2+nvir,ad1+1,ad1+nvir,
197     &                    scra,nvir)
198              ad3=0
199              do a=1,nvir
200                do b=1,nvir
201                  ad3=ad3+1
202                  scra(ad3)=scra(ad3)+t1(i,a)*t1(j,b)
203                enddo
204              enddo
205              call dgemm('n','t',nvir,nbf,nvir,1.0d00,scra,nvir,
206     &                   cmo(1,ncor+nocc+1),nbf,0.0d00,scrb,nvir)
207              call dgemm('n','n',nbf,nbf,nvir,1.0d00,cmo(1,ncor+nocc+1),
208     &                   nbf,scrb,nvir,0.0d00,scra,nbf)
209              ad3=0
210              do ksh=1,nsh
211                if (.not. bas_cn2bfr(basis,ksh,klo,khi))
212     &              call errquit('vvvv: bas_cn2bfr',ksh, BASIS_ERR)
213                do lsh=1,ksh
214                  if (.not. bas_cn2bfr(basis,lsh,llo,lhi))
215     &                call errquit('vvvv: bas_cn2bfr',lsh, BASIS_ERR)
216                  do k=klo,khi
217                    do l=llo,lhi
218                      ad3=ad3+1
219                      scrb(ad3)=scra((k-1)*nbf+l)
220                    enddo
221                  enddo
222                enddo
223              enddo
224              ad3=(i-1)*nocc+j
225              call ga_put(g_st2,ad3,ad3,1,nbfdim,
226     &                    scrb,1)
227      enddo
228      call ga_sync()
229c
230c ------------------------------------------------------------
231c - form t2+/- (see gustavos paper)
232c - t2+ = t_ij^ab + t_ij^ba
233c - t2- = t_ij^ab - t_ij^ba
234c - some prefactors absorbed
235c ------------------------------------------------------------
236      call ga_distribution(g_st2,iam,g_jlo,g_jhi,g_ilo,g_ihi)
237      xy=0
238      do ish=1,nsh
239        if (.not. bas_cn2bfr(basis,ish,ilo,ihi))
240     $      call errquit('vvvv: bas_cn2bfr',ish, BASIS_ERR)
241        do jsh=1,ish
242          if (.not. bas_cn2bfr(basis,jsh,jlo,jhi))
243     $        call errquit('vvvv: bas_cn2bfr',jsh, BASIS_ERR)
244          do x=ilo,ihi
245            do y=jlo,jhi
246              xy=xy+1
247              if (xy.ge.g_ilo.and.xy.le.g_ihi)then
248                call ga_get(g_st2,1,lnoo,xy,xy,scra,lnoo)
249                ipp=0
250                imm=nocc*(nocc+1)/2
251                do ii=1,nocc
252                  do jj=1,ii-1
253                    iijj=(ii-1)*nocc+jj
254                    jjii=(jj-1)*nocc+ii
255                    ipp=ipp+1
256                    imm=imm+1
257                    scrb(ipp)=(scra(iijj)+scra(jjii))
258                    scrb(imm)=(scra(iijj)-scra(jjii))
259                  enddo
260                  iijj=(ii-1)*nocc+ii
261                  ipp=ipp+1
262                  scrb(ipp)=scra(iijj)+scra(iijj)
263                enddo
264                if (x.eq.y)then
265                  call dscal(lnoo,0.25d00,scrb,1)
266                else
267                  call dscal(lnoo,0.5d00,scrb,1)
268                endif
269                call ga_put(g_st2,1,lnoo,xy,xy,scrb,lnoo)
270              endif
271            enddo
272          enddo
273        enddo
274      enddo
275      call ga_sync()
276c
277      if (occsdps) then
278         call pstat_off(ps_t2pm)
279      else
280         call qexit('t2pm',0)
281      endif
282c
283      return
284      end
285      subroutine ccsd_sxy(basis,nsh,ncor,nocc,nvir,nact,nbf,g_st2,
286     &                    g_sht2,g_c,g_x,offsh,snsi,sisn,lssni,scre,
287     &                    mem2,max2e,eri1,eri2,t1,cmo,t1ao,scra,
288     &                    scrb,lscr,kscr,tol2e,iprt,tklst,maxints,
289     &                    ish_idx,shinf,max_sht2_blk,nbfdim,a_st2,
290     &                    use_ccsd_omp)
291      implicit none
292#include "errquit.fh"
293      integer basis,nsh,ncor,nocc,nvir,nact,nbf,g_st2,g_sht2,
294     &        g_c,g_x,lssni,
295     &        mem2,max2e,lscr,kscr,iprt,max_sht2_blk,nbfdim
296      integer offsh(nsh,nsh,2),tklst(nsh*(nsh+1)/2,2),maxints
297      integer ish_idx(nsh*nsh+2),shinf(nsh,3)
298      double precision tol2e,eri1(maxints),eri2(maxints),scre(mem2),
299     &                 t1(nocc*nvir),cmo(nbf,nbf),scra(kscr),
300     &                 scrb(lscr),snsi(lssni),sisn(lssni),
301     &                 t1ao(nocc*nbf),tx(2),a_st2(nocc*nocc,nbfdim)
302      logical schwarz1,schwarz2
303
304#include "mafdecls.fh"
305#include "global.fh"
306#include "tcgmsg.fh"
307#include "bas.fh"
308#include "geom.fh"
309#include "rtdb.fh"
310#include "schwarz.fh"
311#include "eaf.fh"
312*rak-s
313#include "ccsd_time.fh"
314*rak-e
315#include "ccsdps.fh"
316#include "ccsd_data.fh"
317      integer ad1,ish,ilo,ihi,jsh,jlo,jhi,ksh,klo,khi,lsh,llo,lhi,next,
318     &        icnt,nfi,nfk,lnijkl,offjl,offlj,lnoo,nfj,nfl,nodes,iam
319      logical flush, store, use_storage
320      integer nxtask,blocksize,logsize,intsize,ijk(2),cnijkl
321      integer icount(2),jcnt,btki,kcount,jcount,ishb,iblk,bnijkl,ival
322      integer cart_2e4c,ad2
323      double precision store_need,store_used
324      double precision faddr,kaddr
325      external nxtask, cart_2e4c
326      integer trishel,iamoff,iii,jjj,istart
327      logical, optional, intent(in) :: use_ccsd_omp
328      if (.not.present(use_ccsd_omp)) then
329          call errquit('ccsd_pzamp: use_ccsd_omp not present!',0,0)
330      endif
331c
332      if (occsdps) then
333         call pstat_on(ps_sxy)
334      else
335         call qenter('sxy',0)
336      endif
337c
338      lnoo=nocc*nocc
339      iam=ga_nodeid()
340      nodes=ga_nnodes()
341c
342c ------------------------------------------------------------
343c - work out an index array for shell offset
344c ------------------------------------------------------------
345      call ga_sync()
346      ad1=0
347      ad2=0
348      do ish=1,nsh
349        if (.not. bas_cn2bfr(basis,ish,ilo,ihi))
350     $      call errquit('vvvv: bas_cn2bfr',ish, BASIS_ERR)
351        nfi=ihi-ilo+1
352        shinf(ish,1)=nfi
353        shinf(ish,2)=ilo
354        shinf(ish,3)=ihi
355        do jsh=1,nsh
356          if (.not. bas_cn2bfr(basis,jsh,jlo,jhi))
357     $        call errquit('vvvv: bas_cn2bfr',jsh, BASIS_ERR)
358          nfj=jhi-jlo+1
359          offsh(ish,jsh,1)=ad1
360          ad1=ad1+nfi*nfj
361          if (jsh.le.ish) then
362             offsh(ish,jsh,2)=ad2
363             ad2=ad2+nfi*nfj
364          endif
365        enddo
366      enddo
367c ------------------------------------------------------------
368c create t1ao
369c ------------------------------------------------------------
370      call dgemm('n','t',nbf,nocc,nvir,1.0d00,cmo(1,ncor+nocc+1),
371     &           nbf,t1,nocc,0.0d00,t1ao,nbf)
372c
373c ------------------------------------------------------------
374c - loop over the integral generation this is what needs to be
375c - efficient, so gather timings for this
376c - at the moment we are shell blocked, we could go to atom
377c - blocking if this is inefficient
378c - note integrals computed 4 times minimal list
379c ------------------------------------------------------------
380c
381c      call ga_sync()
382      tx(1)=tcgtime()
383      if (iam.eq.0.and.iprt.gt.5)
384     &    print *,' begin parallel integral generation'
385      call ga_zero(g_c)
386c      call ga_sync()
387      faddr=0.d0
388      kaddr=0.d0
389      logsize=ma_sizeof(MT_LOG,2,MT_BYTE)
390      intsize=ma_sizeof(MT_INT,1,MT_BYTE)
391c
392c If this is the second call to ccsd_sxy we may have the integrals and other
393c info on disk. We should read from there....
394c
395      if (repeat.and.use_disk) then
396c
397c Get the number of lsh,lsh blocks this node is doing
398c
399      if (eaf_read(sxy_hl,faddr,icount,intsize*2).ne.0)
400     &   call errquit('ccsd_sxy: read failed',1,DISK_ERR)
401      faddr=faddr+intsize*2
402c
403c Start the main loop, using the integrals that are stored on disk
404c
405      do jcnt=1,icount(1)
406c
407c Get the jsh and lsh shell info
408c
409         if (eaf_read(sxy_hl,faddr,ijk,intsize*2).ne.0)
410     &      call errquit('ccsd_sxy: read failed',2,DISK_ERR)
411         faddr=faddr+intsize*2
412c
413         jsh=(tklst(ijk(1)+1,2)-1)/nsh+1
414         lsh=tklst(ijk(1)+1,2)-(jsh-1)*nsh
415c
416         nfj=shinf(jsh,1)
417         jlo=shinf(jsh,2)
418         jhi=shinf(jsh,3)
419         nfl=shinf(lsh,1)
420         llo=shinf(lsh,2)
421         lhi=shinf(lsh,3)
422c
423         call dcopy(lnoo*nfj*nfl,0.0d00,0,scrb,1)
424         call dcopy(nfj*nfl*nbf*nocc,0.0d00,0,snsi,1)
425         call dcopy(nfj*nfl*nbf*nocc,0.0d00,0,sisn,1)
426c
427c ijk(2) contains number of ksh*(ish-block) blocks stored
428c
429         do iblk=1,ijk(2)
430c
431c Get data for ksh and ish-block. First get number of blocks
432c for the ksh to be calculated, then get whole block of appropriate
433c length
434c
435            call qenter('r_read',0)
436            if (eaf_read(sxy_hl,faddr,ival,intsize).ne.0)
437     &         call errquit('ccsd_sxy: read failed',
438     &         3,DISK_ERR)
439            blocksize=intsize*ival
440            if (eaf_read(sxy_hl,faddr,ish_idx,blocksize).ne.0)
441     &         call errquit('ccsd_sxy: read failed',
442     &         4,DISK_ERR)
443            faddr=faddr+blocksize
444c
445c Get integrals, number of integrals in ish_idx(1,2)
446c
447            blocksize=ma_sizeof(MT_DBL,ish_idx(2),MT_BYTE)
448            if (eaf_read(sxy_hl,faddr,eri1,blocksize).ne.0)
449     &         call errquit('ccsd_sxy: read failed',
450     &         5,DISK_ERR)
451            faddr=faddr+blocksize
452            if (eaf_read(sxy_hl,faddr,eri2,blocksize).ne.0)
453     &         call errquit('ccsd_sxy: read failed',
454     &         6,DISK_ERR)
455            faddr=faddr+blocksize
456            call qexit('r_read',0)
457            call t2eri(ish_idx(1),ish_idx,jlo,jhi,nfj,llo,lhi,nfl,nsh,
458     &                 eri1,eri2,scra,scrb,lnoo,nocc,offsh,nbf,g_st2,
459     &                 shinf,max_sht2_blk,snsi,sisn,t1ao,nbfdim,a_st2,
460     &                 use_ccsd_omp)
461         enddo
462         ad1=offsh(jsh,lsh,1)
463         call ga_put(g_sht2,1,lnoo,ad1+1,ad1+nfj*nfl,scrb,lnoo)
464         offjl=offsh(jsh,lsh,1)
465         offlj=offsh(lsh,jsh,1)
466         if (use_ccsd_omp) then
467           call ccsd_idx2_omp(snsi,sisn,cmo,lscr,
468     &                        nfj,nfl,ncor,nocc,nact,nbf,
469     &                        jlo,jhi,llo,lhi,offjl,offlj,g_x,g_c)
470         else
471           call ccsd_idx2(snsi,sisn,cmo,scra,scrb,lscr,
472     &                    nfj,nfl,ncor,nocc,nact,nbf,
473     &                    jlo,jhi,llo,lhi,offjl,offlj,g_x,g_c)
474         endif
475      enddo
476c
477c     End reading integrals from disk in repeat iteration
478c
479      endif!!!!if (repeat.and.use_disk) then
480c
481      if (.not.repeat.or..not.use_disk.or.icount(2).gt.0) then
482c
483c We got here in one of two cases:
484c
485c 1. First time through, use dynamical (nxtask) distribution to
486c    calculate and store the integrals....
487c    - Start from iam
488c    - use_storage = true
489c    - reset icount array
490c 2. Repeat time through, needing to recalculate some integrals
491c    as they could not be stored to disk...
492c    If no disk was used, case is like 1, but with use_storage = false
493c    - Start from icount(2)
494c    - use_storage = false
495c
496      if (repeat.and.use_disk) then
497         istart = icount(2)
498         use_storage=.false.
499      else
500         istart = iam
501         icount(1)=0
502         icount(2)=0
503         if (use_disk) then
504            if (eaf_write(sxy_hl,faddr,icount,intsize*2).ne.0)
505     &         call errquit('ccsd_sxy: write failed',1,DISK_ERR)
506            faddr=faddr+intsize*2
507         endif
508         use_storage=use_disk
509         store_used=0.d0
510      endif
511      do icnt=istart,nsh*(nsh+1)/2-1,nodes
512c
513c ijk contains jsh and lsh info plus number of ksh and ish blocks
514c
515            ijk(1)=icnt
516            jsh=(tklst(icnt+1,2)-1)/nsh+1
517            lsh=tklst(icnt+1,2)-(jsh-1)*nsh
518c
519            nfj=shinf(jsh,1)
520            jlo=shinf(jsh,2)
521            jhi=shinf(jsh,3)
522            nfl=shinf(lsh,1)
523            llo=shinf(lsh,2)
524            lhi=shinf(lsh,3)
525c
526            call dcopy(lnoo*nfj*nfl,0.0d00,0,scrb,1)
527            call dcopy(nfj*nfl*nbf*nocc,0.0d00,0,snsi,1)
528            call dcopy(nfj*nfl*nbf*nocc,0.0d00,0,sisn,1)
529c
530            bnijkl=1
531            kcount=2
532            flush=.false.
533            jcount=0
534            trishel=nsh*(nsh+1)/2
535            iamoff=iam*(trishel/nodes)
536c
537c Determine if we can still store the integrals to disk
538c If not, calculate and set integral recalculation point icnt
539c
540            if (use_storage) then
541              store_need=dfloat(ma_sizeof(MT_DBL,2*nfj*nfl*nbf*(nbf+1),
542     $                                    MT_BYTE))
543              store_need=store_need+dfloat(ma_sizeof(MT_INT,nsh*nsh+2,
544     $                                               MT_BYTE))
545              if ((store_need+store_used).gt.store_avail) then
546                use_storage=.false.
547                icount(2)=icnt
548              else
549                icount(1)=icount(1)+1
550                kaddr=faddr
551                if (eaf_write(sxy_hl,faddr,ijk,intsize*2).ne.0) call
552     &                 errquit('ccsd_sxy: write failed',2,DISK_ERR)
553                faddr=faddr+intsize*2
554              endif
555            endif
556c
557            do iii=iamoff,trishel+iamoff-1 !begin loop
558               jjj=iii
559               if(iii.gt.trishel-1) jjj=iii-trishel
560               ish=(tklst(jjj+1,2)-1)/nsh+1
561               ksh=tklst(jjj+1,2)-(ish-1)*nsh
562               nfk=shinf(ksh,1)
563                schwarz1=schwarz_shell(ish,jsh)*
564     &                   schwarz_shell(ksh,lsh).ge.tol2e
565                schwarz2=schwarz_shell(ish,lsh)*
566     &                   schwarz_shell(ksh,jsh).ge.tol2e
567                nfi=shinf(ish,1)
568                lnijkl=nfi*nfj*nfk*nfl
569                cnijkl=cart_2e4c(basis,ish,jsh,ksh,lsh)
570c
571c Check if the integral buffers are full or we have done all ish and ksh
572c If so store and process integrals, else add another block
573c
574                store=(bnijkl+cnijkl.gt.maxints)
575  111           continue
576                if ((store.and.(schwarz1.or.schwarz2)).or.flush) then
577c
578c Store the integrals and process them in t2eri and ccsd_idx1
579c
580c First four spots in ish_idx are used to store additional data:
581c  (1)=kcount  -> number of ish+ksh blocks plus 2 for this data (2*2)
582c  (2)=bnijkl  -> number of integrals
583c
584                   ish_idx(1)=kcount
585                   ish_idx(2)=bnijkl
586                   if (use_storage) then
587                     jcount=jcount+1
588                     blocksize=intsize*kcount
589                     call qenter('f_write',0)
590                     if (eaf_write(sxy_hl,faddr,ish_idx,blocksize).ne.0)
591     &                  call errquit('ccsd_sxy: write failed',
592     &                  3,DISK_ERR)
593                     faddr=faddr+blocksize
594                     store_used=store_used+blocksize
595                     blocksize=ma_sizeof(MT_DBL,bnijkl,MT_BYTE)
596                     if (eaf_write(sxy_hl,faddr,eri1,blocksize).ne.0)
597     &                  call errquit('ccsd_sxy: write failed',
598     &                  4,DISK_ERR)
599                     faddr=faddr+blocksize
600                     if (eaf_write(sxy_hl,faddr,eri2,blocksize).ne.0)
601     &                  call errquit('ccsd_sxy: write failed',
602     &                  5,DISK_ERR)
603                     faddr=faddr+blocksize
604                     call qexit('f_write',0)
605                     store_used=store_used+2*blocksize
606                   endif
607c
608c Process the integral block
609c
610                   call t2eri(kcount,ish_idx,jlo,jhi,nfj,llo,lhi,nfl,
611     &                        nsh,eri1,
612     &                        eri2,scra,scrb,lnoo,nocc,offsh,nbf,g_st2,
613     &                        shinf,max_sht2_blk,snsi,sisn,t1ao,
614     &                        nbfdim,a_st2,
615     &                        use_ccsd_omp)
616c
617c Reset some indices
618c
619                   bnijkl=1
620                   kcount=2
621                endif
622c
623c Add next block of integrals
624c
625                if (.not.flush.and.(schwarz1.or.schwarz2)) then
626                   call qenter('ints',0)
627                   kcount=kcount+1
628                   ish_idx(kcount)=(ksh-1)*nsh+ish
629                   if (schwarz1) then
630                      if(bnijkl+max2e.gt.maxints) then
631                         write(6,'(I6,A,I20,A)') ga_nodeid(),
632     W                        ' : increase stack memory ',
633     W                     (trishel+iamoff-iii+1)*2*max2e*8/1024/1024,
634     M                        ' MBytes '
635                         call util_flush(6)
636                         call errquit(' ccsdsxy: not enough MA',0,0)
637                      endif
638                     call int_2e4c(basis, ish, jsh, basis, ksh, lsh,
639     $                             mem2, scre, max2e, eri1(bnijkl) )
640                   else
641                     call dcopy(lnijkl,0.0d00,0,eri1(bnijkl),1)
642                   endif
643                   if (schwarz2) then
644                     if (jsh.eq.lsh) then
645                        call dcopy(lnijkl,eri1(bnijkl),1,eri2(bnijkl),1)
646                     else
647                      if(bnijkl+max2e.gt.maxints) then
648                         write(6,*) ' increase stack memory'
649                         call util_flush(6)
650                         call errquit(' ccsdsxy: not enough MA',0,0)
651                      endif
652                        call int_2e4c(basis,ish,lsh,basis,ksh,jsh,
653     $                                mem2,scre,max2e,eri2(bnijkl))
654                     endif
655                   else
656                     call dcopy(lnijkl,0.0d00,0,eri2(bnijkl),1)
657                   endif
658                   bnijkl=bnijkl+lnijkl
659                   call qexit('ints',0)
660                endif
661                if (iii.eq.(trishel+iamoff-1).and..not.flush) then
662                   flush=.true.
663                   goto 111
664                endif
665                flush=.false.
666            enddo!end loop
667c
668            if (use_storage) then
669              ijk(2)=jcount
670              call qenter('f_write',0)
671              if (eaf_write(sxy_hl,kaddr,ijk,intsize*2).ne.0)
672     &           call errquit('ccsd_sxy: write failed',6,DISK_ERR)
673              call qexit('f_write',0)
674            endif
675c
676            ad1=offsh(jsh,lsh,1)
677            call ga_put(g_sht2,1,lnoo,ad1+1,ad1+nfj*nfl,scrb,lnoo)
678            offjl=offsh(jsh,lsh,1)
679            offlj=offsh(lsh,jsh,1)
680            if (use_ccsd_omp) then
681              call ccsd_idx2_omp(snsi,sisn,cmo,lscr,
682     &                           nfj,nfl,ncor,nocc,nact,nbf,
683     &                           jlo,jhi,llo,lhi,offjl,offlj,g_x,g_c)
684            else
685              call ccsd_idx2(snsi,sisn,cmo,scra,scrb,lscr,
686     &                       nfj,nfl,ncor,nocc,nact,nbf,
687     &                       jlo,jhi,llo,lhi,offjl,offlj,g_x,g_c)
688            endif
689      enddo
690      faddr=0.d0
691      if (.not.repeat.and.use_disk) then
692         if (eaf_write(sxy_hl,faddr,icount,intsize*2).ne.0) call
693     &       errquit('ccsd_sxy: write failed',7,DISK_ERR)
694      endif
695      repeat=.true.
696      endif!!!if (.not.repeat.or..not.use_disk.or.icount(2).gt.0)
697cedo
698      call ga_sync()
699c sync before 3 and 4 index transformation
700      if (use_ccsd_omp) then
701        call ccsd_idx34_omp(basis,cmo,
702     &                      nsh,ncor,nocc,nact,nbf,
703     &                      g_x,g_c)
704      else
705        call ccsd_idx34(basis,cmo,scra,scrb,
706     &                  nsh,ncor,nocc,nact,nbf,
707     &                  g_x,g_c)
708      endif
709c
710      call ga_sync()
711      tx(2)=tcgtime()
712      if (iam.eq.0) then
713*rak     write(6,*)'Time around main block',tx(2)-tx(1)
714         main_block_time = tx(2)-tx(1)
715      endif
716c
717      if (occsdps) then
718         call pstat_off(ps_sxy)
719      else
720         call qexit('sxy',0)
721      endif
722c
723      return
724      end
725