1      SUBROUTINE tce_mo2e_zones_4a_disk_4s_new_patch_m3(rtdb,d_v2,
2     1                                kax_v2_alpha_offset,size_2e)
3C     $Id$
4C     This is a Fortran77 program generated by Tensor Contraction Engine v.1.0
5C     Copyright (c) Battelle & Pacific Northwest National Laboratory (2002)
6C     t ( p1 p2 h3 h4 )_t
7      IMPLICIT NONE
8#include "rtdb.fh"
9#include "global.fh"
10#include "mafdecls.fh"
11#include "stdio.fh"
12#include "util.fh"
13#include "bas.fh"
14#include "schwarz.fh"
15#include "sym.fh"
16#include "sf.fh"
17#include "errquit.fh"
18#include "tce.fh"
19#include "tce_main.fh"
20c
21c
22      integer rtdb                 ! Run-time database
23      integer d_v2                 ! MO integrals
24      integer kax_v2_alpha_offset  ! MO integrals offset
25      integer size_2e              ! 2e file size
26c
27      INTEGER size_2g2a,l_2g2a,k_2g2a
28      INTEGER azone1,azone2,azone3,azone4
29      INTEGER g1b,g2b,g3b,g4b
30      INTEGER igi1,igi2,igi3,igi4
31      INTEGER ii,i,j,k,l,N,ipos1,ipos2,ipos3,ipos4
32      INTEGER del1,del2,p1rel,p2rel
33      INTEGER size_4a,l_4a,k_4a
34c
35      integer mu,nu,rho,sigma
36      integer mu_lo,mu_hi
37      integer nu_lo,nu_hi
38      integer rho_lo,rho_hi
39      integer sigma_lo,sigma_hi
40      integer mu_range
41      integer nu_range
42      integer rho_range
43      integer sigma_range
44      integer mu1,nu1,rho1,sigma1
45      integer shift_mu,shift_nu
46      integer shift_rho,shift_sigma
47      integer work1,work2          ! Work array sizes
48      integer l_work1,k_work1      ! Work array 1
49      integer l_work2,k_work2      ! Work array 2
50      integer imu1,inu1,irho1,isigma1
51c
52      integer l_movecs_orb,k_movecs_orb
53      integer l_gpair,k_gpair
54      integer len_pair,g12_shift
55c ATTENTION,ACHTUNG,UWAGA 2000 - max # of CPU
56c
57      integer size_2g2z,l_2g2z,k_2g2z
58      integer tot_azone1_sh,tot_azone2_sh
59      integer tot_azone3_sh,tot_azone4_sh
60      integer ixi,jxi,point_pair
61c
62      integer iha,ihb !number of corr. alpha, beta holes
63      integer ipa,ipb !number of corr. alpha, beta particles
64c compression
65      integer max_size_temp,size_temp,sumx
66      double precision wall,cpu,wall1,cpu1,wall2,cpu2,wall3,cpu3
67c
68      double precision tot_zone(1000)
69c
70      integer l_3a1m_offset,k_3a1m_offset,size_3a1mf
71      integer l_1a3m_offset,k_1a3m_offset,size_1a3mf
72      integer d_1a3m,d_3a1m
73      integer size_3a1m
74      integer l_3a1m,k_3a1m
75      integer size_amc
76      integer key_3a1m,offset_3a1m
77      integer l_4a_sort,k_4a_sort
78      integer size_2a2m,l_2a2m,k_2a2m
79      integer l_3a1m_sort,k_3a1m_sort
80      integer l_2a2m_aux,k_2a2m_aux
81      integer key_1a3m,offset_1a3m
82      integer l_1a3m,k_1a3m,size_1a3m
83      integer l_1a3m_sort,k_1a3m_sort
84      INTEGER IROW,ICOL,IRES
85      INTEGER INDEX_PAIR
86      integer l_4m,k_4m,size_4m,key_4m,offset_4m
87      integer sf3a1m_chunk,sf1a3m_chunk
88c
89      integer l_integral,l_coeff
90      integer k_integral,k_coeff
91      integer size_ic,size_icc,size_integral,size_coeff,max_na
92c
93      integer l_4af_offset,k_4af_offset,d_4af
94      integer k_2a2m_offset,l_2a2m_offset
95      integer sf_chunk,request
96      integer key_4af,offset_4af,size_4af
97      integer sf2a2m_chunk,key_2a2m,offset_2a2m
98      integer d_2a2m,size_2a2mf
99      character*255 filename,filename_aux
100c
101      integer size_max1,size_max2,ip_max1,ip_max2
102      integer chunk_max1,chunk_max2,i_from_x,i_to_x
103      integer isizef(5),imaxch(5)
104      integer d_max1,d_max2
105c
106      integer l_a4_offset,k_a4_offset
107      integer offset_bl4a,size_bl4a,a4_ini,a4_fin
108      integer offset_ismall,length_a4,key_ini
109      integer l_a3_offset,k_a3_offset
110      integer length_a3,offset_bl3a1m,size_bl3a1m
111      integer size_3a1mi,a3_ini,a3_fin
112c
113      integer l_a2_offset,k_a2_offset,length_a2
114      integer offset_bl2a2m,size_bl2a2m,a2_ini,a2_fin
115      integer size_1a3mi,size_2a2mi
116c
117      integer length_a1
118      integer l_a1_offset,k_a1_offset
119      integer offset_bl1a3m,size_bl1a3m,a1_ini,a1_fin
120c I/O improvements
121      integer handle_4a,handle_3a1m,handle_2a2m,handle_1a3m
122c
123c *** debug ***
124c      double precision xxx,xmax
125c      integer iqx
126      integer ierrcode1,ierrcode2
127c *************
128c
129c multiple files ===============
130       integer k_offset_4a_m(1000),l_offset_4a_m(1000)
131       integer k_offset_3a1m_m(1000),l_offset_3a1m_m(1000)
132       integer k_offset_2a2m_m(1000),l_offset_2a2m_m(1000)
133       integer k_offset_1a3m_m(1000),l_offset_1a3m_m(1000)
134       integer size_4a_m(1000),size_3a1m_m(1000)
135       integer size_2a2m_m(1000),size_1a3m_m(1000)
136       integer n_4a_files,n_3a1m_files
137       integer n_2a2m_files,n_1a3m_files
138       integer handle_4a_m(1000),handle_3a1m_m(1000)
139       integer handle_2a2m_m(1000),handle_1a3m_m(1000)
140       integer my_file_to_write,my_file_to_read
141       integer n_max_files
142c
143       integer size_l(1000),n_str,size_g
144       integer addr
145c
146c *** debug ***
147       double precision xmax
148c ==============================
149      logical parallel
150c
151      INTEGER length
152      INTEGER next
153      INTEGER nprocs
154      INTEGER count
155      integer nxtask
156      external nxtask
157      logical nodezero
158      logical idiskl
159c
160c
161c
162ccx      idisk=.true.
163      if(idisk.eq.0) then
164       idiskl=.false.
165      else
166       idiskl=.true.
167      end if
168c
169      parallel = .true.
170c
171c
172      max_size_temp=imaxsize**4
173c
174c
175      do ii=1,1000
176       tot_zone(ii)=0.0d0
177      enddo
178      if(atpart.gt.1000)
179     &  call errquit('tce_zones: atpart too big',1,MA_ERR)
180      sumx=0
181      do ii=1,atpart
182       tot_zone(ii)=sumx
183       sumx=sumx+nalength(ii)
184      enddo
185c
186      nodezero=(ga_nodeid().eq.0)
187c
188c *** debug ***
189c      if(nodezero) then
190c        write(6,*)'--------- NEW LOOP STRUCTURE --------'
191c        call util_flush(6)
192c      end if
193c *************
194c
195c this module is called only if intorb = .true.
196c N is the number of correlated orbitals
197        N = nmo(1) - nfc(1) - nfv(1)
198        iha = nocc(1)-nfc(1)
199        ihb = nocc(ipol)-nfc(ipol)
200        ipa = nmo(1)-nocc(1)-nfv(1)
201        ipb = nmo(ipol)-nocc(ipol)-nfv(ipol)
202c
203      sf_chunk=(imaxsize+10)**4
204      sf3a1m_chunk=((imaxsize+10)**3)*tile_dim
205      sf2a2m_chunk=((tile_dim)**2)*((imaxsize+10)**2)
206      sf1a3m_chunk=((tile_dim)**3)*(imaxsize+10)
207c
208c l_integral and l_coeff local files are opened here
209c ATTENTION,ACHTUNG,UWAGA - "manually" defined
210cc      size_ic=2*(imaxsize+10)**4
211cc      size_ic=2*(imaxsize)**4
212      max_na=0
213      do ixi=1,atpart
214       if(nalength(ixi).gt.max_na) max_na=nalength(ixi)
215      enddo
216c     it was
217cc      size_ic=2*(max_na)**4+1
218c it is =================================================
219        size_ic=2*((max_na)**4)+1
220        size_icc=tile_dim*max_na
221c
222       if(nodezero) then
223        write(6,*)'2_EL_BATCH AND COEFF MATRIX'
224        write(6,*)'size_ic',size_ic
225        write(6,*)'size_icc',size_icc
226        write(6,*)'tile_dim',tile_dim
227        write(6,*)'max_na',max_na
228        call util_flush(6)
229       end if
230c
231      if (.not.ma_push_get(mt_dbl,size_ic,'l_int',
232     1  l_integral,k_integral))
233     1  call errquit('tce_4s: MA problem l_int',0,MA_ERR)
234c
235      if (.not.ma_push_get(mt_dbl,size_icc,'l_coeff',
236     1  l_coeff,k_coeff))
237     1  call errquit('tce_4s: MA problem l_coeff',0,MA_ERR)
238c
239c
240c
241       chunk_max1=size_ic
242       chunk_max2=size_ic
243c
244c offsets for 2
245c
246c
247c DISK USAGE ONLY !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
248c
249      if(idiskl) then
250        if(.not.parallel)
251     1       call errquit('sf only for parallel runs',1,DISK_ERR)
252c
253c
254c dummy for idiskl true
255c
256c
257      else
258       if(.not.parallel)
259     1       call errquit('sf only for parallel runs',1,DISK_ERR)
260       if(parallel) call ga_sync()
261        if(size_4af.le.size_max2) then
262         call createfile(filename,d_max1,size_max1)
263         call reconcilefile(d_max1,size_max1)
264         call ga_zero(d_max1)
265         call createfile(filename,d_max2,size_max2)
266         call reconcilefile(d_max2,size_max2)
267         call ga_zero(d_max2)
268        else
269         call createfile(filename,d_max2,size_max2)
270         call reconcilefile(d_max2,size_max2)
271         call ga_zero(d_max2)
272         call createfile(filename,d_max1,size_max1)
273         call reconcilefile(d_max1,size_max1)
274         call ga_zero(d_max1)
275        end if
276      end if
277c
278c
279c alpha orbitals only
280c
281      if (.not.ma_push_get(mt_dbl,nbf*(iha+ipa)
282     1  ,"sorted MO coeffs",
283     2  l_movecs_orb,k_movecs_orb))
284     3  call errquit('tce_mo2e_zone: MA problem 1',0,
285     2    BASIS_ERR)
286      call dfill(nbf*(iha+ipa),0.0d0, dbl_mb(k_movecs_orb), 1)
287      do i=1,iha
288      do isigma1=1,nbf
289       dbl_mb(k_movecs_orb+(i-1)*nbf+isigma1-1)=
290     & dbl_mb(k_movecs_sorted+(i-1)*nbf+isigma1-1)
291      enddo
292      enddo
293      do i=iha+1,iha+ipa
294      do isigma1=1,nbf
295       dbl_mb(k_movecs_orb+(i-1)*nbf+isigma1-1)=
296     & dbl_mb(k_movecs_sorted+(i+ihb-1)*nbf+isigma1-1)
297      enddo
298      enddo
299c
300c
301c
302c
303       n_max_files=1000
304       do i=1,n_max_files
305        k_offset_4a_m(i)=0
306        l_offset_4a_m(i)=0
307        k_offset_3a1m_m(i)=0
308        l_offset_3a1m_m(i)=0
309        k_offset_2a2m_m(i)=0
310        l_offset_2a2m_m(i)=0
311        k_offset_1a3m_m(i)=0
312        l_offset_1a3m_m(i)=0
313        size_4a_m(i)=0
314        size_3a1m_m(i)=0
315        size_2a2m_m(i)=0
316        size_1a3m_m(i)=0
317       enddo
318c
319c      write(6,*)'step1'
320c      call util_flush(6)
321c
322cxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
323ccx       call tce_1a3m_offsetyy1(l_offset_1a3m_m,k_offset_1a3m_m,   xx
324ccx     &      size_1a3m_m,n_1a3m_files)                             xx
325cxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
326c
327      do i=1,1000
328       size_l(i)=0
329      enddo
330c
331      i = 0
332      DO g3b = 1,noa+nva   !k
333      i=i+1
334      length=0
335      DO g4b = g3b,noa+nva !l
336      DO g2b = 1,noa+nva   !i
337      DO azone1 = 1,atpart !nu
338      length = length + 1
339      END DO
340      END DO
341      END DO
342      size_l(i)=length
343      END DO
344      n_str=i
345      n_1a3m_files=n_str
346c
347      do i=1,n_str
348      IF (.not.MA_PUSH_GET(mt_int,2*size_l(i)+1,'noname4',
349     & l_offset_1a3m_m(i),k_offset_1a3m_m(i)))
350     & CALL ERRQUIT('tce_t2_offset',0,MA_ERR)
351      int_mb(k_offset_1a3m_m(i)) = size_l(i)
352c
353c      write(6,*)'1a3m yy1 i size_l(i)',i,size_l(i)
354c      call util_flush(6)
355c
356      enddo
357c
358      i = 0
359      DO g3b = 1,noa+nva   !k
360        i=i+1
361        addr=0
362        size_g=0
363        length=size_l(i)
364      DO g4b = g3b,noa+nva !l
365      DO g2b = 1,noa+nva   !i
366      DO azone1 = 1,atpart !nu
367      addr = addr + 1
368      int_mb(k_offset_1a3m_m(i)+addr)=azone1-1+atpart*(g2b-1+(noa+nva)*
369     &  (g4b-1+(noa+nva)*(g3b-1)))
370      int_mb(k_offset_1a3m_m(i)+length+addr) = size_g
371      size_g = size_g +
372     &   int_mb(k_range_alpha+g3b-1) * int_mb(k_range_alpha+g4b-1)
373     &      * int_mb(k_range_alpha+g2b-1) *  nalength(azone1)
374      END DO
375      END DO
376      END DO
377        size_1a3m_m(i)=size_g
378      END DO
379c
380c
381c
382c
383c      write(6,*)'step2 n_1a3m_files',n_1a3m_files
384c      call util_flush(6)
385c
386cxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
387cxxxxx       call tce_2a2m_offsetyy1(l_offset_2a2m_m,k_offset_2a2m_m, xx
388cxxxxx     &      size_2a2m_m,n_2a2m_files)                           xx
389cxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
390c
391      do i=1,1000
392       size_l(i)=0
393      enddo
394c
395      i = 0
396      DO g3b = 1,noa+nva   !k
397      i=i+1
398      length=0
399      DO g4b = g3b,noa+nva !l
400      DO azone1 = 1,atpart !nu
401      DO azone2 = azone1,atpart !mu
402      length = length + 1
403      END DO
404      END DO
405      END DO
406      size_l(i)=length
407      END DO
408      n_str=i
409      n_2a2m_files=n_str
410c
411      do i=1,n_str
412      IF (.not.MA_PUSH_GET(mt_int,2*size_l(i)+1,'noname3',
413     & l_offset_2a2m_m(i),k_offset_2a2m_m(i)))
414     & CALL ERRQUIT('tce_t2_offset',0,MA_ERR)
415      int_mb(k_offset_2a2m_m(i)) = size_l(i)
416c
417c      write(6,*)'2a2m yy1 i size_l(i)',i,size_l(i)
418c      call util_flush(6)
419c
420      enddo
421c
422      i = 0
423      DO g3b = 1,noa+nva   !k
424        i=i+1
425        addr=0
426        size_g=0
427        length=size_l(i)
428      DO g4b = g3b,noa+nva !l
429      DO azone1 = 1,atpart !nu
430      DO azone2 = azone1,atpart !mu
431      addr = addr + 1
432      int_mb(k_offset_2a2m_m(i)+addr)=azone2 - 1 + atpart*( azone1 - 1 +
433     & atpart*(g4b - 1 + (noa+nva) * (g3b-1)))
434      int_mb(k_offset_2a2m_m(i)+length+addr) = size_g
435      size_g = size_g +
436     &   int_mb(k_range_alpha+g3b-1) * int_mb(k_range_alpha+g4b-1)
437     &      *  nalength(azone1) * nalength(azone2)
438      END DO
439      END DO
440      END DO
441        size_2a2m_m(i)=size_g
442      END DO
443c
444c
445c
446c      write(6,*)'step3,n_2a2m_files',n_2a2m_files
447c      call util_flush(6)
448cxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
449ccx       call tce_3a1m_offsetyy1(l_offset_3a1m_m,k_offset_3a1m_m,   xxx
450ccx     &      size_3a1m_m,n_3a1m_files)                             xxx
451cxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
452c      write(6,*)'step4 n_3a1m_files'
453c      call util_flush(6)
454c
455      do i=1,1000
456       size_l(i)=0
457      enddo
458c
459      i=0
460      DO azone1 = 1,atpart      !nu
461      i=i+1
462      length=0
463      DO azone2 = azone1,atpart !mu
464      DO g4b    =  1,noa+nva    !g2b
465      DO azone3 = 1,atpart      !sigma
466      length = length + 1
467      END DO
468      END DO
469      END DO
470      size_l(i)=length
471      END DO
472      n_str=i
473      n_3a1m_files=n_str
474c *** debug ***
475c      write(6,*)'from 3a1m'
476c      do i=1,n_str
477c       write(6,*)'size_l(i) i',size_l(i),i
478c       call util_flush(6)
479c      enddo
480c      write(6,*)
481c *************
482c
483      do i=1,n_str
484      IF (.not.MA_PUSH_GET(mt_int,2*size_l(i)+1,'noname2',
485     & l_offset_3a1m_m(i),k_offset_3a1m_m(i)))
486     & CALL ERRQUIT('tce_t2_offset',0,MA_ERR)
487      int_mb(k_offset_3a1m_m(i)) = size_l(i)
488c
489c      write(6,*)'3a1m yy1 i size_l(i)',i,size_l(i)
490c      call util_flush(6)
491c
492      enddo
493c
494      i=0
495      DO azone1 = 1,atpart      !nu
496        i=i+1
497        addr=0
498        size_g=0
499        length=size_l(i)
500      DO azone2 = azone1,atpart !mu
501      DO g4b    =  1,noa+nva    !g2b
502      DO azone3 =  1,atpart     !sigma
503      addr = addr + 1
504      int_mb(k_offset_3a1m_m(i)+addr)=azone3-1+atpart*(g4b-1+(noa+nva)*
505     &    (azone2 - 1 + atpart * (azone1 - 1)))
506      int_mb(k_offset_3a1m_m(i)+length+addr) = size_g
507      size_g = size_g + nalength(azone1) * nalength(azone2) *
508     &  nalength(azone3) * int_mb(k_range_alpha+g4b-1)
509      END DO
510      END DO
511      END DO
512        size_3a1m_m(i)=size_g
513      END DO
514c
515cxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
516cxx       call tce_4a_offsetyy1(l_offset_4a_m,k_offset_4a_m,   xxxx
517cxx     &      size_4a_m,n_4a_files)                           xxxx
518cxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
519c
520      do i=1,1000
521       size_l(i)=0
522      enddo
523c
524      i = 0
525      DO azone1 = 1,atpart      !nu
526      i=i+1
527      length=0
528      DO azone2 = azone1,atpart !mu
529      DO azone3 = 1,atpart      !sigma
530      DO azone4 = azone3,atpart !rho
531      length = length + 1
532      END DO
533      END DO
534      END DO
535      size_l(i)=length
536      END DO
537      n_str=i
538      n_4a_files=n_str
539c
540      do i=1,n_str
541      IF (.not.MA_PUSH_GET(mt_int,2*size_l(i)+1,'noname1',
542     & l_offset_4a_m(i),k_offset_4a_m(i)))
543     & CALL ERRQUIT('tce_t2_offset',0,MA_ERR)
544      int_mb(k_offset_4a_m(i)) = size_l(i)
545c
546c      write(6,*)'4a yy1 i size_l(i)',i,size_l(i)
547c      call util_flush(6)
548c
549      enddo
550c
551      i=0
552      DO azone1 = 1,atpart      !nu
553        i=i+1
554        addr=0
555        size_g=0
556        length=size_l(i)
557      DO azone2 = azone1,atpart !mu
558      DO azone3 = 1,atpart      !sigma
559      DO azone4 = azone3,atpart !rho
560      addr = addr + 1
561      int_mb(k_offset_4a_m(i)+addr)=azone4 - 1 + atpart * (azone3 - 1 +
562     &  atpart * (azone2 - 1 + atpart * (azone1 - 1)))
563      int_mb(k_offset_4a_m(i)+length+addr) = size_g
564      size_g = size_g + nalength(azone1) * nalength(azone2) *
565     &  nalength(azone3) * nalength(azone4)
566      END DO
567      END DO
568      END DO
569        size_4a_m(i)=size_g
570      END DO
571c
572c
573c
574c
575c
576c      write(6,*)'step5 n_4a_files',n_4a_files
577c      call util_flush(6)
578c debug *******
579c        stop
580c *************
581c
582c *** debug ***
583          if(nodezero) then
584           write(6,*)'n_4a_files',n_4a_files
585           write(6,*)'n_3a1m_files',n_3a1m_files
586           write(6,*)'n_2a2m_files',n_2a2m_files
587           write(6,*)'n_1a3m_files',n_1a3m_files
588           write(6,*)'atpart',atpart
589           write(6,*)'noa+nva',noa+nva
590           write(6,*)'noa',noa
591           write(6,*)'nva',nva
592           call util_flush(6)
593          end if
594c *************
595c
596c
597      call int_mem_2e4c(work1,work2)
598      if (.not.ma_push_get(mt_dbl,work1,'work1',l_work1,k_work1))
599     1  call errquit('tce_ao2e: MA problem work1',0,MA_ERR)
600      if (.not.ma_push_get(mt_dbl,work2,'work2',l_work2,k_work2))
601     1  call errquit('tce_ao2e: MA problem work2',1,MA_ERR)
602c
603c
604c
605
606      if(nodezero) then
607        write(6,*)'before sf_create'
608        write(6,*)'n_4a_files',n_4a_files
609        write(6,*)'sizes:'
610        do i=1,n_4a_files
611         write(6,301) i,size_4a_m(i)
612        enddo
613        write(6,*)'-----------------'
614  301   format('size_4a_m(i)',2x,i4,2x,i10)
615        call util_flush(6)
616      end if
617c
618      if(idiskl) then
619       do i=1,n_4a_files
620ccx        call tce_filenameindexed(i,'a4f',filename_aux)
621ccx        call util_file_name(filename_aux,.false.,.false.,filename)
622        call tce_filename_4ind(i,'a4f',filename)
623c *** debug ***
624c        write(6,*)'filename=',filename,i,ga_nodeid()
625c        write(6,*)'size_4a_m(i)',size_4a_m(i),i,ga_nodeid()
626c
627c *************
628        if (sf_create(filename,dfloat(bytes)*dfloat(size_4a_m(i)),
629     1   dfloat(bytes)*dfloat(size_4a_m(i)),chunk_max1,handle_4a_m(i))
630     1   .ne.0)
631     2   call errquit('4-index: sf problem opening 4af',0,DISK_ERR)
632       enddo
633      end if
634c
635      if(nodezero) then
636        write(6,*)'after sf_create'
637        call util_flush(6)
638      end if
639c
640c ============================
641c
642             cpu = - util_cpusec()
643             wall = - util_wallsec()
644c
645      nprocs = GA_NNODES()
646      count = 0
647      next = NXTASK(nprocs, 1)
648c attention: I have interchanged the orderign
649c from    azone1 azone2 azone3 azone4
650c to      azone3 azone4 azone1 azone2
651c      DO azone1 = 1,atpart      !nu
652c      DO azone2 = azone1,atpart !mu
653      DO azone3 = 1,atpart      !sigma
654      DO azone4 = azone3,atpart !rho
655      DO azone1 = 1,atpart      !nu
656      DO azone2 = azone1,atpart !mu
657      IF (next.eq.count) THEN
658c ---------------------------
659        my_file_to_write=azone1
660c(atpart*(atpart+1))/2-
661c     1   ((atpart-azone1+1)*(atpart-azone1))/2-
662c     1   (atpart-azone2+1)+1
663        size_4a = nalength(azone1)*nalength(azone2)*
664     1            nalength(azone3)*nalength(azone4)
665        if(.not.ma_push_get(mt_dbl,size_4a,'4a',l_4a,k_4a))
666     1     call errquit('tce_4af_zones1: MA problem',0,MA_ERR)
667        call dfill(size_4a, 0.0d0, dbl_mb(k_4a), 1)
668         shift_mu = 0
669         do mu    = a2length(azone2)+1,a2length(azone2+1)
670            if (.not.bas_cn2bfr(ao_bas_han,mu,mu_lo,mu_hi))
671     1      call errquit('tce_ao2e: basis fn range problem 1',0,
672     2      BASIS_ERR)
673            mu_range = mu_hi - mu_lo + 1
674         shift_nu = 0
675         do nu    = a2length(azone1)+1,a2length(azone1+1)
676            if (.not.bas_cn2bfr(ao_bas_han,nu,nu_lo,nu_hi))
677     1      call errquit('tce_ao2e: basis fn range problem 1',0,
678     2      BASIS_ERR)
679            nu_range = nu_hi - nu_lo + 1
680         shift_rho = 0
681         do rho   = a2length(azone4)+1,a2length(azone4+1)
682            if (.not.bas_cn2bfr(ao_bas_han,rho,rho_lo,rho_hi))
683     1      call errquit('tce_ao2e: basis fn range problem 1',0,
684     2      BASIS_ERR)
685            rho_range = rho_hi - rho_lo + 1
686         shift_sigma = 0
687         do sigma = a2length(azone3)+1,a2length(azone3+1)
688            if (.not.bas_cn2bfr(ao_bas_han,sigma,sigma_lo,sigma_hi))
689     1      call errquit('tce_ao2e: basis fn range problem 1',0,
690     2      BASIS_ERR)
691            sigma_range = sigma_hi - sigma_lo + 1
692            if (schwarz_shell(rho,sigma)*schwarz_shell(mu,nu)
693     1          .ge. tol2e) then
694            call int_2e4c(ao_bas_han,mu,nu,ao_bas_han,rho,sigma,
695     1           work2,dbl_mb(k_work2),work1,dbl_mb(k_work1))
696            i=0
697             do mu1     = 1,mu_range
698             do nu1     = 1,nu_range
699             do rho1    = 1,rho_range
700             do sigma1  = 1,sigma_range
701            i=i+1
702            inu1=nu1+shift_nu
703            isigma1=sigma1+shift_sigma
704            imu1=mu1+shift_mu
705            irho1=rho1+shift_rho
706c (isigma1,irho1|inu1, imu1)
707c differnce with other versiosn
708c here
709            ipos1=(((inu1-1)*nalength(azone2)+imu1-1)*
710     1            nalength(azone3)+isigma1-1)*nalength(azone4)
711     2            +irho1
712            dbl_mb(k_4a+ipos1-1)=dbl_mb(k_work1+i-1)
713            enddo
714            enddo
715            enddo
716            enddo
717            end if !schwarz  screening
718         shift_sigma = shift_sigma + sigma_range
719         enddo !sigma
720         shift_rho   = shift_rho + rho_range
721         enddo !rho
722         shift_nu    = shift_nu + nu_range
723         enddo !nu
724         shift_mu    = shift_mu + mu_range
725         enddo !mu
726c
727c fixing offsets and sf_writing
728         key_4af=azone4 - 1 + atpart * (azone3 - 1 +
729     &          atpart * (azone2 - 1 + atpart * (azone1 - 1)))
730         call tce_hash(int_mb(k_offset_4a_m(my_file_to_write)),
731     &          key_4af,offset_4af)
732      if(idiskl) then !----
733c *** debug ***
734c       xmax=0.0d0
735c       do i=1,size_4a
736c        xmax=xmax+dbl_mb(k_4a+i-1)
737c       enddo
738c       write(6,200) azone4,azone3,azone2,azone1,xmax
739c 200   format('4A ',4i5,3x,f19.8)
740c       call util_flush(6)
741c *************
742        ierrcode1=sf_write(handle_4a_m(my_file_to_write),
743     1    dfloat(bytes)*dfloat(offset_4af),
744     1    dfloat(bytes)*dfloat(size_4a),dbl_mb(k_4a),request)
745ccx        if (sf_write(i_to_x,dfloat(bytes)*dfloat(offset_4af),
746ccx     1    dfloat(bytes)*dfloat(size_4a),dbl_mb(k_4a),request).ne.0)
747ccx     1    THEN
748        if(ierrcode1.ne.0) then
749          write(6,8154)ierrcode1
750          write(6,8155)ga_nodeid(),key_4af,offset_4af,size_4a
751          call util_flush(6)
752          call errquit('zones put: sf problem2x1',1,DISK_ERR)
753        end if
754        ierrcode2=sf_wait(request)
755ccx        if (sf_wait(request).ne.0)
756        if(ierrcode2.ne.0) then
757           write(6,8153)ierrcode2
758           call util_flush(6)
759           call errquit('zones put: sf problem3x1',2,DISK_ERR)
760        end if
761      else
762        call ga_put(i_to_x,offset_4af+1,offset_4af+size_4a,1,1,
763     1    dbl_mb(k_4a),1)
764      end if          !----
765c closing l_4a file
766        if (.not.ma_pop_stack(l_4a))
767     1   call errquit('tce_mo2e_4af2: l_4a',15,MA_ERR)
768c ---------------------------
769      next = NXTASK(nprocs, 1)
770      END IF
771      count = count + 1
772      ENDDO !azone4
773      ENDDO !azone3
774      ENDDO !azone2
775      ENDDO !azone1
776c
777c
778       call ga_sync()
779       next = nxtask(-nprocs,1)
780c
781      cpu = cpu + util_cpusec()
782      wall = wall + util_wallsec()
783      if(nodezero) then
784       write(6,*)'4A integrals cpu,wall ',cpu,wall
785       call util_flush(6)
786      end if
787c
788      if (.not.ma_pop_stack(l_work2))
789     1  call errquit('tce_ao2e: MA problem',14,MA_ERR)
790      if (.not.ma_pop_stack(l_work1))
791     1  call errquit('tce_ao2e: MA problem',15,MA_ERR)
792c
793c
794      if(nodezero) then
795        write(6,*)'STEP2 4index'
796        call util_flush(6)
797      end if
798c
799c
800c
801c
802c
803c DISK CHANGES =============== 3A1M
804c closing/reopening handle_4a_m files
805c opening handle_3a1m_m files
806      if(idiskl) then
807c
808      do i=1,n_4a_files
809        if(sf_rwtor(handle_4a_m(i)).ne.0)
810     2    call errquit('4-index: sf_rwtor 4a',0,DISK_ERR)
811      enddo
812c
813       do i=1,n_3a1m_files
814ccx        call tce_filenameindexed(i,'3a1mf',filename_aux)
815ccx        call util_file_name(filename_aux,.false.,.false.,filename)
816        call tce_filename_4ind(i,'3a1mf',filename)
817        if (sf_create(filename,dfloat(bytes)*dfloat(size_3a1m_m(i)),
818     1   dfloat(bytes)*dfloat(size_3a1m_m(i)),chunk_max1,
819     1   handle_3a1m_m(i)).ne.0)
820     2   call errquit('4-index: sf problem opening 3a1mf',0,DISK_ERR)
821       enddo
822c
823      end if ! idiskl
824c ============================
825c
826      cpu = - util_cpusec()
827      wall = - util_wallsec()
828c
829c 3A1M part here
830c
831      nprocs = GA_NNODES()
832      count = 0
833      next = NXTASK(nprocs, 1)
834c
835c changed ordering
836c from:
837c      azone1 azone2 g4b azone3
838c to:
839c      g4b azone3 azone1 azone2
840c
841c      DO azone1 = 1,atpart      !nu
842c      DO azone2 = azone1,atpart !mu
843cuu      DO g4b = 1,noa+nva        !k
844cuu      DO azone3 = 1,atpart      !sigma
845cuu      DO azone1 = 1,atpart      !nu
846cuu      DO azone2 = azone1,atpart !mu
847c --------------------------------------
848cb      DO g4b = 1,noa+nva        !k
849cb      DO azone3 = 1,atpart      !sigma
850cb      DO azone2= 1,atpart
851cb      DO azone1=1,azone2
852cc      DO azone2= 1,atpart
853cc      DO azone1=1,azone2
854cc      DO g4b = 1,noa+nva        !k
855cc      DO azone3 = 1,atpart      !sigma
856
857      DO azone1=1,azone2
858      DO azone2= 1,atpart
859      DO g4b = 1,noa+nva        !k
860      DO azone3 = 1,atpart      !sigma
861c --------------------------------------
862c
863      IF (next.eq.count) THEN
864c
865        my_file_to_read=azone1
866c(atpart*(atpart+1))/2-
867c     1   ((atpart-azone1+1)*(atpart-azone1))/2-
868c     1   (atpart-azone2+1)+1
869        my_file_to_write=azone1
870c(atpart*(atpart+1))/2-
871c     1   ((atpart-azone1+1)*(atpart-azone1))/2-
872c     1   (atpart-azone2+1)+1
873c
874      size_3a1m=int_mb(k_range_alpha+g4b-1)
875     1          *nalength(azone1)*nalength(azone2)*nalength(azone3)
876      if (.not.ma_push_get(mt_dbl,size_3a1m,'3a1m',l_3a1m,k_3a1m))
877     1    call errquit('tce_4ind: step1_1',0,MA_ERR)
878      call dfill(size_3a1m, 0.0d0, dbl_mb(k_3a1m), 1)
879c
880c    routine for offset
881c
882      call a4_offset(azone1,azone2,azone3,size_ic,
883     &           l_a4_offset,k_a4_offset)
884c
885      length_a4=int_mb(k_a4_offset)
886      DO ii=1,length_a4
887       a4_ini      = int_mb(k_a4_offset+ii)
888       a4_fin      = int_mb(k_a4_offset+length_a4+ii)
889       key_ini=a4_ini - 1 + atpart * (azone3 - 1 +
890     &          atpart * (azone2 - 1 + atpart * (azone1 - 1)))
891       call tce_hash(int_mb(k_offset_4a_m(my_file_to_read)),
892     &          key_ini,offset_bl4a)
893c
894       size_bl4a   = int_mb(k_a4_offset+2*length_a4+ii)
895      if(idiskl) then !---------
896        ierrcode1=sf_read(handle_4a_m(my_file_to_read),
897     1    dfloat(bytes)*dfloat(offset_bl4a),
898     1    dfloat(bytes)*dfloat(size_bl4a),dbl_mb(k_integral),
899     1    request)
900ccx        if (sf_read(i_from_x,dfloat(bytes)*dfloat(offset_bl4a),
901ccx     1    dfloat(bytes)*dfloat(size_bl4a),dbl_mb(k_integral),
902ccx     1    request).ne.0) THEN
903         if(ierrcode1.ne.0) THEN
904          write(6,*)'STEP2_step1'
905          write(6,8153)ierrcode1
906          write(6,8155)ga_nodeid(),key_ini,offset_bl4a,size_bl4a
907          write(6,*)'length_a4',length_a4
908          write(6,*)'ii=',ii
909          write(6,*)'i_from_x',i_from_x
910          call util_flush(6)
911          call errquit('zones put: sf problem2x2',1,DISK_ERR)
912        end if
913        ierrcode2=sf_wait(request)
914ccx        if (sf_wait(request).ne.0)
915        if (ierrcode2.ne.0) then
916           write(6,8154)ierrcode2
917           call util_flush(6)
918           call errquit('zones put: sf problem3x2',2,DISK_ERR)
919        end if
920      else
921        call ga_get(i_from_x,offset_bl4a+1,offset_bl4a+size_bl4a,1,1,
922     1    dbl_mb(k_integral),1)
923      end if !--------------
924c
925      offset_ismall=0
926      do i=a4_ini,a4_fin ! over small i ------------------
927      tot_azone4_sh=tot_zone(i)
928      j=0
929      do irho1   =  1,nalength(i)
930      do igi4=1,int_mb(k_range_alpha+g4b-1)
931       j=j+1
932       ipos1=(int_mb(k_offset_alpha+g4b-1)+igi4-1)*nbf+tot_azone4_sh
933     &       +irho1
934       dbl_mb(k_coeff+j-1)=dbl_mb(k_movecs_orb+ipos1-1)
935      enddo
936      enddo
937c C(g4b irho)* v(irho,isigma,imu,inu)
938      call dgemm('N','N',
939     &   int_mb(k_range_alpha+g4b-1),
940     &   nalength(azone1)*nalength(azone2)*nalength(azone3),
941     &   nalength(i),
942     &   1.0d0,dbl_mb(k_coeff),int_mb(k_range_alpha+g4b-1),
943     &   dbl_mb(k_integral+offset_ismall),nalength(i),1.0d0,
944     &   dbl_mb(k_3a1m),int_mb(k_range_alpha+g4b-1))
945c
946       offset_ismall=offset_ismall+nalength(azone1)*nalength(azone2)*
947     1            nalength(azone3)*nalength(i)
948       enddo ! over small i ----------------------
949c
950      ENDDO  ! over ii
951        if (.not.ma_pop_stack(l_a4_offset))
952     1   call errquit('tce_mo2e_trans_zones: uu l_2g2z',15,MA_ERR)
953c
954c
955c
956c
957      DO azone4 = 1, azone3-1  ! second part ----
958c getting piece of 4a
959        size_4a = nalength(azone1)*nalength(azone2)*
960     1            nalength(azone3)*nalength(azone4)
961       if(azone3.le.azone4) then
962         key_4af=azone4 - 1 + atpart * (azone3 - 1 +
963     &          atpart * (azone2 - 1 + atpart * (azone1 - 1)))
964       else
965         key_4af=azone3 - 1 + atpart * (azone4 - 1 +
966     &          atpart * (azone2 - 1 + atpart * (azone1 - 1)))
967       end if
968         call tce_hash(int_mb(k_offset_4a_m(my_file_to_read)),
969     1        key_4af,offset_4af)
970      if(idiskl) then !---------
971        ierrcode1=sf_read(handle_4a_m(my_file_to_read),
972     1 dfloat(bytes)*dfloat(offset_4af),
973     1 dfloat(bytes)*dfloat(size_4a),dbl_mb(k_integral),request)
974ccx        if (sf_read(i_from_x,dfloat(bytes)*dfloat(offset_4af),
975ccx     1 dfloat(bytes)*dfloat(size_4a),dbl_mb(k_integral),request).ne.0)
976ccx     1    THEN
977        if(ierrcode1.ne.0) then
978          write(6,*)'STEP2-step1b'
979          write(6,8153)ierrcode1
980          write(6,8155)ga_nodeid(),key_4af,offset_4af,size_4a
981          call errquit('zones put: sf problem2x3',1,DISK_ERR)
982        end if
983        ierrcode2=sf_wait(request)
984ccx        if (sf_wait(request).ne.0)
985        if(ierrcode2.ne.0) then
986           write(6,8154)ierrcode2
987           call util_flush(6)
988           call errquit('zones put: sf problem3x3',2,DISK_ERR)
989        end if
990      else
991        call ga_get(i_from_x,offset_4af+1,offset_4af+size_4a,1,1,
992     1    dbl_mb(k_integral),1)
993      end if !--------------
994c
995      if(azone4.lt.azone3) then
996        if(.not.ma_push_get(mt_dbl,size_4a,'4a_s',l_4a_sort,k_4a_sort))
997     1     call errquit('tce_4af_zones1: MA problem',0,MA_ERR)
998      CALL TCE_SORT_4KG_(dbl_mb(k_integral),dbl_mb(k_4a_sort),
999     & nalength(azone3),nalength(azone4),
1000     & nalength(azone2),nalength(azone1),
1001     &2,1,3,4,1.0d0)
1002      do i=1,size_4a
1003       dbl_mb(k_integral+i-1)=dbl_mb(k_4a_sort+i-1)
1004      enddo
1005        if (.not.ma_pop_stack(l_4a_sort))
1006     1   call errquit('tce_mo2e_trans_zones: uu l_2g2z',15,MA_ERR)
1007      end if
1008c
1009c  C(g4b irho) ==> l_2g2z
1010c
1011      size_amc=int_mb(k_range_alpha+g4b-1) *
1012     1         nalength(azone4)
1013      i=0
1014      tot_azone4_sh=tot_zone(azone4)
1015      do irho1   =  1,nalength(azone4)
1016      do igi4=1,int_mb(k_range_alpha+g4b-1)
1017       i=i+1
1018       ipos1=(int_mb(k_offset_alpha+g4b-1)+igi4-1)*nbf+tot_azone4_sh
1019     &       +irho1
1020       dbl_mb(k_coeff+i-1)=dbl_mb(k_movecs_orb+ipos1-1)
1021      enddo
1022      enddo
1023c C(g4b irho)* v(irho,isigma,imu,inu)
1024      call dgemm('N','N',
1025     &   int_mb(k_range_alpha+g4b-1),
1026     &   nalength(azone1)*nalength(azone2)*nalength(azone3),
1027     &   nalength(azone4),
1028     &   1.0d0,dbl_mb(k_coeff),int_mb(k_range_alpha+g4b-1),
1029     &   dbl_mb(k_integral),nalength(azone4),1.0d0,
1030     &   dbl_mb(k_3a1m),int_mb(k_range_alpha+g4b-1))
1031c
1032      ENDDO
1033c write to file
1034          key_3a1m =  azone3-1+atpart*(g4b-1+(noa+nva)*
1035     &    (azone2 - 1 + atpart * (azone1 - 1)))
1036         call tce_hash(int_mb(k_offset_3a1m_m(my_file_to_write)),
1037     &     key_3a1m,offset_3a1m)
1038       if(idiskl) then ! --------
1039c *** debug ***
1040c       xmax=0.0d0
1041c       do i=1,size_3a1m
1042c        xmax=xmax+dbl_mb(k_3a1m+i-1)
1043c       enddo
1044c       write(6,201) azone3,g4b,azone2,azone1,xmax
1045c 201   format('3A1M ',4i5,3x,f19.8)
1046c       call util_flush(6)
1047c *************
1048        ierrcode1=sf_write(handle_3a1m_m(my_file_to_write),
1049     1    dfloat(bytes)*dfloat(offset_3a1m),
1050     1    dfloat(bytes)*dfloat(size_3a1m),dbl_mb(k_3a1m),request)
1051        if (ierrcode1.ne.0) then
1052          write(6,*)'WRITE STEP3'
1053          write(6,8156) ierrcode1
1054          write(6,*)'offset_3a1m=',offset_3a1m
1055          write(6,*)'size_3a1m=',size_3a1m
1056          write(6,*)'i_to_x',i_to_x
1057          call util_flush(6)
1058          call errquit('zones put: sf problem21-b',1,DISK_ERR)
1059        end if
1060        ierrcode2=sf_wait(request)
1061ccx        if (sf_wait(request).ne.0)
1062        if(ierrcode2.ne.0) then
1063          write(6,8157)ierrcode2
1064          call errquit('zones put: sf problem31-b',2,DISK_ERR)
1065        end if
1066       else
1067        call ga_put(i_to_x,offset_3a1m+1,offset_3a1m+size_3a1m,1,1,
1068     1    dbl_mb(k_3a1m),1)
1069       end if ! ---------
1070c
1071      if (.not.ma_pop_stack(l_3a1m))
1072     1  call errquit('tce_mo2e_trans_zones: abc-MA problem',15,MA_ERR)
1073c
1074      next = NXTASK(nprocs, 1)
1075      END IF
1076      count = count + 1
1077      ENDDO
1078      ENDDO
1079      ENDDO
1080      ENDDO
1081c 3A1P fully done here
1082      call ga_sync()
1083      next = nxtask(-nprocs,1)
1084c
1085      cpu = cpu + util_cpusec()
1086      wall = wall + util_wallsec()
1087      if(nodezero) then
1088       write(6,*)'3A1M integrals cpu wall',cpu,wall
1089       call util_flush(6)
1090      end if
1091c
1092c Destroying all l_offset_4a_m(n_4a_files)
1093c
1094      if(idiskl) then
1095       do i=n_4a_files,1,-1
1096        if (.not.ma_pop_stack(l_offset_4a_m(i)))
1097     1   call errquit('4-ind n5: l_offset_4a_m err.',15,MA_ERR)
1098       enddo
1099      endif
1100c
1101c
1102      if(nodezero) then
1103        write(6,*)'STEP3 4index'
1104        call util_flush(6)
1105      end if
1106c
1107c
1108c
1109c
1110c
1111c 2A2M part here
1112c
1113c DISK CHANGES =============== 2A2M
1114c destroying handle_4a
1115c closing/reopening handle_3a1m file
1116c opening handle_2a2m file
1117      if(idiskl) then
1118c
1119      do i=1,n_4a_files
1120       if (sf_destroy(handle_4a_m(i)).ne.0)
1121     1   call errquit('tce_sf_destroy: sf problem handle_4a',15,MA_ERR)
1122      enddo
1123c
1124      do i=1,n_3a1m_files
1125        if(sf_rwtor(handle_3a1m_m(i)).ne.0)
1126     2    call errquit('4-index: sf_rwtor 3a1m',0,DISK_ERR)
1127      enddo
1128c
1129       do i=1,n_2a2m_files
1130ccx        call tce_filenameindexed(i,'2a2mf',filename_aux)
1131ccx        call util_file_name(filename_aux,.false.,.false.,filename)
1132        call tce_filename_4ind(i,'2a2mf',filename)
1133        if (sf_create(filename,dfloat(bytes)*dfloat(size_2a2m_m(i)),
1134     1   dfloat(bytes)*dfloat(size_2a2m_m(i)),chunk_max1,
1135     1   handle_2a2m_m(i)).ne.0)
1136     2   call errquit('4-index: sf problem opening 2a2mf',0,DISK_ERR)
1137       enddo
1138c
1139      end if ! idiskl true
1140c ============================
1141c
1142      cpu = - util_cpusec()
1143      wall = - util_wallsec()
1144c
1145      nprocs = GA_NNODES()
1146      count = 0
1147      next = NXTASK(nprocs, 1)
1148c
1149cuu      DO g3b = 1,noa+nva   !k
1150cuu      DO g4b = g3b,noa+nva !l
1151cuu      DO azone1 = 1,atpart      !nu
1152cuu      DO azone2 = azone1,atpart !mu
1153c -----------------------------------------
1154      DO g4b = 1,noa+nva !l
1155      DO g3b = 1,g4b
1156      DO azone2 = 1,atpart !mu
1157      DO azone1 = 1,azone2
1158c -----------------------------------------
1159ccx      DO g3b = 1,noa+nva   !k
1160ccx      DO g4b = g3b,noa+nva !l
1161c
1162      IF (next.eq.count) THEN
1163c
1164        my_file_to_read=azone1
1165c(atpart*(atpart+1))/2-
1166c     1   ((atpart-azone1+1)*(atpart-azone1))/2-
1167c     1   (atpart-azone2+1)+1
1168        my_file_to_write=g3b
1169c((noa+nva)*(noa+nva+1))/2-
1170c     1   ((noa+nva-g3b+1)*(noa+nva-g3b))/2-
1171c     1   (noa+nva-g4b+1)+1
1172c
1173      size_2a2m=int_mb(k_range_alpha+g4b-1)*int_mb(k_range_alpha+g3b-1)
1174     1          *nalength(azone1)*nalength(azone2)
1175      if (.not.ma_push_get(mt_dbl,size_2a2m,'2a2m',l_2a2m,k_2a2m))
1176     1    call errquit('tce_4ind: step1_1',0,MA_ERR)
1177      call dfill(size_2a2m, 0.0d0, dbl_mb(k_2a2m), 1)
1178c
1179      call a3_offset(azone1,azone2,g3b,size_ic,
1180     &           l_a3_offset,k_a3_offset)
1181      length_a3=int_mb(k_a3_offset)
1182c
1183      DO ii=1,length_a3     ! ii---loop
1184
1185       a3_ini    = int_mb(k_a3_offset+ii)
1186       a3_fin    = int_mb(k_a3_offset+length_a3+ii)
1187       key_ini   = a3_ini-1+atpart*(g3b-1+(noa+nva)*
1188     &    (azone2 - 1 + atpart * (azone1 - 1)))
1189       call tce_hash(int_mb(k_offset_3a1m_m(my_file_to_read)),
1190     &   key_ini,offset_bl3a1m)
1191c
1192       size_bl3a1m   = int_mb(k_a3_offset+2*length_a3+ii)
1193c
1194c
1195       if(idiskl) then ! -------
1196        ierrcode1=sf_read(handle_3a1m_m(my_file_to_read),
1197     2    dfloat(bytes)*dfloat(offset_bl3a1m),
1198     2    dfloat(bytes)*dfloat(size_bl3a1m),
1199     2    dbl_mb(k_integral),request)
1200         if(ierrcode1.ne.0) then
1201          write(6,*)'STEP3-step1a'
1202          write(6,8153) ierrcode1
1203          write(6,*)'offset_bl3a1m=',offset_bl3a1m
1204          write(6,*)'size_bl3a1m=',size_bl3a1m
1205          write(6,*)'ii=',ii
1206          write(6,*)'length_a3=',length_a3
1207          call util_flush(6)
1208          call errquit('zones put: sf problem21-b',1,DISK_ERR)
1209         end if
1210         ierrcode2=sf_wait(request)
1211         if(ierrcode2.ne.0) then
1212          write(6,*)'STEP3-step1a'
1213          write(6,8154) ierrcode2
1214          call util_flush(6)
1215          call errquit('zones put: sf problem31-b',2,DISK_ERR)
1216         end if
1217       else
1218        call ga_get(i_from_x,offset_bl3a1m+1,offset_bl3a1m+size_bl3a1m,
1219     1       1,1,dbl_mb(k_integral),1)
1220       end if ! -------
1221      offset_ismall=0
1222      do i=a3_ini,a3_fin ! over small i ------------------
1223      size_3a1mi=nalength(azone1)*nalength(azone2)*
1224     1            nalength(i)*int_mb(k_range_alpha+g3b-1)
1225      if (.not.ma_push_get(mt_dbl,size_3a1mi,'3a1m',
1226     1 l_3a1m_sort,k_3a1m_sort))
1227     1    call errquit('tce_4ind: step1_1',0,MA_ERR)
1228      CALL TCE_SORT_4KG_(dbl_mb(k_integral+offset_ismall),
1229     & dbl_mb(k_3a1m_sort),
1230     & int_mb(k_range_alpha+g3b-1),nalength(i),
1231     & nalength(azone2),nalength(azone1),
1232     &2,1,3,4,1.0d0)
1233      do j=1,size_3a1mi
1234       dbl_mb(k_integral+offset_ismall+j-1)=dbl_mb(k_3a1m_sort+j-1)
1235      enddo
1236      if (.not.ma_pop_stack(l_3a1m_sort))
1237     1  call errquit('tce_mo2e_trans_zones: abc-MA problem',15,MA_ERR)
1238c  C(g4b isigma)
1239      j=0
1240      tot_azone3_sh=tot_zone(i)
1241      do isigma1   =  1,nalength(i)
1242      do igi4=1,int_mb(k_range_alpha+g4b-1)
1243       j=j+1
1244       ipos1=(int_mb(k_offset_alpha+g4b-1)+igi4-1)*nbf+tot_azone3_sh
1245     &       +isigma1
1246       dbl_mb(k_coeff+j-1)=dbl_mb(k_movecs_orb+ipos1-1)
1247      enddo
1248      enddo
1249c C(g4b isigma)* v(isigma,g3b,imu,inu)
1250      call dgemm('N','N',
1251     &   int_mb(k_range_alpha+g4b-1),
1252     &   nalength(azone1)*nalength(azone2)*int_mb(k_range_alpha+g3b-1),
1253     &   nalength(i),
1254     &   1.0d0,dbl_mb(k_coeff),int_mb(k_range_alpha+g4b-1),
1255     &   dbl_mb(k_integral+offset_ismall),nalength(i),1.0d0,
1256     &   dbl_mb(k_2a2m),int_mb(k_range_alpha+g4b-1))
1257c
1258      offset_ismall=offset_ismall+nalength(azone1)*nalength(azone2)*
1259     1            nalength(i)*int_mb(k_range_alpha+g3b-1)
1260      enddo ! over small i ----------------------
1261      ENDDO                 ! ii---loop
1262c
1263        if (.not.ma_pop_stack(l_a3_offset))
1264     1   call errquit('tce_mo2e_trans_zones: uu l_2g2z',15,MA_ERR)
1265c write to file
1266          key_2a2m= azone2 - 1 + atpart*( azone1 - 1 +
1267     & atpart*(g4b - 1 + (noa+nva) * (g3b-1)))
1268         call tce_hash(int_mb(k_offset_2a2m_m(my_file_to_write)),
1269     &     key_2a2m,offset_2a2m)
1270       if(idiskl) then ! -------
1271c *** debug ***
1272c       xmax=0.0d0
1273c       do i=1,size_2a2m
1274c        xmax=xmax+dbl_mb(k_2a2m+i-1)
1275c       enddo
1276c       write(6,202) azone2,azone1,g4b,g3b,xmax,
1277c     &              my_file_to_write,ga_nodeid(),key_2a2m,
1278c     &              offset_2a2m
1279c 202   format('2A2M ',4i5,3x,f19.8,2x,2i6,2x,2i8)
1280c       call util_flush(6)
1281c *************
1282        ierrcode1=sf_write(handle_2a2m_m(my_file_to_write),
1283     2                   dfloat(bytes)*dfloat(offset_2a2m),
1284     2    dfloat(bytes)*dfloat(size_2a2m),dbl_mb(k_2a2m),request)
1285        if(ierrcode1.ne.0) then
1286         write(6,*)'STEP3 write'
1287         write(6,8156) ierrcode1
1288         write(6,*)'offset_2a2m=',offset_2a2m
1289         write(6,*)'size_2a2m=',size_2a2m
1290         write(6,*)'key_2a2m=',key_2a2m
1291         write(6,*)'i_to_x=',i_to_x
1292         call util_flush(6)
1293         call errquit('zones put: sf problem21-b',1,DISK_ERR)
1294        end if
1295        ierrcode2=sf_wait(request)
1296        if(ierrcode2.ne.0) then
1297         call errquit('zones put: sf problem31-b',2,DISK_ERR)
1298        end if
1299       else
1300        call ga_put(i_to_x,offset_2a2m+1,offset_2a2m+size_2a2m,1,1,
1301     1    dbl_mb(k_2a2m),1)
1302       end if  ! --------
1303c
1304      if (.not.ma_pop_stack(l_2a2m))
1305     1  call errquit('tce_mo2e_trans_zones: abc-MA problem',15,MA_ERR)
1306c
1307      next = NXTASK(nprocs, 1)
1308      END IF
1309      count = count + 1
1310      ENDDO
1311      ENDDO
1312      ENDDO
1313      ENDDO
1314c 2A2P fully done here
1315      call ga_sync()
1316      next = nxtask(-nprocs,1)
1317c
1318      cpu = cpu + util_cpusec()
1319      wall = wall + util_wallsec()
1320      if(nodezero) then
1321       write(6,*)'2A2M integrals cpu wall',cpu,wall
1322       call util_flush(6)
1323      end if
1324c
1325c
1326c Destroying all l_offset_3a1m_m(n_3a1m_files)
1327c
1328      if(idiskl) then
1329      do i=n_3a1m_files,1,-1
1330       if (.not.ma_pop_stack(l_offset_3a1m_m(i)))
1331     1  call errquit('4-ind n5: l_offset_3a1m_m err.',15,MA_ERR)
1332      enddo
1333      end if
1334c
1335c
1336c
1337      if(nodezero) then
1338        write(6,*)'STEP4 4index'
1339        call util_flush(6)
1340      end if
1341c
1342c
1343c
1344c
1345c
1346c
1347c 1A3M part here ([g4b][g3b]|[mu][nu]) =>([g4b][g3b]|[g2b][nu])
1348c                            [mu]>[nu] (azone2>=azone1)
1349c
1350c
1351c
1352c DISK CHANGES =============== 1A3M
1353c destroying handle_3a1m
1354c closing/reopening handle_2a2m file
1355c opening handle_1a3m file
1356      if(idiskl) then
1357c
1358      do i=1,n_3a1m_files
1359       if (sf_destroy(handle_3a1m_m(i)).ne.0)
1360     1  call errquit('tce_sf_destroy: sf problem handle_3a1m',15,MA_ERR)
1361      enddo
1362c
1363      do i=1,n_2a2m_files
1364        if(sf_rwtor(handle_2a2m_m(i)).ne.0)
1365     2    call errquit('4-index: sf_rwtor 2a2m',0,DISK_ERR)
1366      enddo
1367c
1368       do i=1,n_1a3m_files
1369ccx        call tce_filenameindexed(i,'1a3mf',filename_aux)
1370ccx        call util_file_name(filename_aux,.false.,.false.,filename)
1371        call tce_filename_4ind(i,'1a3mf',filename)
1372        if (sf_create(filename,dfloat(bytes)*dfloat(size_1a3m_m(i)),
1373     1   dfloat(bytes)*dfloat(size_1a3m_m(i)),chunk_max1,
1374     1   handle_1a3m_m(i)).ne.0)
1375     2   call errquit('4-index: sf problem opening 1a3mf',0,DISK_ERR)
1376       enddo
1377c
1378      end if !idiskl
1379c ============================
1380c
1381      cpu = - util_cpusec()
1382      wall = - util_wallsec()
1383c
1384c
1385      nprocs = GA_NNODES()
1386      count = 0
1387      next = NXTASK(nprocs, 1)
1388c
1389cuu      DO g2b = 1,noa+nva   !
1390cuu      DO azone1 = 1,atpart      !nu
1391cuu      DO g3b = 1,noa+nva   !k
1392cuu      DO g4b = g3b,noa+nva !l
1393c -------------------------------------------------
1394      DO g2b = 1,noa+nva   !
1395      DO azone1 = 1,atpart      !nu
1396      DO g4b = 1,noa+nva !l
1397      DO g3b = 1,g4b
1398c -------------------------------------------------
1399      IF (next.eq.count) THEN
1400        my_file_to_read=g3b
1401c((noa+nva)*(noa+nva+1))/2-
1402c     1   ((noa+nva-g3b+1)*(noa+nva-g3b))/2-
1403c     1   (noa+nva-g4b+1)+1
1404        my_file_to_write=g3b
1405c((noa+nva)*(noa+nva+1))/2-
1406c     1   ((noa+nva-g3b+1)*(noa+nva-g3b))/2-
1407c     1   (noa+nva-g4b+1)+1
1408      size_1a3m=int_mb(k_range_alpha+g4b-1)*int_mb(k_range_alpha+g3b-1)
1409     1         *int_mb(k_range_alpha+g2b-1)*nalength(azone1)
1410      if (.not.ma_push_get(mt_dbl,size_1a3m,'1a3m',l_1a3m,k_1a3m))
1411     1    call errquit('tce_4ind: step1_1',0,MA_ERR)
1412      call dfill(size_1a3m, 0.0d0, dbl_mb(k_1a3m), 1)
1413c
1414      DO azone2 = 1,atpart ! two cases here
1415       size_2a2m=int_mb(k_range_alpha+g4b-1)*int_mb(k_range_alpha+g3b-1)
1416     1           *nalength(azone1)*nalength(azone2)
1417       if (.not.ma_push_get(mt_dbl,size_2a2m,'2a2m',l_2a2m,k_2a2m))
1418     1    call errquit('tce_4ind: step1_1',0,MA_ERR)
1419       if(azone2.ge.azone1) then !-----------------------------------------
1420        key_2a2m=azone2 - 1 + atpart*( azone1 - 1 +
1421     &  atpart*(g4b - 1 + (noa+nva) * (g3b-1)))
1422         call tce_hash(int_mb(k_offset_2a2m_m(my_file_to_read)),
1423     &        key_2a2m,offset_2a2m)
1424        if(idiskl) then !------
1425        if (sf_read(handle_2a2m_m(my_file_to_read),
1426     &    dfloat(bytes)*dfloat(offset_2a2m),
1427     &    dfloat(bytes)*dfloat(size_2a2m),dbl_mb(k_2a2m),request).ne.0)
1428     &    call errquit('zones put: sf problem21-b',1,DISK_ERR)
1429        if (sf_wait(request).ne.0)
1430     1    call errquit('zones put: sf problem31-b',2,DISK_ERR)
1431        else
1432        call ga_get(i_from_x,offset_2a2m+1,offset_2a2m+size_2a2m,1,1,
1433     1    dbl_mb(k_2a2m),1)
1434        end if !-----
1435       else  ! ------------------------------------------------------------
1436        key_2a2m=azone1 - 1 + atpart*( azone2 - 1 +
1437     &  atpart*(g4b - 1 + (noa+nva) * (g3b-1)))
1438         call tce_hash(int_mb(k_offset_2a2m_m(my_file_to_read)),
1439     &                 key_2a2m,offset_2a2m)
1440        if(idiskl) then ! ------
1441        if (sf_read(handle_2a2m_m(my_file_to_read),
1442     &    dfloat(bytes)*dfloat(offset_2a2m),
1443     &    dfloat(bytes)*dfloat(size_2a2m),dbl_mb(k_2a2m),request).ne.0)
1444     &    call errquit('zones put: sf problem21-b',1,DISK_ERR)
1445        if (sf_wait(request).ne.0)
1446     1    call errquit('zones put: sf problem31-b',2,DISK_ERR)
1447        else
1448        call ga_get(i_from_x,offset_2a2m+1,offset_2a2m+size_2a2m,1,1,
1449     1    dbl_mb(k_2a2m),1)
1450        end if ! ------
1451        end if !  ---------------------------------------------------------
1452       if(azone2.ge.azone1) then
1453ccc       TRANSPOZYCJA
1454        if (.not.ma_push_get(mt_dbl,size_2a2m,'2a2m_sort',
1455     1   l_2a2m_aux,k_2a2m_aux))
1456     1     call errquit('tce_4ind: step1_1',0,MA_ERR)
1457        CALL TCE_SORT_4KG_(dbl_mb(k_2a2m),dbl_mb(k_2a2m_aux),
1458     &   int_mb(k_range_alpha+g4b-1),int_mb(k_range_alpha+g3b-1),
1459     &   nalength(azone2),nalength(azone1),
1460     &   1,2,4,3,1.0d0)
1461        do i=1,size_2a2m
1462         dbl_mb(k_2a2m+i-1)=dbl_mb(k_2a2m_aux+i-1)
1463        enddo
1464        if (.not.ma_pop_stack(l_2a2m_aux))
1465     1  call errquit('tce_mo2e_trans_zones: abc-MA problem',15,MA_ERR)
1466       end if
1467c
1468c  C(imu g2b) ==> l_2g2z
1469c
1470      size_amc=int_mb(k_range_alpha+g2b-1) *
1471     1         nalength(azone2)
1472      if (.not.ma_push_get(mt_dbl,size_amc,'2g2z',l_2g2z,k_2g2z))
1473     1    call errquit('tce_r2_divide3: xx-MA problem',0,MA_ERR)
1474      i=0
1475      tot_azone2_sh=tot_zone(azone2)
1476      do igi2=1,int_mb(k_range_alpha+g2b-1)
1477      do imu1   =  1,nalength(azone2)
1478       i=i+1
1479c       ipos1=(int_mb(k_offset_alpha+g2b-1)+igi2-1)*nbf+tot_azone1_sh
1480c     &       +inu1
1481       ipos1=(int_mb(k_offset_alpha+g2b-1)+igi2-1)*nbf+tot_azone2_sh
1482     &       +imu1
1483       dbl_mb(k_2g2z+i-1)=dbl_mb(k_movecs_orb+ipos1-1)
1484      enddo
1485      enddo
1486c  v(g4b,g3b,inu,imu) C(imu,g2b)
1487        if (.not.ma_push_get(mt_dbl,size_1a3m,'1a3m',
1488     1   l_1a3m_sort,k_1a3m_sort))
1489     1     call errquit('tce_4ind: step1_1',0,MA_ERR)
1490c dgemm
1491      call dgemm('N','N',
1492     &   int_mb(k_range_alpha+g4b-1)*int_mb(k_range_alpha+g3b-1)
1493     &   *nalength(azone1),
1494     &   int_mb(k_range_alpha+g2b-1),
1495     &   nalength(azone2),
1496     &   1.0d0,dbl_mb(k_2a2m),
1497     &   int_mb(k_range_alpha+g4b-1)*int_mb(k_range_alpha+g3b-1)
1498     &  *nalength(azone1),
1499     &   dbl_mb(k_2g2z),nalength(azone2),0.0d0,
1500     &   dbl_mb(k_1a3m_sort),
1501     &   int_mb(k_range_alpha+g4b-1)*int_mb(k_range_alpha+g3b-1)
1502     &  *nalength(azone1))
1503c
1504        i=0
1505        do inu1=1,nalength(azone1)
1506        do igi2=1,int_mb(k_range_alpha+g2b-1)
1507        do igi3=1,int_mb(k_range_alpha+g3b-1)
1508        do igi4=1,int_mb(k_range_alpha+g4b-1)
1509        i=i+1
1510        ipos1=igi4+int_mb(k_range_alpha+g4b-1)*(
1511     &   igi3-1+int_mb(k_range_alpha+g3b-1)*(
1512     &   inu1-1+nalength(azone1)*(igi2-1)))
1513        dbl_mb(k_1a3m+i-1)=dbl_mb(k_1a3m+i-1)+
1514     &                     dbl_mb(k_1a3m_sort+ipos1-1)
1515        enddo
1516        enddo
1517        enddo
1518        enddo
1519c
1520        if (.not.ma_pop_stack(l_1a3m_sort))
1521     1  call errquit('tce_mo2e_trans_zones: abc-MA problem',15,MA_ERR)
1522        if (.not.ma_pop_stack(l_2g2z))
1523     1   call errquit('tce_mo2e_trans_zones: uu l_2g2z',15,MA_ERR)
1524        if (.not.ma_pop_stack(l_2a2m))
1525     1   call errquit('tce_mo2e_trans_zones: uu l_2g2z',15,MA_ERR)
1526      ENDDO
1527c write to file
1528        key_1a3m= azone1-1+atpart*(g2b-1+(noa+nva)*
1529     &     (g4b-1+(noa+nva)*(g3b-1)))
1530         call tce_hash(int_mb(k_offset_1a3m_m(my_file_to_write)),
1531     &      key_1a3m,offset_1a3m)
1532c
1533       if(idiskl) then !---------
1534c *** debug ***
1535c       xmax=0.0d0
1536c       do i=1,size_1a3m
1537c        xmax=xmax+dbl_mb(k_1a3m+i-1)
1538c       enddo
1539c       write(6,203) azone1,g2b,g4b,g3b,xmax,my_file_to_write,
1540c     &              ga_nodeid(),key_1a3m,offset_1a3m
1541c 203   format('1A3M ',4i5,3x,f19.8,2x,i4,2x,i3,2x,2i8)
1542c       call util_flush(6)
1543c *************
1544        if (sf_write(handle_1a3m_m(my_file_to_write),
1545     &    dfloat(bytes)*dfloat(offset_1a3m),
1546     &    dfloat(bytes)*dfloat(size_1a3m),dbl_mb(k_1a3m),request).ne.0)
1547     &    call errquit('zones put: sf problem21-b',1,DISK_ERR)
1548        if (sf_wait(request).ne.0)
1549     1    call errquit('zones put: sf problem31-b',2,DISK_ERR)
1550       else
1551        call ga_put(i_to_x,offset_1a3m+1,offset_1a3m+size_1a3m,1,1,
1552     1    dbl_mb(k_1a3m),1)
1553       end if ! -------------
1554c
1555      if (.not.ma_pop_stack(l_1a3m))
1556     1  call errquit('tce_mo2e_trans_zones: abc-MA problem',15,MA_ERR)
1557c
1558      next = NXTASK(nprocs, 1)
1559      END IF
1560      count = count + 1
1561      ENDDO
1562      ENDDO
1563      ENDDO
1564      ENDDO
1565c 1A3P fully done here
1566      call ga_sync()
1567      next = nxtask(-nprocs,1)
1568c
1569      cpu = cpu + util_cpusec()
1570      wall = wall + util_wallsec()
1571      if(nodezero) then
1572       write(6,*)'1A3M integrals cpu wall',cpu,wall
1573       call util_flush(6)
1574      end if
1575c
1576c Destroying all l_offset_2a2m_m(n_2a2m_files)
1577c
1578      if(idiskl) then
1579       do i=n_2a2m_files,1,-1
1580        if (.not.ma_pop_stack(l_offset_2a2m_m(i)))
1581     1   call errquit('4-ind n5: l_offset_2a2m_m err.',15,MA_ERR)
1582       enddo
1583      end if
1584c
1585c
1586      if(nodezero) then
1587        write(6,*)'STEP5 4index'
1588        call util_flush(6)
1589      end if
1590c
1591c
1592c
1593c
1594c 4M starts here
1595c
1596c
1597c DISK CHANGES =============== 4M
1598c destroying handle_2a2m
1599c closing/reopening handle_1a3m file
1600c opening handle_4m file
1601      if(idiskl) then
1602c
1603      do i=1,n_2a2m_files
1604       if (sf_destroy(handle_2a2m_m(i)).ne.0)
1605     1  call errquit('tce_sf_destroy: sf problem handle_2a2m',15,MA_ERR)
1606      enddo
1607c
1608      do i=1,n_1a3m_files
1609        if(sf_rwtor(handle_1a3m_m(i)).ne.0)
1610     2    call errquit('4-index: sf_rwtor 1a3m',0,DISK_ERR)
1611      enddo
1612c
1613      end if
1614c ============================
1615c
1616      cpu = - util_cpusec()
1617      wall = - util_wallsec()
1618c
1619      nprocs = GA_NNODES()
1620      count = 0
1621      next = NXTASK(nprocs, 1)
1622c
1623c
1624c Ordering of the do-loops has been changed from:
1625c g1b g2b g3b g4b
1626c to :
1627c g3b g4b g1b g2b (g1b,g2b are used to define file's number)
1628c
1629cx      DO g1b = 1,noa+nva      !nu
1630cx      DO g2b = g1b,noa+nva   !
1631cuu      DO g3b = 1,noa+nva   !k
1632cuu      DO g4b = g3b,noa+nva !l
1633cuu      DO g1b = 1,noa+nva      !nu
1634cuu      DO g2b = g1b,noa+nva   !
1635c ---------------------------------------------------
1636      DO g3b = 1,noa+nva   !k
1637      DO g4b = g3b,noa+nva !l
1638      DO g2b = 1,noa+nva   !
1639      Do g1b = 1,g2b
1640c ---------------------------------------------------
1641      IF (next.eq.count) THEN
1642        my_file_to_read=g1b
1643c((noa+nva)*(noa+nva+1))/2-
1644c     1   ((noa+nva-g1b+1)*(noa+nva-g1b))/2-
1645c     1   (noa+nva-g2b+1)+1
1646      IF (int_mb(k_spin_alpha+g3b-1)+int_mb(k_spin_alpha+g4b-1).eq.
1647     &int_mb(k_spin_alpha+g1b-1)+int_mb(k_spin_alpha+g2b-1)) THEN
1648      IF (ieor(int_mb(k_sym_alpha+g3b-1),ieor(int_mb(k_sym_alpha+g4b-1),
1649     &    ieor(int_mb(k_sym_alpha+g1b-1),int_mb(k_sym_alpha+g2b-1))))
1650     &    .eq. irrep_v) THEN
1651c reversed order
1652      ICOL=INDEX_PAIR(g4b,g3b)
1653      IROW=INDEX_PAIR(g2b,g1b)
1654      IF(IROW.GE.ICOL) THEN
1655      IRES=INDEX_PAIR(IROW,ICOL)
1656      size_4m=int_mb(k_range_alpha+g4b-1)*int_mb(k_range_alpha+g3b-1)
1657     1         *int_mb(k_range_alpha+g2b-1)*int_mb(k_range_alpha+g1b-1)
1658      if (.not.ma_push_get(mt_dbl,size_4m,'4m',l_4m,k_4m))
1659     1    call errquit('tce_4ind: step1_1',0,MA_ERR)
1660      call dfill(size_4m, 0.0d0, dbl_mb(k_4m), 1)
1661c
1662      call a1_offset(g1b,g2b,g4b,size_ic,
1663     &           l_a1_offset,k_a1_offset)
1664      length_a1=int_mb(k_a1_offset)
1665      DO ii=1,length_a1     ! ii---loop
1666       a1_ini    = int_mb(k_a1_offset+ii)
1667       a1_fin    = int_mb(k_a1_offset+length_a1+ii)
1668       key_ini   = a1_ini-1+atpart*(g4b-1+(noa+nva)*
1669     &  (g2b-1+(noa+nva)*(g1b-1)))
1670       call tce_hash(int_mb(k_offset_1a3m_m(my_file_to_read)),
1671     &    key_ini,offset_bl1a3m)
1672c
1673       size_bl1a3m  = int_mb(k_a1_offset+2*length_a1+ii)
1674       if(idiskl) then ! -------
1675        if (sf_read(handle_1a3m_m(my_file_to_read),
1676     &    dfloat(bytes)*dfloat(offset_bl1a3m),
1677     2    dfloat(bytes)*dfloat(size_bl1a3m),
1678     2    dbl_mb(k_integral),request).ne.0)
1679     2    call errquit('zones put: sf problem21-b',1,DISK_ERR)
1680        if (sf_wait(request).ne.0)
1681     1    call errquit('zones put: sf problem31-b',2,DISK_ERR)
1682       else
1683        call ga_get(i_from_x,offset_bl1a3m+1,offset_bl1a3m+size_bl1a3m,
1684     1       1,1,dbl_mb(k_integral),1)
1685       end if ! -------
1686      offset_ismall=0
1687      do i=a1_ini,a1_fin ! over small i ------------------
1688c  C(azone1 g3b)
1689      j=0
1690      tot_azone1_sh=tot_zone(i)
1691      do igi3=1,int_mb(k_range_alpha+g3b-1)
1692      do inu1   =  1,nalength(i)
1693       j=j+1
1694       ipos1=(int_mb(k_offset_alpha+g3b-1)+igi3-1)*nbf+tot_azone1_sh
1695     &       +inu1
1696       dbl_mb(k_coeff+j-1)=dbl_mb(k_movecs_orb+ipos1-1)
1697      enddo
1698      enddo
1699c v(g2b,g1b,g4b,inu)*C(inu g3b)
1700      call dgemm('N','N',
1701     &   int_mb(k_range_alpha+g2b-1)*int_mb(k_range_alpha+g1b-1)
1702     &   *int_mb(k_range_alpha+g4b-1),
1703     &   int_mb(k_range_alpha+g3b-1),
1704     &   nalength(i),
1705     &   1.0d0,dbl_mb(k_integral+offset_ismall),
1706     &   int_mb(k_range_alpha+g2b-1)*int_mb(k_range_alpha+g1b-1)
1707     &  *int_mb(k_range_alpha+g4b-1),
1708     &   dbl_mb(k_coeff),nalength(i),1.0d0,
1709     &   dbl_mb(k_4m),
1710     &   int_mb(k_range_alpha+g2b-1)*int_mb(k_range_alpha+g1b-1)
1711     &  *int_mb(k_range_alpha+g4b-1))
1712c
1713      offset_ismall=offset_ismall+nalength(i)
1714     1       *int_mb(k_range_alpha+g4b-1)
1715     1       *int_mb(k_range_alpha+g1b-1)*int_mb(k_range_alpha+g2b-1)
1716c
1717      enddo ! over small i ----------------------
1718      ENDDO                 ! ii---loop
1719        if (.not.ma_pop_stack(l_a1_offset))
1720     1   call errquit('tce_mo2e_trans_zones: uu l_2g2z',15,MA_ERR)
1721c
1722c zapis
1723c
1724c *** debug ***
1725c       xmax=0.0d0
1726c       do i=1,size_4m
1727c        xmax=xmax+dbl_mb(k_4m+i-1)
1728c       enddo
1729c       write(6,204) ires,xmax
1730c 204   format('4M ',i5,3x,f19.8)
1731c       call util_flush(6)
1732c *************
1733         key_4m=ires
1734         call tce_hash_n(int_mb(kax_v2_alpha_offset),key_4m,offset_4m)
1735      call put_block(d_v2,dbl_mb(k_4m),size_4m,
1736     &               offset_4m)
1737c
1738      if (.not.ma_pop_stack(l_4m))
1739     1  call errquit('tce_mo2e_trans_zones: abc-MA problem',15,MA_ERR)
1740
1741      END IF
1742      END IF
1743      END IF
1744      next = NXTASK(nprocs, 1)
1745      END IF
1746      count = count + 1
1747      ENDDO
1748      ENDDO
1749      ENDDO
1750      ENDDO
1751c
1752       call ga_sync()
1753       next = nxtask(-nprocs,1)
1754c
1755c
1756      cpu = cpu + util_cpusec()
1757      wall = wall + util_wallsec()
1758      if(nodezero) then
1759       write(6,*)'4M integrals cpu wall',cpu,wall
1760       call util_flush(6)
1761      end if
1762c
1763c
1764c
1765c
1766c
1767c
1768c
1769      if(nodezero) then
1770      write(6,*)'DONE --- DONE ---- DONE ---- DONE'
1771      end if
1772c
1773ccx      if (.not.ma_pop_stack(l_work2))
1774ccx     1  call errquit('tce_ao2e: MA problem',14,MA_ERR)
1775ccx      if (.not.ma_pop_stack(l_work1))
1776ccx     1  call errquit('tce_ao2e: MA problem',15,MA_ERR)
1777c
1778c
1779c Destroying all l_offset_1a3m_m(n_1a3m_files)
1780c
1781c
1782      do i=n_1a3m_files,1,-1
1783       if (.not.ma_pop_stack(l_offset_1a3m_m(i)))
1784     1  call errquit('4-ind n5: l_offset_1a3m_m err.',15,MA_ERR)
1785      enddo
1786c
1787      if (.not.ma_pop_stack(l_movecs_orb))
1788     1  call errquit('tce_ao2e: MA problem',15,MA_ERR)
1789c
1790c
1791c DISK CHANGES =============== AFTER 4M
1792c destroying handle_1a3m_m files
1793c
1794      if(idiskl) then
1795      do i=1,n_1a3m_files
1796       if (sf_destroy(handle_1a3m_m(i)).ne.0)
1797     1  call errquit('tce_sf_destroy: sf problem handle_1a3m',15,MA_ERR)
1798      enddo
1799      end if
1800c ============================
1801c
1802c
1803      if(idiskl) then
1804c
1805c dummy part
1806c
1807      else
1808       call deletefile(i_from_x)
1809       call ga_sync()
1810cccx       call deletefile(d_4af)
1811cccx       call ga_sync()
1812      end if
1813c
1814ccx      if (.not.ma_pop_stack(l_1a3m_offset))
1815ccx     1  call errquit('tce_off_1a3m: MA problem',15,MA_ERR)
1816c
1817ccx      if (.not.ma_pop_stack(l_2a2m_offset))
1818ccx     1  call errquit('tce_off_2a2m: MA problem',15,MA_ERR)
1819c
1820ccx      if (.not.ma_pop_stack(l_3a1m_offset))
1821ccx     1  call errquit('tce_off_3a1m: MA problem',15,MA_ERR)
1822c
1823ccx      if (.not.ma_pop_stack(l_4af_offset))
1824ccx     1  call errquit('tce_off_4a: MA problem',15,MA_ERR)
1825c
1826c
1827c     closing two large local files
1828c
1829      if (.not.ma_pop_stack(l_coeff))
1830     1  call errquit('tce_off_4a: MA problem',15,MA_ERR)
1831c
1832      if (.not.ma_pop_stack(l_integral))
1833     1  call errquit('tce_off_4a: MA problem',15,MA_ERR)
1834c
1835
1836c
1837      call ga_sync()
1838c *** debug ***
1839c 800  format('DGEMM1 MAX',i5,2x,3f15.5)
1840c 801  format('DGEMM2 ',i5,2x,3f15.5)
1841 8153 format('SF_READ ERROR CODE = ',2x,i10)
1842 8154 format('SF_READ_WAIT ERROR CODE = ',2x,i10)
1843 8156 format('SF_WRITE ERROR CODE = ',2x,i10)
1844 8157 format('SF_WRITE_WAIT ERROR CODE = ',2x,i10)
1845 8155 format('FU',i6,2x,3i20)
1846 9000 format('PART1',i4,1x,'Cpu  wall ',2(f17.12,1x),3x,'g4b g3b',2i5)
1847c 9001 format('PART2',i4,1x,'Cpu  wall ',2(f17.12,1x),3x,'g4b g3b',2i5)
1848c 9003 format('PART1-4a',i4,1x,'Cpu  wall ',2(f17.12,1x))
1849c 9004 format('PART1-2g2z',i4,1x,'Cpu  wall ',2(f17.12,1x))
1850c 9005 format('PART1-dgemm',i4,1x,'Cpu  wall ',2(f17.12,1x))
1851c 9010 format('  P1-mnrs',i3,1x,2i5,1x,2i5,1x,'Cpu  wall ',2(f17.12,1x))
1852  555  format('atom loop ',2x,i5,3x,2i5,3x,2i5,i12)
1853  556  format('atom time',2x,i5,3x,2i5,3x,2i5,'Cpu wall ',2(f12.7,1x))
1854  777  format('main do loop ',2x,i5,3x,2i5,3x,2i5)
1855  775  format('main loop step1 ',2x,i5,3x,2i5,3x,2i5)
1856  776  format('main loop step2 ',2x,i5,3x,2i5,3x,2i5)
1857  778  format('PART1',2x,i5,3x,2i5,3x,2i5,2x,'Cpu  wall ',2(f17.12,1x))
1858  779  format('PART2',2x,i5,3x,2i5,3x,2i5,2x,'Cpu  wall ',2(f17.12,1x))
1859  780  format('ADD BLOCK',2x,i5,3x,2i5,3x,2i5,2x,'Cpu  wall ',
1860     &        2(f17.12,1x))
1861c *************
1862c
1863      RETURN
1864      END
1865c
1866c
1867c
1868