1      subroutine uccsdt_makex(g_omega,g_x_big,spini,spinj)
2csh   subroutine uccsdt_makex(g_omega,d_x,spini,spinj)
3c
4c     x(k,l,i,j) = Sum(u,v) t(u,v,i,j)C(u,k,spini)C(v,l,spinj)
5c     where u & v are SO's (no spin), k,b,i,j are MO's with spin labels
6c
7      implicit none
8#include "mafdecls.fh"
9#include "cuccsdtP.fh"
10#include "amplitudes.fh"
11#include "global.fh"
12      integer g_x_big      ! [output] DRA handle for X(k,l,i,j)
13      integer g_omega      ! [input] GA handle for Omega(u,v,i,j)
14      integer spini, spinj ! [input] Spins of i and j
15      integer spink, spinl ! Spins of k and l
16      integer g_tmp        ! Temporary GA x(u,l,i,j) & x(u,v,i,j)
17      integer g_x          ! Temporary GA x(k,l,i,j)
18      integer k_buf,l_buf  ! Temporary MA for tranformation
19      integer k_tmp,l_tmp  ! Temporary MA for tranformation
20      integer k_c,l_c      ! Temporary MA for MO coeffs
21      integer lenij,lenul,lenuv,lenu,lenv,lenl
22      integer lenik,lenjl,leni,lenk,leniu
23      integer maxlenuv
24      integer symij,symj,symi,symul,syml
25      integer symuv,symu,symv,symjl,symik,symk
26      integer joff(nbf,0:7)
27      integer loff(nbf,0:7)
28      integer symvoff(0:7,0:7)
29      integer symloff(0:7,0:7)
30      integer ij,jl,ik
31      integer ijlo,ijhi,jllo,jlhi
32      integer i,j,k,l
33      integer llo,lhi,vlo,klo,ilo,ihi,ulo,uhi
34      integer dummy
35csh
36      logical odebug
37      odebug = .false.
38csh
39c
40csh
41      if (odebug) then
42      write(*,*) 'no = ',no(1),no(2)
43      write(*,*) 'nv = ',nv(1),nv(2)
44      write(*,*) 'no_sym = ',(no_sym(i,1),i=0,7),(no_sym(i,2),i=0,7)
45      write(*,*) 'nv_sym = ',(nv_sym(i,1),i=0,7),(nv_sym(i,2),i=0,7)
46      write(*,*) 'o_sym(a) = ',(o_sym(1,i,1),i=0,7),(o_sym(2,i,1),i=0,7)
47      write(*,*) 'o_sym(b) = ',(o_sym(1,i,2),i=0,7),(o_sym(2,i,2),i=0,7)
48      write(*,*) 'v_sym(a) = ',(v_sym(1,i,1),i=0,7),(v_sym(2,i,1),i=0,7)
49      write(*,*) 'v_sym(b) = ',(v_sym(1,i,2),i=0,7),(v_sym(2,i,2),i=0,7)
50      write(*,*) 'bf_per_ir = ',(bf_per_ir(i),i=0,7)
51      endif
52csh
53      spink = spini
54      spinl = spinj
55csh ***** We have to create an empty g_b_big
56      dummy = 0
57      do symjl = 0,7
58      do syml = 0,7
59      symj = ieor(symjl,syml)
60      do l = 1,no_sym(syml,spinl)
61      do j = 1,no_sym(symj,spinj)
62      do symk = 0,7
63      symi = ieor(symjl,symk)
64      do k = 1,no_sym(symk,spink)
65      do i = 1,no_sym(symi,spini)
66      dummy = dummy + 1
67      enddo
68      enddo
69      enddo
70      enddo
71      enddo
72      enddo
73      enddo
74      if (odebug) write(*,*) 'Size of g_x_big = ',dummy
75      if (.not.ga_create(mt_dbl,dummy,1,'x_big',
76     $ -1,1,g_x_big)) call errquit
77     $ ('uccsdt_makex: room for x_big?',dummy)
78csh
79c
80c     Addressing & sizes
81c
82      lenij = 0
83      do symij = 0,7
84       do symj = 0,7
85        symi = ieor(symij,symj)
86        do j = o_sym(1,symj,spinj),o_sym(2,symj,spinj)
87         joff(j,symij) = lenij
88         lenij = lenij + no_sym(symi,spini)
89        end do
90       end do
91      end do
92      if (odebug) write(*,*) 'lenij = ',lenij
93c
94c     Addressing & sizes
95c
96      do symul = 0,7
97       lenul = 0
98       do syml = 0,7
99        symu = ieor(symul,syml)
100        do l = o_sym(1,syml,spinl),o_sym(2,syml,spinl)
101         loff(l,symul) = lenul
102         lenul = lenul + bf_per_ir(symu)
103        end do
104       end do
105      end do
106c
107c     Addressing & sizes
108c
109      do symul = 0,7
110       lenul = 0
111       do syml = 0,7
112        symu = ieor(symul,syml)
113        symloff(syml,symul) = lenul
114        lenul = bf_per_ir(symu)*no_sym(syml,spinl)
115       end do
116      end do
117c
118c     Addressing & sizes
119c
120      do symuv = 0,7
121       lenuv = 0
122       do symv = 0,7
123        symu = ieor(symuv,symv)
124        symvoff(symv,symuv) = lenuv
125        lenuv = lenuv + bf_per_ir(symu)*bf_per_ir(symv)
126       end do
127       if (lenuv .gt. maxlenuv) maxlenuv = lenuv
128      end do
129      if (odebug) write(*,*) 'maxlenuv = ',maxlenuv
130c
131c     Transformed MO coefficients
132c
133      if (.not. ma_push_get(mt_dbl, nbf*no(spinl), 'c',
134     $ l_c, k_c))
135     $ call errquit('ma? nbf*no',nbf*no(spinl),0)
136      call ga_get(g_part(spinl),1,nbf,1,no(spinl),
137     $ dbl_mb(k_c),nbf)
138c
139c     Allocate temporary GA
140c
141      if (.not.ga_create(mt_dbl,maxlenuv,lenij,'X tmp',
142     $ maxlenuv,-1,g_tmp)) call errquit
143     $ ('uccsdt_makex: room for tmp?',maxlenuv*lenij)
144c
145c     Data parallel transform Omega(u,v,i,j)C(v,l) = tmp(u,l,i,j)
146c     (tmp was allocated big enough to do this).
147c
148      call ga_distribution(g_omega,ga_nodeid(),ijlo,ijhi,dummy,dummy)
149      ij = 0
150      do symij = 0,7
151       symuv = symij
152       symul = symij
153       do symj = 0,7
154        symi = ieor(symij,symj)
155        do j = o_sym(1,symj,spinj),o_sym(2,symj,spinj)
156         do i = o_sym(1,symi,spini),o_sym(2,symi,spini)
157          ij = ij + 1
158          if (ij.ge.ijlo .and. ij.le.ijhi) then
159           do syml = 0,7
160            symu = ieor(symij,syml)
161            symv = syml
162            lenl = no_sym(syml,spinl)
163            lenu = bf_per_ir(symu)
164            lenv = bf_per_ir(symv)
165            lenul = lenu * lenl
166            lenuv = lenu * lenv
167            llo = o_sym(1,syml,spinl)
168            vlo = bf_per_ir_cum(symv) + 1
169csh
170            if (odebug) write(*,*) 'symij,symj,j,i,syml = ',
171     $                  symij,symj,j,i,syml
172csh
173            if (lenuv.gt.0 .and. lenl.gt.0) then
174csh
175             if (odebug) write(*,*) 'lenul,lenuv = ',lenul,lenuv
176csh
177             if (.not. ma_push_get(mt_dbl, lenuv, 'buf',
178     $        l_buf, k_buf))
179     $       call errquit('ma? lenuv',lenuv,0)
180             if (.not. ma_push_get(mt_dbl, lenul, 'tmp',
181     $        l_tmp, k_tmp))
182     $        call errquit('ma? lenul',lenul,0)
183c
184c     (uv) -> (ul) for l (occupied) & v (SO) of same irrep
185c
186             call ga_get(g_omega,1,lenuv,ij,ij,
187     $        dbl_mb(k_buf),lenu)
188csh
189             if (odebug) call ma_print(dbl_mb(k_buf),lenu,lenv,'UV')
190csh
191             call dgemm('n','n',lenu,lenl,lenv,
192     $        1.0d0,dbl_mb(k_buf),lenu,
193     $        dbl_mb(k_c),nbf,
194     $        0.0d0,dbl_mb(k_tmp),lenu)
195csh
196             if (odebug) call ma_print(dbl_mb(k_tmp),lenu,lenl,'UL')
197             if (odebug) write(*,*) 'ga_put in the range ... ',
198     $       symloff(syml,symul)+1,'~',symloff(syml,symul)+lenul,
199     $       ij,'~',ij
200csh
201             call ga_put(g_tmp,
202     $        symloff(syml,symul)+1,symloff(syml,symul)+lenul,
203     $        ij,ij,dbl_mb(k_tmp),lenu)
204             if (.not. ma_pop_stack(l_tmp))
205     $        call errquit('ma_pop?',l_tmp,0)
206             if (.not. ma_pop_stack(l_buf))
207     $        call errquit('ma_pop?',l_buf,0)
208            end if
209           end do
210          end if
211         end do
212        end do
213       end do
214      end do
215csh
216      if (odebug) write(*,*) '**************************'
217      if (odebug) write(*,*) 'transformation v -> l done'
218      if (odebug) write(*,*) '**************************'
219csh
220c
221c     Perform quarter transformation from tmp(u,l,i,j) to t(k,l,i,j)
222c     t(i,k,symk,j,l,syml,symjl)
223c
224      if (.not. ma_pop_stack(l_c))
225     $ call errquit('ma_pop?',l_c,0)
226      if (.not. ma_push_get(mt_dbl, nbf*no(spink), 'c',
227     $ l_c, k_c))
228     $ call errquit('ma? nbf*no',nbf*no(spink),0)
229      call ga_get(g_part(spink),1,nbf,1,no(spink),
230     $ dbl_mb(k_c),nbf)
231      do symjl = 0,7
232       symik = symjl
233       do syml = 0,7
234        symj = ieor(symjl,syml)
235c
236c     Allocate X block
237c
238        llo = o_sym(1,syml,spinl)
239        lhi = o_sym(2,syml,spinl)
240        lenik = 0
241        do symk = 0,7
242         symi = ieor(symik,symk)
243         lenik = lenik + no_sym(symi,spini) * no_sym(symk,spink)
244        end do
245        lenjl = (lhi-llo+1)*no_sym(symj,spinj)
246        if ((lenik.gt.0).and.(lenjl.gt.0)) then
247csh
248         if (odebug) write(*,*) 'lenik,lenjl = ',lenik,lenjl
249csh
250         if (.not.ga_create(mt_dbl,lenik,lenjl,'l',
251     $    lenik,-1,g_x)) call errquit
252     $    ('uccsdt_makex: room for x?',lenik*lenjl)
253c
254c     Data-parallel transformation of k
255c
256         call ga_distribution(g_x,ga_nodeid(),jllo,jlhi,dummy,dummy)
257         jl = 0
258         do l = llo, lhi
259          do j = o_sym(1,symj,spinj),o_sym(2,symj,spinj)
260           jl = jl + 1
261           if (jl.ge.jllo .and. jl.le.jlhi) then
262            ik = 0
263c     Compare the next line with the comments in Amplitude.F
264c     I think the comments are not accurate
265            do symk = 0, 7
266             symi = ieor(symik,symk)
267             symij = ieor(symi,symj)
268             symu = symk
269             symul = ieor(symu,syml)
270             lenk = no_sym(symk,spink)
271             leni = no_sym(symi,spini)
272             lenu = bf_per_ir(symu)
273             lenik = leni*lenk
274             leniu = leni*lenu
275             klo = o_sym(1,symk,spink)
276             ilo = o_sym(1,symi,spini)
277             ihi = o_sym(2,symi,spini)
278             ulo = bf_per_ir_cum(symu) + 1
279             uhi = ulo + lenu - 1
280             if (lenu.gt.0) then
281              if (.not. ma_push_get(mt_dbl, leniu, 'buf',
282     $         l_buf, k_buf))
283     $         call errquit('ma? leniu',leniu,0)
284              if (.not. ma_push_get(mt_dbl, lenik, 'tmp',
285     $         l_tmp, k_tmp))
286     $         call errquit('ma? lenik',lenik,0)
287c
288c     (iu) -> (ik) for k (occupied) & u (SO) of same irrep
289c
290              ij = joff(j,symij)
291              do i = ilo,ihi
292               ij = ij + 1
293csh
294               if (odebug) write(*,*) 'i = ',i
295               if (odebug) write(*,*) 'get from range ... ',
296     $          loff(l,symul)+1,'~',loff(l,symul)+lenu,
297     $          ij,'~',ij,
298     $          ' to k_buf+',lenu*(i-ilo)
299csh
300               call ga_get(g_tmp,
301     $          loff(l,symul)+1,loff(l,symul)+lenu,
302     $          ij,ij,dbl_mb(k_buf+lenu*(i-ilo)),lenu)
303              end do
304csh
305              if (odebug) call ma_print(dbl_mb(k_buf),lenu,leni,'X(UI)')
306csh
307              call dgemm('t','n',leni,lenk,lenu,
308     $         1.0d0,dbl_mb(k_buf),lenu,
309     $         dbl_mb(k_c),nbf,
310     $         0.0d0,dbl_mb(k_tmp),leni)
311csh
312              if (odebug) call ma_print(dbl_mb(k_tmp),leni,lenk,'X(IK)')
313              if (odebug) write(*,*) 'writing to range ... ',
314     $        ik+1,'~',ik+lenik,jl,'~',jl
315csh
316              call ga_put(g_x,ik+1,ik+lenik,jl,jl,
317     $         dbl_mb(k_tmp),leni)
318              ik = ik + lenik
319              if (.not. ma_pop_stack(l_tmp))
320     $         call errquit('ma_pop?',l_tmp,0)
321              if (.not. ma_pop_stack(l_buf))
322     $         call errquit('ma_pop?',l_buf,0)
323             end if
324            end do
325           end if
326          end do
327         end do
328c
329c     Write one symmetry block of X to disk
330c
331         if (spini.eq.spinj) call ga_scale(g_x,2.0d0)
332         if (odebug) call ga_print(g_x)
333         call sh_put_x(g_x,syml,symjl,spini,spinj,g_x_big)
334csh      if (.not. uccsdt_ampfile_write_t2(d_x,
335csh  $    spini, spink, spinj, spinl, symik,
336csh  $    llo, lhi, g_x))
337csh  $    call errquit('write_t2 failed',0)
338         if (.not. ga_destroy(g_x)) call errquit('GA 999',0,0)
339        end if
340       end do
341      end do
342      if (.not. ga_destroy(g_tmp)) call errquit('GA 999',0,0)
343      if (.not. ma_pop_stack(l_c))
344     $ call errquit('ma_pop?',l_c,0)
345c
346      return
347      end
348c $Id$
349