1      SUBROUTINE cr_eomccsd_t(d_x1,k_x1_offset,d_x2,k_x2_offset,
2     1                        d_t1,k_t1_offset,d_t2,k_t2_offset,
3     2                        d_f1,k_f1_offset,d_v2,k_v2_offset,
4     3                        d_e,k_e_offset,
5     4                        d_ex1,k_ex1_offset,size_ex1,
6     5                        d_ex2,k_ex2_offset,size_ex2,
7     6                        d_c1,k_c1_offset,size_c1,
8     7                        d_c2,k_c2_offset,size_c2,
9     8                        excit,energy1,energy2,size_t1,size_x1,
10     9                        r0xx)
11C
12C     $Id$
13C
14c from this point on we assume that corresponding one- and two-
15c body components of R(d_x1,d_x2) and T(d_t1,d_t2) operators,
16c and corresponding excitation energy are available.
17c objects c1 and c2 are created in tce_energy right before
18c calling this procedure. The same structure of x1 and c1 and
19c x2 and c2 is assumed (irrep_c = irrep_x)
20c
21c 1. calculate R0
22c 2. calculate c1 and c2 vectors (c1excit,c2excit)
23c 3. calculate d
24c
25c local representations go only to "triple" procedures
26c see size_t1 size_x1  in the header of this procedure
27c
28      IMPLICIT NONE
29#include "global.fh"
30#include "mafdecls.fh"
31#include "util.fh"
32#include "errquit.fh"
33#include "tce.fh"
34#include "tce_main.fh"
35      integer d_t1
36      integer k_t1_offset
37      integer d_t2
38      integer k_t2_offset
39c --- sizes ---
40      integer size_c1,size_c2
41      integer size_ex1,size_ex2
42      integer size_d2
43c -------------
44      integer d_x1
45      integer k_x1_offset
46      integer d_x2
47      integer k_x2_offset
48      double precision r0xx  ! r0
49      double precision dr0xx ! r0*r0
50      double precision d1xx  ! <(singles)+singles>
51      double precision d2xx  ! <(doubles)+doubles>
52      double precision d12xx ! <(singles+doubles)+(singles+doubles)>
53      double precision d1xxt,d2xxt
54      double precision d1xxr,d2xxr
55      double precision d1xxtr,d2xxtr
56      double precision excit ! eomsd excitation energy
57      logical lr0            ! (true) r0*M3-calculated
58      integer d_ex1
59      integer k_ex1_offset
60      integer d_ex2
61      integer k_ex2_offset
62      integer d_c1
63      integer k_c1_offset
64      integer d_c2
65      integer k_c2_offset
66      integer d_d2
67      integer l_d2_offset
68      integer k_d2_offset
69c ----
70      integer d_f1
71      integer k_f1_offset
72      integer d_v2
73      integer k_v2_offset
74      integer d_e
75      integer k_e_offset
76      integer t_h1b, t_h1
77      integer t_h2b, t_h2
78      integer t_h3b, t_h3
79      integer t_p4b, t_p4
80      integer t_p5b, t_p5
81      integer t_p6b, t_p6
82      integer k_den,l_den
83      integer k_right,l_right  !r0*M3(T1,T2)
84c --- cr_ccsd_t_E ---
85      integer k_left,l_left
86c -------------------
87c - T1/X1 LOCALIZATION -------------------
88      integer l_t1_local,k_t1_local
89      integer l_x1_local,k_x1_local
90      integer size_t1,size_x1
91c ---------------------------------------
92      integer size,i
93      integer g_energy
94      integer nxtask
95      integer next
96      integer nprocs
97      integer count
98c --- new intermediates ---
99c  cr_ccsd_t_N1 or cr_ccsd_t_N
100      integer d_i1_1,d_i1_2,d_i1_3
101      integer k_i1_offset_1,k_i1_offset_2,k_i1_offset_3
102      integer l_i1_offset_1,l_i1_offset_2,l_i1_offset_3
103c  cr_ccsd_t_E
104      integer d_i1_4,k_i1_offset_4,l_i1_offset_4
105c  cr_ccsd_t_N2
106      integer d_i2_1,d_i2_2,d_i2_3,d_i2_4,d_i2_5,d_i2_6
107      integer k_i2_offset_1,k_i2_offset_2,k_i2_offset_3
108      integer k_i2_offset_4,k_i2_offset_5,k_i2_offset_6
109      integer l_i2_offset_1,l_i2_offset_2,l_i2_offset_3
110      integer l_i2_offset_4,l_i2_offset_5,l_i2_offset_6
111c  q3rexpt
112      integer d_i3_1,d_i3_2
113      integer k_i3_offset_1,k_i3_offset_2
114      integer l_i3_offset_1,l_i3_offset_2
115c ----------------------------
116      double precision energy1,energy2
117      double precision factor
118      double precision den1,num1
119      double precision den2,num2
120      double precision denex
121      double precision cpu,wall
122      character*255 filename
123      logical nodezero
124      external nxtask
125c
126c
127c
128      nodezero=(ga_nodeid().eq.0)
129c
130c Getting R0
131c
132c === jaguar ===
133      if(.not.read_in3) then ! --- read_in3 ---
134      r0xx=0.0d0
135      call  nr0(d_f1,d_e,d_t1,d_v2,d_x1,d_x2,k_f1_offset,
136     &k_e_offset,k_t1_offset,k_v2_offset,k_x1_offset,k_x2_offset)
137      call reconcilefile(d_e,1)
138      call get_block(d_e,r0xx,1,0)
139      if(dabs(excit).gt.1.0d-7) then
140         r0xx = r0xx/excit
141      else
142         write(6,1000)
143      end if
144      end if ! --- read_in3
145      dr0xx = r0xx*r0xx
146      lr0 = .true.
147      if(dabs(r0xx).lt.1.0d-7) lr0 = .false.
148c === jaguar ===
149      if(nodezero) write(6,2000)
150      if(nodezero) write(6,2001) r0xx
151      call util_flush(6)
152c ==============
153c
154c Now on ga with handle d_e we store corresponding R0 value
155c
156c
157c
158c Calculating one- and two-body overlaps
159c
160      if(lr0) then !symmetry of the reference
161         call tce_zero(d_ex1,size_ex1)
162         call tce_zero(d_ex2,size_ex2)
163         call c1_c1(d_t1,d_ex1,k_t1_offset,k_ex1_offset)
164         call reconcilefile(d_ex1,1)
165         d1xxt = 0.0d0
166         call get_block(d_ex1,d1xxt,1,0)
167         call t2t12(d_c2,d_t1,d_t2,k_c2_offset,k_t1_offset,k_t2_offset)
168         call c2_c2(d_c2,d_ex2,k_c2_offset,k_ex2_offset)
169         call reconcilefile(d_ex2,1)
170         d2xxt = 0.0d0
171         call get_block(d_ex2,d2xxt,1,0)
172c forming vector R1T1+R2, the d_d2 irrep corresponds to irrep_x
173c (in this case this is fully symmetric situation)
174c on d_c2 we have (T2+1/2T1*T1)|Phi\rangle
175            irrep_d=irrep_c
176            call tce_zero(d_ex1,size_ex1)
177            call tce_x2_offset(l_d2_offset,k_d2_offset,size_d2)
178            call tce_filename('d2',filename)
179            call createfile(filename,d_d2,size_d2)
180            call tce_zero(d_d2,size_d2)
181         call c2excit2(d_d2,d_t1,d_x1,d_x2,k_d2_offset,k_t1_offset,k_
182     &   x1_offset,k_x2_offset)
183c <C2+ D2>  here
184      call c2_d2(d_c2,d_d2,d_ex1,k_c2_offset,k_d2_offset,k_ex1_offset)
185         call reconcilefile(d_ex1,1)
186         d2xxtr = 0.0d0
187         call get_block(d_ex1,d2xxtr,1,0)
188            call deletefile(d_d2)
189            if (.not.ma_pop_stack(l_d2_offset))
190     1        call errquit('cr_eomccsd_t: MA problem',36,MA_ERR)
191c
192         call tce_zero(d_ex1,size_ex1)
193         call tce_zero(d_ex2,size_ex2)
194         call x1_t1(d_ex1,d_t1,d_x1,k_ex1_offset,k_t1_offset,
195     &   k_x1_offset)
196         call reconcilefile(d_ex1,1)
197         d1xxtr=0.0d0
198         call get_block(d_ex1,d1xxtr,1,0)
199c
200         call tce_zero(d_ex1,size_ex1)
201         call tce_zero(d_ex2,size_ex2)
202         call tce_zero(d_c1,size_c1)
203         call tce_zero(d_c2,size_c2)
204         call c1excit2(d_c1,d_x1,k_c1_offset,k_x1_offset)
205         call c2excit2(d_c2,d_t1,d_x1,d_x2,k_c2_offset,k_t1_offset,k_
206     &   x1_offset,k_x2_offset)
207         call c1_c1(d_c1,d_ex1,k_c1_offset,k_ex1_offset)
208         call reconcilefile(d_ex1,1)
209         d1xxr = 0.0d0
210         call get_block(d_ex1,d1xxr,1,0)
211         call c2_c2(d_c2,d_ex2,k_c2_offset,k_ex2_offset)
212         call reconcilefile(d_ex2,1)
213         d2xxr = 0.0d0
214         call get_block(d_ex2,d2xxr,1,0)
215         d12xx = d1xxt*r0xx*r0xx+d2xxt*r0xx*r0xx+d1xxr+d2xxr+
216     &           2.0d0*d1xxtr*r0xx+2.0d0*d2xxtr*r0xx
217      else         !symmetry different form the symmetry of the reference
218c        <(R1+R0*T1)^{\dagger} (R1+R0*T1)>
219         call c1excit2(d_c1,d_x1,k_c1_offset,k_x1_offset)
220         call c2excit2(d_c2,d_t1,d_x1,d_x2,k_c2_offset,k_t1_offset,k_
221     &   x1_offset,k_x2_offset)
222         call c1_c1(d_c1,d_ex1,k_c1_offset,k_ex1_offset)
223         call reconcilefile(d_ex1,1)
224         d1xx = 0.0d0
225         call get_block(d_ex1,d1xx,1,0)
226c        <(R2+R1T1)^{\dagger} (R2+R1T1)>
227         call c2_c2(d_c2,d_ex2,k_c2_offset,k_ex2_offset)
228         call reconcilefile(d_ex2,1)
229         d2xx = 0.0d0
230         call get_block(d_ex2,d2xx,1,0)
231         d12xx = d1xx+d2xx
232      end if
233c - T1/X1 LOCALIZATION ----------
234c    opening l_t1_local and l_x1_local
235        if (.not.MA_PUSH_GET(mt_dbl,size_t1,'t1_local',
236     1      l_t1_local,k_t1_local))
237     1      call errquit('t1_local',1,MA_ERR)
238        if (.not.MA_PUSH_GET(mt_dbl,size_x1,'x1_local',
239     1      l_x1_local,k_x1_local))
240     1      call errquit('x1_local',1,MA_ERR)
241        call ma_zero(dbl_mb(k_t1_local),size_t1)
242        call ma_zero(dbl_mb(k_x1_local),size_x1)
243c    copy d_t1 ==> l_t1_local
244c    copy x1(ivec) ==> l_x1_local
245#if 1
246        call util_mygabcast(d_t1,size_t1,1,dbl_mb(k_t1_local),
247     L       size_t1)
248        call util_mygabcast(d_x1,size_x1,1,dbl_mb(k_x1_local),
249     L       size_x1)
250#else
251        call ga_get(d_t1,1,size_t1,1,1,dbl_mb(k_t1_local),1)
252        call ga_get(d_x1,1,size_x1,1,1,dbl_mb(k_x1_local),1)
253#endif
254c -------------------------------
255c
256c     Caution! k_right & k_den are not even allocated yet
257c     but they will not be used.
258c --------------- initialization -----------------------------------
259      cpu=-util_cpusec()
260      wall=-util_wallsec()
261      if(lr0) then
262      call cr_ccsd_t_N(dbl_mb(k_right),d_f1,d_i1_1,d_i1_2,
263     1  k_t1_local,d_t2,d_v2,k_f1_offset,k_i1_offset_1,k_i1_offset_2,
264     2  k_t1_offset,k_t2_offset,k_v2_offset,l_i1_offset_1,
265     3  l_i1_offset_2,t_h1b,t_h2b,t_h3b,t_p4b,t_p5b,t_p6b,1)
266      end if
267c
268      call creomsd_t_n2_mem(dbl_mb(k_right),d_f1,d_i2_1,d_i2_2,
269     &d_i2_3,d_i2_4,k_t1_local,d_t2,d_v2,k_x1_local,
270     &d_x2,k_f1_offset,k_i2_offset_1,
271     &k_i2_offset_2,k_i2_offset_3,k_i2_offset_4,k_t1_offset,k_t2_offset,
272     &k_v2_offset,k_x1_offset,k_x2_offset,l_i2_offset_1,l_i2_offset_2,
273     &l_i2_offset_3,l_i2_offset_4,t_h1b,t_h2b,t_h3b,t_p4b,t_p5b,t_p6b,1)
274c
275c
276      if(lr0) then
277      call cr_ccsd_t_E(dbl_mb(k_left),d_i1_4,
278     1  k_t1_local,d_t2,k_i1_offset_4,k_t1_offset,k_t2_offset,
279     2  l_i1_offset_4,t_h1b,t_h2b,t_h3b,t_p4b,t_p5b,t_p6b,1)
280      call q3rexpt2(dbl_mb(k_left),d_i3_1,
281     &k_t1_local,d_t2,k_x1_local,d_x2,
282     &k_i3_offset_1,k_t1_offset,k_t2_offset,k_x1_offset,k_x2_offset,
283     &l_i3_offset_1,t_h1b,t_h2b,t_h3b,t_p4b,t_p5b,t_p6b,1)
284      else
285      call q3rexpt2(dbl_mb(k_left),d_i3_1,k_t1_local,d_t2,
286     &k_x1_local,d_x2,
287     &k_i3_offset_1,k_t1_offset,k_t2_offset,k_x1_offset,k_x2_offset,
288     &l_i3_offset_1,t_h1b,t_h2b,t_h3b,t_p4b,t_p5b,t_p6b,1)
289      end if
290      cpu=cpu+util_cpusec()
291      wall=wall+util_wallsec()
292      if(nodezero) then
293       write(6,*)'         '
294       write(6,111) cpu,wall
295       if(lr0) then
296        write(6,*)'cr_N i1_1 nr of boxes ',int_mb(k_i1_offset_1)
297        write(6,*)'cr_N i1_2 nr of boxes ',int_mb(k_i1_offset_2)
298       end if
299       write(6,*)'creom_N i2_1 nr of boxes ',int_mb(k_i2_offset_1)
300       write(6,*)'creom_N i2_2 nr of boxes ',int_mb(k_i2_offset_2)
301       write(6,*)'creom_N i2_3 nr of boxes ',int_mb(k_i2_offset_3)
302       write(6,*)'creom_N i2_4 nr of boxes ',int_mb(k_i2_offset_4)
303       write(6,*)'         '
304       call util_flush(6)
305      end if
306  111 format('CR-EOMCCSD(T) Intermediates cpu wall time ',2f15.1)
307c ------------------------------------------------------------------
308c
309c     Get the numerator
310c
311      num1 = 0.0d0
312      den1 = 0.0d0
313c
314      num2 = 0.0d0
315      den2 = 0.0d0
316c
317      cpu=-util_cpusec()
318      wall=-util_wallsec()
319c
320      if (.not.ga_create(mt_dbl,1,1,'perturbative',1,1,g_energy))
321     1  call errquit('ccsd_t: GA problem',0,GA_ERR)
322      nprocs = GA_NNODES()
323      count = 0
324      next = nxtask(nprocs,1)
325      do t_p4b = noab+1,noab+nvab
326       do t_p5b = t_p4b,noab+nvab
327        do t_p6b = t_p5b,noab+nvab
328         do t_h1b = 1,noab
329          do t_h2b = t_h1b,noab
330           do t_h3b = t_h2b,noab
331            if (int_mb(k_spin+t_p4b-1)
332     1         +int_mb(k_spin+t_p5b-1)
333     2         +int_mb(k_spin+t_p6b-1)
334     3      .eq.int_mb(k_spin+t_h1b-1)
335     4         +int_mb(k_spin+t_h2b-1)
336     5         +int_mb(k_spin+t_h3b-1)) then
337            if ((.not.restricted).or.
338     1         (int_mb(k_spin+t_p4b-1)
339     1         +int_mb(k_spin+t_p5b-1)
340     2         +int_mb(k_spin+t_p6b-1)
341     3         +int_mb(k_spin+t_h1b-1)
342     4         +int_mb(k_spin+t_h2b-1)
343     5         +int_mb(k_spin+t_h3b-1).le.8)) then
344            if (ieor(int_mb(k_sym+t_p4b-1),
345     1          ieor(int_mb(k_sym+t_p5b-1),
346     2          ieor(int_mb(k_sym+t_p6b-1),
347     3          ieor(int_mb(k_sym+t_h1b-1),
348     4          ieor(int_mb(k_sym+t_h2b-1),
349     5               int_mb(k_sym+t_h3b-1)))))).eq.irrep_x) then
350c
351            if (next.eq.count) then
352c
353c Symmetry control (above)
354c
355            size = int_mb(k_range+t_p4b-1)
356     1           * int_mb(k_range+t_p5b-1)
357     2           * int_mb(k_range+t_p6b-1)
358     3           * int_mb(k_range+t_h1b-1)
359     4           * int_mb(k_range+t_h2b-1)
360     5           * int_mb(k_range+t_h3b-1)
361            if (.not.MA_PUSH_GET(mt_dbl,size,'moment 2,3',
362     1        l_right,k_right)) call errquit('eomccsd_t',3,MA_ERR)
363            if (.not.MA_PUSH_GET(mt_dbl,size,'denominator',
364     1        l_left,k_left)) call errquit('ccsd_t',3,MA_ERR)
365ccx            do i = 1, size
366ccx             dbl_mb(k_right+i-1) = 0.0d0
367ccx             dbl_mb(k_left+i-1) = 0.0d0
368ccx            enddo
369c zeroing ---
370        call dfill(size, 0.0d0, dbl_mb(k_right), 1)
371        call dfill(size, 0.0d0, dbl_mb(k_left), 1)
372c -----------
373c
374c Moments are calculated here
375c
376c
377      if(lr0) then
378      call cr_ccsd_t_N(dbl_mb(k_right),d_f1,d_i1_1,d_i1_2,
379     1  k_t1_local,d_t2,d_v2,k_f1_offset,k_i1_offset_1,k_i1_offset_2,
380     2  k_t1_offset,k_t2_offset,k_v2_offset,l_i1_offset_1,
381     3  l_i1_offset_2,t_h1b,t_h2b,t_h3b,t_p4b,t_p5b,t_p6b,2)
382c
383cc            do i = 1, size
384cc             dbl_mb(k_right+i-1) = r0xx*dbl_mb(k_right+i-1)
385cc            enddo
386        call dscal(size,r0xx,dbl_mb(k_right),1)
387c
388      end if
389c
390      call creomsd_t_n2_mem(dbl_mb(k_right),d_f1,d_i2_1,d_i2_2,
391     &d_i2_3,d_i2_4,k_t1_local,d_t2,d_v2,
392     &k_x1_local,d_x2,k_f1_offset,k_i2_offset_1,
393     &k_i2_offset_2,k_i2_offset_3,k_i2_offset_4,k_t1_offset,k_t2_offset,
394     &k_v2_offset,k_x1_offset,k_x2_offset,l_i2_offset_1,l_i2_offset_2,
395     &l_i2_offset_3,l_i2_offset_4,t_h1b,t_h2b,t_h3b,t_p4b,t_p5b,t_p6b,2)
396c
397c
398c Q3(R0+R1+R2)exp(T1+T2)|Ref> calculated here
399c
400      if(lr0) then
401      call cr_ccsd_t_E(dbl_mb(k_left),d_i1_4,
402     1  k_t1_local,d_t2,k_i1_offset_4,k_t1_offset,k_t2_offset,
403     2  l_i1_offset_4,t_h1b,t_h2b,t_h3b,t_p4b,t_p5b,t_p6b,2)
404c
405cc            do i = 1, size
406cc             dbl_mb(k_left+i-1) = r0xx*dbl_mb(k_left+i-1)
407cc            enddo
408        call dscal(size,r0xx,dbl_mb(k_left),1)
409c
410      call q3rexpt2(dbl_mb(k_left),d_i3_1,k_t1_local,d_t2,
411     &k_x1_local,d_x2,
412     &k_i3_offset_1,k_t1_offset,k_t2_offset,k_x1_offset,k_x2_offset,
413     &l_i3_offset_1,t_h1b,t_h2b,t_h3b,t_p4b,t_p5b,t_p6b,2)
414      else
415      call q3rexpt2(dbl_mb(k_left),d_i3_1,k_t1_local,d_t2,
416     &k_x1_local,d_x2,
417     &k_i3_offset_1,k_t1_offset,k_t2_offset,k_x1_offset,k_x2_offset,
418     &l_i3_offset_1,t_h1b,t_h2b,t_h3b,t_p4b,t_p5b,t_p6b,2)
419      end if
420c
421            if (restricted) then
422              factor = 2.0d0
423            else
424              factor = 1.0d0
425            endif
426            if ((t_p4b.eq.t_p5b).and.(t_p5b.eq.t_p6b)) then
427              factor = factor / 6.0d0
428            else if ((t_p4b.eq.t_p5b).or.(t_p5b.eq.t_p6b)) then
429              factor = factor / 2.0d0
430            endif
431            if ((t_h1b.eq.t_h2b).and.(t_h2b.eq.t_h3b)) then
432              factor = factor / 6.0d0
433            else if ((t_h1b.eq.t_h2b).or.(t_h2b.eq.t_h3b)) then
434              factor = factor / 2.0d0
435            endif
436c
437c
438c
439            i = 0
440            do t_p4 = 1, int_mb(k_range+t_p4b-1)
441             do t_p5 = 1, int_mb(k_range+t_p5b-1)
442              do t_p6 = 1, int_mb(k_range+t_p6b-1)
443               do t_h1 = 1, int_mb(k_range+t_h1b-1)
444                do t_h2 = 1, int_mb(k_range+t_h2b-1)
445                 do t_h3 = 1, int_mb(k_range+t_h3b-1)
446                  i = i + 1
447            denex=-dbl_mb(k_evl_sorted+int_mb(k_offset+t_p4b-1)+t_p4-1)
448     4         -dbl_mb(k_evl_sorted+int_mb(k_offset+t_p5b-1)+t_p5-1)
449     5         -dbl_mb(k_evl_sorted+int_mb(k_offset+t_p6b-1)+t_p6-1)
450     6         +dbl_mb(k_evl_sorted+int_mb(k_offset+t_h1b-1)+t_h1-1)
451     7         +dbl_mb(k_evl_sorted+int_mb(k_offset+t_h2b-1)+t_h2-1)
452     8         +dbl_mb(k_evl_sorted+int_mb(k_offset+t_h3b-1)+t_h3-1)
453     9         +excit
454c numerator
455                  num1 = num1 + factor *
456     1            (dbl_mb(k_right+i-1)*dbl_mb(k_right+i-1))/denex
457                  num1 = num1 + factor *
458     1            dbl_mb(k_left+i-1)*dbl_mb(k_right+i-1)
459c denominator
460c
461                  den1 = den1 + factor *
462     1            (dbl_mb(k_left+i-1)*dbl_mb(k_right+i-1))/denex
463                  den1 = den1 + factor *
464     1            (dbl_mb(k_left+i-1)*dbl_mb(k_left+i-1))
465c
466                 enddo
467                enddo
468               enddo
469              enddo
470             enddo
471            enddo
472c
473c
474            if (.not.MA_POP_STACK(l_left))
475     1        call errquit('eomccsd_t',6,MA_ERR)
476            if (.not.MA_POP_STACK(l_right))
477     1        call errquit('eomccsd_t',6,MA_ERR)
478c
479            next = nxtask(nprocs,1)
480            endif
481            count = count + 1
482c
483            endif
484            endif
485            endif
486           enddo
487          enddo
488         enddo
489        enddo
490       enddo
491      enddo
492      next = nxtask(-nprocs,1)
493c
494      cpu=cpu+util_cpusec()
495      wall=wall+util_wallsec()
496      if(nodezero) then
497       write(6,*)'         '
498       write(6,112) cpu,wall
499       write(6,*)'         '
500       call util_flush(6)
501      end if
502  112 format('CR-EOMCCSD(T) triples loop cpu wall ',2f15.1)
503c --- toggle = 3 ---
504      if(lr0) then
505      call q3rexpt2(dbl_mb(k_left),d_i3_1,k_t1_local,d_t2,
506     &k_x1_local,d_x2,
507     &k_i3_offset_1,k_t1_offset,k_t2_offset,k_x1_offset,k_x2_offset,
508     &l_i3_offset_1,t_h1b,t_h2b,t_h3b,t_p4b,t_p5b,t_p6b,3)
509      call cr_ccsd_t_E(dbl_mb(k_left),d_i1_4,
510     1  k_t1_local,d_t2,k_i1_offset_4,k_t1_offset,k_t2_offset,
511     2  l_i1_offset_4,t_h1b,t_h2b,t_h3b,t_p4b,t_p5b,t_p6b,3)
512      else
513      call q3rexpt2(dbl_mb(k_left),d_i3_1,k_t1_local,d_t2,
514     &k_x1_local,d_x2,
515     &k_i3_offset_1,k_t1_offset,k_t2_offset,k_x1_offset,k_x2_offset,
516     &l_i3_offset_1,t_h1b,t_h2b,t_h3b,t_p4b,t_p5b,t_p6b,3)
517      end if
518c
519      call creomsd_t_n2_mem(dbl_mb(k_right),d_f1,d_i2_1,d_i2_2,
520     &d_i2_3,d_i2_4,k_t1_local,d_t2,d_v2,
521     &k_x1_local,d_x2,k_f1_offset,k_i2_offset_1,
522     &k_i2_offset_2,k_i2_offset_3,k_i2_offset_4,k_t1_offset,k_t2_offset,
523     &k_v2_offset,k_x1_offset,k_x2_offset,l_i2_offset_1,l_i2_offset_2,
524     &l_i2_offset_3,l_i2_offset_4,t_h1b,t_h2b,t_h3b,t_p4b,t_p5b,t_p6b,3)
525c
526c
527      if(lr0) then
528      call cr_ccsd_t_N(dbl_mb(k_right),d_f1,d_i1_1,d_i1_2,
529     1  k_t1_local,d_t2,d_v2,k_f1_offset,k_i1_offset_1,k_i1_offset_2,
530     2  k_t1_offset,k_t2_offset,k_v2_offset,l_i1_offset_1,
531     3  l_i1_offset_2,t_h1b,t_h2b,t_h3b,t_p4b,t_p5b,t_p6b,3)
532      end if
533c ------------------
534c
535#if 1
536      call ga_dgop(1975,num1,1,'+')
537      call ga_dgop(1976,num2,1,'+')
538      call ga_dgop(1977,den1,1,'+')
539      call ga_dgop(1978,den2,1,'+')
540#else
541      call ga_zero(g_energy)
542      call ga_acc(g_energy,1,1,1,1,num1,1,1.0d0)
543      call ga_sync()
544      call ga_get(g_energy,1,1,1,1,num1,1)
545c
546      call ga_zero(g_energy)
547      call ga_acc(g_energy,1,1,1,1,den1,1,1.0d0)
548      call ga_sync()
549      call ga_get(g_energy,1,1,1,1,den1,1)
550c
551      call ga_zero(g_energy)
552      call ga_acc(g_energy,1,1,1,1,num2,1,1.0d0)
553      call ga_sync()
554      call ga_get(g_energy,1,1,1,1,num2,1)
555c
556      call ga_zero(g_energy)
557      call ga_acc(g_energy,1,1,1,1,den2,1,1.0d0)
558      call ga_sync()
559      call ga_get(g_energy,1,1,1,1,den2,1)
560#endif
561c
562      if (.not.ga_destroy(g_energy))
563     1  call errquit('creom_t: GA problem',1,GA_ERR)
564       energy1 = num1/(dr0xx+d12xx+den1)
565      call util_flush(6)
566c - T1/X1 LOCALIZATION ------
567         if(.not.MA_POP_STACK(l_x1_local))
568     &      call errquit('l_x1_local',4,MA_ERR)
569         if(.not.MA_POP_STACK(l_t1_local))
570     &      call errquit('l_t1_local',4,MA_ERR)
571c ---------------------------
572 1000    format('corresponding excitation energy = 0')
573 2000    format('======= CR-EOMCCSD(T) calculation ======')
574 2001    format('corresponding R0 value = ',f12.6)
575c
576      return
577      end
578c
579
580