      logical function tce_ccsdtq_response_driver(title,thresh,maxiter,
     &        d_tr1,k_tr1_offset,size_tr1,d_tr2,k_tr2_offset,size_tr2,
     &        d_tr3,k_tr3_offset,size_tr3,d_tr4,k_tr4_offset,size_tr4,
     &        d_t1,k_t1_offset,size_t1,d_t2,k_t2_offset,size_t2,
     &        d_t3,k_t3_offset,size_t3,d_t4,k_t4_offset,size_t4,
     &        d_f1,k_f1_offset,d_v2,k_v2_offset,d_o1,k_o1_offset,omega,
     &        diis_tr4)
c
c $Id$
c 
c Main routine for many-electron theory calculations.
c Some of the subroutines have been generated by 
c operator/tensor contraction engines.
c
c Written by Jeff Hammond, January 2008.
c
      implicit none
#include "mafdecls.fh"
#include "tcgmsg.fh"
#include "global.fh"
#include "bas.fh"
#include "geom.fh"
#include "sym.fh"
#include "util.fh"
#include "msgids.fh"
#include "stdio.fh"
#include "sf.fh"
#include "inp.fh"
#include "errquit.fh"
c#include "tce.fh"
c#include "tce_main.fh"
c#include "tce_prop.fh"
c#include "tce_ints.fh"
c#include "tce_amps.fh"
c#include "tce_diis.fh"
c
c     CI, CC, & MBPT
c
      integer iter,maxiter
      logical nodezero         ! True if node 0
      logical recompf          ! True if recompute Fock
      double precision cpu     ! CPU sec counter
      double precision wall    ! WALL sec counter
      integer irrep
      integer irrep_g
      integer d_tr1,k_tr1_offset,size_tr1       
      integer d_tr2,k_tr2_offset,size_tr2    
      integer d_tr3,k_tr3_offset,size_tr3         
      integer d_tr4,k_tr4_offset,size_tr4      
      integer d_rr1,d_rr2,d_rr3,d_rr4
      integer d_o1,k_o1_offset
      integer d_f1,k_f1_offset     
      integer d_v2,k_v2_offset   
      integer d_t1,k_t1_offset,size_t1          
      integer d_t2,k_t2_offset,size_t2       
      integer d_t3,k_t3_offset,size_t3          
      integer d_t4,k_t4_offset,size_t4       
      integer d_y1,k_y1_offset 
      integer d_y2,k_y2_offset     
      integer d_y3,k_y3_offset   
      integer d_y4,k_y4_offset 
      double precision rr1,rr2,rr3,rr4
      double precision residual! Largest residual
      double precision thresh
      double precision ddotfile
      double precision omega
      external ddotfile
      integer dummy            ! Dummy argument for DIIS
      character*255 filename
      character*20 title
      character*4 irrepname
      character*6 rr1filename,rr2filename,rr3filename,rr4filename
      data rr1filename/'rr1   '/
      data rr2filename/'rr2   '/
      data rr3filename/'rr3   '/
      data rr4filename/'rr4   '/
      logical diis_tr4
#if defined(CCSDTQ)
      nodezero=(ga_nodeid().eq.0)
      call tce_diis_init()
      do iter=1,maxiter
        cpu=-util_cpusec()
        wall=-util_wallsec()
        if (nodezero.and.(iter.eq.1)) write(LuOut,9400) title
        call tce_filename(rr1filename,filename)
        call createfile(filename,d_rr1,size_tr1)
        call tce_zero(d_rr1,size_tr1)
c        write(LuOut,*) "daxpyfile"
        call daxpyfile(1,omega,d_tr1,d_rr1,size_tr1)
c        write(LuOut,*) "eomccsdtq_x1"
        call eomccsdtq_x1(d_f1,d_rr1,d_t1,d_t2,d_v2,
     1       d_tr1,d_tr2,d_tr3,k_f1_offset,k_tr1_offset,
     4       k_t1_offset,k_t2_offset,k_v2_offset,
     5       k_tr1_offset,k_tr2_offset,k_tr3_offset)
c        write(LuOut,*) "ccsdtq1_o1"
        call ccsdtq_o1(d_rr1,d_o1,d_t1,d_t2,k_tr1_offset,k_o1_offset,
     &       k_t1_offset,k_t2_offset)
c
        call tce_filename(rr2filename,filename)
        call createfile(filename,d_rr2,size_tr2)
        call tce_zero(d_rr2,size_tr2)
c        write(LuOut,*) "daxpyfile"
        call daxpyfile(1,omega,d_tr2,d_rr2,size_tr2)
c        write(LuOut,*) "eomccsdtq_x2"
        call eomccsdtq_x2(d_f1,d_rr2,d_t1,d_t2,d_t3,d_v2,
     1       d_tr1,d_tr2,d_tr3,d_tr4,k_f1_offset,k_tr2_offset,
     4       k_t1_offset,k_t2_offset,k_t3_offset,k_v2_offset,
     5       k_tr1_offset,k_tr2_offset,k_tr3_offset,k_tr4_offset)
c        write(LuOut,*) "ccsdtq_o2"
        call ccsdtq_o2(d_rr2,d_o1,d_t1,d_t2,d_t3,k_tr2_offset,
     1       k_o1_offset,k_t1_offset,k_t2_offset,k_t3_offset)
c
        call tce_filename(rr3filename,filename)
        call createfile(filename,d_rr3,size_tr3)
        call tce_zero(d_rr3,size_tr3)
c        write(LuOut,*) "daxpyfile"
        call daxpyfile(1,omega,d_tr3,d_rr3,size_tr3)
