1      logical function jantest(rtdb)
2*
3*     $$
4*
5      implicit none
6#include "errquit.fh"
7#include "global.fh"
8#include "mafdecls.fh"
9#include "bas.fh"
10#include "geom.fh"
11#include "rtdb.fh"
12#include "inp.fh"
13      integer rtdb
14c
15      integer basis, geom, nbf
16      integer l_aoint, k_aoint        ! nbf**4 array of AO integrals
17      integer l_moint, k_moint        ! nbf**4 array of AO integrals
18      character*255 movecs      ! Name of movector file
19      character*80 title, name_of_basis, scftype
20      integer nbf_file, nsets, nmo_file(2)
21      logical movecs_read, movecs_read_header
22      external movecs_read, movecs_read_header
23c
24      integer nmo, g_tmp, l_occ, k_occ,
25     $     l_eval, k_eval, l_mos, k_mos, l_most, k_most
26      integer nocc, nvirt, nopen, nclosed, nso, noso
27      integer l_t, k_t, l_t2, k_t2, l_f, k_f
28      integer l_t1, k_t1, l_r1, k_r1
29      integer l_tg, k_tg, l_tf, k_tf
30      integer l_w, k_w, l_v, k_v
31      logical int_normalize
32      external int_normalize
33c
34c     load the geometry/basis set and get info
35c
36      if (.not. geom_create(geom, 'geometry'))
37     $     call errquit('scf_init: geom_create?', 0, GEOM_ERR)
38      if (.not. geom_rtdb_load(rtdb, geom, 'geometry'))
39     $     call errquit('scf_init: no geometry ', 0, RTDB_ERR)
40      if (.not. bas_create(basis, 'ao basis'))
41     $     call errquit('scf_init: bas_create?', 0, BASIS_ERR)
42      if (.not. bas_rtdb_load(rtdb, geom, basis, 'ao basis'))
43     $        call errquit('scf_init: no ao basis set', 0, RTDB_ERR)
44      if (.not.int_normalize(rtdb,basis))
45     $           call errquit('scf:int_normalize failed', 0, INT_ERR)
46      if (.not. bas_numbf(basis, nbf)) call errquit
47     $     ('scf_init: basis info',0, BASIS_ERR)
48c
49c     Read the MO vectors and evals from a RHF calculation
50c
51      call util_file_name('movecs',.false.,.false.,movecs)
52      if (.not. movecs_read_header(movecs, title, name_of_basis,
53     $     scftype, nbf_file, nsets, nmo_file, 2)) call errquit
54     $     ('jantest: failed to read movecs header',911, DISK_ERR)
55      write(6,*) ' Read movecs header from ', movecs
56      write(6,*) ' Job title :                ',
57     $     title(1:inp_strlen(title))
58      write(6,*) ' Basis name:                ',
59     $     name_of_basis(1:inp_strlen(name_of_basis))
60      nmo = nmo_file(1)
61      if (.not. rtdb_get(rtdb, 'scf:nclosed', mt_int, 1, nocc))
62     $     call errquit('nocc?',0, RTDB_ERR)
63      if (.not. rtdb_get(rtdb, 'scf:nopen', mt_int, 1, nopen))
64     $     call errquit('nopen?',0, RTDB_ERR)
65      if (nopen .ne. 0) call errquit('asjdlfkadjsl',0, UNKNOWN_ERR)
66      nvirt= nmo - nocc
67      write(6,*) ' No. of closed shells       ', nocc
68      write(6,*) ' No. of molecular orbitals: ', nmo
69      write(6,*) ' No. of basis functions:    ', nbf
70c
71*ga:1:0
72      if (.not. ga_create(mt_dbl, nbf, nmo, 'tmp', 0, 0, g_tmp))
73     &     call errquit('scf_v_g: tmp', 0, GA_ERR)
74      if (.not. ma_push_get(mt_dbl, nbf,'occ',l_occ, k_occ))
75     $     call errquit('ma occ', nbf, MA_ERR)
76      if (.not. ma_push_get(mt_dbl, nbf,'eval',l_eval, k_eval))
77     $     call errquit('ma eval', nbf, MA_ERR)
78      if (.not. ma_push_get(mt_dbl, nbf*nbf,'mos', l_mos, k_mos))
79     $     call errquit('ma mos', nbf*nbf, MA_ERR)
80      if (.not. ma_push_get(mt_dbl, nbf*nbf,'mos', l_most, k_most))
81     $     call errquit('ma mos', nbf*nbf, MA_ERR)
82c
83      if (.not. movecs_read(movecs, 1, dbl_mb(k_occ), dbl_mb(k_eval),
84     $     g_tmp)) call errquit('movecs_read of amos failed ',0,
85     &       DISK_ERR)
86      call ga_get(g_tmp, 1, nbf, 1, nmo, dbl_mb(k_mos), nbf)
87      call util_transpose(dbl_mb(k_mos),nbf,dbl_mb(k_most),nmo,
88     $     nbf,nmo)
89c
90      write(6,*) ' Orbital eigenvalues '
91      call output(dbl_mb(k_eval),1,nmo,1,1,nmo,1,1)
92      write(6,*) ' MOs'
93      call output(dbl_mb(k_mos),1,nbf,1,nmo,nbf,nmo,1)
94      write(6,*) ' MOs T'
95      call output(dbl_mb(k_most),1,nmo,1,nbf,nmo,nbf,1)
96c
97      if (.not. ga_destroy(g_tmp)) call errquit(' ga bad?',0, GA_ERR)
98c
99c     Make all AO integrals
100c
101      if (.not. ma_push_get(mt_dbl,nbf**4,'aoint',l_aoint,k_aoint))
102     $     call errquit('allocation of AO integrals failed',nbf**4,
103     &       MA_ERR)
104      call jan_all_ao_integrals(rtdb,basis,nbf,'dirac',dbl_mb(k_aoint))
105c     call jan_debug_print('AOINTS',dbl_mb(k_aoint), nbf,  nbf,  nbf,
106c    $     nbf)
107c
108c     Make all MO integrals in Dirac order
109c
110      if (.not. ma_push_get(mt_dbl,nmo**4,'moint',l_moint,k_moint))
111     $     call errquit('allocation of MO integrals failed',nbf**4,
112     &       MA_ERR)
113      call jan_full_transform(
114     $     rtdb, basis,
115     $     nmo, nmo, nmo, nmo,
116     $     nmo, nmo, nmo, nmo,
117     $     dbl_mb(k_most),dbl_mb(k_most),dbl_mb(k_most),dbl_mb(k_most),
118     $     dbl_mb(k_moint), 'Dirac')
119c      call jan_debug_print('MOINTS',dbl_mb(k_moint), nmo,  nmo,  nmo,
120c     $     nmo)
121c
122c     do some incore cc
123c     set nso=2*nmo, noso=2*nocc
124c
125      nso=2*nmo
126      noso=2*nocc
127      if (.not. ma_push_get(mt_dbl,nso**4,'t amps',l_t,k_t))
128     $     call errquit('allocation of t amplitudes failed',nso**4,
129     &       MA_ERR)
130      if (.not. ma_push_get(mt_dbl,nso**4,'t2 amps',l_t2,k_t2))
131     $     call errquit('allocation of t2 amplitudes failed',nso**4,
132     &       MA_ERR)
133      if (.not. ma_push_get(mt_dbl,nmo**2,'Fock Matrix',l_f,k_f))
134     $     call errquit('allocation of t2 amplitudes failed',nmo**2,
135     &       MA_ERR)
136      if (.not. ma_push_get(mt_dbl,nso**2,'t1 amps',l_t1,k_t1))
137     $     call errquit('allocation of t1 amplitudes failed',nso**2,
138     &       MA_ERR)
139      if (.not. ma_push_get(mt_dbl,nso**2,'t1 resi',l_r1,k_r1))
140     $     call errquit('allocation of t1 residual failed',nso**2,
141     &       MA_ERR)
142c
143      call ccsd_incore(rtdb, basis, dbl_mb(k_moint),
144     &                 dbl_mb(k_eval), dbl_mb(k_t), dbl_mb(k_t2),
145     &                 dbl_mb(k_f), dbl_mb(k_t1), dbl_mb(k_r1),
146     &                 nbf, nmo, nocc, nso, noso)
147c
148      if (.not. ma_push_get(mt_dbl,noso**3*nso**3,'t3 tg',l_tg,k_tg))
149     $     call errquit('allocation of t3 tg failed',nso**6,
150     &       MA_ERR)
151      if (.not. ma_push_get(mt_dbl,noso**3*nso**3,'t3 tf',l_tg,k_tf))
152     $     call errquit('allocation of t3 tf failed',nso**6, MA_ERR)
153      if (.not. ma_push_get(mt_dbl,noso**3*nso**3,'t3 w',l_tg,k_w))
154     $     call errquit('allocation of t3 w failed',nso**6, MA_ERR)
155      if (.not. ma_push_get(mt_dbl,noso**3*nso**3,'t3 v',l_tg,k_v))
156     $     call errquit('allocation of t3 v failed',nso**6, MA_ERR)
157c
158      if (.not. ga_create(mt_dbl, nbf, nmo, 'tmp', 0, 0, g_tmp))
159     &     call errquit('scf_v_g: tmp', 0, GA_ERR)
160      if (.not. movecs_read(movecs, 1, dbl_mb(k_occ), dbl_mb(k_eval),
161     $     g_tmp)) call errquit('movecs_read of amos failed ',0,
162     &       DISK_ERR)
163      call ga_get(g_tmp, 1, nbf, 1, nmo, dbl_mb(k_mos), nbf)
164      call util_transpose(dbl_mb(k_mos),nbf,dbl_mb(k_most),nmo,
165     $     nbf,nmo)
166      if (.not. ga_destroy(g_tmp)) call errquit(' ga bad?',0, GA_ERR)
167c
168      call jan_full_transform(
169     $     rtdb, basis,
170     $     nmo, nmo, nmo, nmo,
171     $     nmo, nmo, nmo, nmo,
172     $     dbl_mb(k_most),dbl_mb(k_most),dbl_mb(k_most),dbl_mb(k_most),
173     $     dbl_mb(k_moint), 'Dirac')
174c
175      call triples_incore(rtdb, basis,
176     &                 dbl_mb(k_moint),
177     &                 dbl_mb(k_eval), dbl_mb(k_t2),
178     &                 dbl_mb(k_t1),
179     &                 dbl_mb(k_tg), dbl_mb(k_tf),
180     &                 dbl_mb(k_w), dbl_mb(k_v),
181     &                 nbf, nmo, nocc, nso, noso)
182c
183c     Tidy up
184c
185      if (.not. ma_chop_stack(l_occ)) call errquit(' ma chop?', 0,
186     &       MA_ERR)
187      if (.not. bas_destroy(basis)) call errquit(' bas ?',0, BASIS_ERR)
188      if (.not. geom_destroy(geom)) call errquit(' geom ?',0, GEOM_ERR)
189      jantest = .true.
190c
191      end
192      subroutine jan_all_ao_integrals(rtdb, basis, nbf, order, ao)
193      implicit none
194#include "errquit.fh"
195#include "mafdecls.fh"
196#include "bas.fh"
197      integer rtdb, basis, nbf
198      double precision ao(nbf,nbf,nbf,nbf)
199      character*(*) order
200c
201      integer nsh, k_i, k_j, k_k, k_l, l_i, l_j, l_k, l_l,
202     $     maxg2, maxs2, k_buf, l_buf, k_scr, l_scr
203c
204      call int_init(rtdb, 1, basis)
205      if ( .not. bas_numcont(basis, nsh) ) call errquit(
206     $     'ao_fock_2e: problem with call to bas_numcont', basis,
207     &       BASIS_ERR)
208      call int_mem_2e4c(maxg2,maxs2)
209      if (.not. ma_push_get(mt_dbl,maxs2,'scr',l_scr, k_scr))
210     $   call errquit('ma scr',maxg2, MA_ERR)
211      if (.not. ma_push_get(mt_dbl,maxg2,'buf',l_buf, k_buf))
212     $   call errquit('ma buf',maxg2, MA_ERR)
213      if (.not. ma_push_get(mt_int,maxg2,'i',l_i, k_i))
214     $   call errquit('ma ibuf',maxg2, MA_ERR)
215      if (.not. ma_push_get(mt_int,maxg2,'j',l_j, k_j))
216     $   call errquit('ma jbuf',maxg2, MA_ERR)
217      if (.not. ma_push_get(mt_int,maxg2,'k',l_k, k_k))
218     $   call errquit('ma kbuf',maxg2, MA_ERR)
219      if (.not. ma_push_get(mt_int,maxg2,'l',l_l, k_l))
220     $   call errquit('ma lbuf',maxg2, MA_ERR)
221c
222      call jan_do_all_ao_integrals(basis, dbl_mb(k_buf), dbl_mb(k_scr),
223     $     int_mb(k_i), int_mb(k_j), int_mb(k_k), int_mb(k_l),
224     $     maxg2, maxs2, nbf, nsh, order, ao)
225c
226      if (.not. ma_chop_stack(l_scr)) call errquit('janallao: ma?',0,
227     &       MA_ERR)
228c
229      call int_terminate
230c
231      end
232      subroutine jan_do_all_ao_integrals(
233     $     basis, buf, scr, ilab, jlab, klab,
234     $  llab, maxg2, maxs2, nbf, nsh, order, ao)
235      implicit none
236#include "errquit.fh"
237c
238      integer basis, nbf, nsh, maxg2, maxs2
239      double precision buf(maxg2), scr(maxs2)
240      integer ilab(maxg2), jlab(maxs2), klab(maxs2), llab(maxs2)
241      integer i, j, k ,l, ish, jsh, ksh, lsh, ijkl, nint
242      character*(*) order
243      double precision ao(nbf,nbf,nbf,nbf)
244      double precision zerotol
245      logical omulliken
246c
247      omulliken = .false.   ! avoids compiler warning
248      if (order .eq. 'mulliken') then
249         omulliken = .true.
250      else if (order .eq. 'dirac') then
251         omulliken = .false.
252      else
253         call errquit(' unknown order',0, UNKNOWN_ERR)
254      end if
255c
256      call dfill(nbf**4, 0.0d0, ao, 1)
257      zerotol = 1d-12
258c
259      do ish = 1, nsh
260         do jsh = 1, nsh
261            do ksh = 1, nsh
262               do lsh = 1, nsh
263                  call int_l2e4c(basis, ish, jsh, basis, ksh, lsh,
264     &                 zerotol, .false., maxg2, buf, nint,
265     $                 ilab, jlab, klab, llab, maxs2, scr)
266                  do ijkl = 1, nint
267                     i = ilab(ijkl)
268                     j = jlab(ijkl)
269                     k = klab(ijkl)
270                     l = llab(ijkl)
271                     if (omulliken) then
272                        ao(i,j,k,l) = buf(ijkl)
273                     else
274                        ao(i,k,j,l) = buf(ijkl)
275                     end if
276                  end do
277               end do
278            end do
279         end do
280      end do
281c
282c$$$      write(6,*)
283c$$$      write(6,*) ' AO integrals '
284c$$$      write(6,*)
285c$$$      do i = 1, nbf
286c$$$         do j = 1, nbf
287c$$$            do k = 1, nbf
288c$$$               do l = 1, nbf
289c$$$                  if ( abs(ao(i,j,k,l)) .gt. 1e-6 )
290c$$$     $                 write(6,7) i,j,k,l,ao(i,j,k,l)
291c$$$ 7                format(1x,4i5,2x,f12.6)
292c$$$               end do
293c$$$            end do
294c$$$         end do
295c$$$      end do
296c
297      end
298      subroutine jan_full_transform(
299     $     rtdb, basis,
300     $     n1, n2, n3, n4,
301     $     ld1, ld2, ld3, ld4,
302     $     c1t, c2t, c3t, c4t,
303     $     full, order)
304      implicit none
305#include "errquit.fh"
306#include "schwarz.fh"
307#include "bas.fh"
308#include "mafdecls.fh"
309#include "inp.fh"
310c
311      integer rtdb
312      integer basis             ! AO basis handle
313      integer n1, n2, n3, n4    ! Dimension of each MO set
314      integer ld1, ld2, ld3, ld4
315      double precision c1t(ld1,*), c2t(ld2,*), ! Transposed MO coeffs
316     $     c3t(ld3,*), c4t(ld4,*)
317      double precision full(n1,n2,n3,n4)
318      character*(*) order
319c
320c     Generate the specified block of MO integrals with
321c     no assumptions of equivalence between the sets of coefficients.
322c
323c     Order can be either
324c     .    ChargeCloud -> full(p,q,r,s) = (pq|rs)
325c     or
326c     .          Dirac -> full(p,q,r,s) = <pq|rs>
327c     or
328c     .  LeftAsymDirac -> full(p,q,r,s) = <pq|rs>-<qp|rs>
329c     .  (must have c1t=c2t, n1=n2)
330c     or
331c     . RightAsymDirac -> full(p,q,r,s) = <pq|rs>-<pq|sr>
332c     .  (must have c3t=c4t, n3=n4)
333c
334c     Presently the antisymmetrization is done at the top level
335c     and the storage of full is not reduced to use the symmetry.
336c
337c     Memory requirements are
338c     .  n1*n2*n3*n4 + S*n2*n3*n4 + S*S*n3*n4 + S*S*S*n4 +
339c     .  maxs2 + maxg2*(1 + 4*integer)
340c
341c     Index 4 is the first transformed so there is advantage in
342c     making it the smallest range.
343c
344      double precision tol2e
345      parameter (tol2e = 1d-12)
346C
347      integer nsh, maxbfsh, lenhalf, lenthird, geom
348      integer l_half, k_half, l_third, k_third
349      integer lsh, ksh, llo, lhi, klo, khi
350      integer p, q, r, s
351      logical ochargecloud, oasym
352      character*8 side
353c
354      ochargecloud = .false.   ! avoids compiler warning
355      oasym = .false.          ! avoids compiler warning
356      if (inp_compare(.false.,order,'chargecloud')) then
357         ochargecloud = .true.
358         oasym = .false.
359         side  = ' '
360      else if (inp_compare(.false.,order,'dirac')) then
361         ochargecloud = .false.
362         oasym = .false.
363         side  = ' '
364      else if (inp_compare(.false.,order,'leftasymdirac')) then
365         ochargecloud = .false.
366         oasym = .true.
367         side  = 'left'
368      else if (inp_compare(.false.,order,'rightasymdirac')) then
369         ochargecloud = .false.
370         oasym = .true.
371         side  = 'right'
372      else
373         call errquit('jan_full_trans: unkown integral option',0,
374     &       INT_ERR)
375      endif
376c
377c     Initialize integrals and Schwarz screening
378c
379      if (.not. bas_geom(basis, geom))
380     $     call errquit('jan_transform: basis ', basis, BASIS_ERR)
381      call int_init(rtdb, 1, basis)
382      call schwarz_init(geom, basis)
383c
384      if (.not. bas_numcont(basis, nsh)) call errquit(
385     $     'jan_transform: bas_numcont', basis, BASIS_ERR)
386      if (.not. bas_nbf_cn_max(basis,maxbfsh)) call errquit(
387     $     'jan_transform: bas_nbf_cn_max', basis, BASIS_ERR)
388c
389      lenhalf = n3*n4*maxbfsh**2
390      lenthird= n2*n3*n4*maxbfsh
391c
392      if (.not. ma_push_get(mt_dbl,lenhalf,'half',l_half, k_half))
393     $     call errquit('ma half', lenhalf, MA_ERR)
394      if (.not. ma_push_get(mt_dbl,lenthird,'third',l_third, k_third))
395     $     call errquit('ma third', lenthird, MA_ERR)
396c
397      call dfill(n1*n2*n3*n4, 0.0d0, full, 1)
398      do ksh = 1, nsh
399         if (.not. bas_cn2bfr(basis, ksh, klo, khi))
400     $        call errquit('jan_transform: bas_cn2bfr',basis, BASIS_ERR)
401         call dfill(n2*n3*n4*(khi-klo+1), 0.0d0, dbl_mb(k_third), 1)
402         do lsh = 1, nsh
403            if (.not. bas_cn2bfr(basis, lsh, llo, lhi))
404     $           call errquit('jan_transform: bas_cn2bfr',basis,
405     &       BASIS_ERR)
406            if (schwarz_shell(ksh,lsh)*schwarz_max()
407     $           .gt. tol2e) then
408c
409c     Make (rs|kl) all rs (indices 3 and 4) given shells k and l
410c
411               call jan_half_transform(basis, ksh, lsh, n3, n4,
412     $              c3t, c4t, ld3, ld4,
413     $              dbl_mb(k_half), ochargecloud, tol2e)
414
415*               write(6,*) ' ksh, lsh ', ksh, lsh
416*               call jan_debug_print('half',
417*     $              dbl_mb(k_half), n3, n4, khi-klo+1, lhi-llo+1)
418
419c
420               call jan_third_transform(llo, lhi, klo, khi,
421     $              n2, n3, n4, dbl_mb(k_half), dbl_mb(k_third),
422     $              c2t, ld2, tol2e)
423            end if
424         end do
425*         write(6,*) ' lsh ', lsh
426*         call jan_debug_print('third',
427*     $        dbl_mb(k_third), n2, n3, n4, lhi-llo+1)
428         call jan_final_transform(klo, khi, n1, n2, n3, n4,
429     $        dbl_mb(k_third), full, c1t, ld1, tol2e)
430      end do
431c
432      if (oasym) call jan_asym_trans(full,n1,n2,n3,n4,side)
433c
434      do s = 1, n4
435         do r = 1, n3
436            do q = 1, n2
437               do p = 1, n1
438                  if (abs(full(p,q,r,s)).lt.1d-10)
439     $                 full(p,q,r,s) = 0.0d0
440               end do
441            end do
442         end do
443      end do
444c
445      if (.not. ma_pop_stack(l_third)) call errquit('ma third',0,
446     &       MA_ERR)
447      if (.not. ma_pop_stack(l_half)) call errquit('ma half',0,
448     &       MA_ERR)
449c
450      call schwarz_tidy()
451      call int_terminate
452c
453      end
454      subroutine jan_full_transform_noinit(
455     $     rtdb, basis,
456     $     n1, n2, n3, n4,
457     $     ld1, ld2, ld3, ld4,
458     $     c1t, c2t, c3t, c4t,
459     $     full, order)
460      implicit none
461#include "errquit.fh"
462#include "schwarz.fh"
463#include "bas.fh"
464#include "mafdecls.fh"
465#include "inp.fh"
466c
467      integer rtdb
468      integer basis             ! AO basis handle
469      integer n1, n2, n3, n4    ! Dimension of each MO set
470      integer ld1, ld2, ld3, ld4
471      double precision c1t(ld1,*), c2t(ld2,*), ! Transposed MO coeffs
472     $     c3t(ld3,*), c4t(ld4,*)
473      double precision full(n1,n2,n3,n4)
474      character*(*) order
475c
476c     Generate the specified block of MO integrals with
477c     no assumptions of equivalence between the sets of coefficients.
478c
479c     Order can be either
480c     .    ChargeCloud -> full(p,q,r,s) = (pq|rs)
481c     or
482c     .          Dirac -> full(p,q,r,s) = <pq|rs>
483c     or
484c     .  LeftAsymDirac -> full(p,q,r,s) = <pq|rs>-<qp|rs>
485c     .  (must have c1t=c2t, n1=n2)
486c     or
487c     . RightAsymDirac -> full(p,q,r,s) = <pq|rs>-<pq|sr>
488c     .  (must have c3t=c4t, n3=n4)
489c
490c     Presently the antisymmetrization is done at the top level
491c     and the storage of full is not reduced to use the symmetry.
492c
493c     Memory requirements are
494c     .  n1*n2*n3*n4 + S*n2*n3*n4 + S*S*n3*n4 + S*S*S*n4 +
495c     .  maxs2 + maxg2*(1 + 4*integer)
496c
497c     Index 4 is the first transformed so there is advantage in
498c     making it the smallest range.
499c
500      double precision tol2e
501      parameter (tol2e = 1d-12)
502C
503      integer nsh, maxbfsh, lenhalf, lenthird, geom
504      integer l_half, k_half, l_third, k_third
505      integer lsh, ksh, llo, lhi, klo, khi
506      integer p, q, r, s
507      logical ochargecloud, oasym
508      character*8 side
509c
510      ochargecloud = .false.   ! avoids compiler warning
511      oasym = .false.          ! avoids compiler warning
512      if (inp_compare(.false.,order,'chargecloud')) then
513         ochargecloud = .true.
514         oasym = .false.
515         side  = ' '
516      else if (inp_compare(.false.,order,'dirac')) then
517         ochargecloud = .false.
518         oasym = .false.
519         side  = ' '
520      else if (inp_compare(.false.,order,'leftasymdirac')) then
521         ochargecloud = .false.
522         oasym = .true.
523         side  = 'left'
524      else if (inp_compare(.false.,order,'rightasymdirac')) then
525         ochargecloud = .false.
526         oasym = .true.
527         side  = 'right'
528      else
529         call errquit('jan_full_trans: unkown integral option',0,
530     &       INT_ERR)
531      endif
532c
533c     Initialize integrals and Schwarz screening
534c
535      if (.not. bas_geom(basis, geom))
536     $     call errquit('jan_transform: basis ', basis, BASIS_ERR)
537*      call int_init(rtdb, 1, basis)
538*      call schwarz_init(geom, basis)
539c
540      if (.not. bas_numcont(basis, nsh)) call errquit(
541     $     'jan_transform: bas_numcont', basis, BASIS_ERR)
542      if (.not. bas_nbf_cn_max(basis,maxbfsh)) call errquit(
543     $     'jan_transform: bas_nbf_cn_max', basis, BASIS_ERR)
544c
545      lenhalf = n3*n4*maxbfsh**2
546      lenthird= n2*n3*n4*maxbfsh
547c
548      if (.not. ma_push_get(mt_dbl,lenhalf,'half',l_half, k_half))
549     $     call errquit('ma half', lenhalf, MA_ERR)
550      if (.not. ma_push_get(mt_dbl,lenthird,'third',l_third, k_third))
551     $     call errquit('ma third', lenthird, MA_ERR)
552c
553      call dfill(n1*n2*n3*n4, 0.0d0, full, 1)
554      do ksh = 1, nsh
555         if (.not. bas_cn2bfr(basis, ksh, klo, khi))
556     $        call errquit('jan_transform: bas_cn2bfr',basis, BASIS_ERR)
557         call dfill(n2*n3*n4*(khi-klo+1), 0.0d0, dbl_mb(k_third), 1)
558         do lsh = 1, nsh
559            if (.not. bas_cn2bfr(basis, lsh, llo, lhi))
560     $           call errquit('jan_transform: bas_cn2bfr',basis,
561     &       BASIS_ERR)
562            if (schwarz_shell(ksh,lsh)*schwarz_max()
563     $           .gt. tol2e) then
564c
565c     Make (rs|kl) all rs (indices 3 and 4) given shells k and l
566c
567               call jan_half_transform(basis, ksh, lsh, n3, n4,
568     $              c3t, c4t, ld3, ld4,
569     $              dbl_mb(k_half), ochargecloud, tol2e)
570
571*               write(6,*) ' ksh, lsh ', ksh, lsh
572*               call jan_debug_print('half',
573*     $              dbl_mb(k_half), n3, n4, khi-klo+1, lhi-llo+1)
574
575c
576               call jan_third_transform(llo, lhi, klo, khi,
577     $              n2, n3, n4, dbl_mb(k_half), dbl_mb(k_third),
578     $              c2t, ld2, tol2e)
579            end if
580         end do
581*         write(6,*) ' lsh ', lsh
582*         call jan_debug_print('third',
583*     $        dbl_mb(k_third), n2, n3, n4, lhi-llo+1)
584         call jan_final_transform(klo, khi, n1, n2, n3, n4,
585     $        dbl_mb(k_third), full, c1t, ld1, tol2e)
586      end do
587c
588      if (oasym) call jan_asym_trans(full,n1,n2,n3,n4,side)
589c
590      do s = 1, n4
591         do r = 1, n3
592            do q = 1, n2
593               do p = 1, n1
594                  if (abs(full(p,q,r,s)).lt.1d-10)
595     $                 full(p,q,r,s) = 0.0d0
596               end do
597            end do
598         end do
599      end do
600c
601      if (.not. ma_pop_stack(l_third)) call errquit('ma third',0,
602     &       MA_ERR)
603      if (.not. ma_pop_stack(l_half)) call errquit('ma half',0,
604     &       MA_ERR)
605c
606*      call schwarz_tidy()
607*      call int_terminate
608c
609      end
610      subroutine jan_asym_trans(full,n1,n2,n3,n4,side)
611      implicit none
612#include "errquit.fh"
613c
614      integer n1, n2, n3, n4
615      double precision full(n1,n2,n3,n4)
616      character*(*) side
617c
618      integer p, q, r, s
619      double precision tmp
620c
621      if (side .eq. 'left') then
622         if (n1 .ne. n2) call errquit('jan_asym_trans: left', n1,
623     &       UNKNOWN_ERR)
624         do s = 1, n4
625            do r = 1, n3
626               do q = 1, n2
627                  do p = 1, q
628                     tmp = full(p,q,r,s) - full(q,p,r,s)
629                     full(p,q,r,s) = tmp
630                     full(q,p,r,s) =-tmp
631                  enddo
632               enddo
633            enddo
634         enddo
635      else
636         if (n3 .ne. n4) call errquit('jan_asym_trans: right', n3,
637     &       UNKNOWN_ERR)
638         do s = 1, n4
639            do r = 1, s
640               do q = 1, n2
641                  do p = 1, n1
642                     full(p,q,r,s) = full(p,q,r,s) - full(p,q,s,r)
643                  enddo
644               enddo
645               if (r .ne. s) then
646                  do q = 1, n2
647                     do p = 1, n1
648                        full(p,q,s,r) = -full(p,q,r,s)
649                     enddo
650                  enddo
651               else
652                  do q = 1, n2
653                     do p = 1, n1
654                        full(p,q,s,r) = 0.0d0
655                     enddo
656                  enddo
657               endif
658            enddo
659         enddo
660      endif
661c
662      end
663      subroutine jan_final_transform(klo, khi,
664     $     n1, n2, n3, n4, third, full, c1t, ld1, tol2e)
665      implicit none
666c
667      integer klo, khi, n1, n2, n3, n4, ld1
668      double precision third(n2,n3,n4,klo:khi)
669      double precision full(n1,n2,n3,n4)
670      double precision c1t(ld1,*)
671      double precision tol2e
672c
673      integer k, s, r, q, p
674      double precision g
675c
676      do k = klo, khi
677         do s = 1, n4
678            do r = 1, n3
679               do q = 1, n2
680                  g = third(q,r,s,k)
681                  if (abs(g) .gt. tol2e) then
682                     do p = 1, n1
683                        full(p,q,r,s) = full(p,q,r,s) + g*c1t(p,k)
684                     end do
685                  end if
686               end do
687            end do
688         end do
689      end do
690c
691      end
692      subroutine jan_third_transform(llo, lhi, klo, khi,
693     $     n2, n3, n4, half, third, c2t, ld2, tol2e)
694      implicit none
695c
696      integer llo, lhi, klo, khi, n2, n3, n4, ld2
697      double precision half(n3,n4,klo:khi,llo:lhi)
698      double precision third(n2,n3,n4,klo:khi)
699      double precision c2t(ld2,*)
700      double precision tol2e
701c
702      integer k, l, s, r, q
703      double precision g
704c
705      do l = llo, lhi
706         do k = klo, khi
707            do s = 1, n4
708               do r = 1, n3
709                  g = half(r,s,k,l)
710                  if (abs(g) .gt. tol2e) then
711                     do q  = 1, n2
712                        third(q,r,s,k) = third(q,r,s,k) + g*c2t(q,l)
713                     end do
714                  end if
715               end do
716            end do
717         end do
718      end do
719c
720      end
721      subroutine jan_half_transform(basis, ksh, lsh, n1, n2,
722     $     c1t, c2t, ld1, ld2, half, ochargecloud, tol2e)
723      implicit none
724#include "errquit.fh"
725#include "bas.fh"
726#include "mafdecls.fh"
727      integer basis, ksh, lsh, n1, n2, ld1, ld2
728      double precision c1t(*), c2t(*) ! Transposed MO coeffs
729      double precision half(*)  ! n1*n2*kdim*ldim (pq|kl)
730      logical ochargecloud
731      double precision tol2e
732c
733c     For a pair of shells k and l fill
734c
735c     .   half(p,q,k,l) = (pq|kl)
736c
737c     for all p, q, and k, l within their respective shells
738c     where p and q are transformed into the new bases and
739c     k and l are AO indices.
740c
741c     Eventually exploiting sparsity and abelian symmetry
742c     ... should also eventually use the texas integrals
743c
744c     Assumes that integrals and schwarz have been initialized.
745c
746      integer nbf, nsh, k_i, k_j, k_k, k_l, l_i, l_j, l_k, l_l,
747     $     maxg2, maxs2, k_buf, l_buf, k_scr, l_scr, maxbfsh
748      integer llo, lhi, klo, khi, lenaobuf, l_aobuf, k_aobuf
749c
750c     Get dimensions and required scratch space info
751c
752      if (.not. bas_numcont(basis, nsh)) call errquit(
753     $     'jan_transform: bas_numcont', basis, BASIS_ERR)
754      if (.not. bas_nbf_cn_max(basis,maxbfsh)) call errquit(
755     $     'jan_transform: bas_nbf_cn_max', basis, BASIS_ERR)
756      if (.not. bas_numbf(basis, nbf))
757     $     call errquit('jan_transform: nbf',basis, BASIS_ERR)
758      if (.not. bas_cn2bfr(basis, ksh, klo, khi))
759     $     call errquit('jan_transform: bas_cn2bfr',basis, BASIS_ERR)
760      if (.not. bas_cn2bfr(basis, lsh, llo, lhi))
761     $     call errquit('jan_transform: bas_cn2bfr',basis, BASIS_ERR)
762      call int_mem_2e4c(maxg2,maxs2)
763      lenaobuf = (khi-klo+1)*(lhi-llo+1)*maxbfsh*n2 ! (iq|kl)
764c
765c     Allocate scratch space for integrals and buffers for
766c     transformation
767c
768      if (.not. ma_push_get(mt_dbl,maxs2,'scr',l_scr, k_scr))
769     $   call errquit('ma scr',maxs2, MA_ERR)
770      if (.not. ma_push_get(mt_dbl,maxg2,'buf',l_buf, k_buf))
771     $   call errquit('ma buf',maxg2, MA_ERR)
772      if (.not. ma_push_get(mt_int,maxg2,'i',l_i, k_i))
773     $   call errquit('ma ibuf',maxg2, MA_ERR)
774      if (.not. ma_push_get(mt_int,maxg2,'j',l_j, k_j))
775     $   call errquit('ma jbuf',maxg2, MA_ERR)
776      if (.not. ma_push_get(mt_int,maxg2,'k',l_k, k_k))
777     $   call errquit('ma kbuf',maxg2, MA_ERR)
778      if (.not. ma_push_get(mt_int,maxg2,'l',l_l, k_l))
779     $   call errquit('ma lbuf',maxg2, MA_ERR)
780      if (.not. ma_push_get(mt_dbl,lenaobuf,'aobuf',l_aobuf, k_aobuf))
781     $   call errquit('ma aobuf',lenaobuf, MA_ERR)
782c
783      call jan_do_half_transform(
784     $     basis,
785     $     dbl_mb(k_buf), dbl_mb(k_scr),
786     $     int_mb(k_i), int_mb(k_j), int_mb(k_k), int_mb(k_l),
787     $     maxg2, maxs2,
788     $     nbf, nsh, maxbfsh, n1, n2, tol2e,
789     $     ksh, lsh, klo, khi, llo, lhi,
790     $     half, c1t, c2t, ld1, ld2, dbl_mb(k_aobuf), ochargecloud)
791c
792      if (.not. ma_chop_stack(l_scr)) call errquit
793     $     ('jan_transform: chopping stack', 0, MA_ERR)
794c
795      end
796      subroutine jan_do_half_transform(
797     $     basis,
798     $     buf, scr, ilab, jlab, klab, llab, maxg2, maxs2,
799     $     nbf, nsh, maxbfsh, n1, n2, tol2e,
800     $     ksh, lsh, klo, khi, llo, lhi,
801     $     half, c1t, c2t, ld1, ld2, aobuf, ochargecloud)
802      implicit none
803#include "errquit.fh"
804#include "schwarz.fh"
805#include "bas.fh"
806      integer basis
807      integer maxg2, maxs2
808      integer nbf, nsh, n1, n2, maxbfsh, ld1, ld2
809      double precision buf(maxg2), scr(maxs2)
810      double precision c1t(ld1,nbf), c2t(ld2,nbf)
811      integer ilab(maxg2), jlab(maxg2), klab(maxg2), llab(maxg2)
812      integer ksh, lsh, klo, khi, llo, lhi
813      double precision aobuf(n2,maxbfsh,klo:khi,llo:lhi)
814      double precision half(n1,n2,klo:khi,llo:lhi)
815      double precision tol2e
816      logical ochargecloud
817c
818c     For a pair of shells k and l, fill aobuf with integrals (ij|kl)
819c     for all i>=j ... eventually exploiting sparsity and abelian symmetry
820c     ... should also use the texas integrals
821c
822      double precision skl, g
823      integer ish, jsh, i, j, k, l, p, q, ijkl, nint, ilo, ihi,
824     $     idim, kdim, ldim
825c
826      kdim = khi - klo + 1
827      ldim = lhi - llo + 1
828      skl = schwarz_shell(ksh,lsh)
829      call dfill(n1*n2*kdim*ldim,0.0d0,half,1)
830c
831      do ish = 1, nsh
832         if (.not. bas_cn2bfr(basis,ish,ilo,ihi))
833     $        call errquit('jan_do_half_transform',ish, BASIS_ERR)
834         idim = ihi-ilo+1
835         do l = llo,lhi
836            do k = klo, khi
837               do i = 1, idim
838                  do q = 1, n2
839                     aobuf(q,i,k,l) = 0.0d0
840                  end do
841               end do
842            end do
843         end do
844c
845         do jsh = 1, nsh
846            if (ochargecloud) then ! (ij|kl)
847               if (schwarz_shell(ish,jsh)*skl .gt. tol2e) then
848                  call int_l2e4c(basis, ish, jsh, basis, ksh, lsh,
849     &                 tol2e, .false., maxg2, buf, nint,
850     $                 ilab, jlab, klab, llab, maxs2, scr)
851                  do ijkl = 1, nint
852                     i = ilab(ijkl)-ilo+1
853                     j = jlab(ijkl)
854                     k = klab(ijkl)
855                     l = llab(ijkl)
856                     g = buf(ijkl)
857                     if (abs(g) .gt. tol2e) then
858                        do q = 1, n2
859                           aobuf(q,i,k,l) = aobuf(q,i,k,l) + g*c2t(q,j)
860                        end do
861                     end if
862                  end do
863               end if
864            else                ! <ij|kl> = (ik|jl)
865               if (schwarz_shell(ish,ksh)*schwarz_shell(jsh,lsh)
866     $              .gt. tol2e) then
867                  call int_l2e4c(basis, ish, ksh, basis, jsh, lsh,
868     &                 tol2e, .false., maxg2, buf, nint,
869     $                 ilab, klab, jlab, llab, maxs2, scr)
870                  do ijkl = 1, nint
871                     i = ilab(ijkl)-ilo+1
872                     j = jlab(ijkl)
873                     k = klab(ijkl)
874                     l = llab(ijkl)
875                     g = buf(ijkl)
876                     if (abs(g) .gt. tol2e) then
877                        do q = 1, n2
878                           aobuf(q,i,k,l) = aobuf(q,i,k,l) + g*c2t(q,j)
879                        end do
880                     end if
881                  end do
882               end if
883            endif
884         end do
885         do l = llo, lhi
886            do k = klo, khi
887               do i = ilo, ihi
888                  do q = 1, n2
889                     g = aobuf(q,i-ilo+1,k,l)
890                     if (abs(g) .gt. tol2e) then
891                        do p = 1, n1
892                           half(p,q,k,l) = half(p,q,k,l) + g*c1t(p,i)
893                        end do
894                     end if
895                  end do
896               end do
897            end do
898         end do
899      end do
900c
901      end
902      subroutine jan_debug_print(string,full, n1, n2, n3, n4)
903      implicit none
904      character*(*) string
905      integer n1, n2, n3, n4
906      double precision full(n1, n2, n3, n4)
907c
908      integer p, q, r, s
909      write(6,*) ' DEBUG FOR ', string, n1, n2, n3, n4
910      do s = 1, n4
911         do r = 1, n3
912            do q = 1, n2
913               do p = 1, n1
914                  if (abs(full(p,q,r,s)).gt.1e-6) then
915                     write(6,1) p,q,r,s,full(p,q,r,s)
916 1                   format(4i5,2x,f14.8)
917                  end if
918               end do
919            end do
920         end do
921      end do
922c
923      end
924      subroutine ccsd_incore(rtdb, basis, g_mo, e, r2, t2, fock,
925     &                       t1, r1, nbf, nmo, nocc, nso, noso)
926      implicit none
927      integer rtdb, basis, nbf, nmo, nocc, nso, noso
928c
929c     nbf = number of basis functions
930c     nmo = number of molecular orbitals
931c     nocc = number of occupied orbitals
932c     nso = number of spin orbitals (for now set to 2*nmo)
933c     noso = number of occupied spin orbitals (for now set to 2*nocc)
934c
935c     toy coupled cluster program based on the green book
936c
937      double precision g_mo(nmo,nmo,nmo,nmo)
938      double precision e(nmo), fock(nmo,nmo)
939      double precision r2(nso,nso,nso,nso)
940      double precision t2(nso,nso,nso,nso)
941      double precision t1(nso,nso)
942      double precision r1(nso,nso)
943      integer n
944c
945c     cc stuff (greek indices denote occupied spin-orbitals)
946c     cc stuff (r, s, t, u, ... denote unoccupied spin-orbitals)
947c
948      write(6,*)' nbf, nmo, nocc, nso, noso: ',
949     &            nbf, nmo, nocc, nso, noso
950c
951c     Initially fill the Fock matrix with diagonal elements only
952c
953      call dfill(nmo**2, 0.0d0, fock, 1)
954      do n = 1, nmo
955        fock(n,n) = e(n)
956      enddo
957c      write(6,*) ' Fock Matrix'
958c      call output(fock,1,nmo,1,nmo,nmo,nmo,1)
959c
960c     First generate an initial guess for t2
961c
962      call dfill(nso**4, 0.0d0, r2, 1)
963      call dfill(nso**4, 0.0d0, t2, 1)
964      call dfill(nso**2, 0.0d0, r1, 1)
965      call dfill(nso**2, 0.0d0, t1, 1)
966c
967c      call t2_init(g_mo, e, t, nmo, nso, noso)
968c      write(6,*)' T2 initial guess. '
969c      call writet2(t,nso)
970c      call correlation(g_mo,t,nmo,nso,noso)
971c
972c
973c     next put this guess into t2 expression keeping all terms linear in t2
974c
975c      call t2_l_gb(g_mo, e, t, t2, nmo, nso, noso)
976c      do n = 1, 3
977c      call dfill(nso**4, 0.0d0, t2, 1)
978c      call t2_l_rjh(g_mo, t, t2, fock, nmo, nso, noso)
979c      call t2_update_rjh(g_mo, t, t2, fock, nmo, nso, noso)
980c      write(6,*)' T2 with linear terms. '
981c      call writet2(t,nso)
982c      call correlation(g_mo,t,nmo,nso,noso)
983c      enddo
984c
985c     next put this into t2 expression keeping all terms
986c
987      do n = 1, 50
988         call dfill(nso**4, 0.0d0, r2, 1)
989         call dfill(nso**2, 0.0d0, r1, 1)
990c         call t2_l_q_gb(g_mo, e, t, t2, nmo, nso, noso)
991c         call t2_l_rjh(g_mo, t, t2, fock, nmo, nso, noso)
992c
993c        Generate delta T1
994c
995         call t1_rjh(g_mo, r1, t2, fock, nmo, nso, noso)
996c
997c        Generate delta T2
998c
999         call t2_l_q_rjh(g_mo, t2, r2, fock, nmo, nso, noso)
1000c
1001c        Update T1
1002c
1003         call t1_update_rjh(g_mo, t1, r1, fock, nmo, nso, noso)
1004c         write(6,*)' T1: '
1005c         call writet1(t1,nso)
1006c
1007c        Update T2
1008c
1009         call t2_update_rjh(g_mo, t2, r2, fock, nmo, nso, noso)
1010c         write(6,*)' T2 with quadratic terms. '
1011c         call writet2(t2,nso)
1012c
1013c        Compute Energy expression
1014c
1015         call correlation(g_mo,t1,t2,nmo,nso,noso)
1016c
1017c        Use T1 to transform 1e- and 2e- ints
1018c
1019         call get_new_f_g(rtdb, basis, fock, g_mo, t1, nso, nmo)
1020      enddo
1021c
1022      return
1023      end
1024      subroutine triples_incore(rtdb, basis, g_mo, ener, t2,
1025     &                          t1, tg, tf, w, v,
1026     &                          nbf, nmo, nocc, nso, noso)
1027      implicit none
1028      integer rtdb, basis, nbf, nmo, nocc, nso, noso
1029c
1030c     nbf = number of basis functions
1031c     nmo = number of molecular orbitals
1032c     nocc = number of occupied orbitals
1033c     nso = number of spin orbitals (for now set to 2*nmo)
1034c     noso = number of occupied spin orbitals (for now set to 2*nocc)
1035c
1036      double precision g_mo(nmo,nmo,nmo,nmo)
1037      double precision ener(nmo)
1038      double precision t2(nso,nso,nso,nso)
1039      double precision t1(nso,nso)
1040      double precision tg(nso,nso,nso,noso,noso,noso)
1041      double precision tf(nso,nso,nso,noso,noso,noso)
1042      double precision w(nso,nso,nso,noso,noso,noso)
1043      double precision v(nso,nso,nso,noso,noso,noso)
1044      double precision g, et
1045c occupied orbitals
1046      integer i, j, k, l, m, n
1047      integer i_orb, j_orb, k_orb, l_orb, m_orb, n_orb
1048      integer i_spin, j_spin, k_spin, l_spin, m_spin, n_spin
1049c unoccupied orbitals
1050      integer a, b, c, d, e, f
1051      integer a_orb, b_orb, c_orb, d_orb, e_orb, f_orb
1052      integer a_spin, b_spin, c_spin, d_spin, e_spin, f_spin
1053c
1054      write(6,*)' nbf, nmo, nocc, nso, noso: ',
1055     &            nbf, nmo, nocc, nso, noso
1056c
1057      call dfill(noso**3*nso**3, 0.0d0, tg, 1)
1058      call dfill(noso**3*nso**3, 0.0d0, tf, 1)
1059      call dfill(noso**3*nso**3, 0.0d0, w, 1)
1060      call dfill(noso**3*nso**3, 0.0d0, v, 1)
1061c
1062      do a = noso+1, nso
1063         a_orb = (a-1)/2 + 1
1064         a_spin = 0
1065         if ((a-2*a_orb).ne.0)a_spin = 1
1066         do b = noso+1, nso
1067            b_orb = (b-1)/2 + 1
1068            b_spin = 0
1069            if ((b-2*b_orb).ne.0)b_spin = 1
1070            do c = noso+1, nso
1071               c_orb = (c-1)/2 + 1
1072               c_spin = 0
1073               if ((c-2*c_orb).ne.0)c_spin = 1
1074c
1075               do i = 1, noso
1076                  i_orb = (i-1)/2 + 1
1077                  i_spin = 0
1078                  if ((i-2*i_orb).ne.0)i_spin = 1
1079                  do j = 1, noso
1080                     j_orb = (j-1)/2 + 1
1081                     j_spin = 0
1082                     if ((j-2*j_orb).ne.0)j_spin = 1
1083                     do k = 1, noso
1084                        k_orb = (k-1)/2 + 1
1085                        k_spin = 0
1086                        if ((k-2*k_orb).ne.0)k_spin = 1
1087c
1088                        do e = noso+1, nso
1089                           e_orb = (e-1)/2 + 1
1090                           e_spin = 0
1091                           if ((e-2*e_orb).ne.0)e_spin = 1
1092c
1093                           g=0.0d0
1094                           if (a_spin.eq.e_spin.and.
1095     &                         b_spin.eq.k_spin)then
1096                              g = g_mo(a_orb,b_orb,e_orb,k_orb)
1097                           endif
1098                           if (a_spin.eq.k_spin.and.
1099     &                         b_spin.eq.e_spin)then
1100                              g = g - g_mo(a_orb,b_orb,k_orb,e_orb)
1101                           endif
1102c                           write(6,*)' g:',g
1103                           tg(a,b,c,i,j,k) = tg(a,b,c,i,j,k) +
1104     &                        g*t2(c,e,i,j)
1105                        enddo
1106c
1107                        do m = 1, noso
1108                           m_orb = (m-1)/2 + 1
1109                           m_spin = 0
1110                           if ((m-2*m_orb).ne.0)m_spin = 1
1111c
1112                           g=0.0d0
1113                           if (c_spin.eq.i_spin.and.
1114     &                         m_spin.eq.j_spin)then
1115                              g = g_mo(c_orb,m_orb,i_orb,j_orb)
1116                           endif
1117                           if (c_spin.eq.j_spin.and.
1118     &                         m_spin.eq.i_spin)then
1119                              g = g - g_mo(c_orb,m_orb,j_orb,i_orb)
1120                           endif
1121c                           write(6,*)' g:',g
1122                           tg(a,b,c,i,j,k) = tg(a,b,c,i,j,k) -
1123     &                        g*t2(a,b,m,k)
1124                        enddo
1125c
1126                        g=0.0d0
1127                        if (i_spin.eq.a_spin.and.j_spin.eq.b_spin)then
1128                           g = g_mo(i_orb,j_orb,a_orb,b_orb)
1129                        endif
1130                        if (i_spin.eq.b_spin.and.j_spin.eq.a_spin)then
1131                           g = g - g_mo(i_orb,j_orb,b_orb,a_orb)
1132                        endif
1133c                        write(6,*)' g:',g
1134                        tf(a,b,c,i,j,k) = tg(a,b,c,i,j,k) +
1135     &                     g*t1(c,k)
1136                     enddo
1137                  enddo
1138               enddo
1139            enddo
1140         enddo
1141      enddo
1142c
1143      do a = noso+1, nso
1144         a_orb = (a-1)/2 + 1
1145         a_spin = 0
1146         if ((a-2*a_orb).ne.0)a_spin = 1
1147         do b = noso+1, nso
1148            b_orb = (b-1)/2 + 1
1149            b_spin = 0
1150            if ((b-2*b_orb).ne.0)b_spin = 1
1151            do c = noso+1, nso
1152               c_orb = (c-1)/2 + 1
1153               c_spin = 0
1154               if ((c-2*c_orb).ne.0)c_spin = 1
1155c
1156               do i = 1, noso
1157                  i_orb = (i-1)/2 + 1
1158                  i_spin = 0
1159                  if ((i-2*i_orb).ne.0)i_spin = 1
1160                  do j = 1, noso
1161                     j_orb = (j-1)/2 + 1
1162                     j_spin = 0
1163                     if ((j-2*j_orb).ne.0)j_spin = 1
1164                     do k = 1, noso
1165                        k_orb = (k-1)/2 + 1
1166                        k_spin = 0
1167                        if ((k-2*k_orb).ne.0)k_spin = 1
1168c
1169                        w(a,b,c,i,j,k) = w(a,b,c,i,j,k) +
1170     &                              tg(a,b,c,i,j,k) -
1171     &                              tg(a,b,c,k,j,i) -
1172     &                              tg(a,b,c,i,k,j) -
1173     &                              tg(c,b,a,i,j,k) -
1174     &                              tg(a,c,b,i,j,k) +
1175     &                              tg(c,b,a,k,j,i) +
1176     &                              tg(a,c,b,k,j,i) +
1177     &                              tg(c,b,a,i,k,j) +
1178     &                              tg(a,c,b,i,k,j)
1179c
1180                        v(a,b,c,i,j,k) = v(a,b,c,i,j,k) +
1181     &                              tf(a,b,c,i,j,k) -
1182     &                              tf(a,b,c,k,j,i) -
1183     &                              tf(a,b,c,i,k,j) -
1184     &                              tf(c,b,a,i,j,k) -
1185     &                              tf(a,c,b,i,j,k) +
1186     &                              tf(c,b,a,k,j,i) +
1187     &                              tf(a,c,b,k,j,i) +
1188     &                              tf(c,b,a,i,k,j) +
1189     &                              tf(a,c,b,i,k,j)
1190c
1191                     enddo
1192                  enddo
1193               enddo
1194            enddo
1195         enddo
1196      enddo
1197c
1198      Et = 0.0d0
1199c
1200      do a = noso+1, nso
1201         a_orb = (a-1)/2 + 1
1202         a_spin = 0
1203         if ((a-2*a_orb).ne.0)a_spin = 1
1204         do b = noso+1, nso
1205            b_orb = (b-1)/2 + 1
1206            b_spin = 0
1207            if ((b-2*b_orb).ne.0)b_spin = 1
1208            do c = noso+1, nso
1209               c_orb = (c-1)/2 + 1
1210               c_spin = 0
1211               if ((c-2*c_orb).ne.0)c_spin = 1
1212c
1213               do i = 1, noso
1214                  i_orb = (i-1)/2 + 1
1215                  i_spin = 0
1216                  if ((i-2*i_orb).ne.0)i_spin = 1
1217                  do j = 1, noso
1218                     j_orb = (j-1)/2 + 1
1219                     j_spin = 0
1220                     if ((j-2*j_orb).ne.0)j_spin = 1
1221                     do k = 1, noso
1222                        k_orb = (k-1)/2 + 1
1223                        k_spin = 0
1224                        if ((k-2*k_orb).ne.0)k_spin = 1
1225c
1226                        Et = Et - 1.0d0/36.0d0*w(a,b,c,i,j,k)*
1227     &                                         v(a,b,c,i,j,k)/
1228     &                        ( ener(a_orb) +
1229     &                          ener(b_orb) +
1230     &                          ener(c_orb) -
1231     &                          ener(i_orb) -
1232     &                          ener(j_orb) -
1233     &                          ener(k_orb) )
1234
1235
1236                     enddo
1237                  enddo
1238               enddo
1239            enddo
1240         enddo
1241      enddo
1242c
1243      write(6,*)' (t) energy = ',et
1244c
1245c      do a = noso+1, nso
1246c         a_orb = (a-1)/2 + 1
1247c         a_spin = 0
1248c         if ((a-2*a_orb).ne.0)a_spin = 1
1249c         do b = noso+1, nso
1250c            b_orb = (b-1)/2 + 1
1251c            b_spin = 0
1252c            if ((b-2*b_orb).ne.0)b_spin = 1
1253c            do c = noso+1, nso
1254c               c_orb = (c-1)/2 + 1
1255c               c_spin = 0
1256c               if ((c-2*c_orb).ne.0)c_spin = 1
1257cc
1258c               do i = 1, noso
1259c                  i_orb = (i-1)/2 + 1
1260c                  i_spin = 0
1261c                  if ((i-2*i_orb).ne.0)i_spin = 1
1262c                  do j = 1, noso
1263c                     j_orb = (j-1)/2 + 1
1264c                     j_spin = 0
1265c                     if ((j-2*j_orb).ne.0)j_spin = 1
1266c                     do k = 1, noso
1267c                        k_orb = (k-1)/2 + 1
1268c                        k_spin = 0
1269c                        if ((k-2*k_orb).ne.0)k_spin = 1
1270cc
1271c                        if(dabs(w(a,b,c,i,j,k)).gt.1d-10)then
1272c                          write(6,*)' i, j, k, a, b, c, w, v: ',
1273c     &                                i, j, k, a, b, c,
1274c     &                                w(a,b,c,i,j,k),
1275c     &                                v(a,b,c,i,j,k)
1276c                        endif
1277cc
1278c                     enddo
1279c                  enddo
1280c               enddo
1281c            enddo
1282c         enddo
1283c      enddo
1284c
1285      return
1286      end
1287      subroutine t2_init(g_mo, e, t, nmo, nso, noso)
1288      implicit none
1289      integer nmo, nso, noso
1290      double precision g_mo(nmo,nmo,nmo,nmo)
1291      double precision e(nmo)
1292      double precision t(nso,nso,nso,nso)
1293      double precision g, denom
1294c
1295c     cc stuff (greek indices denote occupied spin-orbitals)
1296c     cc stuff (r, s, t, u, ... denote unoccupied spin-orbitals)
1297c
1298c occupied orbitals
1299      integer alpha, beta
1300      integer alpha_orb, beta_orb
1301      integer alpha_spin, beta_spin
1302c unoccupied orbitals
1303      integer m, n
1304      integer m_orb, n_orb
1305      integer m_spin, n_spin
1306c
1307c     First generate an initial guess for t2
1308c
1309      call dfill(nso**4, 0.0d0, t, 1)
1310      do m = noso+1, nso
1311c get spatial orbital and whether alpha or beta spin
1312         m_orb = (m-1)/2 + 1
1313         m_spin = 0
1314         if ((m-2*m_orb).ne.0)m_spin = 1
1315c         write(6,*)' m, m_orb, m_spin: ',
1316c     &               m, m_orb, m_spin
1317         do n = noso+1, nso
1318            n_orb = (n-1)/2 + 1
1319            n_spin = 0
1320            if ((n-2*n_orb).ne.0)n_spin = 1
1321c            write(6,*)' n, n_orb, n_spin: ',
1322c     &                  n, n_orb, n_spin
1323            do alpha = 1, noso
1324               alpha_orb = (alpha-1)/2 + 1
1325               alpha_spin = 0
1326               if ((alpha-2*alpha_orb).ne.0)alpha_spin = 1
1327c               write(6,*)' alpha, alpha_orb, alpha_spin: ',
1328c     &                     alpha, alpha_orb, alpha_spin
1329               do beta = 1, noso
1330                  beta_orb = (beta-1)/2 + 1
1331                  beta_spin = 0
1332                  if ((beta-2*beta_orb).ne.0)beta_spin = 1
1333c                  write(6,*)' beta, beta_orb, beta_spin: ',
1334c     &                        beta, beta_orb, beta_spin
1335c
1336                  denom = -e(m_orb)-e(n_orb)+e(alpha_orb)+e(beta_orb)
1337c                  write(6,*)' denom: ', denom
1338                  g=0.0d0
1339                  if (m_spin.eq.alpha_spin.and.n_spin.eq.beta_spin)then
1340                     g = g_mo(m_orb,n_orb,alpha_orb,beta_orb)
1341                  endif
1342                  if (m_spin.eq.beta_spin.and.n_spin.eq.alpha_spin)then
1343                     g = g - g_mo(m_orb,n_orb,beta_orb,alpha_orb)
1344                  endif
1345c                  write(6,*)' g:',g
1346                  t(m,n,alpha,beta) = g/denom
1347c                  write(6,*)' m,n,alpha,beta,t(m,n,alpha,beta):',
1348c     &                        m,n,alpha,beta,t(m,n,alpha,beta)
1349               enddo
1350            enddo
1351         enddo
1352      enddo
1353      return
1354      end
1355      subroutine t2_l_gb(g_mo, e, t, t2, nmo, nso, noso)
1356      implicit none
1357      integer nmo, nso, noso
1358      double precision g_mo(nmo,nmo,nmo,nmo)
1359      double precision e(nmo)
1360      double precision t(nso,nso,nso,nso)
1361      double precision t2(nso,nso,nso,nso)
1362      double precision g, denom
1363c
1364c     cc stuff (greek indices denote occupied spin-orbitals)
1365c     cc stuff (r, s, t, u, ... denote unoccupied spin-orbitals)
1366c
1367c occupied orbitals
1368      integer alpha, beta, delta, gamma
1369      integer alpha_orb, beta_orb, delta_orb, gamma_orb
1370      integer alpha_spin, beta_spin, delta_spin, gamma_spin
1371c unoccupied orbitals
1372      integer m, n, p, q
1373      integer m_orb, n_orb, p_orb, q_orb
1374      integer m_spin, n_spin, p_spin, q_spin
1375c
1376c     put t into t2 expression keeping all terms linear in t2
1377c
1378      do m = noso+1, nso
1379c get spatial orbital and whether alpha or beta spin
1380         m_orb = (m-1)/2 + 1
1381         m_spin = 0
1382         if ((m-2*m_orb).ne.0)m_spin = 1
1383         do n = noso+1, nso
1384            n_orb = (n-1)/2 + 1
1385            n_spin = 0
1386            if ((n-2*n_orb).ne.0)n_spin = 1
1387            do alpha = 1, noso
1388               alpha_orb = (alpha-1)/2 + 1
1389               alpha_spin = 0
1390               if ((alpha-2*alpha_orb).ne.0)alpha_spin = 1
1391               do beta = 1, noso
1392                  beta_orb = (beta-1)/2 + 1
1393                  beta_spin = 0
1394                  if ((beta-2*beta_orb).ne.0)beta_spin = 1
1395                  denom = e(m_orb)+e(n_orb)-e(alpha_orb)-e(beta_orb)
1396c                  write(6,*)' denom: ', denom
1397                  g=0.0d0
1398                  if (m_spin.eq.alpha_spin.and.n_spin.eq.beta_spin)then
1399                     g = g_mo(m_orb,n_orb,alpha_orb,beta_orb)
1400                  endif
1401                  if (m_spin.eq.beta_spin.and.n_spin.eq.alpha_spin)then
1402                     g = g - g_mo(m_orb,n_orb,beta_orb,alpha_orb)
1403                  endif
1404c                  write(6,*)' g:',g
1405                  t2(m,n,alpha,beta) = g
1406c
1407                  do p = noso+1, nso
1408                     p_orb = (p-1)/2 + 1
1409                     p_spin = 0
1410                     if ((p-2*p_orb).ne.0)p_spin = 1
1411                     do q = noso+1, p-1
1412                        q_orb = (q-1)/2 + 1
1413                        q_spin = 0
1414                        if ((q-2*q_orb).ne.0)q_spin = 1
1415c
1416                        g=0.0d0
1417                        if (m_spin.eq.p_spin.and.n_spin.eq.q_spin)then
1418                           g = g_mo(m_orb,n_orb,p_orb,q_orb)
1419                        endif
1420                        if (m_spin.eq.q_spin.and.n_spin.eq.p_spin)then
1421                           g = g - g_mo(m_orb,n_orb,q_orb,p_orb)
1422                        endif
1423c                        write(6,*)' g:',g
1424                        t2(m,n,alpha,beta) = t2(m,n,alpha,beta) -
1425     &                                       g*t(p,q,alpha,beta)
1426c
1427                     enddo
1428                  enddo
1429c
1430                  do gamma = 1, noso
1431                     gamma_orb = (gamma-1)/2 + 1
1432                     gamma_spin = 0
1433                     if ((gamma-2*gamma_orb).ne.0)gamma_spin = 1
1434                     do delta = 1, gamma-1
1435                        delta_orb = (delta-1)/2 + 1
1436                        delta_spin = 0
1437                        if ((delta-2*delta_orb).ne.0)delta_spin = 1
1438c
1439                        g=0.0d0
1440                        if (gamma_spin.eq.alpha_spin.and.
1441     &                      delta_spin.eq.beta_spin)then
1442                           g =
1443     &                     g_mo(gamma_orb,delta_orb,alpha_orb,beta_orb)
1444                        endif
1445                        if (gamma_spin.eq.beta_spin.and.
1446     &                      delta_spin.eq.alpha_spin)then
1447                           g = g -
1448     &                     g_mo(gamma_orb,delta_orb,beta_orb,alpha_orb)
1449                        endif
1450c                        write(6,*)' g:',g
1451                        t2(m,n,alpha,beta) = t2(m,n,alpha,beta) -
1452     &                                       g*t(m,n,gamma,delta)
1453c
1454                     enddo
1455                  enddo
1456c
1457                  do p = noso+1, nso
1458                     p_orb = (p-1)/2 + 1
1459                     p_spin = 0
1460                     if ((p-2*p_orb).ne.0)p_spin = 1
1461                     do gamma = 1, noso
1462                        gamma_orb = (gamma-1)/2 + 1
1463                        gamma_spin = 0
1464                        if ((gamma-2*gamma_orb).ne.0)gamma_spin = 1
1465c
1466                        g=0.0d0
1467                        if (gamma_spin.eq.beta_spin.and.
1468     &                      n_spin.eq.p_spin)then
1469                           g = g_mo(gamma_orb,n_orb,beta_orb,p_orb)
1470                        endif
1471                        if (gamma_spin.eq.p_spin.and.
1472     &                      n_spin.eq.beta_spin)then
1473                           g = g - g_mo(gamma_orb,n_orb,p_orb,beta_orb)
1474                        endif
1475c                        write(6,*)' g:',g
1476                        t2(m,n,alpha,beta) = t2(m,n,alpha,beta) +
1477     &                                       g*t(m,p,alpha,gamma)
1478c
1479                        g=0.0d0
1480                        if (gamma_spin.eq.beta_spin.and.
1481     &                      m_spin.eq.p_spin)then
1482                           g = g_mo(gamma_orb,m_orb,beta_orb,p_orb)
1483                        endif
1484                        if (gamma_spin.eq.p_spin.and.
1485     &                      m_spin.eq.beta_spin)then
1486                           g = g - g_mo(gamma_orb,m_orb,p_orb,beta_orb)
1487                        endif
1488c                        write(6,*)' g:',g
1489                        t2(m,n,alpha,beta) = t2(m,n,alpha,beta) -
1490     &                                       g*t(n,p,alpha,gamma)
1491c
1492                        g=0.0d0
1493                        if (gamma_spin.eq.alpha_spin.and.
1494     &                      n_spin.eq.p_spin)then
1495                           g = g_mo(gamma_orb,n_orb,alpha_orb,p_orb)
1496                        endif
1497                        if (gamma_spin.eq.p_spin.and.
1498     &                      n_spin.eq.alpha_spin)then
1499                           g = g - g_mo(gamma_orb,n_orb,p_orb,alpha_orb)
1500                        endif
1501c                        write(6,*)' g:',g
1502                        t2(m,n,alpha,beta) = t2(m,n,alpha,beta) -
1503     &                                       g*t(m,p,beta,gamma)
1504c
1505                        g=0.0d0
1506                        if (gamma_spin.eq.alpha_spin.and.
1507     &                      m_spin.eq.p_spin)then
1508                           g = g_mo(gamma_orb,m_orb,alpha_orb,p_orb)
1509                        endif
1510                        if (gamma_spin.eq.p_spin.and.
1511     &                      m_spin.eq.alpha_spin)then
1512                           g = g - g_mo(gamma_orb,m_orb,p_orb,alpha_orb)
1513                        endif
1514c                        write(6,*)' g:',g
1515                        t2(m,n,alpha,beta) = t2(m,n,alpha,beta) +
1516     &                                       g*t(n,p,beta,gamma)
1517                     enddo
1518                  enddo
1519                  t2(m,n,alpha,beta) = t2(m,n,alpha,beta)/denom
1520c                  write(6,*)' t2(m,n,alpha,beta):',
1521c     &                        t2(m,n,alpha,beta)
1522               enddo
1523            enddo
1524         enddo
1525      enddo
1526c      write(6,*)' T2 with linear terms. '
1527c      call writet2(t2,nso)
1528c      call correlation(g_mo,t2,nmo,nso,noso)
1529      return
1530      end
1531      subroutine t2_l_q_gb(g_mo, e, t, t2, nmo, nso, noso)
1532      implicit none
1533      integer nmo, nso, noso
1534      double precision g_mo(nmo,nmo,nmo,nmo)
1535      double precision e(nmo)
1536      double precision t(nso,nso,nso,nso)
1537      double precision t2(nso,nso,nso,nso)
1538      double precision g, denom
1539c
1540c     cc stuff (greek indices denote occupied spin-orbitals)
1541c     cc stuff (r, s, t, u, ... denote unoccupied spin-orbitals)
1542c
1543c occupied orbitals
1544      integer alpha, beta, delta, gamma
1545      integer alpha_orb, beta_orb, delta_orb, gamma_orb
1546      integer alpha_spin, beta_spin, delta_spin, gamma_spin
1547c unoccupied orbitals
1548      integer m, n, p, q
1549      integer m_orb, n_orb, p_orb, q_orb
1550      integer m_spin, n_spin, p_spin, q_spin
1551c
1552c     put t into t2 expression keeping all terms (linear and quadratic)
1553c
1554      do m = noso+1, nso
1555c get spatial orbital and whether alpha or beta spin
1556         m_orb = (m-1)/2 + 1
1557         m_spin = 0
1558         if ((m-2*m_orb).ne.0)m_spin = 1
1559         do n = noso+1, nso
1560            n_orb = (n-1)/2 + 1
1561            n_spin = 0
1562            if ((n-2*n_orb).ne.0)n_spin = 1
1563            do alpha = 1, noso
1564               alpha_orb = (alpha-1)/2 + 1
1565               alpha_spin = 0
1566               if ((alpha-2*alpha_orb).ne.0)alpha_spin = 1
1567               do beta = 1, noso
1568                  beta_orb = (beta-1)/2 + 1
1569                  beta_spin = 0
1570                  if ((beta-2*beta_orb).ne.0)beta_spin = 1
1571                  denom = e(m_orb)+e(n_orb)-e(alpha_orb)-e(beta_orb)
1572c                  write(6,*)' denom: ', denom
1573                  g=0.0d0
1574                  if (m_spin.eq.alpha_spin.and.n_spin.eq.beta_spin)then
1575                     g = g_mo(m_orb,n_orb,alpha_orb,beta_orb)
1576                  endif
1577                  if (m_spin.eq.beta_spin.and.n_spin.eq.alpha_spin)then
1578                     g = g - g_mo(m_orb,n_orb,beta_orb,alpha_orb)
1579                  endif
1580c                  write(6,*)' g:',g
1581                  t2(m,n,alpha,beta) = g
1582c
1583                  do p = noso+1, nso
1584                     p_orb = (p-1)/2 + 1
1585                     p_spin = 0
1586                     if ((p-2*p_orb).ne.0)p_spin = 1
1587                     do q = noso+1, p-1
1588                        q_orb = (q-1)/2 + 1
1589                        q_spin = 0
1590                        if ((q-2*q_orb).ne.0)q_spin = 1
1591c
1592                        g=0.0d0
1593                        if (m_spin.eq.p_spin.and.n_spin.eq.q_spin)then
1594                           g = g_mo(m_orb,n_orb,p_orb,q_orb)
1595                        endif
1596                        if (m_spin.eq.q_spin.and.n_spin.eq.p_spin)then
1597                           g = g - g_mo(m_orb,n_orb,q_orb,p_orb)
1598                        endif
1599c                        write(6,*)' g:',g
1600                        t2(m,n,alpha,beta) = t2(m,n,alpha,beta) -
1601     &                                       g*t(p,q,alpha,beta)
1602c
1603                     enddo
1604                  enddo
1605c
1606                  do gamma = 1, noso
1607                     gamma_orb = (gamma-1)/2 + 1
1608                     gamma_spin = 0
1609                     if ((gamma-2*gamma_orb).ne.0)gamma_spin = 1
1610                     do delta = 1, gamma-1
1611                        delta_orb = (delta-1)/2 + 1
1612                        delta_spin = 0
1613                        if ((delta-2*delta_orb).ne.0)delta_spin = 1
1614c
1615                        g=0.0d0
1616                        if (gamma_spin.eq.alpha_spin.and.
1617     &                      delta_spin.eq.beta_spin)then
1618                           g =
1619     &                     g_mo(gamma_orb,delta_orb,alpha_orb,beta_orb)
1620                        endif
1621                        if (gamma_spin.eq.beta_spin.and.
1622     &                      delta_spin.eq.alpha_spin)then
1623                           g = g -
1624     &                     g_mo(gamma_orb,delta_orb,beta_orb,alpha_orb)
1625                        endif
1626c                        write(6,*)' g:',g
1627                        t2(m,n,alpha,beta) = t2(m,n,alpha,beta) -
1628     &                                       g*t(m,n,gamma,delta)
1629c
1630                     enddo
1631                  enddo
1632c
1633                  do p = noso+1, nso
1634                     p_orb = (p-1)/2 + 1
1635                     p_spin = 0
1636                     if ((p-2*p_orb).ne.0)p_spin = 1
1637                     do gamma = 1, noso
1638                        gamma_orb = (gamma-1)/2 + 1
1639                        gamma_spin = 0
1640                        if ((gamma-2*gamma_orb).ne.0)gamma_spin = 1
1641c
1642                        g=0.0d0
1643                        if (gamma_spin.eq.beta_spin.and.
1644     &                      n_spin.eq.p_spin)then
1645                           g = g_mo(gamma_orb,n_orb,beta_orb,p_orb)
1646                        endif
1647                        if (gamma_spin.eq.p_spin.and.
1648     &                      n_spin.eq.beta_spin)then
1649                           g = g - g_mo(gamma_orb,n_orb,p_orb,beta_orb)
1650                        endif
1651c                        write(6,*)' g:',g
1652                        t2(m,n,alpha,beta) = t2(m,n,alpha,beta) +
1653     &                                       g*t(m,p,alpha,gamma)
1654c
1655                        g=0.0d0
1656                        if (gamma_spin.eq.beta_spin.and.
1657     &                      m_spin.eq.p_spin)then
1658                           g = g_mo(gamma_orb,m_orb,beta_orb,p_orb)
1659                        endif
1660                        if (gamma_spin.eq.p_spin.and.
1661     &                      m_spin.eq.beta_spin)then
1662                           g = g - g_mo(gamma_orb,m_orb,p_orb,beta_orb)
1663                        endif
1664c                        write(6,*)' g:',g
1665                        t2(m,n,alpha,beta) = t2(m,n,alpha,beta) -
1666     &                                       g*t(n,p,alpha,gamma)
1667c
1668                        g=0.0d0
1669                        if (gamma_spin.eq.alpha_spin.and.
1670     &                      n_spin.eq.p_spin)then
1671                           g = g_mo(gamma_orb,n_orb,alpha_orb,p_orb)
1672                        endif
1673                        if (gamma_spin.eq.p_spin.and.
1674     &                      n_spin.eq.alpha_spin)then
1675                           g = g - g_mo(gamma_orb,n_orb,p_orb,alpha_orb)
1676                        endif
1677c                        write(6,*)' g:',g
1678                        t2(m,n,alpha,beta) = t2(m,n,alpha,beta) -
1679     &                                       g*t(m,p,beta,gamma)
1680c
1681                        g=0.0d0
1682                        if (gamma_spin.eq.alpha_spin.and.
1683     &                      m_spin.eq.p_spin)then
1684                           g = g_mo(gamma_orb,m_orb,alpha_orb,p_orb)
1685                        endif
1686                        if (gamma_spin.eq.p_spin.and.
1687     &                      m_spin.eq.alpha_spin)then
1688                           g = g - g_mo(gamma_orb,m_orb,p_orb,alpha_orb)
1689                        endif
1690c                        write(6,*)' g:',g
1691                        t2(m,n,alpha,beta) = t2(m,n,alpha,beta) +
1692     &                                       g*t(n,p,beta,gamma)
1693                     enddo
1694                  enddo
1695c
1696                  do gamma = 1, noso
1697                     gamma_orb = (gamma-1)/2 + 1
1698                     gamma_spin = 0
1699                     if ((gamma-2*gamma_orb).ne.0)gamma_spin = 1
1700                     do delta = 1, gamma-1
1701                        delta_orb = (delta-1)/2 + 1
1702                        delta_spin = 0
1703                        if ((delta-2*delta_orb).ne.0)delta_spin = 1
1704                        do p = noso+1, nso
1705                           p_orb = (p-1)/2 + 1
1706                           p_spin = 0
1707                           if ((p-2*p_orb).ne.0)p_spin = 1
1708                           do q = noso+1, p-1
1709                              q_orb = (q-1)/2 + 1
1710                              q_spin = 0
1711                              if ((q-2*q_orb).ne.0)q_spin = 1
1712c
1713                              g=0.0d0
1714                              if (gamma_spin.eq.p_spin.and.
1715     &                           delta_spin.eq.q_spin)then
1716                              g = g_mo(gamma_orb,delta_orb,p_orb,q_orb)
1717                              endif
1718                              if (gamma_spin.eq.q_spin.and.
1719     &                            delta_spin.eq.p_spin)then
1720                          g = g - g_mo(gamma_orb,delta_orb,q_orb,p_orb)
1721                              endif
1722c                              write(6,*)' g:',g
1723c
1724                              t2(m,n,alpha,beta) = t2(m,n,alpha,beta) +
1725     &                      g*(t(p,q,alpha,beta)*t(m,n,gamma,delta) -
1726     &                  2.0d0*(t(m,p,alpha,beta)*t(n,q,gamma,delta) +
1727     &                         t(n,q,alpha,beta)*t(m,p,gamma,delta))-
1728     &                  2.0d0*(t(m,n,alpha,gamma)*t(p,q,beta,delta) +
1729     &                         t(p,q,alpha,gamma)*t(m,n,beta,delta))+
1730     &                  4.0d0*(t(m,p,alpha,gamma)*t(n,q,beta,delta) +
1731     &                         t(n,q,alpha,gamma)*t(m,p,beta,delta)))
1732                           enddo
1733                        enddo
1734                     enddo
1735                  enddo
1736                  t2(m,n,alpha,beta) = t2(m,n,alpha,beta)/denom
1737c                  write(6,*)' t2(m,n,alpha,beta):',
1738c     &                        t2(m,n,alpha,beta)
1739               enddo
1740            enddo
1741         enddo
1742      enddo
1743c      write(6,*)' T2 with all terms. '
1744c      call writet2(t2,nso)
1745c      call correlation(g_mo,t2,nmo,nso,noso)
1746      return
1747      end
1748      subroutine writet2(t,len)
1749      implicit none
1750      integer len, i, j, k, l
1751      double precision t(len,len,len,len)
1752      do i = 1, len
1753         do j = 1, len
1754            do k = 1, len
1755               do l = 1, len
1756                  if (abs(t(i,j,k,l)).gt.1d-6)
1757     &               write(6,*)i,j,k,l,t(i,j,k,l)
1758               enddo
1759            enddo
1760         enddo
1761      enddo
1762      return
1763      end
1764      subroutine writet1(t,len)
1765      implicit none
1766      integer len, i, j
1767      double precision t(len,len)
1768      do i = 1, len
1769         do j = 1, len
1770            if (abs(t(i,j)).gt.1d-6)
1771     &         write(6,*)i,j,t(i,j)
1772         enddo
1773      enddo
1774      return
1775      end
1776      subroutine correlation(g_mo,t1,t2,nmo,nso,noso)
1777      implicit none
1778      integer nmo, nso, noso
1779      double precision g_mo(nmo,nmo,nmo,nmo)
1780      double precision t2(nso,nso,nso,nso)
1781      double precision t1(nso,nso)
1782      double precision g, e2
1783c occupied orbitals
1784      integer alpha, beta
1785      integer alpha_orb, beta_orb
1786      integer alpha_spin, beta_spin
1787c unoccupied orbitals
1788      integer m, n
1789      integer m_orb, n_orb
1790      integer m_spin, n_spin
1791c
1792      e2 = 0.0d0
1793      do m = noso+1, nso
1794c get spatial orbital and whether alpha or beta spin
1795         m_orb = (m-1)/2 + 1
1796         m_spin = 0
1797         if ((m-2*m_orb).ne.0)m_spin = 1
1798         do n = noso+1, nso
1799            n_orb = (n-1)/2 + 1
1800            n_spin = 0
1801            if ((n-2*n_orb).ne.0)n_spin = 1
1802            do alpha = 1, noso
1803               alpha_orb = (alpha-1)/2 + 1
1804               alpha_spin = 0
1805               if ((alpha-2*alpha_orb).ne.0)alpha_spin = 1
1806               do beta = 1, noso
1807                  beta_orb = (beta-1)/2 + 1
1808                  beta_spin = 0
1809                  if ((beta-2*beta_orb).ne.0)beta_spin = 1
1810c
1811                  g=0.0d0
1812                  if (m_spin.eq.beta_spin.and.n_spin.eq.alpha_spin)then
1813                     g = g_mo(beta_orb,alpha_orb,m_orb,n_orb)
1814                  endif
1815                  if (m_spin.eq.alpha_spin.and.n_spin.eq.beta_spin)then
1816                     g = g - g_mo(beta_orb,alpha_orb,n_orb,m_orb)
1817                  endif
1818c                 write(6,*)' g:',g
1819                  e2 = e2 - 0.25d0*g*(t2(m,n,alpha,beta)+
1820     &                                t1(m,alpha)*t1(n,beta) -
1821     &                                t1(n,alpha)*t1(m,beta))
1822c
1823               enddo
1824            enddo
1825         enddo
1826      enddo
1827      write(6,*)' e2 = ',e2
1828c
1829      return
1830      end
1831      subroutine t2_l_rjh(g_mo, t, t2, fock, nmo, nso, noso)
1832      implicit none
1833      integer nmo, nso, noso
1834      double precision g_mo(nmo,nmo,nmo,nmo)
1835      double precision fock(nmo,nmo)
1836      double precision t(nso,nso,nso,nso)
1837      double precision t2(nso,nso,nso,nso)
1838      double precision g, f_t, denom
1839c occupied orbitals
1840      integer i, j, k, l, m, n
1841      integer i_orb, j_orb, k_orb, l_orb, m_orb, n_orb
1842      integer i_spin, j_spin, k_spin, l_spin, m_spin, n_spin
1843c unoccupied orbitals
1844      integer a, b, c, d, e, f
1845      integer a_orb, b_orb, c_orb, d_orb, e_orb, f_orb
1846      integer a_spin, b_spin, c_spin, d_spin, e_spin, f_spin
1847c
1848c     put t into t2 expression keeping all terms linear in t2
1849c
1850      do e = noso+1, nso
1851c get spatial orbital and whether alpha or beta spin
1852         e_orb = (e-1)/2 + 1
1853         e_spin = 0
1854         if ((e-2*e_orb).ne.0)e_spin = 1
1855         do f = noso+1, nso
1856            f_orb = (f-1)/2 + 1
1857            f_spin = 0
1858            if ((f-2*f_orb).ne.0)f_spin = 1
1859            do m = 1, noso
1860               m_orb = (m-1)/2 + 1
1861               m_spin = 0
1862               if ((m-2*m_orb).ne.0)m_spin = 1
1863               do n = 1, noso
1864                  n_orb = (n-1)/2 + 1
1865                  n_spin = 0
1866                  if ((n-2*n_orb).ne.0)n_spin = 1
1867c
1868                  g=0.0d0
1869                  if (e_spin.eq.m_spin.and.f_spin.eq.n_spin)then
1870                     g = g_mo(e_orb,f_orb,m_orb,n_orb)
1871                  endif
1872                  if (e_spin.eq.n_spin.and.f_spin.eq.m_spin)then
1873                     g = g - g_mo(e_orb,f_orb,n_orb,m_orb)
1874                  endif
1875c                  write(6,*)' g:',g
1876                  t2(e,f,m,n) = g
1877c
1878                  do a = noso+1, nso
1879                     a_orb = (a-1)/2 + 1
1880                     a_spin = 0
1881                     if ((a-2*a_orb).ne.0)a_spin = 1
1882                     f_t = 0.0d0
1883                     if (f_spin.eq.a_spin)then
1884                        f_t = fock(f_orb,a_orb)
1885                     endif
1886                     t2(e,f,m,n) = t2(e,f,m,n) +
1887     &                          f_t*t(e,a,m,n)
1888                     f_t = 0.0d0
1889                     if (e_spin.eq.a_spin)then
1890                        f_t = fock(e_orb,a_orb)
1891                     endif
1892                     t2(e,f,m,n) = t2(e,f,m,n) -
1893     &                          f_t*t(f,a,m,n)
1894                  enddo
1895c
1896                  do i = 1, noso
1897                     i_orb = (i-1)/2 + 1
1898                     i_spin = 0
1899                     if ((i-2*i_orb).ne.0)i_spin = 1
1900                     f_t = 0.0d0
1901                     if (i_spin.eq.n_spin)then
1902                        f_t = fock(i_orb,n_orb)
1903                     endif
1904                     t2(e,f,m,n) = t2(e,f,m,n) -
1905     &                          f_t*t(e,f,m,i)
1906                     f_t = 0.0d0
1907                     if (i_spin.eq.m_spin)then
1908                        f_t = fock(i_orb,m_orb)
1909                     endif
1910                     t2(e,f,m,n) = t2(e,f,m,n) +
1911     &                          f_t*t(e,f,n,i)
1912                  enddo
1913c
1914                  do a = noso+1, nso
1915                     a_orb = (a-1)/2 + 1
1916                     a_spin = 0
1917                     if ((a-2*a_orb).ne.0)a_spin = 1
1918                     do i = 1, noso
1919                        i_orb = (i-1)/2 + 1
1920                        i_spin = 0
1921                        if ((i-2*i_orb).ne.0)i_spin = 1
1922                        g=0.0d0
1923                        if (i_spin.eq.a_spin.and.f_spin.eq.n_spin)then
1924                           g = g_mo(i_orb,f_orb,a_orb,n_orb)
1925                        endif
1926                        if (i_spin.eq.n_spin.and.a_spin.eq.f_spin)then
1927                           g = g - g_mo(i_orb,f_orb,n_orb,a_orb)
1928                        endif
1929c                        write(6,*)' g:',g
1930                        t2(e,f,m,n) = t2(e,f,m,n) +
1931     &                             g*t(e,a,m,i)
1932c ef flip
1933                        g=0.0d0
1934                        if (i_spin.eq.a_spin.and.e_spin.eq.n_spin)then
1935                           g = g_mo(i_orb,e_orb,a_orb,n_orb)
1936                        endif
1937                        if (i_spin.eq.n_spin.and.a_spin.eq.e_spin)then
1938                           g = g - g_mo(i_orb,e_orb,n_orb,a_orb)
1939                        endif
1940c                        write(6,*)' g:',g
1941                        t2(e,f,m,n) = t2(e,f,m,n) -
1942     &                             g*t(f,a,m,i)
1943c mn flip
1944                        g=0.0d0
1945                        if (i_spin.eq.a_spin.and.f_spin.eq.m_spin)then
1946                           g = g_mo(i_orb,f_orb,a_orb,m_orb)
1947                        endif
1948                        if (i_spin.eq.m_spin.and.a_spin.eq.f_spin)then
1949                           g = g - g_mo(i_orb,f_orb,m_orb,a_orb)
1950                        endif
1951c                        write(6,*)' g:',g
1952                        t2(e,f,m,n) = t2(e,f,m,n) -
1953     &                             g*t(e,a,n,i)
1954c ef+mn flip
1955                        g=0.0d0
1956                        if (i_spin.eq.a_spin.and.e_spin.eq.m_spin)then
1957                           g = g_mo(i_orb,e_orb,a_orb,m_orb)
1958                        endif
1959                        if (i_spin.eq.m_spin.and.a_spin.eq.e_spin)then
1960                           g = g - g_mo(i_orb,e_orb,m_orb,a_orb)
1961                        endif
1962c                        write(6,*)' g:',g
1963                        t2(e,f,m,n) = t2(e,f,m,n) +
1964     &                             g*t(f,a,n,i)
1965                     enddo
1966                  enddo
1967c
1968                  do i = 1, noso
1969                     i_orb = (i-1)/2 + 1
1970                     i_spin = 0
1971                     if ((i-2*i_orb).ne.0)i_spin = 1
1972                     do j = 1, noso
1973                        j_orb = (j-1)/2 + 1
1974                        j_spin = 0
1975                        if ((j-2*j_orb).ne.0)j_spin = 1
1976                        g=0.0d0
1977                        if (i_spin.eq.m_spin.and.j_spin.eq.n_spin)then
1978                           g = g_mo(i_orb,j_orb,m_orb,n_orb)
1979                        endif
1980                        if (i_spin.eq.n_spin.and.j_spin.eq.m_spin)then
1981                           g = g - g_mo(i_orb,j_orb,n_orb,m_orb)
1982                        endif
1983c                        write(6,*)' g:',g
1984                        t2(e,f,m,n) = t2(e,f,m,n) +
1985     &                             0.5d0*g*t(e,f,i,j)
1986                     enddo
1987                  enddo
1988                  do a = noso+1, nso
1989                     a_orb = (a-1)/2 + 1
1990                     a_spin = 0
1991                     if ((a-2*a_orb).ne.0)a_spin = 1
1992                     do b = noso+1, nso
1993                        b_orb = (b-1)/2 + 1
1994                        b_spin = 0
1995                        if ((b-2*b_orb).ne.0)b_spin = 1
1996                        g=0.0d0
1997                        if (e_spin.eq.a_spin.and.f_spin.eq.b_spin)then
1998                           g = g_mo(e_orb,f_orb,a_orb,b_orb)
1999                        endif
2000                        if (e_spin.eq.b_spin.and.f_spin.eq.a_spin)then
2001                           g = g - g_mo(e_orb,f_orb,b_orb,a_orb)
2002                        endif
2003c                        write(6,*)' g:',g
2004                        t2(e,f,m,n) = t2(e,f,m,n) +
2005     &                             0.5d0*g*t(a,b,m,n)
2006                     enddo
2007                  enddo
2008               enddo
2009            enddo
2010         enddo
2011      enddo
2012c      write(6,*)' t2(e,f,m,n) L residual:'
2013c      call writet2(t2,nso)
2014      return
2015      end
2016      subroutine t2_l_q_rjh(g_mo, t2, r2, fock, nmo, nso, noso)
2017      implicit none
2018      integer nmo, nso, noso
2019      double precision g_mo(nmo,nmo,nmo,nmo)
2020      double precision fock(nmo,nmo)
2021      double precision r2(nso,nso,nso,nso)
2022      double precision t2(nso,nso,nso,nso)
2023      double precision g, f_t, denom
2024c occupied orbitals
2025      integer i, j, k, l, m, n
2026      integer i_orb, j_orb, k_orb, l_orb, m_orb, n_orb
2027      integer i_spin, j_spin, k_spin, l_spin, m_spin, n_spin
2028c unoccupied orbitals
2029      integer a, b, c, d, e, f
2030      integer a_orb, b_orb, c_orb, d_orb, e_orb, f_orb
2031      integer a_spin, b_spin, c_spin, d_spin, e_spin, f_spin
2032c
2033c     put t2 into r2 expression keeping all terms quadratic in t2
2034c
2035      do e = noso+1, nso
2036c get spatial orbital and whether alpha or beta spin
2037         e_orb = (e-1)/2 + 1
2038         e_spin = 0
2039         if ((e-2*e_orb).ne.0)e_spin = 1
2040         do f = noso+1, nso
2041            f_orb = (f-1)/2 + 1
2042            f_spin = 0
2043            if ((f-2*f_orb).ne.0)f_spin = 1
2044            do m = 1, noso
2045               m_orb = (m-1)/2 + 1
2046               m_spin = 0
2047               if ((m-2*m_orb).ne.0)m_spin = 1
2048               do n = 1, noso
2049                  n_orb = (n-1)/2 + 1
2050                  n_spin = 0
2051                  if ((n-2*n_orb).ne.0)n_spin = 1
2052c
2053                  g=0.0d0
2054                  if (e_spin.eq.m_spin.and.f_spin.eq.n_spin)then
2055                     g = g_mo(e_orb,f_orb,m_orb,n_orb)
2056                  endif
2057                  if (e_spin.eq.n_spin.and.f_spin.eq.m_spin)then
2058                     g = g - g_mo(e_orb,f_orb,n_orb,m_orb)
2059                  endif
2060c                  write(6,*)' g:',g
2061                  r2(e,f,m,n) = g
2062c
2063                  do a = noso+1, nso
2064                     a_orb = (a-1)/2 + 1
2065                     a_spin = 0
2066                     if ((a-2*a_orb).ne.0)a_spin = 1
2067                     f_t = 0.0d0
2068                     if (f_spin.eq.a_spin)then
2069                        f_t = fock(f_orb,a_orb)
2070                     endif
2071                     r2(e,f,m,n) = r2(e,f,m,n) +
2072     &                          f_t*t2(e,a,m,n)
2073                     f_t = 0.0d0
2074                     if (e_spin.eq.a_spin)then
2075                        f_t = fock(e_orb,a_orb)
2076                     endif
2077                     r2(e,f,m,n) = r2(e,f,m,n) -
2078     &                          f_t*t2(f,a,m,n)
2079                  enddo
2080c
2081                  do i = 1, noso
2082                     i_orb = (i-1)/2 + 1
2083                     i_spin = 0
2084                     if ((i-2*i_orb).ne.0)i_spin = 1
2085                     f_t = 0.0d0
2086                     if (i_spin.eq.n_spin)then
2087                        f_t = fock(i_orb,n_orb)
2088                     endif
2089                     r2(e,f,m,n) = r2(e,f,m,n) -
2090     &                          f_t*t2(e,f,m,i)
2091                     f_t = 0.0d0
2092                     if (i_spin.eq.m_spin)then
2093                        f_t = fock(i_orb,m_orb)
2094                     endif
2095                     r2(e,f,m,n) = r2(e,f,m,n) +
2096     &                          f_t*t2(e,f,n,i)
2097                  enddo
2098c
2099                  do a = noso+1, nso
2100                     a_orb = (a-1)/2 + 1
2101                     a_spin = 0
2102                     if ((a-2*a_orb).ne.0)a_spin = 1
2103                     do i = 1, noso
2104                        i_orb = (i-1)/2 + 1
2105                        i_spin = 0
2106                        if ((i-2*i_orb).ne.0)i_spin = 1
2107                        g=0.0d0
2108                        if (i_spin.eq.a_spin.and.f_spin.eq.n_spin)then
2109                           g = g_mo(i_orb,f_orb,a_orb,n_orb)
2110                        endif
2111                        if (i_spin.eq.n_spin.and.a_spin.eq.f_spin)then
2112                           g = g - g_mo(i_orb,f_orb,n_orb,a_orb)
2113                        endif
2114c                        write(6,*)' g:',g
2115                        r2(e,f,m,n) = r2(e,f,m,n) +
2116     &                             g*t2(e,a,m,i)
2117c ef flip
2118                        g=0.0d0
2119                        if (i_spin.eq.a_spin.and.e_spin.eq.n_spin)then
2120                           g = g_mo(i_orb,e_orb,a_orb,n_orb)
2121                        endif
2122                        if (i_spin.eq.n_spin.and.a_spin.eq.e_spin)then
2123                           g = g - g_mo(i_orb,e_orb,n_orb,a_orb)
2124                        endif
2125c                        write(6,*)' g:',g
2126                        r2(e,f,m,n) = r2(e,f,m,n) -
2127     &                             g*t2(f,a,m,i)
2128c mn flip
2129                        g=0.0d0
2130                        if (i_spin.eq.a_spin.and.f_spin.eq.m_spin)then
2131                           g = g_mo(i_orb,f_orb,a_orb,m_orb)
2132                        endif
2133                        if (i_spin.eq.m_spin.and.a_spin.eq.f_spin)then
2134                           g = g - g_mo(i_orb,f_orb,m_orb,a_orb)
2135                        endif
2136c                        write(6,*)' g:',g
2137                        r2(e,f,m,n) = r2(e,f,m,n) -
2138     &                             g*t2(e,a,n,i)
2139c ef+mn flip
2140                        g=0.0d0
2141                        if (i_spin.eq.a_spin.and.e_spin.eq.m_spin)then
2142                           g = g_mo(i_orb,e_orb,a_orb,m_orb)
2143                        endif
2144                        if (i_spin.eq.m_spin.and.a_spin.eq.e_spin)then
2145                           g = g - g_mo(i_orb,e_orb,m_orb,a_orb)
2146                        endif
2147c                        write(6,*)' g:',g
2148                        r2(e,f,m,n) = r2(e,f,m,n) +
2149     &                             g*t2(f,a,n,i)
2150                     enddo
2151                  enddo
2152c
2153                  do i = 1, noso
2154                     i_orb = (i-1)/2 + 1
2155                     i_spin = 0
2156                     if ((i-2*i_orb).ne.0)i_spin = 1
2157                     do j = 1, noso
2158                        j_orb = (j-1)/2 + 1
2159                        j_spin = 0
2160                        if ((j-2*j_orb).ne.0)j_spin = 1
2161                        g=0.0d0
2162                        if (i_spin.eq.m_spin.and.j_spin.eq.n_spin)then
2163                           g = g_mo(i_orb,j_orb,m_orb,n_orb)
2164                        endif
2165                        if (i_spin.eq.n_spin.and.j_spin.eq.m_spin)then
2166                           g = g - g_mo(i_orb,j_orb,n_orb,m_orb)
2167                        endif
2168c                        write(6,*)' g:',g
2169                        r2(e,f,m,n) = r2(e,f,m,n) +
2170     &                             0.5d0*g*t2(e,f,i,j)
2171                     enddo
2172                  enddo
2173                  do a = noso+1, nso
2174                     a_orb = (a-1)/2 + 1
2175                     a_spin = 0
2176                     if ((a-2*a_orb).ne.0)a_spin = 1
2177                     do b = noso+1, nso
2178                        b_orb = (b-1)/2 + 1
2179                        b_spin = 0
2180                        if ((b-2*b_orb).ne.0)b_spin = 1
2181                        g=0.0d0
2182                        if (e_spin.eq.a_spin.and.f_spin.eq.b_spin)then
2183                           g = g_mo(e_orb,f_orb,a_orb,b_orb)
2184                        endif
2185                        if (e_spin.eq.b_spin.and.f_spin.eq.a_spin)then
2186                           g = g - g_mo(e_orb,f_orb,b_orb,a_orb)
2187                        endif
2188c                        write(6,*)' g:',g
2189                        r2(e,f,m,n) = r2(e,f,m,n) +
2190     &                             0.5d0*g*t2(a,b,m,n)
2191                     enddo
2192                  enddo
2193c
2194                  do i = 1, noso
2195                     i_orb = (i-1)/2 + 1
2196                     i_spin = 0
2197                     if ((i-2*i_orb).ne.0)i_spin = 1
2198                     do j = 1, noso
2199                        j_orb = (j-1)/2 + 1
2200                        j_spin = 0
2201                        if ((j-2*j_orb).ne.0)j_spin = 1
2202                        do a = noso+1, nso
2203                           a_orb = (a-1)/2 + 1
2204                           a_spin = 0
2205                           if ((a-2*a_orb).ne.0)a_spin = 1
2206                           do b = noso+1, nso
2207                              b_orb = (b-1)/2 + 1
2208                              b_spin = 0
2209                              if ((b-2*b_orb).ne.0)b_spin = 1
2210                              g=0.0d0
2211                              if (i_spin.eq.a_spin.and.
2212     &                            j_spin.eq.b_spin)then
2213                                 g = g_mo(i_orb,j_orb,a_orb,b_orb)
2214                              endif
2215                              if (i_spin.eq.b_spin.and.
2216     &                            j_spin.eq.a_spin)then
2217                                 g = g - g_mo(i_orb,j_orb,b_orb,a_orb)
2218                              endif
2219c                              write(6,*)' g:',g
2220                              r2(e,f,m,n) = r2(e,f,m,n) +
2221     &                             0.5d0*g* ( t2(f,b,j,n)*t2(e,a,i,m) -
2222     &                                        t2(e,b,j,n)*t2(f,a,i,m) -
2223     &                                        t2(f,b,j,m)*t2(e,a,i,n) +
2224     &                                        t2(e,b,j,m)*t2(f,a,i,n) -
2225     &                                        t2(f,b,m,n)*t2(e,a,i,j) +
2226     &                                        t2(e,b,m,n)*t2(f,a,i,j) -
2227     &                                        t2(e,f,m,i)*t2(a,b,n,j) +
2228     &                                        t2(e,f,n,i)*t2(a,b,m,j) +
2229     &                                  0.5d0*t2(e,f,i,j)*t2(a,b,m,n) )
2230
2231                           enddo
2232                        enddo
2233                     enddo
2234                  enddo
2235               enddo
2236            enddo
2237         enddo
2238      enddo
2239c      write(6,*)' t2(e,f,m,n) L+Q residual:'
2240c      call writet2(r2,nso)
2241      return
2242      end
2243      subroutine t2_update_rjh(g_mo, t2, r2, fock, nmo, nso, noso)
2244      implicit none
2245      integer nmo, nso, noso
2246      double precision g_mo(nmo,nmo,nmo,nmo)
2247      double precision fock(nmo,nmo)
2248      double precision t2(nso,nso,nso,nso)
2249      double precision r2(nso,nso,nso,nso)
2250      double precision g, denom, r2norm
2251c occupied orbitals
2252      integer i, j, k, l, m, n
2253      integer i_orb, j_orb, k_orb, l_orb, m_orb, n_orb
2254      integer i_spin, j_spin, k_spin, l_spin, m_spin, n_spin
2255c unoccupied orbitals
2256      integer a, b, c, d, e, f
2257      integer a_orb, b_orb, c_orb, d_orb, e_orb, f_orb
2258      integer a_spin, b_spin, c_spin, d_spin, e_spin, f_spin
2259      double precision ddot
2260      external ddot
2261c
2262c     update t with delta t (in t2)
2263c
2264      do e = noso+1, nso
2265c get spatial orbital and whether alpha or beta spin
2266         e_orb = (e-1)/2 + 1
2267         e_spin = 0
2268         if ((e-2*e_orb).ne.0)e_spin = 1
2269         do f = noso+1, nso
2270            f_orb = (f-1)/2 + 1
2271            f_spin = 0
2272            if ((f-2*f_orb).ne.0)f_spin = 1
2273            do m = 1, noso
2274               m_orb = (m-1)/2 + 1
2275               m_spin = 0
2276               if ((m-2*m_orb).ne.0)m_spin = 1
2277               do n = 1, noso
2278                  n_orb = (n-1)/2 + 1
2279                  n_spin = 0
2280                  if ((n-2*n_orb).ne.0)n_spin = 1
2281c
2282                  denom = fock(e_orb,e_orb) + fock(f_orb,f_orb) -
2283     &                    fock(m_orb,m_orb) - fock(n_orb,n_orb)
2284c                 write(6,*)' denom: ', denom
2285c
2286                  t2(e,f,m,n) = t2(e,f,m,n) - r2(e,f,m,n)/denom
2287               enddo
2288            enddo
2289         enddo
2290      enddo
2291      r2norm =  dsqrt(ddot(nso**4,r2,1,r2,1))
2292      write(6,*)' r2norm = ', r2norm
2293      return
2294      end
2295      subroutine t1_rjh(g_mo, r1, t2, fock, nmo, nso, noso)
2296      implicit none
2297      integer nmo, nso, noso
2298      double precision g_mo(nmo,nmo,nmo,nmo)
2299      double precision fock(nmo,nmo)
2300      double precision r1(nso,nso)
2301      double precision t2(nso,nso,nso,nso)
2302      double precision g, f_t, denom
2303c occupied orbitals
2304      integer i, j, k, l, m, n
2305      integer i_orb, j_orb, k_orb, l_orb, m_orb, n_orb
2306      integer i_spin, j_spin, k_spin, l_spin, m_spin, n_spin
2307c unoccupied orbitals
2308      integer a, b, c, d, e, f
2309      integer a_orb, b_orb, c_orb, d_orb, e_orb, f_orb
2310      integer a_spin, b_spin, c_spin, d_spin, e_spin, f_spin
2311c
2312c     put t into t2 expression keeping all terms linear in t2
2313c
2314c      write(6,*)' debug inside t1_rjh '
2315c      call jan_debug_print('MOINTS', g_mo, nmo, nmo,
2316c     $                      nmo, nmo)
2317c      write(6,*)' T2 with quadratic terms. '
2318c      call writet2(t2,nso)
2319c      write(6,*) ' Fock Matrix'
2320c      call output(fock,1,nmo,1,nmo,nmo,nmo,1)
2321c
2322      do e = noso+1, nso
2323c get spatial orbital and whether alpha or beta spin
2324         e_orb = (e-1)/2 + 1
2325         e_spin = 0
2326         if ((e-2*e_orb).ne.0)e_spin = 1
2327         do m = 1, noso
2328            m_orb = (m-1)/2 + 1
2329            m_spin = 0
2330            if ((m-2*m_orb).ne.0)m_spin = 1
2331c
2332            r1(e,m) = 0.0d0
2333            if (e_spin.eq.m_spin)then
2334               r1(e,m) = fock(e_orb,m_orb)
2335            endif
2336c
2337            do j = 1, noso
2338               j_orb = (j-1)/2 + 1
2339               j_spin = 0
2340               if ((j-2*j_orb).ne.0)j_spin = 1
2341               do b = noso+1, nso
2342                  b_orb = (b-1)/2 + 1
2343                  b_spin = 0
2344                  if ((b-2*b_orb).ne.0)b_spin = 1
2345c
2346                  f_t = 0.0d0
2347                  if (j_spin.eq.b_spin)then
2348                     f_t = fock(j_orb,b_orb)
2349                  endif
2350c
2351                  r1(e,m) = r1(e,m) + f_t*t2(e,b,m,j)
2352c
2353                  do i = 1, noso
2354                     i_orb = (i-1)/2 + 1
2355                     i_spin = 0
2356                     if ((i-2*i_orb).ne.0)i_spin = 1
2357c
2358                     g=0.0d0
2359                     if (i_spin.eq.m_spin.and.j_spin.eq.b_spin)then
2360                        g = g_mo(i_orb,j_orb,m_orb,b_orb)
2361                     endif
2362                     if (i_spin.eq.b_spin.and.j_spin.eq.m_spin)then
2363                        g = g - g_mo(i_orb,j_orb,b_orb,m_orb)
2364                     endif
2365
2366c                     write(6,*)' e,m,j,b,i,g:',e,m,j,b,i,g
2367c
2368                     r1(e,m) = r1(e,m) - 0.5d0*g*t2(e,b,i,j)
2369c
2370                  enddo
2371               enddo
2372            enddo
2373            do b = noso+1, nso
2374               b_orb = (b-1)/2 + 1
2375               b_spin = 0
2376               if ((b-2*b_orb).ne.0)b_spin = 1
2377c
2378               do i = 1, noso
2379                  i_orb = (i-1)/2 + 1
2380                  i_spin = 0
2381                  if ((i-2*i_orb).ne.0)i_spin = 1
2382c
2383                  do a = noso+1, nso
2384                     a_orb = (a-1)/2 + 1
2385                     a_spin = 0
2386                     if ((a-2*a_orb).ne.0)a_spin = 1
2387                     g=0.0d0
2388                     if (i_spin.eq.a_spin.and.e_spin.eq.b_spin)then
2389                        g = g_mo(i_orb,e_orb,a_orb,b_orb)
2390                     endif
2391                     if (i_spin.eq.b_spin.and.e_spin.eq.a_spin)then
2392                        g = g - g_mo(i_orb,e_orb,b_orb,a_orb)
2393                     endif
2394c
2395c                     write(6,*)' g:',g
2396c
2397                     r1(e,m) = r1(e,m) + 0.5d0*g*t2(a,b,i,m)
2398c
2399                  enddo
2400               enddo
2401            enddo
2402         enddo
2403      enddo
2404c      write(6,*)' t1(e,m) residual:'
2405c      call writet1(r1,nso)
2406      return
2407      end
2408      subroutine t1_update_rjh(g_mo, t1, r1, fock, nmo, nso, noso)
2409      implicit none
2410      integer nmo, nso, noso
2411      double precision g_mo(nmo,nmo,nmo,nmo)
2412      double precision fock(nmo,nmo)
2413      double precision r1(nso,nso)
2414      double precision t1(nso,nso)
2415      double precision g, denom, r1norm
2416c occupied orbitals
2417      integer i, j, k, l, m, n
2418      integer i_orb, j_orb, k_orb, l_orb, m_orb, n_orb
2419      integer i_spin, j_spin, k_spin, l_spin, m_spin, n_spin
2420c unoccupied orbitals
2421      integer a, b, c, d, e, f
2422      integer a_orb, b_orb, c_orb, d_orb, e_orb, f_orb
2423      integer a_spin, b_spin, c_spin, d_spin, e_spin, f_spin
2424      double precision ddot
2425      external ddot
2426c
2427c     update t with delta t (in t2)
2428c
2429      do e = noso+1, nso
2430c get spatial orbital and whether alpha or beta spin
2431         e_orb = (e-1)/2 + 1
2432         e_spin = 0
2433         if ((e-2*e_orb).ne.0)e_spin = 1
2434         do m = 1, noso
2435            m_orb = (m-1)/2 + 1
2436            m_spin = 0
2437            if ((m-2*m_orb).ne.0)m_spin = 1
2438c
2439            denom = fock(e_orb,e_orb) - fock(m_orb,m_orb)
2440c            write(6,*)' denom: ', denom
2441c
2442            t1(e,m) = t1(e,m) - r1(e,m)/denom
2443         enddo
2444      enddo
2445      r1norm =  dsqrt(ddot(nso**2,r1,1,r1,1))
2446      write(6,*)' r1norm = ', r1norm
2447      return
2448      end
2449      subroutine jan_h(
2450     $     rtdb, basis,
2451     $     n1, n2,
2452     $     lda1, lda2,
2453     $     c1t, c2t,
2454     $     h)
2455      implicit none
2456#include "errquit.fh"
2457#include "bas.fh"
2458#include "global.fh"
2459#include "mafdecls.fh"
2460      integer rtdb, basis
2461      integer lda1, lda2, n1, n2
2462      double precision c1t(lda1, n1), c2t(lda2, n2)
2463      double precision h(n1, n2)
2464c
2465c     Return hij = sum(kl) c1t(i,k) hAO(k,l) c2t(j,l)
2466c
2467c     Note that the transposed MO coeffs are passed in
2468c     to be compatible with jan_full_transform
2469c
2470      integer geom, nbf, g_tmp
2471      integer l_tmp1, k_tmp1, l_tmp2, k_tmp2
2472c
2473      call int_init(rtdb, 1, basis)
2474      if (.not. bas_geom(basis, geom))
2475     $     call errquit('jan_transform: basis ', basis, BASIS_ERR)
2476      call schwarz_init(geom, basis)
2477      if (.not. bas_numbf(basis, nbf))
2478     $     call errquit('jan_h: nbf',basis, BASIS_ERR)
2479c
2480      if (.not. ma_push_get(mt_dbl, nbf*nbf,'tmp1', l_tmp1, k_tmp1))
2481     $     call errquit('tmp1', nbf*nbf, MA_ERR)
2482      if (.not. ma_push_get(mt_dbl, nbf*nbf,'tmp2', l_tmp2, k_tmp2))
2483     $     call errquit('tmp2', nbf*nbf, MA_ERR)
2484c
2485      if (.not. ga_create(mt_dbl, nbf, nbf, 'tmp', 1, 1, g_tmp))
2486     &     call errquit('scf_v_g: tmp', 0, GA_ERR)
2487      call ga_zero(g_tmp)
2488      call int_1e_ga(basis, basis, g_tmp, 'kinetic', .false.)
2489      call int_1e_ga(basis, basis, g_tmp, 'potential', .false.)
2490      call ga_get(g_tmp, 1, nbf, 1, nbf, dbl_mb(k_tmp1), nbf)
2491      if (.not. ga_destroy(g_tmp)) call errquit('ga?',0, GA_ERR)
2492c
2493      call dgemm('n', 'n', n1, nbf, nbf,
2494     $     1.0d0, c1t, lda1, dbl_mb(k_tmp1), nbf,
2495     $     0.0d0, dbl_mb(k_tmp2), n1)
2496      call dgemm('n', 't', n1, n2, nbf,
2497     $     1.0d0, dbl_mb(k_tmp2), n1, c2t, lda2,
2498     $     0.0d0, h, n1)
2499c
2500c      write(6,*) ' Transformed H'
2501c      call output(h, 1, n1, 1, n2, n1, n2, 1)
2502c
2503      call schwarz_tidy()
2504      call int_terminate
2505c
2506      if (.not. ma_chop_stack(l_tmp1)) call errquit('ma',0, MA_ERR)
2507c
2508      return
2509      end
2510      subroutine uccsdtest_lambda(nbf, nmo, no, nv,
2511     $        t1, ct, pt, ht, work)
2512      implicit none
2513c
2514      integer nbf, nmo, no, nv
2515      double precision t1(nv, no)
2516      double precision ct(nmo, nbf)
2517      double precision pt(nmo, nbf)
2518      double precision ht(nmo, nbf)
2519      double precision work(nmo, nmo)
2520c
2521      integer i, a
2522c
2523c     Particle transformation
2524c
2525      call dfill(nmo*nmo, 0.0d0, work, 1)
2526      call dfill(nmo, 1.0d0, work, nmo+1)
2527      do i = 1, no
2528         do a = 1, nv
2529            work(i,a+no) = work(i,a+no) - t1(a,i)
2530         end do
2531      end do
2532      call dgemm('t','n',nmo,nbf,nmo,
2533     $     1.0d0, work, nmo, ct, nmo,
2534     $     0.0d0, pt, nmo)
2535c
2536c     Hole transformation
2537c
2538      call dfill(nmo*nmo, 0.0d0, work, 1)
2539      call dfill(nmo, 1.0d0, work, nmo+1)
2540      do i = 1, no
2541         do a = 1, nv
2542            work(a+no,i) = work(a+no,i) + t1(a,i)
2543         end do
2544      end do
2545      call dgemm('t','n',nmo,nbf,nmo,
2546     $     1.0d0, work, nmo, ct, nmo,
2547     $     0.0d0, ht, nmo)
2548c
2549      end
2550      subroutine get_new_f_g(rtdb, basis, fock, g_mo, t1, nso, nmo)
2551      implicit none
2552#include "errquit.fh"
2553#include "global.fh"
2554#include "mafdecls.fh"
2555#include "bas.fh"
2556#include "geom.fh"
2557#include "rtdb.fh"
2558#include "inp.fh"
2559      integer rtdb
2560c
2561      integer nmo, g_tmp, l_occ, k_occ,
2562     $     l_eval, k_eval, l_mos, k_mos, l_most, k_most
2563      integer nocc, nvirt, nopen, nclosed, nso, noso
2564      integer l_t, k_t, l_t2, k_t2, l_f, k_f
2565      integer l_t1, k_t1, l_r1, k_r1
2566      integer l_pmost, k_pmost, l_hmost, k_hmost
2567      integer l_work, k_work, l_h, k_h
2568      integer l_t1_vomos, k_t1_vomos
2569c
2570      double precision g_mo(nmo,nmo,nmo,nmo)
2571      double precision fock(nmo,nmo)
2572      double precision t1(nso,nso)
2573      integer basis, nbf
2574      character*255 movecs      ! Name of movector file
2575      character*80 title, name_of_basis, scftype
2576      integer nbf_file, nsets, nmo_file(2)
2577      logical movecs_read, movecs_read_header
2578      external movecs_read, movecs_read_header
2579c
2580      logical int_normalize
2581      external int_normalize
2582c
2583c     Read the MO vectors and evals from a RHF calculation
2584c
2585      call util_file_name('movecs',.false.,.false.,movecs)
2586      if (.not. movecs_read_header(movecs, title, name_of_basis,
2587     $     scftype, nbf_file, nsets, nmo_file, 2)) call errquit
2588     $     ('jantest: failed to read movecs header',911,
2589     &       DISK_ERR)
2590c      write(6,*) ' Read movecs header from ', movecs
2591c      write(6,*) ' Job title :                ',
2592c     $     title(1:inp_strlen(title))
2593c      write(6,*) ' Basis name:                ',
2594c     $     name_of_basis(1:inp_strlen(name_of_basis))
2595      nmo = nmo_file(1)
2596      if (.not. rtdb_get(rtdb, 'scf:nclosed', mt_int, 1, nocc))
2597     $     call errquit('nocc?',0, RTDB_ERR)
2598      if (.not. rtdb_get(rtdb, 'scf:nopen', mt_int, 1, nopen))
2599     $     call errquit('nopen?',0, RTDB_ERR)
2600      if (nopen .ne. 0) call errquit('asjdlfkadjsl',0, UNKNOWN_ERR)
2601      if (.not. bas_numbf(basis, nbf)) call errquit
2602     $     ('scf_init: basis info',0, BASIS_ERR)
2603      nvirt= nmo - nocc
2604c      write(6,*) ' No. of closed shells       ', nocc
2605c      write(6,*) ' No. of molecular orbitals: ', nmo
2606c      write(6,*) ' No. of basis functions:    ', nbf
2607c
2608      if (.not. ga_create(mt_dbl, nbf, nmo, 'tmp', 0, 0, g_tmp))
2609     &     call errquit('scf_v_g: tmp', 0, GA_ERR)
2610      if (.not. ma_push_get(mt_dbl, nbf,'occ',l_occ, k_occ))
2611     $     call errquit('ma occ', nbf, MA_ERR)
2612      if (.not. ma_push_get(mt_dbl, nbf,'eval',l_eval, k_eval))
2613     $     call errquit('ma eval', nbf, MA_ERR)
2614      if (.not. ma_push_get(mt_dbl, nbf*nbf,'mos', l_mos, k_mos))
2615     $     call errquit('ma mos', nbf*nbf, MA_ERR)
2616      if (.not. ma_push_get(mt_dbl, nbf*nbf,'mos', l_most, k_most))
2617     $     call errquit('ma mos', nbf*nbf, MA_ERR)
2618c
2619      if (.not. movecs_read(movecs, 1, dbl_mb(k_occ), dbl_mb(k_eval),
2620     $     g_tmp)) call errquit('movecs_read of amos failed ',0,
2621     &       DISK_ERR)
2622      call ga_get(g_tmp, 1, nbf, 1, nmo, dbl_mb(k_mos), nbf)
2623      call util_transpose(dbl_mb(k_mos),nbf,dbl_mb(k_most),nmo,
2624     $     nbf,nmo)
2625c
2626c      write(6,*) ' Orbital eigenvalues '
2627c      call output(dbl_mb(k_eval),1,nmo,1,1,nmo,1,1)
2628c      write(6,*) ' MOs'
2629c      call output(dbl_mb(k_mos),1,nbf,1,nmo,nbf,nmo,1)
2630c      write(6,*) ' MOs T'
2631c      call output(dbl_mb(k_most),1,nmo,1,nbf,nmo,nbf,1)
2632c
2633      if (.not. ga_destroy(g_tmp)) call errquit(' ga bad?',0, GA_ERR)
2634c
2635c     Transform MO integrals in Dirac order to new T1 bases
2636c
2637      if (.not. ma_push_get(mt_dbl, nbf*nbf,'pmost', l_pmost, k_pmost))
2638     $     call errquit('ma pmost', nbf*nbf, MA_ERR)
2639      if (.not. ma_push_get(mt_dbl, nbf*nbf,'hmost', l_hmost, k_hmost))
2640     $     call errquit('ma hmost', nbf*nbf, MA_ERR)
2641      if (.not. ma_push_get(mt_dbl, nbf*nbf,'work', l_work, k_work))
2642     $     call errquit('ma work', nbf*nbf, MA_ERR)
2643      if (.not. ma_push_get(mt_dbl, nbf*nbf,'1e-', l_h, k_h))
2644     $     call errquit('ma h', nbf*nbf, MA_ERR)
2645      if (.not. ma_push_get(mt_dbl, nvirt*nocc,'t1_vomos',
2646     $     l_t1_vomos, k_t1_vomos))
2647     $     call errquit('ma t1_vomos', nvirt*nocc, MA_ERR)
2648      call pack_t1(t1, dbl_mb(k_t1_vomos), nso, nmo, nocc, nvirt)
2649      call uccsdtest_lambda(nbf, nmo, nocc, nvirt,
2650     $                      dbl_mb(k_t1_vomos), dbl_mb(k_most),
2651     $                      dbl_mb(k_pmost),
2652     $                      dbl_mb(k_hmost), dbl_mb(k_work))
2653c
2654c      write(6,*) ' P mos'
2655c      call output(dbl_mb(k_pmost), 1, nmo, 1, nmo, nmo, nmo, 1)
2656c
2657c      write(6,*) ' H mos'
2658c      call output(dbl_mb(k_hmost), 1, nmo, 1, nmo, nmo, nmo, 1)
2659c
2660      call jan_h(rtdb, basis, nmo, nmo, nmo, nmo,
2661     $           dbl_mb(k_pmost), dbl_mb(k_hmost), dbl_mb(k_h))
2662      call jan_full_transform(
2663     $     rtdb, basis,
2664     $     nmo, nmo, nmo, nmo,
2665     $     nmo, nmo, nmo, nmo,
2666     $     dbl_mb(k_pmost), dbl_mb(k_pmost), dbl_mb(k_hmost),
2667     $     dbl_mb(k_hmost), g_mo, 'Dirac')
2668c      call jan_debug_print('MOINTS', g_mo, nmo, nmo,
2669c     $                      nmo, nmo)
2670      call build_f(fock, g_mo, dbl_mb(k_h), nmo, nocc)
2671c
2672c     Tidy up
2673c
2674      if (.not. ma_chop_stack(l_occ)) call errquit(' ma chop?', 0,
2675     &       MA_ERR)
2676c
2677      return
2678c
2679      end
2680      subroutine pack_t1(t1, t1_vomos, nso, nmo, nocc, nvirt)
2681      implicit none
2682c
2683      integer nso, nmo, nocc, nvirt
2684      integer i, j, i_orb, j_orb, i_so, j_so
2685      double precision t1_vomos(nvirt, nocc)
2686      double precision t1(nso,nso)
2687c
2688c      write(6,*)' t1 from pack_t1 '
2689c      do i = 1, nso
2690c         do j = 1, nso
2691c            if (abs(t1(i,j)).gt.1d-6)
2692c     &         write(6,*)i,j,t1(i,j)
2693c         enddo
2694c      enddo
2695c
2696      call dfill(nvirt*nocc, 0.0d0, t1_vomos, 1)
2697c
2698      do i = 1, nvirt
2699         do j = 1, nocc
2700            i_orb = nocc + i
2701            j_orb = j
2702            i_so = 2*(i_orb-1)+1
2703            j_so = 2*(j_orb-1)+1
2704            t1_vomos(i,j) = t1(i_so,j_so)
2705         enddo
2706      enddo
2707c      write(6,*)' t1_vomos '
2708c      do i = 1, nvirt
2709c         do j = 1, nocc
2710c            if (abs(t1_vomos(i,j)).gt.1d-6)
2711c     &         write(6,*)i,j,t1_vomos(i,j)
2712c         enddo
2713c      enddo
2714      return
2715      end
2716      subroutine build_f(fock, g_mo, h_mo, nmo, nocc)
2717      implicit none
2718c
2719      integer nmo, nocc
2720      double precision g_mo(nmo,nmo,nmo,nmo)
2721      double precision fock(nmo,nmo), h_mo(nmo,nmo)
2722c
2723c     fij = hij + 2*<ik|jk> - <ik|kj>
2724c
2725      integer i, j, k
2726c
2727      call dfill(nmo*nmo, 0.0d0, fock, 1)
2728c
2729      do j = 1, nmo
2730         do i = 1, nmo
2731            fock(i,j) = h_mo(i,j)
2732            do k = 1, nocc
2733               fock(i,j) = fock(i,j) + 2.0d0*g_mo(i,k,j,k)
2734     &                               - g_mo(i,k,k,j)
2735            end do
2736         end do
2737      end do
2738c
2739c      write(6,*) ' FOCK'
2740c      call output(fock, 1, nmo, 1, nmo, nmo, nmo, 1)
2741c
2742      return
2743      end
2744
2745c $Id$
2746