c        write(LuOut,*) "eomccsdtq_x3"
        call eomccsdtq_x3(d_f1,d_rr3,d_t1,d_t2,d_t3,d_t4,
     1       d_v2,d_tr1,d_tr2,d_tr3,d_tr4,k_f1_offset,k_tr3_offset,
     4       k_t1_offset,k_t2_offset,k_t3_offset,k_t4_offset,
     5       k_v2_offset,k_tr1_offset,k_tr2_offset,
     7       k_tr3_offset,k_tr4_offset)
c        write(LuOut,*) "ccsdtq_o3"
        call ccsdtq_o3(d_rr3,d_o1,d_t1,d_t2,d_t3,d_t4,
     &       k_tr3_offset,k_o1_offset,
     &       k_t1_offset,k_t2_offset,k_t3_offset,k_t4_offset)
c
        call tce_filename(rr4filename,filename)
        call createfile(filename,d_rr4,size_tr4)
        call tce_zero(d_rr4,size_tr4)
c        write(LuOut,*) "daxpyfile"
        call daxpyfile(1,omega,d_tr4,d_rr4,size_tr4)
c        write(LuOut,*) "eomccsdtq_x4"
        call eomccsdtq_x4(d_f1,d_rr4,d_t1,d_t2,d_t3,d_t4,
     1       d_v2,d_tr1,d_tr2,d_tr3,d_tr4,k_f1_offset,k_tr4_offset,
     4       k_t1_offset,k_t2_offset,k_t3_offset,k_t4_offset,
     5       k_v2_offset,k_tr1_offset,k_tr2_offset,
     7       k_tr3_offset,k_tr4_offset)
c        write(LuOut,*) "ccsdtq_o4"
        call ccsdtq_o4(d_rr4,d_o1,d_t1,d_t2,d_t3,d_t4,
     &       k_tr4_offset,k_o1_offset,
     &       k_t1_offset,k_t2_offset,k_t3_offset,k_t4_offset)
c
c        write(LuOut,*) "reconcilefile"
        call reconcilefile(d_rr1,size_tr1)
        call reconcilefile(d_rr2,size_tr2)
        call reconcilefile(d_rr3,size_tr3)
        call reconcilefile(d_rr4,size_tr4)
c        write(LuOut,*) "tce_residual_tr1"
        call tce_residual_tr1(d_rr1,k_tr1_offset,rr1)
c        write(LuOut,*) "tce_residual_tr2"
        call tce_residual_tr2(d_rr2,k_tr2_offset,rr2)
c        write(LuOut,*) "tce_residual_tr3"
        call tce_residual_tr3(d_rr3,k_tr3_offset,rr3)
c        write(LuOut,*) "tce_residual_tr4"
        call tce_residual_tr4(d_rr4,k_tr4_offset,rr4)
        residual = max(rr1,rr2,rr3,rr4)
        cpu=cpu+util_cpusec()
        wall=wall+util_wallsec()
        if (nodezero) write(LuOut,9420) iter,residual,cpu,wall
        if (residual .lt. thresh) then
          if (nodezero) write(LuOut,9410)
          call deletefile(d_rr4)
          call deletefile(d_rr3)
          call deletefile(d_rr2)
          call deletefile(d_rr1)
          call tce_diis_tidy()
          tce_ccsdtq_response_driver=.true.
        endif
        call tce_diis2(.false.,iter,.true.,.true.,.true.,diis_tr4,
     1       d_rr1,d_tr1,k_tr1_offset,size_tr1,
     3       d_rr2,d_tr2,k_tr2_offset,size_tr2,
     5       d_rr3,d_tr3,k_tr3_offset,size_tr3,
     7       d_rr4,d_tr4,k_tr4_offset,size_tr4)
        call deletefile(d_rr4)
        call deletefile(d_rr3)
        call deletefile(d_rr2)
        call deletefile(d_rr1)
        if (nodezero) call util_flush(LuOut)
      enddo
      call tce_diis_tidy()
#endif
      tce_ccsdtq_response_driver=.false.
      return  
c
c     ======
c     Format
c     ======
c
 9000 format(1x,A,' file size   = ',i16)
 9010 format(1x,A,' file name   = ',A)
 9090 format(1x,A,' file handle = ',i10)
 9020 format(1x,'Cpu & wall time / sec',2f15.1)
 9480 format(1x,'Cpu & wall time / sec for ',A,2f15.1)
 9050 format(/,1x,A,' iterations',/,
     1  1x,'--------------------------------------------------------',/
     2  1x,'Iter          Residuum       Correlation     Cpu    Wall',/
     3  1x,'--------------------------------------------------------')
 9060 format(
     1  1x,'--------------------------------------------------------',/
     2  1x,'Iterations converged')
 9070 format(1x,A,' correlation energy / hartree = ',f25.15)
 9080 format(1x,A,' total energy / hartree       = ',f25.15)
 9100 format(1x,i4,2f18.13,2f8.1)
 9120 format(1x,A)
 9250 format(1x,'Ground-state symmetry is ',A4)
 9210 format(/,1x,'Iteration ',i3,' using ',i4,' trial vectors')
 9230 format(1x,f17.13,f18.13,f11.5,2f8.1)
 9240 format(1x,
     1'--------------------------------------------------------------'
     2,/,1x,'Iterations converged')
 9310 format(1x,A,' ground state energy / hartree  =',f25.15)
 9400 format(/,1x,A,' iterations',/,
     1  1x,'--------------------------------------',/
     2  1x,'Iter          Residuum     Cpu    Wall',/
     3  1x,'--------------------------------------')
 9410 format(
     1  1x,'--------------------------------------',/
     2  1x,'Iterations converged')
 9420 format(1x,i4,f18.13,2f8.1)
 9440 format(1x,A3,' axis ( ',A4,'symmetry)')
      end
