C 19 Sep 05 - mws - adjustment for h or i function runs
C  5 Jul 05 - MWS - SELECT NEW ATOM,BASIS,EFP,PCM,DAF DIMENSIONS
C 13 feb 05 - mws - pad common block nshel
C  5 Feb 05 - OQ  - Iterative procedures for extended TDHF properties
C
C*MODULE TDXITR  *DECK CompUsum
      subroutine CompUsum(Usum,Ux,Uy,Ty,norb,nocc,mult)
      implicit none
C     Arguments
      integer norb,nocc
      double precision mult
      double precision Usum(norb*norb)
      double precision Ux(norb*norb)
      double precision Uy(norb*norb)
      character*1 Ty
C     Variables
      double precision dummy
      double precision ONE
      parameter (ONE=1.0D+00)
      character*1 sN,sT,sY
      parameter (sN='N',sT='T')
C     Body of the subprogram
      if(Ty.eq.'a') then
         dummy = mult
         sY = sT
      else if(Ty.eq.'f') then
         dummy = -mult
         sY = sN
      else if(Ty.eq.'o') then
         dummy = mult
         sY = sN
      else
         write(*,*) 'Wrong type for Ty [',Ty,'] in CompUsum'
         call abrtx("COMPUSUM")
      end if
      call DGEMM(sN,sY,norb,norb,nocc,
     &           dummy,Ux,norb,Uy,norb,ONE,Usum,norb)
      end
C
C*MODULE TDXITR  *DECK CompDDterm
      subroutine CompDDterm(DD,C0,Usum,norb,nparams)
      implicit none
C     Arguments
      integer norb,nparams
      double precision DD(norb*norb,nparams)
      double precision C0(norb*norb)
      double precision Usum(norb*norb,nparams)
C     Used common block
      DOUBLE PRECISION XX
      COMMON /FMCOM / XX(1)
C     Variables
      integer lWork
      integer iparam
      double precision ZERO,ONE,TWO
      parameter (ZERO=0.0D+00,ONE=1.0D+00,TWO=2.0D+00)
      character*1 sN,sT
      parameter (sN='N',sT='T')
C     Body of the subprogram
      call NEWMEM(lWork,norb*norb)
      do iparam=1,nparams
         call DGEMM(sN,sN,norb,norb,norb,
     &              ONE,C0,norb,Usum(1,iparam),norb,ZERO,XX(lWork),norb)
         call DGEMM(sN,sT,norb,norb,norb,
     &              TWO,XX(lWork),norb,C0,norb,ZERO,DD(1,iparam),norb)
      end do
      call DELETEMEM(lWork,norb*norb)
      end
C
C*MODULE TDXITR  *DECK CompQterm
      subroutine CompQterm(Q,Gx,Ux,Tx,Uy,Ey,norb,mult)
      implicit none
C     Arguments
      integer norb
      double precision mult
      double precision Q(norb*norb)
      double precision Gx(norb*norb),Ux(norb*norb)
      double precision Uy(norb*norb),Ey(norb*norb)
      character*1 Tx
C     Variables and constants
      double precision dummy
      double precision ONE
      parameter (ONE=1.0D+00)
      character*1 sT,sN,sX
      parameter (sT='T',sN='N')
C     Body of the subprogram
      dummy = mult
      call DGEMM(sN,sN,norb,norb,norb,
     &           dummy,Gx,norb,Uy,norb,ONE,Q,norb)
      if(Tx.eq.'a') then
         dummy = mult
         sX = sT
      else if(Tx.eq.'o') then
         dummy = -mult
         sX = sN
      else
         write(*,*) 'Wrong type for Tx [',Tx,'] in CompQterm'
         call abrtx("COMPQTERM")
      end if
      call DGEMM(sX,sN,norb,norb,norb,dummy,Ux,norb,Ey,norb,ONE,Q,norb)
      end
C
C*MODULE TDXITR  *DECK CompQterm2
      subroutine CompQterm2(U,w,E0,norb)
      implicit none
C     Arguments
      integer norb
      double precision w
      double precision U(norb,norb),E0(norb)
C     Variables
      integer j
      double precision dummy
C     Body of the subprogram
      do j=1,norb
         dummy = w-E0(j)
         call DSCAL(norb,dummy,U(1,j),1)
      end do
      end
C
C*MODULE TDXITR  *DECK CompleteGUXE
      subroutine CompleteGUXE(C0,E0,omega,G,U,X,E,norb,nocc,nxyz)
      implicit none
C     Arguments
      integer norb,nocc,nxyz
      double precision omega
      double precision G(norb,norb,nxyz)
      double precision U(norb,norb,nxyz)
      double precision X(norb,norb,nxyz)
      double precision E(norb,norb,nxyz)
      double precision C0(norb,norb)
      double precision E0(norb)
C     Variables and constants
      logical bX,bF,bQ
      double precision ONE
      parameter (ONE=1.0D+00)
C     Called subprogram
      logical LoadFQX
C     Body of the subprogram
      call CheckFQX(bF,bQ,bX)
      if(bF) then
         if(.not.LoadFQX(U,.true.,U,.false.,U,.false.,norb,nxyz)) then
            call abrtx("CompleteGUXE: not loadfqx:1")
         end if
         call DAXPY(norb*norb*nxyz,ONE,U,1,G,1)
      end if
      call MakeGij(G,C0,G,norb,nxyz)
      if(bQ) then
         if(.not.LoadFQX(U,.false.,U,.true.,U,.false.,norb,nxyz)) then
            call abrtx("CompleteGUXE: not loadfqx:2")
         end if
         call DAXPY(norb*norb*nxyz,ONE,G,1,U,1)
      else
         call DCOPY(norb*norb*nxyz,G,1,U,1)
      end if
      call MakeEStep1(E,U,norb,nocc,nxyz)
      call MakeUStep1(U,E0,omega,U,norb,nocc,nxyz)
      if(bX) then
         if(.not.LoadFQX(X,.false.,X,.false.,X,.true.,norb,nxyz)) then
            call abrtx("CompleteGUXE: not loadfqx:3")
         end if
      end if
      call MakeUStep2(U,X,bX,norb,nocc,nxyz)
      call MakeEStep2(E,X,bX,E0,omega,norb,nocc,nxyz)
C
C      call PrintMM(G,norb,nxyz,'G')
C      call PrintMM(U,norb,nxyz,'U')
C      if(bX) call PrintMM(X,norb,nxyz,'X')
C      call PrintMM(E,norb,nxyz,'E')
      end
C
C*MODULE TDXITR  *DECK MakeGij
      SUBROUTINE MakeGij(G,C,F,norb,nxyz)
      IMPLICIT NONE
C     Arguments
      integer norb,nxyz
      DOUBLE PRECISION G(norb*norb,nxyz),F(norb*norb,nxyz)
      DOUBLE PRECISION C(norb*norb)
C     Used common block
      DOUBLE PRECISION XX
      COMMON /FMCOM / XX(1)
C     Variables and constants
      integer ixyz,ialloc
      double precision ZERO,ONE
      parameter (ZERO=0.0D+00,ONE=1.0D+00)
      character*1 sT,sN
      parameter (sT='T',sN='N')
C     Body of the subprogram
      call NEWMEM(ialloc,norb*norb)
      do ixyz=1,nxyz
         call DGEMM(sT,sN,norb,norb,norb,
     &              ONE,C,norb,F(1,ixyz),norb,ZERO,XX(ialloc),norb)
         call DGEMM(sN,sN,norb,norb,norb,
     &              ONE,XX(ialloc),norb,C,norb,ZERO,G(1,ixyz),norb)
      end do
      call DELETEMEM(ialloc,norb*norb)
      END
C
C*MODULE TDXITR  *DECK MakeUStep1
      subroutine MakeUStep1(U,E0,omega,G,norb,nocc,nxyz)
      implicit none
C     Arguments
      integer norb,nocc,nxyz
      double precision omega
      double precision U(norb,norb,nxyz),G(norb,norb,nxyz)
      double precision E0(norb)
C     Variables
      integer i,j,ixyz
      double precision dummy
C     Body of the subprogram
      do ixyz = 1,nxyz
         do j=1,nocc
            do i=nocc+1,norb
               dummy = G(i,j,ixyz)/(E0(j)-E0(i)-omega)
               U(i,j,ixyz)=dummy
            end do
         end do
         do j=nocc+1,norb
            do i=1,nocc
               dummy = G(i,j,ixyz)/(E0(j)-E0(i)-omega)
               U(i,j,ixyz)=dummy
            end do
         end do
      end do
      end
C
C*MODULE TDXITR  *DECK MakeUStep2
      subroutine MakeUStep2(U,X,bX,norb,nocc,nxyz)
      implicit none
C     Arguments
      integer norb,nocc,nxyz
      logical bX
      double precision U(norb,norb,nxyz),X(norb,norb,nxyz)
C     Variables and constants
      integer i,j,ixyz
      double precision ZERO,HALF
      parameter (ZERO=0.0D+00,HALF = 0.5D+00)
C     Body of the subprogram
      if(bX) then
         do ixyz = 1,nxyz
            do j=1,nocc
               do i=1,nocc
                  U(i,j,ixyz)=X(i,j,ixyz)*HALF
               end do
            end do
            do j=nocc+1,norb
               do i=nocc+1,norb
                  U(i,j,ixyz)=X(i,j,ixyz)*HALF
               end do
            end do
         end do
      else
         do ixyz = 1,nxyz
            do j=1,nocc
               do i=1,nocc
                  U(i,j,ixyz)=ZERO
               end do
            end do
            do j=nocc+1,norb
               do i=nocc+1,norb
                  U(i,j,ixyz)=ZERO
               end do
            end do
         end do
      end if
      end
C
C*MODULE TDXITR  *DECK MakeEStep1
      subroutine MakeEStep1(E,G,norb,nocc,nxyz)
      implicit none
C     Arguments
      integer norb,nocc,nxyz
      double precision E(norb,norb,nxyz)
      double precision G(norb,norb,nxyz)
C     Variables and constants
      integer ixyz,i,j
      double precision ZERO
      parameter (ZERO=0.0D+00)
C     Body of the subprogram
      do ixyz = 1,nxyz
         do j = 1,nocc
            do i = 1,nocc
               E(i,j,ixyz) = G(i,j,ixyz)
            end do
            do i = nocc+1,norb
               E(i,j,ixyz) = ZERO
            end do
         end do
         do j = nocc+1,norb
            do i = 1,nocc
               E(i,j,ixyz) = ZERO
            end do
            do i = nocc+1,norb
               E(i,j,ixyz) = G(i,j,ixyz)
            end do
         end do
      end do
      end
C
C*MODULE TDXITR  *DECK MakeEStep2
      subroutine MakeEStep2(E,X,bX,E0,omega,norb,nocc,nxyz)
      implicit none
C     Arguments
      integer norb,nocc,nxyz
      logical bX
      double precision omega
      double precision E(norb,norb,nxyz)
      double precision X(norb,norb,nxyz)
      double precision E0(norb)
C     Variables and constants
      integer i,j,ixyz
      double precision dummy
      double precision HALF
      parameter (HALF=0.5D+00)
C     Body of the subprogram
      if(bX) then
         do ixyz = 1,nxyz
            do j=1,nocc
               do i=1,nocc
                  dummy=X(i,j,ixyz)*(E0(i)-E0(j)+omega)
                  E(i,j,ixyz)=E(i,j,ixyz)+HALF*dummy
               end do
            end do
            do j=nocc+1,norb
               do i=nocc+1,norb
                  dummy=X(i,j,ixyz)*(E0(i)-E0(j)+omega)
                  E(i,j,ixyz)=E(i,j,ixyz)+HALF*dummy
               end do
            end do
         end do
      end if
      end
C
C*MODULE TDXITR  *DECK MakeDDStep0
      subroutine MakeDDStep0
      implicit none
      call FreeFQX
      end
C
C*MODULE TDXITR  *DECK MakeDDStep1
      subroutine MakeDDStep1(DD,C0,E0,omega,F,bF,norb,nocc,nvir,nxyz)
      implicit none
C     Arguments
      integer norb,nocc,nvir,nxyz
      logical bF
      double precision omega
      double precision F(norb,norb,nxyz)
      double precision DD(norb,norb,nxyz)
      double precision C0(norb,norb)
      double precision E0(norb)
C     Used common block
      DOUBLE PRECISION XX
      COMMON /FMCOM / XX(1)
C     Variables
      integer ialloc
      double precision conv(1)
C     Body of the subprogram
C
      conv(1) = 1.0d+00  ! suppresses FTNCHEK warning
C
      if(bF) then
         call SaveFQX(F,.true.,F,.false.,F,.false.,norb,nxyz)
         call NEWMEM(ialloc,2*nvir*nocc*nxyz)
         call cycMakeCdFC(C0,F,xx(ialloc),conv,.false.,
     &                    norb,nocc,nvir,nxyz)
         call cycMakeU(E0,omega,xx(ialloc),conv,.false.,
     &                 norb,nocc,nvir,nxyz)
C         call SetU1(xx(ialloc),norb,nocc,nvir,nxyz)
         call cycMakeD(DD,C0,xx(ialloc),norb,nocc,nvir,nxyz)
         call DELETEMEM(ialloc,2*nvir*nocc*nxyz)
C      else
C         call NEWMEM(ialloc,2*nvir*nocc*nxyz)
C         call VCLR(XX(ialloc),1,2*nvir*nocc*nxyz)
C         call SetU1(xx(ialloc),norb,nocc,nvir,nxyz)
C         call DELETEMEM(ialloc,2*nvir*nocc*nxyz)
      end if
      end
C
C*MODULE TDXITR  *DECK MakeDDStep2
      subroutine MakeDDStep2(DD,C0,E0,omega,Q,bQ,norb,nocc,nvir,nxyz)
      implicit none
C     Arguments
      integer norb,nocc,nvir,nxyz
      logical bQ
      double precision omega
      double precision Q(norb,norb,nxyz)
      double precision DD(norb*norb,nxyz)
      double precision C0(norb,norb)
      double precision E0(norb)
C     Used common block
      DOUBLE PRECISION XX
      COMMON /FMCOM / XX(1)
C     Variables
      integer ixyz,i,j,ij,ialloc,iwork
C     Body of the subprogram
      if(bQ) then
         call SaveFQX(Q,.false.,Q,.true.,Q,.false.,norb,nxyz)
         call NEWMEM(ialloc,nocc*nvir*2)
         call NEWMEM(iwork,norb*norb)
         do ixyz=1,nxyz
            ij = 0
            do j=nocc+1,norb
               do i=1,nocc
                  XX(ialloc+ij) = Q(i,j,ixyz)
                  ij=ij+1
               end do
            end do
            do i=nocc+1,norb
               do j=1,nocc
                  XX(ialloc+ij) = Q(i,j,ixyz)
                  ij=ij+1
               end do
            end do
            call cycMakeUx(E0,omega,XX(ialloc),norb,nocc,nvir)
C            call SetU2(xx(ialloc),ixyz,norb,nocc,nvir,nxyz)
            call cycMakeDx(DD(1,ixyz),XX(iwork),C0,XX(ialloc),
     &                     norb,nocc,nvir)
         end do
         call DELETEMEM(iwork,norb*norb)
         call DELETEMEM(ialloc,nocc*nvir*2)
      end if
      end
C
C*MODULE TDXITR  *DECK MakeDDStep3
      subroutine MakeDDStep3(DD,C0,X,bX,norb,nocc,nxyz)
      implicit none
C     Arguments
      integer norb,nocc,nxyz
      logical bX
      double precision X(norb*norb,nxyz)
      double precision DD(norb*norb,nxyz)
      double precision C0(norb*norb)
C     Used common block
      DOUBLE PRECISION XX
      COMMON /FMCOM / XX(1)
C     Variables and constants
      integer ixyz,ialloc
      double precision ZERO,ONE,TWO
      parameter (ZERO=0.0D+00,ONE=1.0D+00,TWO=2.0D+00)
      character*1 sT,sN
      parameter (sT='T',sN='N')
C     Body of the subprogram
      if(bX) then
         call SaveFQX(X,.false.,X,.false.,X,.true.,norb,nxyz)
         call NEWMEM(ialloc,norb*norb)
         do ixyz=1,nxyz
            call DGEMM(sN,sN,norb,norb,nocc,
     &                 ONE,C0(1),norb,X(1,ixyz),norb,
     &                 ZERO,XX(ialloc),norb)
            call DGEMM(sN,sT,norb,norb,norb,
     &                 TWO,XX(ialloc),norb,C0,norb,
     &                 ONE,DD(1,ixyz),norb)
         end do
         call DELETEMEM(ialloc,norb*norb)
      end if
      end
C
C*MODULE TDXITR  *DECK cycMakeCdFC
      subroutine cycMakeCdFC(C0,F,G,conv,bc,norb,nocc,nvir,nxyz)
      implicit none
C     Arguments
      integer norb,nocc,nvir,nxyz
      logical bc
      double precision C0(norb,norb)
      double precision F(norb*norb,nxyz)
      double precision G(nocc*nvir,2,nxyz)
      double precision conv(nxyz)
C     Used common block
      DOUBLE PRECISION XX
      COMMON /FMCOM / XX(1)
C     Variables and constants
      logical test
      integer ixyz,ialloc
      double precision ZERO,ONE
      parameter (ZERO=0.0D+00,ONE=1.0D+00)
      character*1 sT,sN
      parameter (sT='T',sN='N')
C     Body of the subprogram
      call NEWMEM(ialloc,norb*norb)
      do ixyz = 1,nxyz
C-----        how not to touch CONV during initialization pass
C-----   if((.not.bc).or.(conv(ixyz).ne.ZERO)) then
         test = .not.bc
         if(.not.test) test = conv(ixyz).ne.ZERO
         if(test) then
            call DGEMM(sT,sN,nocc,norb,norb,
     &                 ONE,C0,norb,F(1,ixyz),norb,
     &                 ZERO,XX(ialloc),norb)
            call DGEMM(sN,sN,nocc,nvir,norb,
     &                 ONE,XX(ialloc),norb,C0(1,nocc+1),norb,
     &                 ZERO,G(1,1,ixyz),nocc)
C
            call DGEMM(sT,sT,nocc,norb,norb,
     &                 ONE,C0,norb,F(1,ixyz),norb,
     &                 ZERO,XX(ialloc),norb)
            call DGEMM(sN,sN,nocc,nvir,norb,
     &                 ONE,XX(ialloc),norb,C0(1,nocc+1),norb,
     &                 ZERO,G(1,2,ixyz),nocc)
         end if
      end do
      call DELETEMEM(ialloc,norb*norb)
      end
C
C*MODULE TDXITR  *DECK cycMakeU
      subroutine cycMakeU(E0,omega,U,conv,bc,norb,nocc,nvir,nxyz)
      implicit none
C     Arguments
      integer norb,nocc,nvir,nxyz
      logical bc
      double precision omega
      double precision E0(norb)
      double precision U(nocc*nvir*2,nxyz)
      double precision conv(nxyz)
C     Variables and constants
      logical test
      integer ixyz
      double precision ZERO
      parameter (ZERO=0.0D+00)
C     Body of the subprogram
      do ixyz=1,nxyz
C-----        how not to touch CONV during initialization pass
C-----   if((.not.bc).or.(conv(ixyz).ne.ZERO)) then
         test = .not.bc
         if(.not.test) test = conv(ixyz).ne.ZERO
         if(test) then
            call cycMakeUx(E0,omega,U(1,ixyz),norb,nocc,nvir)
         end if
      end do
      end
C
C*MODULE TDXITR  *DECK cycMakeUx
      subroutine cycMakeUx(E0,omega,U,norb,nocc,nvir)
      implicit none
C     Arguments
      integer norb,nocc,nvir
      double precision omega
      double precision E0(norb)
      double precision U(nocc,nvir,2)
C     Variables and constants
      integer i,j
C     Body of the subprogram
      do j=1,nvir
         do i=1,nocc
            U(i,j,1)=U(i,j,1)/(E0(nocc+j)-E0(i)-omega)
         end do
      end do
      do j=1,nvir
         do i=1,nocc
            U(i,j,2)=U(i,j,2)/(E0(i)-E0(nocc+j)-omega)
         end do
      end do
      end
C
C*MODULE TDXITR  *DECK cycMakeD
      subroutine cycMakeD(D,C0,U,norb,nocc,nvir,nxyz)
      implicit none
C     Arguments
      integer norb,nocc,nvir,nxyz
      double precision U(nocc*nvir*2,nxyz)
      double precision D(norb*norb,nxyz)
      double precision C0(norb,norb)
C     Used common block
      DOUBLE PRECISION XX
      COMMON /FMCOM / XX(1)
C     Variables and constants
      integer ixyz,ialloc
C     Body of the subprogram
      call NEWMEM(ialloc,norb*norb)
      do ixyz=1,nxyz
         call cycMakeDx(D(1,ixyz),XX(ialloc),C0,U(1,ixyz),
     &                  norb,nocc,nvir)
      end do
      call DELETEMEM(ialloc,norb*norb)
      end
C
C*MODULE TDXITR  *DECK cycMakeDx
      subroutine cycMakeDx(D,work,C0,U,norb,nocc,nvir)
      implicit none
C     Arguments
      integer norb,nocc,nvir
      double precision U(nocc*nvir,2)
      double precision D(norb*norb),work(norb*norb)
      double precision C0(norb,norb)
C     Variables and constants
      double precision ZERO,ONE,TWO,MINUSTWO
      parameter (ZERO=0.0D+00,ONE=1.0D+00,TWO=2.0D+00,MINUSTWO=-2.0D+00)
      character*1 sT,sN
      parameter (sT='T',sN='N')
C     Body of the subprogram
      call DGEMM(sN,sT,norb,nocc,nvir,
     &           ONE,C0(1,nocc+1),norb,U(1,2),nocc,ZERO,work,norb)
      call DGEMM(sN,sT,norb,norb,nocc,
     &           TWO,work,norb,C0,norb,ONE,D,norb)
      call DGEMM(sN,sN,norb,nvir,nocc,
     &           ONE,C0,norb,U(1,1),nocc,ZERO,work,norb)
      call DGEMM(sN,sT,norb,norb,nvir,
     &           MINUSTWO,work,norb,C0(1,nocc+1),norb,ONE,D,norb)
      end
C
C*MODULE TDXITR  *DECK cycCopyU
      subroutine cycCopyU(U,Unew,conv,nocc,nvir,nxyz)
      implicit none
C     Arguments
      integer nocc,nvir,nxyz
      double precision U(nocc,nvir,2,nxyz)
      double precision Unew(nocc,nvir,2,nxyz)
      double precision conv(nxyz)
C     Variables and constants
      integer i,j,ixyz
      double precision temp,dum,diff
      double precision ZERO,SIGMA
      parameter (ZERO=0.0D+00,SIGMA=-1.0D+00)
C     Called subprogram
      double precision ABS
C     intrinsic ABS
C     Body of the subprogram
      do ixyz = 1,nxyz
         if(conv(ixyz).ne.ZERO) then
            diff = 0.0d+00
            do j = 1,nvir
               do i = 1,nocc
                  temp = Unew(i,j,1,ixyz)-U(i,j,1,ixyz)
                  dum = ABS(temp)
                  if(dum.gt.diff) diff=dum
                  U(i,j,1,ixyz) = Unew(i,j,1,ixyz)
               end do
            end do
            do i = 1,nvir
               do j = 1,nocc
                  temp = Unew(j,i,2,ixyz)-U(j,i,2,ixyz)
                  dum = ABS(temp)
                  if(dum.gt.diff) diff=dum
                  U(j,i,2,ixyz) = Unew(j,i,2,ixyz)
               end do
            end do
            if(diff.eq.ZERO) diff = SIGMA
            conv(ixyz) = diff
         end if
      end do
      end
C
C*MODULE TDXITR  *DECK cycInitU
      subroutine cycInitU(U,X,bX,norb,nocc,nvir,nxyz)
      implicit none
C     Arguments
      integer norb,nocc,nvir,nxyz
      logical bX
      double precision U(nocc,nvir,2,nxyz)
      double precision X(norb,norb,nxyz)
C     Variables
      integer ixyz,i,j
      double precision HALF
      parameter (HALF=0.5D+00)
C     Body of the subprogram
      if(bX) then
C         call PrintMM(X,norb,nxyz,'X')
         do ixyz=1,nxyz
            do j=1,nvir
               do i=1,nocc
                  U(i,j,1,ixyz) = X(i,nocc+j,ixyz)*HALF
               end do
            end do
            do j=1,nvir
               do i=1,nocc
                  U(i,j,2,ixyz) = X(nocc+j,i,ixyz)*HALF
               end do
            end do
         end do
      else
         call vclr(U,1,nocc*nvir*2*nxyz)
      end if
C      call SetU3(U,norb,nocc,nvir,nxyz)
      end
C
C*MODULE TDXITR  *DECK TDHFCycle
      subroutine TDHFCycle(subGetDJK,subPrintConvergedParam,
     &                     C0,E0,omega,D,DD,norb,nocc,nvir,nxyz)
      implicit none
C     Arguments
      external subGetDJK,subPrintConvergedParam
      double precision omega
      integer norb,nocc,nvir,nxyz
      double precision C0(norb*norb)
      double precision E0(norb)
      double precision D(norb*norb,nxyz)
      double precision DD(norb*norb,nxyz)
C     Used common blocks
      integer IR, IW, IP, IS,IPK,IDAF,NAV,IODA
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      character*8 method
      integer itermax
      double precision convergence
      common/OLI_CYCLE/convergence,method,itermax
C     Variables and constants
      external NoAccel,DIISAccel
      character*8 sNONE,sDIIS
      parameter (sNONE='NONE',sDIIS='DIIS')
C     Body of the subprogram
C     Fortran 90 solution
C      select case(method)
C      case (sNONE)
C     Ugly fortran 77 solution
      if(method.eq.sNONE) then
         call TDHFCycleX(convergence,itermax,sNONE,
     &                   subGetDJK,NoAccel,subPrintConvergedParam,
     &                   C0,E0,omega,D,DD,norb,nocc,nvir,nxyz)
C     Fortran 90 solution
C      case (sDIIS)
C     Ugly fortran 77 solution
      else if(method.eq.sDIIS) then
         call TDHFCycleX(convergence,itermax,sDIIS,
     &                   subGetDJK,DIISAccel,subPrintConvergedParam,
     &                   C0,E0,omega,D,DD,norb,nocc,nvir,nxyz)
C     Fortran 90 solution
C      case default
C     Ugly fortran 77 solution
      else
         write(iw,*) 'Wrong KEYWORD [',method,'] in TDHFCycle'
C     Fortran 90 solution
C      end select
C     Ugly fortran 77 solution
      end if
      end
C
C
C*MODULE TDXITR  *DECK TDHFCycleX
      subroutine TDHFCycleX(atol,maxit,method,
     &                      subGetDJK,subAccel,subPrintConvergedParam,
     &                      C0,E0,omega,D,DD,norb,nocc,nvir,nxyz)
      implicit none
C     Arguments
      external subGetDJK,subAccel,subPrintConvergedParam
      character*(*) method
      double precision atol,omega
      integer maxit,norb,nocc,nvir,nxyz
      double precision C0(norb*norb)
      double precision E0(norb)
      double precision D(norb*norb,nxyz)
      double precision DD(norb*norb,nxyz)
C     Used common blocks
      integer IR, IW, IP, IS,IPK,IDAF,NAV,IODA
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      DOUBLE PRECISION XX
      COMMON /FMCOM / XX(1)
C      double precision SOGTOL,ETHRSH
C      integer MAXDII,IRAF
C      COMMON /DIISSO/ SOGTOL,ETHRSH,MAXDII,IRAF
C     Variables and constants
      logical bF,bQ,bX
      logical converged
      integer ixyz,iresponses,icycle,nrot
      integer lUnew,lU,lF,lCONV
      double precision maxerr
      double precision ZERO,ONE
      parameter (ZERO=0.0D+00,ONE=1.0D+00)
      character*8 sCREATE,sITER,sDELETE
      parameter (sCREATE="CREATE",sITER="ITER",sDELETE="DELETE")
C     Called subprogram
      logical LoadFQX
C     Body of the subprogram
      icycle    = 0
      converged = .false.
      nrot      = 2*nocc*nvir
      write(iw,3000) maxit,atol,method,nxyz,nrot,omega
C     Allocation of memory for the Fock Matrices
      call NEWMEM(lF,norb*norb*nxyz)
      call NEWMEM(lU,nrot*nxyz)
      call NEWMEM(lUnew,nrot*nxyz)
      call VCLR(XX(lU),1,nrot*nxyz)
C     Allocation and initialization of memory for the Convergence List
      call NEWMEM(lCONV,nxyz)
      do ixyz=1,nxyz
         xx(lCONV+ixyz-1) = ONE
      end do
C     Initialize Acceleration stuff
      call  subAccel(sCREATE,XX(lUnew),XX(lU),XX(lCONV),maxit,nrot,
     &               nocc,nvir,nxyz)
C
C     ---------
C     prepare D
C     ---------
      call CheckFQX(bF,bQ,bX)
      if(bX) bX = LoadFQX(xx(lF),.false.,xx(lF),.false.,xx(lF),.true.,
     &                    norb,nxyz)
      call cycInitU(XX(lU),XX(lF),bX,norb,nocc,nvir,nxyz)
      call DCOPY(norb*norb*nxyz,DD,1,D,1)
      call cycMakeD(D,C0,XX(lU),norb,nocc,nvir,nxyz)
C     ------
C     Step 1
C     ------
C     Fortran 90 solution
C      do while((.not.converged).and.(icycle.lt.maxit))
C     Ugly Fortran 77 solution
 1    continue
      icycle=icycle+1
      call subGetDJK(D,XX(lF),norb,nxyz)
C      call PrintMM(D,norb,nxyz,'OLIX:D')
C      call PrintMM(XX(lF),norb,nxyz,'OLIX:DJK')
C     ------
C     Step 2
C     ------
      call cycMakeCdFC(C0,XX(lF),XX(lUnew),XX(lCONV),.true.,
     &                 norb,nocc,nvir,nxyz)
C     ------
C     Step 3
C     ------
      call cycMakeU(E0,omega,XX(lUnew),XX(lCONV),.true.,
     &              norb,nocc,nvir,nxyz)
      call subAccel(sITER,XX(lUnew),XX(lU),XX(lCONV),maxit,nrot,
     &              nocc,nvir,nxyz)
      iresponses=0
      maxerr=ZERO
      do ixyz=1,nxyz
         if(XX(lCONV+ixyz-1).gt.maxerr) maxerr=XX(lCONV+ixyz-1)
         if(XX(lCONV+ixyz-1).gt.atol) then
            iresponses=iresponses+1
         else
            if(XX(lCONV+ixyz-1).ne.ZERO) then
               if(XX(lCONV+ixyz-1).le.ZERO) XX(lCONV+ixyz-1)=ZERO
               call subPrintConvergedParam(iw,ixyz,XX(lCONV+ixyz-1))
            end if
            XX(lCONV+ixyz-1)=ZERO
         end if
      end do
      converged=(iresponses.eq.0)
C     ------
C     Step 4
C     ------
      call DCOPY(norb*norb*nxyz,DD,1,D,1)
      call cycMakeD(D,C0,XX(lU),norb,nocc,nvir,nxyz)
C      call PrintMM(XX(lU),norb,nxyz,'OLIX:U')
C     ------
C     Step 5
C     ------
C      call PrintU(XX(lU),nocc,nvir,nxyz)
      write(iw,1001) icycle,iresponses,maxerr
C     Fortran 90 solution
C      end do
C     Ugly Fortran 77 solution
      if((.not.converged).and.(icycle.lt.maxit)) goto 1
C     ---------
C     Converged
C     ---------
C     Copy D[2J-K] into DD
      call DCOPY(norb*norb*nxyz,xx(lf),1,DD,1)
C     Clean up Acceleration stuff
      call subAccel(sDELETE,XX(lUnew),XX(lU),XX(lCONV),maxit,nrot,
     &              nocc,nvir,nxyz)
C     Deallocation of memory for the Convergence List
      call DELETEMEM(lCONV,nxyz)
C     Deallocation of memory for the Fock Matrices
      call DELETEMEM(lUnew,nrot*nxyz)
      call DELETEMEM(lU,nrot*nxyz)
      call DELETEMEM(lF,norb*norb*nxyz)
      return
C     Format
 1001 FORMAT('OLIX ITER',I3,':',I4,' RESPONSES NOT CONVERGED, ',
     &       'MAXIMUM RESPONSE ERROR=',E13.6)
 3000 FORMAT(/,2X,'OLIX:TDHFCycleX',
     &       /,
     &       /,4X,'MAXITER=',I4,
     &       /,4X,'TOLERANCE=',E13.6,
     &       /,4X,'ACCELERATION METHOD=',A,
     &       /,4X,I4,' parameters to solve of size ',I10,
     &       /,4X,'OMEGA=',E13.6,
     &       /)
      end
C
C*MODULE TDXITR  *DECK NoAccel
      subroutine NoAccel(sStep,U,Uold,conv,maxit,nrot,
     &                   nocc,nvir,nxyz)
      implicit none
C     Arguments
      character*8 sStep
      integer maxit,nrot,nocc,nvir,nxyz
      double precision U(nrot,nxyz),Uold(nrot,nxyz),conv(nxyz)
C     Used common blocks
      integer IR, IW, IP, IS,IPK,IDAF,NAV,IODA
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
C     Variables and constants
      character*8 sCREATE,sITER,sDELETE
      parameter (sCREATE="CREATE",sITER="ITER",sDELETE="DELETE")
C     Body of the subprogram
C
C     Fortran 90 solution
C      select case(sStep)
C      case (sITER)
C     Ugly fortran 77 solution
      if(sStep.eq.sITER) then
         call cycCopyU(Uold,U,conv,nocc,nvir,nxyz)
C     Fortran 90 solution
C      case (sCREATE)
C     Ugly fortran 77 solution
      else if(sStep.eq.sCREATE) then
C     Fortran 90 solution
C      case (sDELETE)
C     Ugly fortran 77 solution
      else if(sStep.eq.sDELETE) then
C     Fortran 90 solution
C      case default
C     Ugly fortran 77 solution
      else
         write(iw,*) 'Wrong KEYWORD in NoAccel [',sSTEP,']',maxit
         call abrtx("NOACCEL")
C     Fortran 90 solution
C      end select
C     Ugly fortran 77 solution
      end if
      end
C
C--c*MODULE TDXITR  *DECK PrintConvergedParam
C--      logical function PrintConvergedParam(iw,ixyz,conv)
C--      implicit none
C--c     Arguments
C--      integer iw,ixyz
C--      double precision conv
C--c     Body of the subprogram
C--      write(iw,1001) ixyz,conv
C--      PrintConvergedParam = .true.
C--c     Format
C-- 1001 FORMAT('  Parameter ',I4,' converged within ',E13.6)
C--      end
C
C
C*MODULE TDXITR  *DECK GetDJK
      subroutine GetDJK(D,DJK,norb,nxyz)
      implicit none
C     Arguments
      integer norb,nxyz
      double precision D(norb*norb,nxyz)
      double precision DJK(norb*norb,nxyz)
C     Used common blocks
      DOUBLE PRECISION XX
      COMMON /FMCOM / XX(1)
      integer NPRINT,ITOL,ICUT,NORMF,NORMP,NOPK
      COMMON /OUTPUT/ NPRINT,ITOL,ICUT,NORMF,NORMP,NOPK
      integer MXGTOT,MXSH,MXGSH,MXAO
      PARAMETER (MXGTOT=20000,MXSH=5000,MXAO=8192,MXGSH=30)
      double precision EX,CS,CP,CD,CF,CG,ch,ci
      integer KSTART,KATOM,KTYPE,KNG,KLOC,KMIN,KMAX
      integer NSHELL
      COMMON /NSHEL / EX(MXGTOT),CS(MXGTOT),CP(MXGTOT),CD(MXGTOT),
     &                CF(MXGTOT),CG(MXGTOT),ch(mxgtot),ci(mxgtot),
     &                KSTART(MXSH),KATOM(MXSH),KTYPE(MXSH),KNG(MXSH),
     &                KLOC(MXSH),KMIN(MXSH),KMAX(MXSH),NSHELL
      logical DIRSCF,FDIFF
      COMMON /OPTSCF/ DIRSCF,FDIFF
      integer NINTMX,NHEX,NTUPL,inttyp,igrdtyp
      logical PACK2E
      COMMON /INTFIL/ NINTMX,NHEX,NTUPL,PACK2E,INTtyp,igrdtyp
      integer IA
      COMMON /IJPAIR/ IA(MXAO)
C     Variables and constants
      logical SCHWRZ
      integer NSCHWZ,NINT,NSQR,L2,NSH2,LMAX,NANGM,MAXG
      integer iDUMMY
      double precision DUMMY
      integer LCLINT,LGHOND,LDDIJ,LDSH
      integer LBUF,LIBUF
      DOUBLE PRECISION RHF
      DATA RHF/8HRHF     /
C     Body of the subprogram
C
      NSCHWZ = 0
      SCHWRZ = .FALSE.
      NINT = 0
      NSQR = NORB*NORB
      L2   = (NORB*NORB + NORB)/2
      NSH2 = (NSHELL*NSHELL + NSHELL)/2
C
      CALL BASCHK(LMAX)
                    NANGM =  4
      IF(LMAX.EQ.2) NANGM =  6
      IF(LMAX.EQ.3) NANGM = 10
      IF(LMAX.EQ.4) NANGM = 15
      IF(LMAX.EQ.5) NANGM = 21
      IF(LMAX.EQ.6) NANGM = 28
      MAXG = NANGM**4
C
      call VCLR(DJK,1,nsqr*nxyz)
C
      IF (DIRSCF) THEN
         call NEWMEM(LCLINT,NSH2)
         call NEWMEM(LGHOND,MAXG)
         call NEWMEM(LDDIJ,49*MXGSH*MXGSH)
         call NEWMEM(LDSH,NSH2)
         call TWOEI(RHF,.FALSE.,.TRUE.,.FALSE.,.false.,
     &        inttyp,SCHWRZ,NINT,NSCHWZ,norb,L2,
     &        DUMMY,DUMMY,IDUMMY,1,
     &        XX(LCLINT),NSH2,XX(LGHOND),MAXG,XX(LDDIJ),
     &        IA,DUMMY,DUMMY,DUMMY,DUMMY,XX(LDSH),
     &        D,DJK,nxyz)
         call DELETEMEM(LDSH,NSH2)
         call DELETEMEM(LDDIJ,49*MXGSH*MXGSH)
         call DELETEMEM(LGHOND,MAXG)
         call DELETEMEM(LCLINT,NSH2)
      ELSE
         call NEWMEM(LBUF,NINTMX)
         call NEWMEM(LIBUF,NINTMX)
         CALL F2NONS(D,DJK,XX(LBUF),XX(LIBUF),NINTMX,NOPK,nxyz)
         call DELETEMEM(LIBUF,NINTMX)
         call DELETEMEM(LBUF,NINTMX)
      END IF
      end
C
C*MODULE TDXITR  *DECK DIISAccel
      subroutine DIISAccel(sStep,U,Uold,conv,maxit,nrot,
     &                     nocc,nvir,nxyz)
      implicit none
C     Arguments
      character*8 sStep
      integer maxit,nrot,nocc,nvir,nxyz
      double precision U(nrot,nxyz),Uold(nrot,nxyz),conv(nxyz)
C     Used common blocks
      integer IR, IW, IP, IS,IPK,IDAF,NAV,IODA
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      DOUBLE PRECISION XX
      COMMON /FMCOM / XX(1)
C     Variables and constants
      integer itdiis
      double precision errdiis
      integer lERR,lISIZES,lIODIIS,lA,lX,lB
      save itdiis,lERR,lISIZES,lIODIIS,lA,lX,lB
      double precision ZERO
      parameter (ZERO=0.0D+00)
      character*8 sCREATE,sITER,sDELETE
      parameter (sCREATE="CREATE",sITER="ITER",sDELETE="DELETE")
      integer IRAF,MAXIO
      parameter (IRAF=20,MAXIO=20)
C
C     Body of the subprogram
C
C     Fortran 90 solution
C      select case(sStep)
C      case (sCREATE)
C     Ugly fortran 77 solution
      if(sStep.eq.sCREATE) then
C     ----------------------------
C     Initialization of DIIS stuff
C     ----------------------------
         errdiis = ZERO
         itdiis = 1
         call NEWMEM(lERR,nrot*nxyz)
         call NEWMEM(lISIZES,nxyz)
         call NEWMEM(lIODIIS,MAXIO*2)
         call NEWMEM(lA,MAXIO*MAXIO)
         call NEWMEM(lX,MAXIO*nxyz)
         call NEWMEM(lB,(maxit*maxit+maxit)/2*nxyz)
         CALL RAOPEN(IRAF,XX(lIODIIS),0,2*MAXIO,nrot*nxyz,-5)
         CALL VCLR(XX(lISIZES),1,nxyz)
         CALL VCLR(XX(lB),1,(maxit*maxit+maxit)/2*nxyz)
C     Fortran 90 solution
C      case (sITER)
C     Ugly fortran 77 solution
      else if(sStep.eq.sITER) then
C     -------
C     Iterate
C     -------
         call VSUB(Uold,1,U,1,XX(lERR),1,nrot*nxyz)
         call OLIXDIIS(IRAF,IW,itdiis,U,XX(lERR),conv,
     &        XX(lA),XX(lX),XX(lB),XX(lIODIIS),XX(lISIZES),
     &        maxit,MAXIO,errdiis,nrot,nxyz)
         call cycCopyU(Uold,U,conv,nocc,nvir,nxyz)
C     Fortran 90 solution
C      case (sDELETE)
C     Ugly fortran 77 solution
      else if(sStep.eq.sDELETE) then
C     ----------------------
C     Cleaning up DIIS stuff
C     ----------------------
         call RACLOS(IRAF,'DELETE')
         call DELETEMEM(lB,(maxit*maxit+maxit)/2*nxyz)
         call DELETEMEM(lX,MAXIO*nxyz)
         call DELETEMEM(lA,MAXIO*MAXIO)
         call DELETEMEM(lIODIIS,MAXIO*2)
         call DELETEMEM(lISIZES,nxyz)
         call DELETEMEM(lERR,nrot*nxyz)
C     Fortran 90 solution
C      case default
C     Ugly fortran 77 solution
      else
         write(iw,*) 'Wrong KEYWORD [',sSTEP,'] in DIISAccel'
         call abrtx("DIISACCEL")
C     Fortran 90 solution
C      end select
C     Ugly fortran 77 solution
      end if
      end
C
C*MODULE TDXITR  *DECK OLIXDIIS
      SUBROUTINE OLIXDIIS(IRAF,IW,ITDIIS,VEC,ERR,conv,A,X,B,
     &                    IODIIS,ISIZES,MAXIT,MAXIO,ERDIIS,
     &                    NROT,NXYZ)
      implicit none
C     Arguments
      integer IRAF,IW,ITDIIS,MAXIT,MAXIO,NROT,NXYZ
      double precision ERDIIS
      double precision VEC(nrot,nxyz)
      double precision ERR(nrot,nxyz)
      double precision conv(nxyz)
      double precision A(MAXIO*MAXIO)
      double precision X(MAXIO,nxyz)
      double precision B((maxit*maxit+maxit)/2,nxyz)
      integer ISIZES(NXYZ),IODIIS(MAXIO*2)
C     Uas common blocks
      integer ME,MASTER,NPROC,IBTYP,IPTIM
      logical GOPARR,DSKWRK,MASWRK
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
C     variables and constants
      integer LRAFERR,LRAFVEC,LERR,LVEC,nax
      integer IX,IJ,NA,NAMAX,IERR,ixyz,iBsize
      double precision BIJ,BJJ
      logical GPSAVE,DSKSAV
      double precision ZERO,ONE
      PARAMETER (ZERO=0.0D+00,ONE=1.0D+00)
C     Called subprograms
C     intrinsic ABS,MOD,MIN
      integer MIN,MOD
      double precision DDOT
C     Body of the subprogram
C
C     DIIS DATA IS SAVED ON DIRECT FILE -IRAF- IN CYCLIC ORDER,
C                1,2,3,...,MAXIO = LAST MAXIT ERROR MATRICES
C        MAXIO + 1,2,3,...,MAXIO = LAST MAXIT CPHF SOLUTIONS
      LRAFERR = 0
      LRAFVEC= MAXIO
      DSKSAV = DSKWRK
      DSKWRK = .FALSE.
C     -----------------
C     COMPUTE MAX ERROR
C     -----------------
      erdiis=ZERO
      do ixyz=1,nxyz
         if(conv(ixyz).gt.erdiis) erdiis = conv(ixyz)
      end do
C     ------------------------------
C     SAVE THE CURRENT CPHF SOLUTION
C     ------------------------------
      LVEC = LRAFVEC + MOD(ITDIIS-1,MAXIO) + 1
      CALL RAWRIT(IRAF,IODIIS,VEC,NROT*NXYZ,LVEC,0)
C     -----------------------------
C     SAVE THE CURRENT ERROR MATRIX
C     -----------------------------
      LERR = LRAFERR + MOD(ITDIIS-1,MAXIO) + 1
      CALL RAWRIT(IRAF,IODIIS,ERR,NROT*NXYZ,LERR,0)
C     ----------------------------------------------
C     SETUP INVERSION SPACE SIZE FOR EACH PARAMETERS
C     ----------------------------------------------
      ITDIIS=ITDIIS+1
      NAMAX=MIN(ITDIIS,MAXIO)
      iBsize = (ITDIIS*ITDIIS+ITDIIS)/2
      do ixyz=1,nxyz
         if(conv(ixyz).eq.ZERO) then
            isizes(ixyz) = 0
         else
            isizes(ixyz) = NAMAX
         endif
      end do
      IF (ITDIIS.EQ.2.AND.MASWRK) WRITE(IW,3000)
C        WHEN RUNNING IN PARALLEL, THE O(N**2) COMMUNICATIONS
C        COSTS IN -RAREAD- BELOW ARE THE SAME MAGNITUDE AS
C        THE COMPUTATION IN -TRAREC- AND -DAXPY-.
C        IT SEEMS WISEST JUST TO LET THE MASTER DO THIS I/O,
C        THEN BROADCAST THE INTERPOLATION RESULTS AT THE END.
      GPSAVE = GOPARR
      GOPARR = .FALSE.
C     -------------------------------------------------
C     FORM LAST ROW (ROW ITDIIS) OF TRIANGULAR MATRIX B
C     -------------------------------------------------
C     ---------------------------------------------------------
C     FIRST DOING THE DIAGONAL AND THE FIRST ELEMENT OF THE ROW
C     ---------------------------------------------------------
      do 101 ixyz=1,nxyz
         if (ISIZES(ixyz).eq.0) go to 101
         BJJ=DDOT(nrot,ERR(1,ixyz),1,ERR(1,ixyz),1)
         B(iBsize-ITDIIS+1,ixyz)= -ONE
         B(iBsize         ,ixyz)= BJJ
  101 continue
      if(ITDIIS.EQ.2) then
         GOPARR = GPSAVE
         DSKWRK = DSKSAV
         RETURN
      end if
C     --------------------------
C     THEN THE REST OF THE BIJ'S
C     --------------------------
C     NOTE THAT WE ONLY COMPUTE THE NAMAX-2 VALUES NEXT TO THE
C     DIAGONAL, SO THAT THE OTHER ELEMENTS (EXCEPT THE FIRST
C     COLUMN) ARE RANDOM VALUES.
      IJ=iBsize-NAMAX+1
      do IX=2,NAMAX-1
         IJ=IJ+1
         LERR = LRAFERR + MOD(ITDIIS-NAMAX+IX-2,MAXIO) + 1
         CALL RAREAD(IRAF,IODIIS,VEC,nrot*nxyz,LERR,0)
         do 102 ixyz=1,nxyz
            if (ISIZES(ixyz).eq.0) go to 102
            BIJ=DDOT(nrot,ERR(1,ixyz),1,VEC(1,ixyz),1)
            B(IJ,IXYZ)=BIJ
  102    continue
      end do
C     -------------------------------------------------
C     SET UP AND SOLVE THE SET OF DIIS LINEAR EQUATIONS
C     -------------------------------------------------
      do 1011 ixyz=1,nxyz
         if(ISIZES(ixyz).EQ.0) go to 1011
         do NA=ISIZES(ixyz),1,-1
            CALL DIISEQ(A,NA,B(1,ixyz),iBsize,ITDIIS)
C            A(1) = B(1,ixyz)
C            IX=2
C            do j=ITDIIS-NA+2,ITDIIS
C               IJ    = (j*(j-1))/2
C               A(IX) = B(IJ+1,ixyz)
C               IX=IX+1
C               do i=ITDIIS-NA+2,j
C                  A(IX) = B(IJ+i,ixyz)
C                  IX=IX+1
C               end do
C            end do
C
            CALL VCLR(X(1,ixyz),1,ISIZES(ixyz))
            X(1,ixyz)=-ONE
            IERR=0
            nax = NA
            CALL SLVLEQsvd(A,X(1,ixyz),nax,nax,IERR)
C            ipvt would be a work array of dimension NA
C            call DSPTRF('U',NA,A,IPVT,IERR)
            if(IERR.EQ.0) then
C               call DSPTRS('U',NA,1,A,IPVT,X,NA,IERR)
               ISIZES(ixyz)=NA
               go to 1001
            end if
            IF (MASWRK) WRITE(IW,3001) ixyz,NA
            if(NA.eq.1) then
               IF (MASWRK) WRITE(IW,3002) ixyz
               CALL ABRTX("OLIXDIIS")
            end if
         end do
 1001    continue
 1011 continue
C     -----------------------------------
C     FORM THE INTERPOLATED CPHF SOLUTION
C     -----------------------------------
      CALL VCLR(VEC,1,NROT*NXYZ)
      do IX=2,NAMAX
         LVEC = LRAFVEC + MOD(ITDIIS-NAMAX+IX-2,MAXIO) + 1
         CALL RAREAD(IRAF,IODIIS,ERR,NROT*NXYZ,LVEC,0)
         do 103 ixyz=1,nxyz
            if(ISIZES(IXYZ).EQ.0) go to 103
            if(X(IX,IXYZ).EQ.ZERO) go to 103
            CALL DAXPY(NROT,X(IX,IXYZ),ERR(1,ixyz),1,VEC(1,ixyz),1)
  103    continue
      end do
      GOPARR=GPSAVE
      DSKWRK = DSKSAV
      IF(GOPARR) THEN
         CALL DDI_BCAST(1005,'F',VEC,NXYZ*NROT,MASTER)
      END IF
      RETURN
C     Formats
 3000 FORMAT(' ..... INITIATING DIIS PROCEDURE .....')
 3001 FORMAT('TRUNCATING RESPONSE VECTOR ',I4,
     &       ' TO SOLVE DIIS EQUATIONS (SIZE=',I4,')')
 3002 FORMAT('ERROR: RESPONSE VECTOR ',I4,
     &       ' : DIIS EQUATIONS ARE SINGULAR')
      END
C
C*MODULE TDXITR  *DECK MakeAlpha
      subroutine MakeAlpha(w1)
      implicit none
C     Argument
      double precision w1
C     Used common blocks
      integer soli,ngeom,nocc,nvir,norb,nxyz,nmodes,nfld,ibfld,iefld
      common/OLI/soli,ngeom,nocc,nvir,norb,nxyz,nmodes,nfld,ibfld,iefld
      integer IR, IW, IP, IS,IPK,IDAF,NAV,IODA
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      DOUBLE PRECISION XX
      COMMON /FMCOM / XX(1)
C     Variables and constants
      double precision alpha(9)
      double precision w(2)
      integer lE0,lC0
      double precision dummy(1)
C     Called subprograms
      logical GetHf
C     Body of the subprogram
      call VCLR(alpha,1,9)
C
      w(1) = -w1
      w(2) =  w1
C
      dummy(1) = 0.0d+00
      if(.not.GetHf(dummy,.false.)) call MakeHf()
C
      call NEWMEM(lE0,norb)
      call NEWMEM(lC0,norb*norb)
      call GetEnergy(XX(lE0),norb)
      call GetLCAO(XX(lC0),norb)
C
      call CompAlpha(alpha,w1,XX(lE0),XX(lC0),
     &               norb,nocc,nvir,ibfld,nfld)
C
      call DELETEMEM(lC0,norb*norb)
      call DELETEMEM(lE0,norb)
C
      call SetAlpha(w1,alpha)
      call PrintAlphaX(iw,alpha,w,'Alpha tensor [in au]')
      CALL TIMIT(1)
      end
C
C*MODULE TDXITR  *DECK CompAlpha
      subroutine CompAlpha(alpha,w,E0,C0,
     &                     norb,nocc,nvir,ibfld,nfld)
      implicit none
C     Arguments
      integer norb,nocc,nvir,ibfld,nfld
      double precision w
      double precision alpha(9)
      double precision E0(norb),C0(norb,norb)
C     Used common blocks
      integer IR, IW, IP, IS,IPK,IDAF,NAV,IODA
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      DOUBLE PRECISION XX
      COMMON /FMCOM / XX(1)
C     Variables
      integer nparams
      integer lD,lDD,lHf,lG,lU,lE
      double precision dummyX(1),dummyQ(1)
C     Called subprogram
      logical GetHf
C     External subprogram
      external GetDJK,PrintAlphaConvergedParam
C     Body of the subprogram
      nparams = nfld
      write(iw,3000) -w,w
C     ----------
C     prepare DD
C     ----------
      call NEWMEM(lDD,norb*norb*nparams)
      call VCLR(XX(lDD),1,norb*norb*nparams)
      call MakeDDStep0
      call NEWMEM(lHf,norb*norb*nparams)
      if(.not.GetHf(XX(lHf),.true.)) then
         call abrtx("CompAlpha: not GetHF")
      end if
      call MakeDDStep1(XX(lDD),C0,E0,w,XX(lHf),.true.,
     &                 norb,nocc,nvir,nparams)
      call DELETEMEM(lHf,norb*norb*nparams)
      call MakeDDStep2(XX(lDD),C0,E0,w,dummyQ,.false.,
     &                 norb,nocc,nvir,nparams)
      call MakeDDStep3(XX(lDD),C0,dummyX,.false.,
     &                 norb,nocc,nparams)
C     -----
C     Solve
C     -----
      call NEWMEM(lD,norb*norb*nparams)
      call TDHFCycle(GetDJK,PrintAlphaConvergedParam,
     &               C0,E0,w,XX(lD),XX(lDD),
     &               norb,nocc,nvir,nparams)
C     -----------
C     complete GU
C     -----------
      lG = lDD
      call NEWMEM(lU,norb*norb*nparams)
      call NEWMEM(lE,norb*norb*nparams)
      call CompleteGUXE(C0,E0,w,XX(lG),XX(lU),dummyX,XX(lE),
     &                  norb,nocc,nparams)
C      call PrintMM(XX(lG),norb,nfld,'OLIX:Gf')
C      call PrintMM(XX(lU),norb,nfld,'OLIX:Uf')
C      call PrintMM(XX(lD),norb,nfld,'OLIX:Df')
C      call PrintMM(XX(lE),norb,nfld,'OLIX:Ef')
      call FreeGUDEf(w)
      call SetGUDEf(w,XX(lG),XX(lU),XX(lD),XX(lE),norb,nfld)
      call DELETEMEM(lE,norb*norb*nparams)
      call DELETEMEM(lU,norb*norb*nparams)
C     --------------------
C     compute alpha tensor
C     --------------------
      lHf = lDD
      if(.not.GetHf(XX(lHF),.true.)) then
         call abrtx("CompAlpha: not getHF")
      end if
      call CompAlphaTensor(alpha,XX(lHf),XX(lD),norb,ibfld,nfld)
C     ----------
C     deallocate
C     ----------
      call DELETEMEM(lD,norb*norb*nparams)
      call DELETEMEM(lDD,norb*norb*nparams)
C     Format
 3000 FORMAT('Iterative procedure to compute Alpha(',
     &       F13.6,';',F13.6,')')
      end
C
C*MODULE TDXITR  *DECK CompAlphaTensor
      subroutine CompAlphaTensor(alpha,Hf,Dg,norb,ibfld,nfld)
      implicit none
C     Arguments
      integer norb,ibfld,nfld
      double precision alpha(9)
      double precision Hf(norb*norb,nfld),Dg(norb*norb,nfld)
C     Used common block
      integer NELEMA,IELEMA,INDEXA
      common/NLO_ALPHA/NELEMA,IELEMA(2,9),INDEXA(3,3)
C     Variables
      integer i,idummy,ifld(2),icomp
C     Called subprograms
      double precision dtrace
C     Body of the subprogram
      do i = 1,NELEMA
         do icomp=1,2
            idummy = IELEMA(icomp,i)-ibfld+1
            if((idummy.lt.1).or.(idummy.gt.nfld)) go to 1001
            ifld(icomp) = idummy
            if(icomp.eq.2) alpha(i)=-dtrace('S','N',Hf(1,ifld(1)),norb,
     &                                     Dg(1,ifld(2)),norb,norb,norb)
         end do
 1001    continue
      end do
C
      return
      end
C
C*MODULE TDXITR  *DECK PrintAlphaConvergedParam
      subroutine PrintAlphaConvergedParam(iw,ixyz,conv)
      implicit none
C     Arguments
      integer iw,ixyz
      double precision conv
C     Used common block
      integer soli,ngeom,nocc,nvir,norb,nxyz,nmodes,nfld,ibfld,iefld
      common/OLI/soli,ngeom,nocc,nvir,norb,nxyz,nmodes,nfld,ibfld,iefld
C     Variables and constants
      integer ifld
      character*1 FIELDDIR(3)
C     Body of the subprogram
      DATA FIELDDIR /'x','y','z'/
      ifld = ixyz+ibfld-1
      write(iw,1001) FIELDDIR(ifld),conv
C     Format
 1001 FORMAT('  Field direction ',A1,' converged within ',E13.6)
      end
C
C*MODULE TDXITR  *DECK MakeBeta
      subroutine MakeBeta(w1,w2)
      implicit none
C     Argument
      double precision w1,w2
C     Used common blocks
      integer soli,ngeom,nocc,nvir,norb,nxyz,nmodes,nfld,ibfld,iefld
      common/OLI/soli,ngeom,nocc,nvir,norb,nxyz,nmodes,nfld,ibfld,iefld
      integer IR, IW, IP, IS,IPK,IDAF,NAV,IODA
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      DOUBLE PRECISION XX
      COMMON /FMCOM / XX(1)
C     Variables and constants
      double precision mu(3)
      double precision beta(27)
      double precision w(3),lw(2)
      integer lE0,lC0
      double precision dummy(1)
C     Called subprograms
      logical GetHf,GotDip,GetGUDEf
C     Body of the subprogram
      call VCLR(beta,1,27)
C
      w(1)=-w1-w2
      w(2)=w1
      w(3)=w2
      lw(1)=w1
      lw(2)=w2
C
      dummy(1) = 0.0d+00
      if(.not.GetHf(dummy,.false.)) call MakeHf()
      call GetDipole(dummy,.false.,GotDip)
      if(.not.GotDip) call MakeDipole()
      if(.not.GetGUDEf(w1,dummy,.false.,dummy,.false.,dummy,.false.,
     &     dummy,.false.,norb,nfld)) call MakeAlpha(w1)
      if(.not.GetGUDEf(w2,dummy,.false.,dummy,.false.,dummy,.false.,
     &     dummy,.false.,norb,nfld)) call MakeAlpha(w2)
C
      call NEWMEM(lE0,norb)
      call NEWMEM(lC0,norb*norb)
      call GetEnergy(XX(lE0),norb)
      call GetLCAO(XX(lC0),norb)
C
      call CompBeta(beta,lw,XX(lE0),XX(lC0),
     &              norb,nocc,nvir,ibfld,nfld)
C
      call DELETEMEM(lC0,norb*norb)
      call DELETEMEM(lE0,norb)
C
      call SetBeta(w1,w2,beta)
      call GetDipole(mu,.true.,GotDip)
      if(.not.GotDip) then
         call abrtx("Make Beta:not GotDip")
      end if
      call PrintBetaX(iw,beta,mu,w,'Beta tensor [in au]')
      CALL TIMIT(1)
      end
C
C*MODULE TDXITR  *DECK CompBeta
      subroutine CompBeta(beta,w,E0,C0,
     &                    norb,nocc,nvir,ibfld,nfld)
      implicit none
C     Arguments
      integer norb,nocc,nvir,ibfld,nfld
      double precision w(2)
      double precision beta(27)
      double precision E0(norb),C0(norb,norb)
C     Used common blocks
      integer IR, IW, IP, IS,IPK,IDAF,NAV,IODA
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      DOUBLE PRECISION XX
      COMMON /FMCOM / XX(1)
C     Variables
      integer nparams
      double precision ws
      integer lD,lDD,lHf,lG,lU,lE,lX,lQ
      double precision dummyF(1)
C     Called subprogram
      logical GetHf
C     External subprogram
      external GetDJK,PrintBetaConvergedParam
C     Body of the subprogram
      nparams = nfld*nfld
      ws=w(1)+w(2)
      write(iw,3000) (-ws),w(1),w(2)
C     ----------
C     prepare DD
C     ----------
      call NEWMEM(lDD,norb*norb*nparams)
      call MakeDDfg(XX(lDD),w,C0,norb,nocc,nfld)
      call MakeDDStep0
      call MakeDDStep1(XX(lDD),C0,E0,ws,dummyF,.false.,
     &                 norb,nocc,nvir,nparams)
      call NEWMEM(lQ,norb*norb*nparams)
      call MakeQfg(XX(lQ),w,norb,nfld)
      call MakeDDStep2(XX(lDD),C0,E0,ws,XX(lQ),.true.,
     &                 norb,nocc,nvir,nparams)
C      call PrintMM(XX(lDD),norb,nparams,'DDfg')
      call DELETEMEM(lQ,norb*norb*nparams)
      call NEWMEM(lX,norb*norb*nparams)
      call MakeXfg(XX(lX),w,norb,nfld)
      call MakeDDStep3(XX(lDD),C0,XX(lX),.true.,
     &                 norb,nocc,nparams)
C      call PrintMM(XX(lDD),norb,nparams,'DDfg')
      call DELETEMEM(lX,norb*norb*nparams)
C      call PrintMM(XX(lDD),norb,nparams,'DDfg')
C     -----
C     Solve
C     -----
      call NEWMEM(lD,norb*norb*nparams)
      call TDHFCycle(GetDJK,PrintBetaConvergedParam,
     &               C0,E0,ws,XX(lD),XX(lDD),
     &               norb,nocc,nvir,nparams)
C     -----------
C     complete GU
C     -----------
      lG = lDD
      call NEWMEM(lU,norb*norb*nparams)
      call NEWMEM(lE,norb*norb*nparams)
      call NEWMEM(lX,norb*norb*nparams)
      call CompleteGUXE(C0,E0,ws,XX(lG),XX(lU),XX(lX),XX(lE),
     &                  norb,nocc,nparams)
      call FreeGUXDEfg(w(1),w(2))
      call SetGUXDEfg(w(1),w(2),XX(lG),XX(lU),XX(lX),XX(lD),XX(lE),
     &                norb,nfld)
      call DELETEMEM(lX,norb*norb*nparams)
      call DELETEMEM(lE,norb*norb*nparams)
      call DELETEMEM(lU,norb*norb*nparams)
C     --------------------
C     compute beta tensor
C     --------------------
      lHf = lDD
      if(.not.GetHf(XX(lHF),.true.)) then
         call abrtx("CompBeta: not GetHF")
      end if
      call CompBetaTensor(beta,XX(lHf),XX(lD),norb,ibfld,nfld)
C     ----------
C     deallocate
C     ----------
      call DELETEMEM(lD,norb*norb*nparams)
      call DELETEMEM(lDD,norb*norb*nparams)
C     Format
 3000 FORMAT('Iterative procedure to compute Beta(',
     &       F13.6,';',F13.6,',',F13.6,')')
      end
C
C*MODULE TDXITR  *DECK CompBetaTensor
      subroutine CompBetaTensor(beta,Hf,Dgh,norb,ibfld,nfld)
      implicit none
C     Arguments
      integer norb,ibfld,nfld
      double precision beta(27)
      double precision Hf(norb*norb,nfld),Dgh(norb*norb,nfld,nfld)
C     Used common block
      integer NELEMB,IELEMB,INDEXB
      common/NLO_BETA/NELEMB,IELEMB(3,27),INDEXB(3,3,3)
C     Variables
      integer i,idummy,ifld(3),icomp
C     Called subprograms
      double precision DTrace
C     Body of the subprogram
      do i = 1,NELEMB
         do icomp=1,3
            idummy = IELEMB(icomp,i)-ibfld+1
            if((idummy.lt.1).or.(idummy.gt.nfld)) go to 1001
            ifld(icomp) = idummy
            if(icomp.eq.3) beta(i)=-DTrace('S','N',Hf(1,ifld(1)),norb,
     &                                     Dgh(1,ifld(2),ifld(3)),norb,
     &                                     norb,norb)
         end do
 1001    continue
      end do
C
      return
      end
C
C*MODULE TDXITR  *DECK PrintBetaConvergedParam
      subroutine PrintBetaConvergedParam(iw,ixyz,conv)
      implicit none
C     Arguments
      integer iw,ixyz
      double precision conv
C     Used common block
      integer soli,ngeom,nocc,nvir,norb,nxyz,nmodes,nfld,ibfld,iefld
      common/OLI/soli,ngeom,nocc,nvir,norb,nxyz,nmodes,nfld,ibfld,iefld
C     Variables and constants
      integer ifld,jfld
      character*1 FIELDDIR(3)
C     Body of the subprogram
      DATA FIELDDIR /'x','y','z'/
      ifld = (ixyz-1)/nfld
      jfld = (ixyz-1-ifld*nfld)
      ifld = ifld+ibfld
      jfld = jfld+ibfld
      write(iw,1001) FIELDDIR(ifld),FIELDDIR(jfld),conv
C     Format
 1001 FORMAT('  Field directions ',A1,A1,' converged within ',E13.6)
      end
C
C*MODULE TDXITR  *DECK MakeBetaPerms
      subroutine MakeBetaPerms(w,nfld,nterms,iflds,mults)
      implicit none
C     Arguments
      double precision w(2)
      integer nfld
      integer nterms(2)
      integer iflds(4,*)
      double precision mults(*)
C     Variables and constants
      logical bw12,b12
      integer iterm,ifld,jfld,iifld,jjfld
      double precision mult,ONE,TWO
      parameter (ONE=1.0D+00,TWO=2.0D+00)
      integer iperm,nperm
      parameter (nperm=2)
C     Body of the subprogram
C     ------------------
C     Smart permutations
C     ------------------
      bw12=w(1).eq.w(2)
      iterm = 0
      do iperm=1,nperm
         do jjfld=1,nfld
            do 1011 iifld=1,nfld
               b12=(iifld.eq.jjfld).and.bw12
               if(iperm.eq.1) then
                  ifld = iifld
                  jfld = jjfld
                  if(b12) then
                     mult = TWO
                  else
                     mult = ONE
                  end if
               else if(iperm.eq.2) then
                  if(b12) go to 1011
                  ifld = jjfld
                  jfld = iifld
                  mult = ONE
               end if
               iterm=iterm+1
               iflds(1,iterm) = ifld
               iflds(2,iterm) = jfld
               iflds(3,iterm) = iifld
               iflds(4,iterm) = jjfld
               mults(iterm)   = mult
 1011       continue
         end do
         nterms(iperm) = iterm
      end do
      end
C
C*MODULE TDXITR  *DECK MakeQfg
      subroutine MakeQfg(Qfg,w,norb,nfld)
      implicit none
C     Arguments
      integer norb,nfld
      double precision w(2)
      double precision Qfg(norb*norb,nfld,nfld)
C     Used common block
      DOUBLE PRECISION XX
      COMMON /FMCOM / XX(1)
C     Variables
      logical b12
      integer lGf,lUf,lEf,lGg,lUg,lEg
      double precision dummy
C     Called subprogram
      logical GetGUDEf
C     Body of the subprogram
      b12 = w(1).eq.w(2)
      call NEWMEM(lGf,norb*norb*nfld)
      call NEWMEM(lUf,norb*norb*nfld)
      call NEWMEM(lEf,norb*norb*nfld)
      dummy = 0.0d+00
      if(.not.GetGUDEf(w(1),XX(lGf),.true.,XX(lUf),.true.,
     &     dummy,.false.,XX(lEf),.true.,norb,nfld)) then
         call abrtx("MakeQfg: not GetGUDEf:1")
      end if
      if(.not.b12) then
         call NEWMEM(lGg,norb*norb*nfld)
         call NEWMEM(lUg,norb*norb*nfld)
         call NEWMEM(lEg,norb*norb*nfld)
         if(.not.GetGUDEf(w(2),XX(lGg),.true.,XX(lUg),.true.,
     &        dummy,.false.,XX(lEg),.true.,norb,nfld)) then
            call abrtx("MakeQfg: not GetGUDEf:2")
         end if
      else
         lGg = lGf
         lUg = lUf
         lEg = lEf
      end if
      call CompQfg(w,Qfg,XX(lGf),XX(lUf),XX(lEf),
     &             XX(lGg),XX(lUg),XX(lEg),
     &             norb,nfld)
      if(.not.b12) then
         call DELETEMEM(lEg,norb*norb*nfld)
         call DELETEMEM(lUg,norb*norb*nfld)
         call DELETEMEM(lGg,norb*norb*nfld)
      end if
      call DELETEMEM(lEf,norb*norb*nfld)
      call DELETEMEM(lUf,norb*norb*nfld)
      call DELETEMEM(lGf,norb*norb*nfld)
      end
C
C*MODULE TDXITR  *DECK CompQfg
      subroutine CompQfg(w,Qfg,Gf,Uf,Ef,Gg,Ug,Eg,norb,nfld)
      implicit none
C     Arguments
      double precision w(2)
      integer norb,nfld
      double precision Qfg(norb*norb,nfld,nfld)
      double precision Gf(norb*norb,nfld)
      double precision Uf(norb*norb,nfld)
      double precision Ef(norb*norb,nfld)
      double precision Gg(norb*norb,nfld)
      double precision Ug(norb*norb,nfld)
      double precision Eg(norb*norb,nfld)
C     Variables
      integer iterm,ifld,jfld,iifld,jjfld
      double precision mult
      integer nterms(2)
      integer iflds(4,2*9)
      double precision mults(2*9)
C     Body of the subprogram
      call MakeBetaPerms(w,nfld,nterms,iflds,mults)
      call VCLR(Qfg,1,norb*norb*nfld*nfld)
      do iterm=1,nterms(1)
         ifld = iflds(1,iterm)
         jfld = iflds(2,iterm)
         iifld = iflds(3,iterm)
         jjfld = iflds(4,iterm)
         mult = mults(iterm)
         call CompQterm(Qfg(1,iifld,jjfld),
     &                  Gf(1,ifld),Uf(1,ifld),'o',
     &                  Ug(1,jfld),Eg(1,jfld),
     &                  norb,mult)
      end do
      do iterm=nterms(1)+1,nterms(2)
         ifld = iflds(1,iterm)
         jfld = iflds(2,iterm)
         iifld = iflds(3,iterm)
         jjfld = iflds(4,iterm)
         mult = mults(iterm)
         call CompQterm(Qfg(1,iifld,jjfld),
     &                  Gg(1,ifld),Ug(1,ifld),'o',
     &                  Uf(1,jfld),Ef(1,jfld),
     &                  norb,mult)
      end do
      end
C
C*MODULE TDXITR  *DECK MakeXfg
      subroutine MakeXfg(Xfg,w,norb,nfld)
      implicit none
C     Arguments
      integer norb,nfld
      double precision w(2)
      double precision Xfg(norb*norb,nfld,nfld)
C     Used common block
      DOUBLE PRECISION XX
      COMMON /FMCOM / XX(1)
C     Variables
      logical b12
      integer lUf,lUg
      double precision dummy
C     Called subprogram
      logical GetGUDEf
C     Body of the subprogram
      b12 = w(1).eq.w(2)
      call NEWMEM(lUf,norb*norb*nfld)
      dummy = 0.0d+00
      if(.not.GetGUDEf(w(1),dummy,.false.,XX(lUf),.true.,dummy,.false.,
     &                    dummy,.false.,norb,nfld)) then
         call abrtx("MakeXfg: not GetGUDEf:1")
      end if
      if(.not.b12) then
         call NEWMEM(lUg,norb*norb*nfld)
         if(.not.GetGUDEf(w(2),dummy,.false.,XX(lUg),.true.,
     &        dummy,.false.,dummy,.false.,norb,nfld)) then
            call abrtx("MakeXfg: not GetGUDEf:2")
         end if
      else
         lUg = lUf
      end if
      call CompXfg(w,Xfg,XX(lUf),XX(lUg),norb,nfld)
      if(.not.b12) call DELETEMEM(lUg,norb*norb*nfld)
      call DELETEMEM(lUf,norb*norb*nfld)
C
C      call PrintMM(Xfg,norb,nfld*nfld,"TDX:Xfg")
      end
C
C*MODULE TDXITR  *DECK CompXfg
      subroutine CompXfg(w,Xfg,Uf,Ug,norb,nfld)
      implicit none
C     Arguments
      double precision w(2)
      integer norb,nfld
      double precision Xfg(norb*norb,nfld,nfld)
      double precision Uf(norb*norb,nfld)
      double precision Ug(norb*norb,nfld)
C     Variables
      integer iterm,ifld,jfld,iifld,jjfld
      double precision mult
      integer nterms(2)
      integer iflds(4,2*9)
      double precision mults(2*9)
      double precision ONE
      parameter (ONE=1.0D+00)
      character*1 sN
      parameter (sN='N')
C     Body of the subprogram
      call MakeBetaPerms(w,nfld,nterms,iflds,mults)
      call VCLR(Xfg,1,norb*norb*nfld*nfld)
      do iterm=1,nterms(1)
         ifld = iflds(1,iterm)
         jfld = iflds(2,iterm)
         iifld = iflds(3,iterm)
         jjfld = iflds(4,iterm)
         mult = mults(iterm)
         call DGEMM(sN,sN,norb,norb,norb,
     &              mult,Uf(1,ifld),norb,Ug(1,jfld),norb,
     &              ONE,Xfg(1,iifld,jjfld),norb)
      end do
      do iterm=nterms(1)+1,nterms(2)
         ifld = iflds(1,iterm)
         jfld = iflds(2,iterm)
         iifld = iflds(3,iterm)
         jjfld = iflds(4,iterm)
         mult = mults(iterm)
         call DGEMM(sN,sN,norb,norb,norb,
     &              mult,Ug(1,ifld),norb,Uf(1,jfld),norb,
     &              ONE,Xfg(1,iifld,jjfld),norb)
      end do
      end
C
C*MODULE TDXITR  *DECK MakeDDfg
      subroutine MakeDDfg(DDfg,w,C0,norb,nocc,nfld)
      implicit none
C     Arguments
      integer norb,nocc,nfld
      double precision w(2)
      double precision DDfg(norb*norb,nfld,nfld)
      double precision C0(norb*norb)
C     Used common block
      DOUBLE PRECISION XX
      COMMON /FMCOM / XX(1)
C     Variables
      logical b12
      integer lUf,lUg
      double precision dummy
C     Called subprogram
      logical GetGUDEf
C     Body of the subprogram
      b12 = w(1).eq.w(2)
      call NEWMEM(lUf,norb*norb*nfld)
      dummy = 0.0d+00
      if(.not.GetGUDEf(w(1),dummy,.false.,XX(lUf),.true.,dummy,.false.,
     &                 dummy,.false.,norb,nfld)) then
         call abrtx("MakeDDfg: not GetGUDEf:1")
      end if
      if(.not.b12) then
         call NEWMEM(lUg,norb*norb*nfld)
         if(.not.GetGUDEf(w(2),dummy,.false.,XX(lUg),.true.,
     &        dummy,.false.,dummy,.false.,norb,nfld)) then
            call abrtx("MakeDDfg: not GetGUDEf:1")
         end if
      else
         lUg = lUf
      end if
      call CompDDfg(w,DDfg,XX(lUf),XX(lUg),C0,norb,nocc,nfld)
      if(.not.b12) call DELETEMEM(lUg,norb*norb*nfld)
      call DELETEMEM(lUf,norb*norb*nfld)
      end
C
C*MODULE TDXITR  *DECK CompDDfg
      subroutine CompDDfg(w,DDfg,Uf,Ug,C0,norb,nocc,nfld)
      implicit none
C     Arguments
      double precision w(2)
      integer norb,nocc,nfld
      double precision DDfg(norb*norb,nfld,nfld)
      double precision Uf(norb*norb,nfld)
      double precision Ug(norb*norb,nfld)
      double precision C0(norb*norb)
C     Variables
      integer iterm,ifld,jfld,iifld,jjfld
      double precision mult
      integer nterms(2)
      integer iflds(4,2*9)
      double precision mults(2*9)
C     Body of the subprogram
      call MakeBetaPerms(w,nfld,nterms,iflds,mults)
      call VCLR(DDfg,1,norb*norb*nfld*nfld)
      do iterm=1,nterms(1)
         ifld = iflds(1,iterm)
         jfld = iflds(2,iterm)
         iifld = iflds(3,iterm)
         jjfld = iflds(4,iterm)
         mult = mults(iterm)
         call CompUsum(DDfg(1,iifld,jjfld),
     &                 Uf(1,ifld),Ug(1,jfld),'f',
     &                 norb,nocc,mult)
      end do
      do iterm=nterms(1)+1,nterms(2)
         ifld = iflds(1,iterm)
         jfld = iflds(2,iterm)
         iifld = iflds(3,iterm)
         jjfld = iflds(4,iterm)
         mult = mults(iterm)
         call CompUsum(DDfg(1,iifld,jjfld),
     &                 Ug(1,ifld),Uf(1,jfld),'f',
     &                 norb,nocc,mult)
      end do
      call CompDDterm(DDfg,C0,DDfg,norb,nfld*nfld)
      end
C*MODULE TDXITR  *DECK MakeGamma
      subroutine MakeGamma(w1,w2,w3)
      implicit none
C     Arguments
      double precision w1,w2,w3
C     Used common blocks
      integer soli,ngeom,nocc,nvir,norb,nxyz,nmodes,nfld,ibfld,iefld
      common/OLI/soli,ngeom,nocc,nvir,norb,nxyz,nmodes,nfld,ibfld,iefld
      integer IR, IW, IP, IS,IPK,IDAF,NAV,IODA
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      DOUBLE PRECISION XX
      COMMON /FMCOM / XX(1)
C     Variables and constants
      double precision gamma(21)
      double precision w(4),lw(3)
      integer lE0,lC0
      double precision dummy
C     Called subprograms
      logical GetHf,GetGUDEf,GetGUXDEfg
C     Body of the subprogram
      call VCLR(gamma,1,21)
C
      w(1)=-w1-w2-w3
      w(2)=w1
      w(3)=w2
      w(4)=w3
      lw(1)=w1
      lw(2)=w2
      lw(3)=w3
C
      dummy = 0.0d+00
      if(.not.GetHf(dummy,.false.)) call MakeHf()
      if(.not.GetGUDEf(w1,dummy,.false.,dummy,.false.,dummy,.false.,
     &     dummy,.false.,norb,nfld)) call MakeAlpha(w1)
      if(.not.GetGUDEf(w2,dummy,.false.,dummy,.false.,dummy,.false.,
     &     dummy,.false.,norb,nfld)) call MakeAlpha(w2)
      if(.not.GetGUDEf(w3,dummy,.false.,dummy,.false.,dummy,.false.,
     &     dummy,.false.,norb,nfld)) call MakeAlpha(w3)
      if(.not.GetGUXDEfg(w1,w2,dummy,.false.,dummy,.false.,
     &     dummy,.false.,dummy,.false.,dummy,.false.,norb,nfld))
     &     call MakeBeta(w1,w2)
      if(.not.GetGUXDEfg(w1,w3,dummy,.false.,dummy,.false.,
     &     dummy,.false.,dummy,.false.,dummy,.false.,norb,nfld))
     &     call MakeBeta(w1,w3)
      if(.not.GetGUXDEfg(w2,w3,dummy,.false.,dummy,.false.,
     &     dummy,.false.,dummy,.false.,dummy,.false.,norb,nfld))
     &     call MakeBeta(w2,w3)
C
      call NEWMEM(lE0,norb)
      call NEWMEM(lC0,norb*norb)
      call GetEnergy(XX(lE0),norb)
      call GetLCAO(XX(lC0),norb)
C
      call CompGamma(gamma,lw,XX(lE0),XX(lC0),
     &               norb,nocc,nvir,ibfld,nfld)
C
      call DELETEMEM(lC0,norb*norb)
      call DELETEMEM(lE0,norb)
C
      call SetGamma(w1,w2,w3,gamma)
      call PrintGammaX(iw,gamma,w,'Gamma tensor [in au]')
      CALL TIMIT(1)
      end
C
C*MODULE TDXITR  *DECK CompGamma
      subroutine CompGamma(gamma,w,E0,C0,
     &                     norb,nocc,nvir,ibfld,nfld)
      implicit none
C     Arguments
      integer norb,nocc,nvir,ibfld,nfld
      double precision w(3)
      double precision gamma(21)
      double precision E0(norb),C0(norb,norb)
C     Used common blocks
      integer IR, IW, IP, IS,IPK,IDAF,NAV,IODA
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      DOUBLE PRECISION XX
      COMMON /FMCOM / XX(1)
C     Variables
      integer nparams
      double precision ws
      integer lD,lDD,lHf,lG,lU,lE,lX,lQ
      double precision dummyF(1)
C     Called subprogram
      logical GetHf
C     External subprogram
      external GetDJK,PrintGammaConvergedParam
C     Body of the subprogram
      nparams = nfld*nfld*nfld
      ws=w(1)+w(2)+w(3)
      write(iw,3000) (-ws),w(1),w(2),w(3)
C     ----------
C     prepare DD
C     ----------
      call NEWMEM(lDD,norb*norb*nparams)
      call MakeDDfgh(XX(lDD),w,C0,norb,nocc,nfld)
      call MakeDDStep0
      call MakeDDStep1(XX(lDD),C0,E0,ws,dummyF,.false.,
     &                 norb,nocc,nvir,nparams)
      call NEWMEM(lQ,norb*norb*nparams)
      call MakeQfgh(XX(lQ),w,norb,nfld)
      call MakeDDStep2(XX(lDD),C0,E0,ws,XX(lQ),.true.,
     &                 norb,nocc,nvir,nparams)
C      call PrintMM(XX(lDD),norb,nparams,'DDfgh')
      call DELETEMEM(lQ,norb*norb*nparams)
      call NEWMEM(lX,norb*norb*nparams)
      call MakeXfgh(XX(lX),w,norb,nfld)
      call MakeDDStep3(XX(lDD),C0,XX(lX),.true.,
     &                 norb,nocc,nparams)
C      call PrintMM(XX(lDD),norb,nparams,'DDfgh')
      call DELETEMEM(lX,norb*norb*nparams)
C      call PrintMM(XX(lDD),norb,nparams,'DDfgh')
C     -----
C     Solve
C     -----
      call NEWMEM(lD,norb*norb*nparams)
      call TDHFCycle(GetDJK,PrintGammaConvergedParam,
     &               C0,E0,ws,XX(lD),XX(lDD),
     &               norb,nocc,nvir,nparams)
C     -----------
C     complete GU
C     -----------
      lG = lDD
      call NEWMEM(lU,norb*norb*nparams)
      call NEWMEM(lE,norb*norb*nparams)
      call NEWMEM(lX,norb*norb*nparams)
      call CompleteGUXE(C0,E0,ws,XX(lG),XX(lU),XX(lX),XX(lE),
     &                  norb,nocc,nparams)
      call FreeGUXDEfgh(w(1),w(2),w(3))
      call SetGUXDEfgh(w(1),w(2),w(3),
     &                 XX(lG),XX(lU),XX(lX),XX(lD),XX(lE),
     &                 norb,nfld)
      call DELETEMEM(lX,norb*norb*nparams)
      call DELETEMEM(lE,norb*norb*nparams)
      call DELETEMEM(lU,norb*norb*nparams)
C     --------------------
C     compute gamma tensor
C     --------------------
      lHf = lDD
      if(.not.GetHf(XX(lHF),.true.)) then
         call abrtx("CompGamma: not GetHF")
      end if
      call CompGammaTensor(gamma,XX(lHf),XX(lD),norb,ibfld,nfld)
C     ----------
C     deallocate
C     ----------
      call DELETEMEM(lD,norb*norb*nparams)
      call DELETEMEM(lDD,norb*norb*nparams)
C     Format
 3000 FORMAT('Iterative procedure to compute Gamma(',
     &       F13.6,';',F13.6,',',F13.6,',',F13.6,')')
      end
C
C*MODULE TDXITR  *DECK CompGammaTensor
      subroutine CompGammaTensor(gamma,Hf,Dghi,norb,ibfld,nfld)
      implicit none
C     Arguments
      integer norb,ibfld,nfld
      double precision gamma(21)
      double precision Hf(norb*norb,nfld),Dghi(norb*norb,nfld,nfld,nfld)
C     Used common block
      integer NELEMG,IELEMG,INDEXG
      common/NLO_GAMMA/NELEMG,IELEMG(4,21),INDEXG(3,3,3,3)
C     Variables
      integer i,idummy,ifld(4),icomp
C     Called subprograms
      double precision DTrace
C     Body of the subprogram
      do i = 1,NELEMG
         do icomp=1,4
            idummy = IELEMG(icomp,i)-ibfld+1
            if((idummy.lt.1).or.(idummy.gt.nfld)) go to 1001
            ifld(icomp) = idummy
            if(icomp.eq.4)
     &           gamma(i)=-DTrace('S','N',Hf(1,ifld(1)),norb,
     &                            Dghi(1,ifld(2),ifld(3),ifld(4)),norb,
     &                            norb,norb)
         end do
 1001    continue
      end do
      end
C
C*MODULE TDXITR  *DECK PrintGammaConvergedParam
      subroutine PrintGammaConvergedParam(iw,ixyz,conv)
      implicit none
C     Arguments
      integer iw,ixyz
      double precision conv
C     Used common block
      integer soli,ngeom,nocc,nvir,norb,nxyz,nmodes,nfld,ibfld,iefld
      common/OLI/soli,ngeom,nocc,nvir,norb,nxyz,nmodes,nfld,ibfld,iefld
C     Variables and constants
      integer ifld,jfld,kfld
      character*1 FIELDDIR(3)
C     Body of the subprogram
      DATA FIELDDIR /'x','y','z'/
      ifld = (ixyz-1)/(nfld*nfld)
      jfld = (ixyz-1-ifld*nfld*nfld)/nfld
      kfld = (ixyz-1-ifld*nfld*nfld-jfld*nfld)
      ifld = ifld+ibfld
      jfld = jfld+ibfld
      kfld = kfld+ibfld
      write(iw,1001) FIELDDIR(ifld),FIELDDIR(jfld),FIELDDIR(kfld),conv
C     Format
 1001 FORMAT('  Field directions ',3A1,' converged within ',E13.6)
      end
C
C*MODULE TDXITR  *DECK MakeGammaPerms
      subroutine MakeGammaPerms(w,nfld,nterms,iflds,mults,ws,bgws)
      implicit none
C     Arguments
      double precision w(3)
      integer nfld
      integer nterms(3)
      integer iflds(6,*)
      double precision mults(*)
      double precision ws(3,3)
      logical bgws(3,3)
C     Variables and constants
C--      logical bw12,bw13,bw23,b12,b13,b23
C--      double precision mult,TWO,THREE
C--      parameter (TWO=2.0D+00,THREE=3.0D+00)
C--      integer i,permlist(3)
      integer iterm,ifld,jfld,kfld,iifld,jjfld,kkfld
      double precision ONE
      integer iperm,nperm
      parameter (ONE=1.0D+00)
      parameter (nperm=3)
C     Body of the subprogram
      if(.true.) then
         ws(1,1) = w(1)
         ws(2,1) = w(2)
         ws(3,1) = w(3)
         ws(1,2) = w(2)
         ws(2,2) = w(1)
         ws(3,2) = w(3)
         ws(1,3) = w(3)
         ws(2,3) = w(1)
         ws(3,3) = w(2)
         bgws(1,1) = .true.
         bgws(2,1) = .true.
         bgws(3,1) = .true.
         bgws(1,2) = .true.
         bgws(2,2) = .true.
         bgws(3,2) = .true.
         bgws(1,3) = .true.
         bgws(2,3) = .true.
         bgws(3,3) = .true.
         iterm=0
         do iperm=1,nperm
            do kkfld=1,nfld
               do jjfld=1,nfld
                  do iifld=1,nfld
                     iterm=iterm+1
                     mults(iterm) = ONE
                     if(iperm.eq.1) then
                        ifld = iifld
                        jfld = jjfld
                        kfld = kkfld
                     else if(iperm.eq.2) then
                        ifld = jjfld
                        jfld = iifld
                        kfld = kkfld
                     else if(iperm.eq.3) then
                        ifld = kkfld
                        jfld = iifld
                        kfld = jjfld
                     end if
                     iflds(1,iterm) = ifld
                     iflds(2,iterm) = jfld
                     iflds(3,iterm) = kfld
                     iflds(4,iterm) = iifld
                     iflds(5,iterm) = jjfld
                     iflds(6,iterm) = kkfld
                  end do
               end do
            end do
            nterms(iperm)=iterm
         end do
         return
      end if
C--c     ------------------
C--c     Smart permutations
C--c     ------------------
C--      bw12=w(1).eq.w(2)
C--      bw13=w(1).eq.w(3)
C--      bw23=w(2).eq.w(3)
C--      permlist(1) = 1
C--      permlist(2) = 2
C--      permlist(3) = 3
C--      iterm = 0
C--      do i=1,nperm
C--         iperm = permlist(i)
C--         if(iperm.eq.1) then
C--            ws(1,i) = w(1)
C--            ws(2,i) = w(2)
C--            ws(3,i) = w(3)
C--         else if(iperm.eq.2) then
C--            ws(1,i) = w(2)
C--            ws(2,i) = w(1)
C--            ws(3,i) = w(3)
C--         else if(iperm.eq.3) then
C--            ws(1,i) = w(3)
C--            ws(2,i) = w(1)
C--            ws(3,i) = w(2)
C--         end if
C--         do kkfld=1,nfld
C--            do jjfld=1,nfld
C--               do 1011 iifld=1,nfld
C--                  b12=(iifld.eq.jjfld).and.bw12
C--                  b13=(iifld.eq.kkfld).and.bw13
C--                  b23=(jjfld.eq.kkfld).and.bw23
C--c
C--                  if(iperm.eq.1) then
C--                     ifld = iifld
C--                     jfld = jjfld
C--                     kfld = kkfld
C--                     if(b12.and.b13) then
C--                        mult = THREE
C--                     else if(b12) then
C--                        mult = TWO
C--                     else
C--                        mult = ONE
C--                     end if
C--                  else if(iperm.eq.2) then
C--                     if((b12.and.b13).or.b12) go to 1011
C--                     ifld = jjfld
C--                     jfld = iifld
C--                     kfld = kkfld
C--                     if(b23) then
C--                        mult = TWO
C--                     else
C--                        mult = ONE
C--                     end if
C--                  else if(iperm.eq.3) then
C--                     if((b12.and.b13).or.b23) go to 1011
C--                     ifld = kkfld
C--                     jfld = iifld
C--                     kfld = jjfld
C--                     mult = ONE
C--                  end if
C--                  iterm=iterm+1
C--                  iflds(1,iterm) = ifld
C--                  iflds(2,iterm) = jfld
C--                  iflds(3,iterm) = kfld
C--                  iflds(4,iterm) = iifld
C--                  iflds(5,iterm) = jjfld
C--                  iflds(6,iterm) = kkfld
C--                  mults(iterm)   = mult
C-- 1011          continue
C--            end do
C--         end do
C--         nterms(i) = iterm
C--      end do
C--      do i = 1,3
C--         bgws(i,1) = .true.
C--         bgws(i,2) = ws(i,1).ne.ws(i,2)
C--         bgws(i,3) = ws(i,2).ne.ws(i,3)
C--      end do
      end
C
C*MODULE TDXITR  *DECK MakeQfgh
      subroutine MakeQfgh(Qfgh,w,norb,nfld)
      implicit none
C     Arguments
      integer norb,nfld
      double precision w(3)
      double precision Qfgh(norb*norb,nfld,nfld,nfld)
C     Used common block
      DOUBLE PRECISION XX
      COMMON /FMCOM / XX(1)
C     Variables
      integer lGf,lUf,lEf,lGgh,lUgh,lEgh
C     Body of the subprogram
      call NEWMEM(lGf,norb*norb*nfld)
      call NEWMEM(lUf,norb*norb*nfld)
      call NEWMEM(lEf,norb*norb*nfld)
      call NEWMEM(lGgh,norb*norb*nfld*nfld)
      call NEWMEM(lUgh,norb*norb*nfld*nfld)
      call NEWMEM(lEgh,norb*norb*nfld*nfld)
      call CompQfgh(w,Qfgh,XX(lGf),XX(lUf),XX(lEf),
     &              XX(lGgh),XX(lUgh),XX(lEgh),
     &              norb,nfld)
      call DELETEMEM(lEgh,norb*norb*nfld*nfld)
      call DELETEMEM(lUgh,norb*norb*nfld*nfld)
      call DELETEMEM(lGgh,norb*norb*nfld*nfld)
      call DELETEMEM(lEf,norb*norb*nfld)
      call DELETEMEM(lUf,norb*norb*nfld)
      call DELETEMEM(lGf,norb*norb*nfld)
      end
C
C*MODULE TDXITR  *DECK CompQfgh
      subroutine CompQfgh(w,Qfgh,Gf,Uf,Ef,Ggh,Ugh,Egh,
     &                    norb,nfld)
      implicit none
C     Arguments
      double precision w(3)
      integer norb,nfld
      double precision Qfgh(norb*norb,nfld,nfld,nfld)
      double precision Gf(norb*norb,nfld)
      double precision Uf(norb*norb,nfld)
      double precision Ef(norb*norb,nfld)
      double precision Ggh(norb*norb,nfld,nfld)
      double precision Ugh(norb*norb,nfld,nfld)
      double precision Egh(norb*norb,nfld,nfld)
C     Variables
      integer iterm,ibterm,ieterm,ifld,jfld,kfld,iifld,jjfld,kkfld
      double precision mult,dummy
      integer iperm,nperm
      parameter (nperm=3)
      integer nterms(nperm)
      integer iflds(6,nperm*27)
      double precision mults(nperm*27)
      double precision ws(3,nperm)
      logical bgws(3,nperm)
C     Called subprogram
      logical GetGUDEf,GetGUXDEfg
C     Body of the subprogram
      call VCLR(Qfgh,1,norb*norb*nfld*nfld*nfld)
      call MakeGammaPerms(w,nfld,nterms,iflds,mults,ws,bgws)
C
      ieterm = 0
      dummy = 0.0d+00
      do iperm = 1,nperm
         if(bgws(1,iperm)) then
            if(.not.GetGUDEf(ws(1,iperm),
     &                       Gf,.true.,Uf,.true.,
     &                       dummy,.false.,Ef,.true.,
     &                       norb,nfld)) then
               call abrtx("CompQfgh: not GetGUDEf:1")
            end if
         end if
         if(bgws(2,iperm).or.bgws(3,iperm)) then
            if(.not.GetGUXDEfg(ws(2,iperm),ws(3,iperm),
     &                         Ggh,.true.,Ugh,.true.,dummy,.false.,
     &                         dummy,.false.,Egh,.true.,
     &                         norb,nfld)) then
               call abrtx("CompQfgh: not GetGUDEf:2")
            end if
         end if
         ibterm=ieterm+1
         ieterm=nterms(iperm)
         do iterm = ibterm,ieterm
            ifld = iflds(1,iterm)
            jfld = iflds(2,iterm)
            kfld = iflds(3,iterm)
            iifld = iflds(4,iterm)
            jjfld = iflds(5,iterm)
            kkfld = iflds(6,iterm)
            mult = mults(iterm)
            call CompQterm(Qfgh(1,iifld,jjfld,kkfld),
     &                     Gf(1,ifld),Uf(1,ifld),'o',
     &                     Ugh(1,jfld,kfld),Egh(1,jfld,kfld),
     &                     norb,mult)
            call CompQterm(Qfgh(1,iifld,jjfld,kkfld),
     &                     Ggh(1,jfld,kfld),Ugh(1,jfld,kfld),'o',
     &                     Uf(1,ifld),Ef(1,ifld),
     &                     norb,mult)
         end do
      end do
      end
C
C*MODULE TDXITR  *DECK MakeXfgh
      subroutine MakeXfgh(Xfgh,w,norb,nfld)
      implicit none
C     Arguments
      integer norb,nfld
      double precision w(3)
      double precision Xfgh(norb*norb,nfld,nfld,nfld)
C     Used common block
      DOUBLE PRECISION XX
      COMMON /FMCOM / XX(1)
C     Variables
      integer lUf,lUgh,lXgh
C     Body of the subprogram
      call NEWMEM(lUf,norb*norb*nfld)
      call NEWMEM(lUgh,norb*norb*nfld*nfld)
      call NEWMEM(lXgh,norb*norb*nfld*nfld)
      call CompXfgh(w,Xfgh,XX(lUf),XX(lUgh),XX(lXgh),norb,nfld)
      call DELETEMEM(lXgh,norb*norb*nfld*nfld)
      call DELETEMEM(lUgh,norb*norb*nfld*nfld)
      call DELETEMEM(lUf,norb*norb*nfld)
      end
C
C*MODULE TDXITR  *DECK CompXfgh
      subroutine CompXfgh(w,Xfgh,Uf,Ugh,Udgh,norb,nfld)
      implicit none
C     Arguments
      double precision w(3)
      integer norb,nfld
      double precision Xfgh(norb*norb,nfld,nfld,nfld)
      double precision Uf(norb*norb,nfld)
      double precision Ugh(norb*norb,nfld,nfld)
      double precision Udgh(norb*norb,nfld,nfld)
C     Variables
      integer iterm,ibterm,ieterm,ifld,jfld,kfld,iifld,jjfld,kkfld
      double precision mult,dummy
      integer iperm,nperm
      parameter (nperm=3)
      integer nterms(nperm)
      integer iflds(6,nperm*27)
      double precision mults(nperm*27)
      double precision ws(3,nperm)
      logical bgws(3,nperm)
      double precision ONE,MINUSONE
      parameter (ONE=1.0D+00,MINUSONE=-1.0D+00)
      character*1 sN
      parameter (sN='N')
C     Called subprogram
      logical GetGUDEf,GetGUXDEfg
C     Body of the subprogram
      call VCLR(Xfgh,1,norb*norb*nfld*nfld*nfld)
      call MakeGammaPerms(w,nfld,nterms,iflds,mults,ws,bgws)
C
      ieterm = 0
      dummy = 0.0d+00
      do iperm = 1,nperm
         if(bgws(1,iperm)) then
            if(.not.GetGUDEf(ws(1,iperm),
     &                       dummy,.false.,Uf,.true.,
     &                       dummy,.false.,dummy,.false.,
     &                       norb,nfld)) then
               call abrtx("CompXfgh: not GetGUDEf:1")
            end if
         end if
         if(bgws(2,iperm).or.bgws(3,iperm)) then
            if(.not.GetGUXDEfg(ws(2,iperm),ws(3,iperm),
     &                         dummy,.false.,Ugh,.true.,Udgh,.true.,
     &                         dummy,.false.,dummy,.false.,
     &                         norb,nfld)) then
               call abrtx("CompXfgh: not GetGUDEf:2")
            end if
C           form Udagger(-w2,-w3)=X(w2,w3)-U(w2,w3)
C            call VSUB(Ugh,1,Udgh,1,Udgh,1,norb*norb*nfld*nfld)
            call DAXPY(norb*norb*nfld*nfld,MINUSONE,Ugh,1,Udgh,1)
         end if
         ibterm=ieterm+1
         ieterm=nterms(iperm)
         do iterm = ibterm,ieterm
            ifld = iflds(1,iterm)
            jfld = iflds(2,iterm)
            kfld = iflds(3,iterm)
            iifld = iflds(4,iterm)
            jjfld = iflds(5,iterm)
            kkfld = iflds(6,iterm)
            mult = mults(iterm)
            call DGEMM(sN,sN,norb,norb,norb,
     &                 mult,Uf(1,ifld),norb,Ugh(1,jfld,kfld),norb,
     &                 ONE,Xfgh(1,iifld,jjfld,kkfld),norb)
            call DGEMM(sN,sN,norb,norb,norb,
     &                 -mult,Udgh(1,jfld,kfld),norb,Uf(1,ifld),norb,
     &                 ONE,Xfgh(1,iifld,jjfld,kkfld),norb)
         end do
      end do
      end
C
C*MODULE TDXITR  *DECK MakeDDfgh
      subroutine MakeDDfgh(DDfgh,w,C0,norb,nocc,nfld)
      implicit none
C     Arguments
      integer norb,nocc,nfld
      double precision w(3)
      double precision DDfgh(norb*norb,nfld,nfld,nfld)
      double precision C0(norb*norb)
C     Used common block
      DOUBLE PRECISION XX
      COMMON /FMCOM / XX(1)
C     Variables
      integer lUf,lUgh,lXgh
C     Body of the subprogram
      call NEWMEM(lUf,norb*norb*nfld)
      call NEWMEM(lUgh,norb*norb*nfld*nfld)
      call NEWMEM(lXgh,norb*norb*nfld*nfld)
      call CompDDfgh(w,DDfgh,XX(lUf),XX(lUgh),XX(lXgh),
     &               C0,norb,nocc,nfld)
      call DELETEMEM(lXgh,norb*norb*nfld*nfld)
      call DELETEMEM(lUgh,norb*norb*nfld*nfld)
      call DELETEMEM(lUf,norb*norb*nfld)
      end
C
C*MODULE TDXITR  *DECK CompDDfgh
      subroutine CompDDfgh(w,DDfgh,Uf,Ugh,Udgh,C0,norb,nocc,nfld)
      implicit none
C     Arguments
      double precision w(3)
      integer norb,nocc,nfld
      double precision DDfgh(norb*norb,nfld,nfld,nfld)
      double precision Uf(norb*norb,nfld)
      double precision Ugh(norb*norb,nfld,nfld)
      double precision Udgh(norb*norb,nfld,nfld)
      double precision C0(norb*norb)
C     Variables and constants
      integer iterm,ibterm,ieterm,ifld,jfld,kfld,iifld,jjfld,kkfld
      double precision mult,dummy
      integer iperm,nperm
      parameter (nperm=3)
      integer nterms(nperm)
      integer iflds(6,nperm*27)
      double precision mults(nperm*27)
      double precision ws(3,nperm)
      logical bgws(3,nperm)
      double precision MINUSONE
      parameter (MINUSONE=-1.0D+00)
C     Called subprogram
      logical GetGUDEf,GetGUXDEfg
C     Body of the subprogram
      call VCLR(DDfgh,1,norb*norb*nfld*nfld*nfld)
      call MakeGammaPerms(w,nfld,nterms,iflds,mults,ws,bgws)
C
      ieterm = 0
      dummy = 0.0d+00
      do iperm = 1,nperm
         if(bgws(1,iperm)) then
            if(.not.GetGUDEf(ws(1,iperm),
     &                       dummy,.false.,Uf,.true.,
     &                       dummy,.false.,dummy,.false.,
     &                       norb,nfld)) then
               call abrtx("CompDDfgh: not GetGUDEf:1")
            end if
         end if
         if(bgws(2,iperm).or.bgws(3,iperm)) then
            if(.not.GetGUXDEfg(ws(2,iperm),ws(3,iperm),
     &                         dummy,.false.,Ugh,.true.,Udgh,.true.,
     &                         dummy,.false.,dummy,.false.,
     &                         norb,nfld)) then
               call abrtx("CompDDfgh: not GetGUDEf:2")
            end if
C           form Udagger(-w2,-w3)=X(w2,w3)-U(w2,w3)
C            call VSUB(Ugh,1,Udgh,1,Udgh,1,norb*norb*nfld*nfld)
            call DAXPY(norb*norb*nfld*nfld,MINUSONE,Ugh,1,Udgh,1)
         end if
         ibterm=ieterm+1
         ieterm=nterms(iperm)
         do iterm = ibterm,ieterm
            ifld = iflds(1,iterm)
            jfld = iflds(2,iterm)
            kfld = iflds(3,iterm)
            iifld = iflds(4,iterm)
            jjfld = iflds(5,iterm)
            kkfld = iflds(6,iterm)
            mult = mults(iterm)
            call CompUsum(DDfgh(1,iifld,jjfld,kkfld),
     &                    Uf(1,ifld),Udgh(1,jfld,kfld),'o',
     &                    norb,nocc,mult)
            call CompUsum(DDfgh(1,iifld,jjfld,kkfld),
     &                    Ugh(1,jfld,kfld),Uf(1,ifld),'f',
     &                    norb,nocc,mult)
         end do
      end do
      call CompDDterm(DDfgh,C0,DDfgh,norb,nfld*nfld*nfld)
      end
C*MODULE TDXITR  *DECK MakedAlphadX
      subroutine MakedAlphadX(w1)
      implicit none
C     Arguments
      double precision w1
C     Used common blocks
      integer soli,ngeom,nocc,nvir,norb,nxyz,nmodes,nfld,ibfld,iefld
      common/OLI/soli,ngeom,nocc,nvir,norb,nxyz,nmodes,nfld,ibfld,iefld
      DOUBLE PRECISION XX
      COMMON /FMCOM / XX(1)
C     Variables and constants
      integer lE0,lC0,ldAdX
      double precision w(2)
      double precision dummy(1)
C     Called subprograms
      logical GetHf,GetGUIDEa,GetGUDEf,GetHfa,GetFFfa
C     Body of the subprogram
      w(1)=-w1
      w(2)=w1
C
      dummy(1) = 0.0d+00
      if(.not.GetGUIDEa(dummy,.false.,dummy,.false.,dummy,.false.,
     &     dummy,.false.,dummy,.false.,norb,ngeom)) then
         call abrtx("MakedAlphadX: not GetGUIDEa")
      end if
      if(.not.GetGUDEf(w1,dummy,.false.,dummy,.false.,dummy,.false.,
     &                 dummy,.false.,norb,nfld)) call MakeAlpha(w1)
      if(.not.GetHf(dummy,.false.)) call MakeHf()
      if(.not.GetHfa(dummy,.false.)) call MakeHfa()
      if(.not.GetFFfa(w1,dummy,.false.)) call MakeFFfa(w1)
C
      call NEWMEM(ldAdX,9*ngeom)
      call VCLR(XX(ldAdX),1,9*ngeom)
C
      call NEWMEM(lE0,norb)
      call NEWMEM(lC0,norb*norb)
      call GetEnergy(XX(lE0),norb)
      call GetLCAO(XX(lC0),norb)
C
      call CompdAlphadX(XX(ldAdX),w1,XX(lE0),XX(lC0),
     &                  norb,nocc,nvir,ibfld,nfld,ngeom)
C
      call DELETEMEM(lC0,norb*norb)
      call DELETEMEM(lE0,norb)
C
      call PostdAlphadX(w,XX(ldAdX))
C
      call DELETEMEM(ldAdX,9*ngeom)
      CALL TIMIT(1)
C
      end
C
C*MODULE TDXITR  *DECK CompdAlphadX
      subroutine CompdAlphadX(dAdX,w,E0,C0,
     &                        norb,nocc,nvir,ibfld,nfld,ngeom)
      implicit none
C     Arguments
      integer norb,nocc,nvir,ibfld,nfld,ngeom
      double precision w
      double precision dAdX(9,ngeom)
      double precision E0(norb),C0(norb,norb)
C     Used common blocks
      integer IR, IW, IP, IS,IPK,IDAF,NAV,IODA
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      DOUBLE PRECISION XX
      COMMON /FMCOM / XX(1)
C     Variables
      double precision dummy
      integer nparams
      integer lD,lDD,lFF,lG,lU,lX,lE,lQ,lHf,lHfa,lDg
C     Called subprograms
      logical GetFFfa,GetHf,GetHfa,GetGUDEf
C     External subprogram
      external GetDJK,PrintdAlphadXConvergedParam
C     Body of the subprogram
      nparams = nfld*ngeom
      write(iw,3000) -w,w
C     ----------
C     prepare DD
C     ----------
      call NEWMEM(lDD,norb*norb*nparams)
      call MakeDDfa(XX(lDD),w,C0,norb,nocc,nfld,ngeom)
      call MakeDDStep0
      call NEWMEM(lFF,norb*norb*nparams)
      if(.not.GetFFfa(w,XX(lFF),.true.)) then
         call abrtx("CompdAlphadX: not GetFFfa")
      end if
      call MakeDDStep1(XX(lDD),C0,E0,w,XX(lFF),.true.,
     &                 norb,nocc,nvir,nparams)
      call DELETEMEM(lFF,norb*norb*nparams)
      call NEWMEM(lQ,norb*norb*nparams)
      call MakeQfa(XX(lQ),w,E0,norb,nfld,ngeom)
      call MakeDDStep2(XX(lDD),C0,E0,w,XX(lQ),.true.,
     &                 norb,nocc,nvir,nparams)
      call DELETEMEM(lQ,norb*norb*nparams)
      call NEWMEM(lX,norb*norb*nparams)
      call MakeXfa(XX(lX),w,norb,nfld,ngeom)
      call MakeDDStep3(XX(lDD),C0,XX(lX),.true.,
     &                 norb,nocc,nparams)
      call DELETEMEM(lX,norb*norb*nparams)
C     -----
C     Solve
C     -----
      call NEWMEM(lD,norb*norb*nparams)
      call TDHFCycle(GetDJK,PrintdAlphadXConvergedParam,
     &               C0,E0,w,XX(lD),XX(lDD),
     &               norb,nocc,nvir,nparams)
C     -----------
C     complete GU
C     -----------
      lG = lDD
      call NEWMEM(lU,norb*norb*nparams)
      call NEWMEM(lE,norb*norb*nparams)
      call NEWMEM(lX,norb*norb*nparams)
      call CompleteGUXE(C0,E0,w,XX(lG),XX(lU),XX(lX),XX(lE),
     &                  norb,nocc,nparams)
C      call PrintMM(XX(lG),norb,nparams,'OLIX:Gfa')
C      call PrintMM(XX(lU),norb,nparams,'OLIX:Ufa')
C      call PrintMM(XX(lX),norb,nparams,'OLIX:Xfa')
C      call PrintMM(XX(lD),norb,nparams,'OLIX:Dfa')
C      call PrintMM(XX(lE),norb,nparams,'OLIX:Efa')
      call FreeGUXDEfa(w)
      call SetGUXDEfa(w,XX(lG),XX(lU),XX(lX),XX(lD),XX(lE),
     &                norb,nfld,ngeom)
      call DELETEMEM(lX,norb*norb*nparams)
      call DELETEMEM(lE,norb*norb*nparams)
      call DELETEMEM(lU,norb*norb*nparams)
C     --------------------
C     compute alpha tensor
C     --------------------
      lHf = lDD
      if(.not.GetHf(XX(lHF),.true.)) then
         call abrtx("CompdAlphadX: not GetHf")
      end if
      call NEWMEM(lHfa,norb*norb*nparams)
      if(.not.GetHfa(XX(lHfa),.true.)) then
         call abrtx("CompdAlphadX: not GetHfa")
      end if
      call NEWMEM(lDg,norb*norb*nfld)
      dummy = 0.0d+00
      if(.not.GetGUDEf(w,dummy,.false.,dummy,.false.,XX(lDg),.true.,
     &                 dummy,.false.,norb,nfld)) then
         call abrtx("CompdAlphadX: not GetGUDEf")
      end if
      call CompdAlphadXTensor(dadx,XX(lHf),XX(lD),XX(lHfa),XX(lDg),
     &                        norb,ibfld,nfld,ngeom)
      call DELETEMEM(lDg,norb*norb*nfld)
      call DELETEMEM(lHfa,norb*norb*nparams)
C     ----------
C     deallocate
C     ----------
      call DELETEMEM(lD,norb*norb*nparams)
      call DELETEMEM(lDD,norb*norb*nparams)
C     Format
 3000 FORMAT('Iterative procedure to compute dAlphadX(',
     &       F13.6,';',F13.6,')')
      end
C
C*MODULE TDXITR  *DECK CompdAlphadXTensor
      subroutine CompdAlphadXTensor(dadx,Hf,Dga,Hfa,Dg,
     &                              norb,ibfld,nfld,ngeom)
      implicit none
C     Arguments
      integer norb,ibfld,nfld,ngeom
      double precision dadx(9,ngeom)
      double precision Hf(norb*norb,nfld),Dga(norb*norb,nfld,ngeom)
      double precision Hfa(norb*norb,nfld,ngeom),Dg(norb*norb,nfld)
C     Used common block
      integer NELEMA,IELEMA,INDEXA
      common/NLO_ALPHA/NELEMA,IELEMA(2,9),INDEXA(3,3)
C     Variables
      integer igeom,index,idummy,ifld(2),icomp
C     Called subprograms
      double precision dTrace
C     Body of the subprogram
      do index = 1,NELEMA
         do icomp=1,2
            idummy = IELEMA(icomp,index)-ibfld+1
            if((idummy.lt.1).or.(idummy.gt.nfld)) go to 1001
            ifld(icomp) = idummy
            if(icomp.eq.2) then
               do igeom=1,ngeom
                  dadx(index,igeom)=
     &                 -dTrace('S','N',Hf(1,ifld(1)),norb,
     &                         Dga(1,ifld(2),igeom),norb,norb,norb)
     &                 -dTrace('S','N',Hfa(1,ifld(1),igeom),norb,
     &                         Dg(1,ifld(2)),norb,norb,norb)
               end do
            end if
         end do
 1001    continue
      end do
      end
C
C*MODULE TDXITR  *DECK PrintdAlphadXConvergedParam
      subroutine PrintdAlphadXConvergedParam(iw,ixyz,conv)
      implicit none
C     Arguments
      integer iw,ixyz
      double precision conv
C     Used common block
      integer soli,ngeom,nocc,nvir,norb,nxyz,nmodes,nfld,ibfld,iefld
      common/OLI/soli,ngeom,nocc,nvir,norb,nxyz,nmodes,nfld,ibfld,iefld
C     Variables and constants
      integer ifld,ipar,iat,icoord
      character*1 DIR(3)
      logical bUseQ
C     Called subprogram
      integer IGetOLI
C     Body of the subprogram
      DATA DIR /'x','y','z'/
C
      bUseQ = IGetOLI(4).ne.0
C
      ipar = (ixyz-1)/nfld
      ifld = (ixyz-1-ipar*nfld)
      ipar = ipar+1
      ifld = ifld+ibfld
      if(bUseQ) then
         write(iw,1001) DIR(ifld),ipar,conv
      else
         iat = (ipar-1)/3
         icoord = (ipar-iat*3)
         iat = iat+1
         write(iw,1002) DIR(ifld),iat,DIR(icoord),conv
      end if
C     Format
 1001 FORMAT('  Field direction ',A1,' mode ',I4,
     &       ' converged within ',  E13.6)
 1002 FORMAT('  Field direction ',A1,' coord. ',I3,A1,
     &       ' converged within ',  E13.6)
      end
C
C*MODULE TDXITR  *DECK MakeQfa
      subroutine MakeQfa(Q,w,E0,norb,nfld,ngeom)
      implicit none
C     Arguments
      integer norb,nfld,ngeom
      double precision w
      double precision Q(norb*norb,nfld,ngeom),E0(norb)
C     Used common block
      DOUBLE PRECISION XX
      COMMON /FMCOM / XX(1)
C     Variables
      integer lUa,lGa,lEa,lIa,lUf,lGf,lEf
      double precision dummy
C     Called subprogram
      logical GetGUDEf,GetGUIDEa
C     Body of the subprogram
      call NEWMEM(lUa,norb*norb*ngeom)
      call NEWMEM(lGa,norb*norb*ngeom)
      call NEWMEM(lEa,norb*norb*ngeom)
      call NEWMEM(lIa,norb*norb*ngeom)
      call NEWMEM(lUf,norb*norb*nfld)
      call NEWMEM(lGf,norb*norb*nfld)
      call NEWMEM(lEf,norb*norb*nfld)
      dummy = 0.0d+00
      if(.not.GetGUDEf(w,XX(lGf),.true.,XX(lUf),.true.,dummy,.false.,
     &                 XX(lEf),.true.,norb,nfld)) then
         call abrtx("MakeQfa: not GetGUDEf")
      end if
      if(.not.GetGUIDEa(XX(lGa),.true.,XX(lUa),.true.,XX(lIa),.true.,
     &     dummy,.false.,XX(lEa),.true.,norb,ngeom)) then
         call abrtx("MakeQfa: not GetGUIDEa")
      end if
      call CompQfa(Q,w,E0,XX(lUa),XX(lIa),XX(lGa),XX(lEa),
     &             XX(lUf),XX(lGf),XX(lEf),norb,nfld,ngeom)
      call DELETEMEM(lEf,norb*norb*nfld)
      call DELETEMEM(lGf,norb*norb*nfld)
      call DELETEMEM(lUf,norb*norb*nfld)
      call DELETEMEM(lIa,norb*norb*ngeom)
      call DELETEMEM(lEa,norb*norb*ngeom)
      call DELETEMEM(lGa,norb*norb*ngeom)
      call DELETEMEM(lUa,norb*norb*ngeom)
      end
C
C*MODULE TDXITR  *DECK CompQfa
      subroutine CompQfa(Q,w,E0,Ua,Ia,Ga,Ea,Uf,Gf,Ef,norb,nfld,ngeom)
      implicit none
C     Arguments
      integer norb,nfld,ngeom
      double precision w
      double precision Q(norb*norb,nfld,ngeom),E0(norb)
      double precision Ua(norb*norb,ngeom),Ia(norb*norb,ngeom)
      double precision Ga(norb*norb,ngeom),Ea(norb*norb,ngeom)
      double precision Uf(norb*norb,nfld),Gf(norb*norb,nfld)
      double precision Ef(norb*norb,nfld)
C     Variables and constants
      integer ifld,igeom
      double precision ONE
      parameter (ONE=1.0D+00)
      character*1 sN
      parameter (sN='N')
C     Body of the subprogram
      call VCLR(Q,1,norb*norb*nfld*ngeom)
      do igeom=1,ngeom
         do ifld=1,nfld
            call CompQterm(Q(1,ifld,igeom),Gf(1,ifld),Uf(1,ifld),'o',
     &                     Ua(1,igeom),Ea(1,igeom),norb,ONE)
            call CompQterm(Q(1,ifld,igeom),Ga(1,igeom),Ua(1,igeom),'a',
     &                     Uf(1,ifld),Ef(1,ifld),norb,ONE)
        end do
      end do
      do ifld=1,nfld
         call CompQterm2(Uf(1,ifld),w,E0,norb)
         do igeom=1,ngeom
            call DGEMM(sN,sN,norb,norb,norb,ONE,Ia(1,igeom),norb,
     &                 Uf(1,ifld),norb,ONE,Q(1,ifld,igeom),norb)
         end do
      end do
      end
C
C*MODULE TDXITR  *DECK MakeXfa
      subroutine MakeXfa(X,w,norb,nfld,ngeom)
      implicit none
C     Arguments
      integer norb,nfld,ngeom
      double precision w
      double precision X(norb*norb,nfld,ngeom)
C     Used common block
      DOUBLE PRECISION XX
      COMMON /FMCOM / XX(1)
C     Variables
      integer lUa,lUf
      double precision dummy
C     Called subprogram
      logical GetGUDEf,GetGUIDEa
C     Body of the subprogram
      call NEWMEM(lUa,norb*norb*ngeom)
      call NEWMEM(lUf,norb*norb*nfld)
      dummy = 0.0d+00
      if(.not.GetGUDEf(w,dummy,.false.,XX(lUf),.true.,dummy,.false.,
     &                 dummy,.false.,norb,nfld)) then
         call abrtx("MakeXfa: not GetGUDEf")
      end if
      if(.not.GetGUIDEa(dummy,.false.,XX(lUa),.true.,dummy,.false.,
     &     dummy,.false.,dummy,.false.,norb,ngeom)) then
         call abrtx("MakeXfa: not GetGUIDEa")
      end if
      call CompXfa(X,XX(lUa),XX(lUf),norb,nfld,ngeom)
      call DELETEMEM(lUf,norb*norb*nfld)
      call DELETEMEM(lUa,norb*norb*ngeom)
      end
C
C*MODULE TDXITR  *DECK CompXfa
      subroutine CompXfa(X,Ua,Uf,norb,nfld,ngeom)
      implicit none
C     Arguments
      integer norb,nfld,ngeom
      double precision X(norb*norb,nfld,ngeom)
      double precision Ua(norb*norb,ngeom),Uf(norb*norb,nfld)
C     Variables
      integer ifld,igeom
      double precision ZERO,ONE,MINUSONE
      parameter (ZERO=0.0D+00,ONE=1.0D+00,MINUSONE=-1.0D+00)
      character*1 sN,sT
      parameter (sN='N',sT='T')
C     Body of the subprogram
      do igeom=1,ngeom
         do ifld=1,nfld
            call DGEMM(sN,sN,norb,norb,norb,
     &                 ONE,Ua(1,igeom),norb,Uf(1,ifld),norb,
     &                 ZERO,X(1,ifld,igeom),norb)
            call DGEMM(sN,sT,norb,norb,norb,
     &                 MINUSONE,Uf(1,ifld),norb,Ua(1,igeom),norb,
     &                 ONE,X(1,ifld,igeom),norb)
         end do
      end do
      end
C
C*MODULE TDXITR  *DECK MakeDDfa
      subroutine MakeDDfa(DD,w,C0,norb,nocc,nfld,ngeom)
      implicit none
C     Arguments
      integer norb,nocc,nfld,ngeom
      double precision w
      double precision DD(norb*norb,nfld,ngeom),C0(norb*norb)
C     Used common block
      DOUBLE PRECISION XX
      COMMON /FMCOM / XX(1)
C     Variables
      integer lUa,lUf
      double precision dummy
C     Called subprogram
      logical GetGUDEf,GetGUIDEa
C     Body of the subprogram
      call NEWMEM(lUa,norb*norb*ngeom)
      call NEWMEM(lUf,norb*norb*nfld)
      dummy = 0.0d+00
      dummy = 0.0d+00
      if(.not.GetGUDEf(w,dummy,.false.,XX(lUf),.true.,dummy,.false.,
     &                 dummy,.false.,norb,nfld)) then
         call abrtx("MakeDDfa: not GetGUDEf")
      end if
      if(.not.GetGUIDEa(dummy,.false.,XX(lUa),.true.,dummy,.false.,
     &     dummy,.false.,dummy,.false.,norb,ngeom)) then
         call abrtx("MakeDDfa: not GetGUIDEa")
      end if
      call CompDDfa(DD,C0,XX(lUa),XX(lUf),norb,nocc,nfld,ngeom)
      call DELETEMEM(lUf,norb*norb*nfld)
      call DELETEMEM(lUa,norb*norb*ngeom)
      end
C
C*MODULE TDXITR  *DECK CompDDfa
      subroutine CompDDfa(DD,C0,Ua,Uf,norb,nocc,nfld,ngeom)
      implicit none
C     Arguments
      integer norb,nocc,nfld,ngeom
      double precision DD(norb*norb,nfld,ngeom),C0(norb*norb)
      double precision Ua(norb*norb,ngeom),Uf(norb*norb,nfld)
C     Variables and constant
      integer ifld,igeom
      double precision ONE
      parameter (ONE=1.0D+00)
C     Body of the subprogram
      call VCLR(DD,1,norb*norb*nfld*ngeom)
      do igeom=1,ngeom
         do ifld=1,nfld
            call CompUsum(DD(1,ifld,igeom),
     &                    Uf(1,ifld),Ua(1,igeom),'a',
     &                    norb,nocc,ONE)
            call CompUsum(DD(1,ifld,igeom),
     &                    Ua(1,igeom),Uf(1,ifld),'f',
     &                    norb,nocc,ONE)
         end do
      end do
      call CompDDterm(DD,C0,DD,norb,nfld*ngeom)
      end
C*MODULE TDXITR  *DECK MakedBetadX
      subroutine MakedBetadX(w1,w2)
      implicit none
C     Arguments
      double precision w1,w2
C     Used common blocks
      integer soli,ngeom,nocc,nvir,norb,nxyz,nmodes,nfld,ibfld,iefld
      common/OLI/soli,ngeom,nocc,nvir,norb,nxyz,nmodes,nfld,ibfld,iefld
      DOUBLE PRECISION XX
      COMMON /FMCOM / XX(1)
C     Variables and constants
      integer lE0,lC0,ldBdX
      double precision w(3),ww(2)
      double precision dummy(1)
C     Called subprograms
      logical GetHf,GetGUIDEa,GetGUDEf,GetGUXDEfg,GetHfa
      logical GetGUXDEfa
C     Body of the subprogram
      w(1)=-w1-w2
      w(2)=w1
      w(3)=w2
      ww(1)=w1
      ww(2)=w2
C
      dummy(1) = 0.0d+00
      if(.not.GetGUIDEa(dummy,.false.,dummy,.false.,dummy,.false.,
     &     dummy,.false.,dummy,.false.,norb,ngeom)) then
         call abrtx("MakedBetadX: not GetGUIDEa")
      end if
      if(.not.GetGUDEf(w1,dummy,.false.,dummy,.false.,dummy,.false.,
     &                 dummy,.false.,norb,nfld)) call MakeAlpha(w1)
      if(.not.GetGUDEf(w2,dummy,.false.,dummy,.false.,dummy,.false.,
     &                 dummy,.false.,norb,nfld)) call MakeAlpha(w2)
      if(.not.GetGUXDEfg(w1,w2,dummy,.false.,dummy,.false.,dummy,
     &                   .false.,dummy,.false.,dummy,.false.,norb,nfld))
     &     call MakeBeta(w1,w2)
      if(.not.GetHf(dummy,.false.)) call MakeHf()
      if(.not.GetHfa(dummy,.false.)) call MakeHfa()
      if(.not.GetGUXDEfa(w1,dummy,.false.,dummy,.false.,dummy,.false.,
     &                   dummy,.false.,dummy,.false.,norb,nfld,ngeom))
     &                   call MakedAlphadX(w1)
      if(.not.GetGUXDEfa(w2,dummy,.false.,dummy,.false.,dummy,.false.,
     &                   dummy,.false.,dummy,.false.,norb,nfld,ngeom))
     &                   call MakedAlphadX(w2)
C
      call NEWMEM(ldBdX,27*ngeom)
      call VCLR(XX(ldBdX),1,27*ngeom)
C
      call NEWMEM(lE0,norb)
      call NEWMEM(lC0,norb*norb)
      call GetEnergy(XX(lE0),norb)
      call GetLCAO(XX(lC0),norb)
C
      call CompdBetadX(XX(ldBdX),ww,XX(lE0),XX(lC0),
     &                 norb,nocc,nvir,ibfld,nfld,ngeom)
C
      call DELETEMEM(lC0,norb*norb)
      call DELETEMEM(lE0,norb)
C
      call PostdBetadX(w,XX(ldBdX))
C
      call DELETEMEM(ldBdX,27*ngeom)
      CALL TIMIT(1)
C
      end
C
C*MODULE TDXITR  *DECK CompdBetadX
      subroutine CompdBetadX(dBdX,w,E0,C0,
     &                       norb,nocc,nvir,ibfld,nfld,ngeom)
      implicit none
C     Arguments
      integer norb,nocc,nvir,ibfld,nfld,ngeom
      double precision w(2)
      double precision dBdX(27,ngeom)
      double precision E0(norb),C0(norb,norb)
C     Used common blocks
      integer IR, IW, IP, IS,IPK,IDAF,NAV,IODA
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      DOUBLE PRECISION XX
      COMMON /FMCOM / XX(1)
C     Variables
      double precision ws,dummy
      integer nparams
      integer lD,lDD,lFF,lG,lU,lX,lE,lQ,lHf,lHfa,lDgh
C     Called subprograms
      logical GetHf,GetHfa,GetGUXDEfg
C     External subprogram
      external GetDJK,PrintdBetadXConvergedParam
C     Body of the subprogram
      nparams = nfld*nfld*ngeom
      ws=w(1)+w(2)
      write(iw,3000) (-ws),w(1),w(2)
C     ----------
C     prepare DD
C     ----------
      call NEWMEM(lDD,norb*norb*nparams)
      call MakeDDfga(XX(lDD),w,C0,norb,nocc,nfld,ngeom)
C      call PrintMM(XX(lDD),norb,nparams,'OLIX:DDfga(0)')
      call MakeDDStep0
      call NEWMEM(lFF,norb*norb*nparams)
      call MakeFFfga(XX(lFF),w,norb,nfld,ngeom)
C      call PrintMM(XX(lFF),norb,nparams,'OLIX:MakeFFfga')
      call MakeDDStep1(XX(lDD),C0,E0,ws,XX(lFF),.true.,
     &                 norb,nocc,nvir,nparams)
C      call PrintMM(XX(lDD),norb,nparams,'OLIX:DDfga(1)')
      call DELETEMEM(lFF,norb*norb*nparams)
      call NEWMEM(lQ,norb*norb*nparams)
      call MakeQfga(XX(lQ),w,E0,norb,nfld,ngeom)
C      call PrintMM(XX(lQ),norb,nparams,'OLIX:MakeQfga')
      call MakeDDStep2(XX(lDD),C0,E0,ws,XX(lQ),.true.,
     &                 norb,nocc,nvir,nparams)
C      call PrintMM(XX(lDD),norb,nparams,'OLIX:DDfga(2)')
      call DELETEMEM(lQ,norb*norb*nparams)
      call NEWMEM(lX,norb*norb*nparams)
      call MakeXfga(XX(lX),w,norb,nfld,ngeom)
C      call PrintMM(XX(lX),norb,nparams,'OLIX:MakeXfga')
      call MakeDDStep3(XX(lDD),C0,XX(lX),.true.,
     &                 norb,nocc,nparams)
C      call PrintMM(XX(lDD),norb,nparams,'OLIX:DDfga(3)')
      call DELETEMEM(lX,norb*norb*nparams)
C     -----
C     Solve
C     -----
      call NEWMEM(lD,norb*norb*nparams)
      call TDHFCycle(GetDJK,PrintdBetadXConvergedParam,
     &               C0,E0,ws,XX(lD),XX(lDD),
     &               norb,nocc,nvir,nparams)
C     -----------
C     complete GU
C     -----------
      lG = lDD
      call NEWMEM(lU,norb*norb*nparams)
      call NEWMEM(lE,norb*norb*nparams)
      call NEWMEM(lX,norb*norb*nparams)
      call CompleteGUXE(C0,E0,ws,XX(lG),XX(lU),XX(lX),XX(lE),
     &                  norb,nocc,nparams)
C      call PrintMM(XX(lG),norb,nparams,'OLIX:Gfha')
C      call PrintMM(XX(lU),norb,nparams,'OLIX:Ufha')
C      call PrintMM(XX(lX),norb,nparams,'OLIX:Xfha')
C      call PrintMM(XX(lD),norb,nparams,'OLIX:Dfha')
C      call PrintMM(XX(lE),norb,nparams,'OLIX:Efha')
      call FreeGUXDEfga(w(1),w(2))
      call SetGUXDEfga(w(1),w(2),XX(lG),XX(lU),XX(lX),XX(lD),XX(lE),
     &                norb,nfld,ngeom)
      call DELETEMEM(lX,norb*norb*nparams)
      call DELETEMEM(lE,norb*norb*nparams)
      call DELETEMEM(lU,norb*norb*nparams)
C     -------------------
C     compute beta tensor
C     -------------------
      lHfa = lDD
      if(.not.GetHfa(XX(lHfa),.true.)) then
         call abrtx("CompdBetadX: not GetHfa")
      end if
      call NEWMEM(lHf,norb*norb*nfld)
      if(.not.GetHf(XX(lHf),.true.)) then
         call abrtx("CompdBetadX: not GetHf")
      end if
      call NEWMEM(lDgh,norb*norb*nfld*nfld)
      dummy = 0.0d+00
      if(.not.GetGUXDEfg(w(1),w(2),dummy,.false.,dummy,.false.,
     &                   dummy,.false.,XX(lDgh),.true.,dummy,.false.,
     &                   norb,nfld)) then
         call abrtx("CompdBetadX: not GetGUXDEfg")
      end if
      call CompdBetadXTensor(dbdx,XX(lHf),XX(lD),XX(lHfa),XX(lDgh),
     &                       norb,ibfld,nfld,ngeom)
      call DELETEMEM(lDgh,norb*norb*nfld*nfld)
      call DELETEMEM(lHf,norb*norb*nfld)
C     ----------
C     deallocate
C     ----------
      call DELETEMEM(lD,norb*norb*nparams)
      call DELETEMEM(lDD,norb*norb*nparams)
C     Format
 3000 FORMAT('Iterative procedure to compute dBetadX(',
     &       F13.6,';',F13.6,',',F13.6,')')
      end
C
C*MODULE TDXITR  *DECK CompdBetadXTensor
      subroutine CompdBetadXTensor(dbdx,Hf,Dgha,Hfa,Dgh,
     &                             norb,ibfld,nfld,ngeom)
      implicit none
C     Arguments
      integer norb,ibfld,nfld,ngeom
      double precision dbdx(27,ngeom)
      double precision Hf(norb*norb,nfld)
      double precision Dgha(norb*norb,nfld,nfld,ngeom)
      double precision Hfa(norb*norb,nfld,ngeom)
      double precision Dgh(norb*norb,nfld,nfld)
C     Used common block
      integer NELEMB,IELEMB,INDEXB
      common/NLO_BETA/NELEMB,IELEMB(3,27),INDEXB(3,3,3)
C     Variables
      integer igeom,index,idummy,ifld(3),icomp
C     Called subprograms
      double precision dTrace
C     Body of the subprogram
      do index = 1,NELEMB
         do icomp=1,3
            idummy = IELEMB(icomp,index)-ibfld+1
            if((idummy.lt.1).or.(idummy.gt.nfld)) go to 1001
            ifld(icomp) = idummy
            if(icomp.eq.3) then
               do igeom=1,ngeom
                  dbdx(index,igeom)=
     &                 -dTrace('S','N',Hf(1,ifld(1)),norb,
     &                         Dgha(1,ifld(2),ifld(3),igeom),norb,
     &                         norb,norb)
     &                 -dTrace('S','N',Hfa(1,ifld(1),igeom),norb,
     &                         Dgh(1,ifld(2),ifld(3)),norb,norb,norb)
               end do
            end if
         end do
 1001    continue
      end do
      end
C
C*MODULE TDXITR  *DECK PrintdBetadXConvergedParam
      subroutine PrintdBetadXConvergedParam(iw,ixyz,conv)
      implicit none
C     Arguments
      integer iw,ixyz
      double precision conv
C     Used common block
      integer soli,ngeom,nocc,nvir,norb,nxyz,nmodes,nfld,ibfld,iefld
      common/OLI/soli,ngeom,nocc,nvir,norb,nxyz,nmodes,nfld,ibfld,iefld
C     Variables and constants
      integer ifld,jfld,ipar,iat,icoord
      character*1 DIR(3)
      logical bUseQ
C     Called subprogram
      integer IGetOLI
C     Body of the subprogram
      DATA DIR /'x','y','z'/
C
      bUseQ = IGetOLI(4).ne.0
C
      ipar = (ixyz-1)/(nfld*nfld)
      ifld = (ixyz-1-ipar*nfld*nfld)/nfld
      jfld = (ixyz-1-ipar*nfld*nfld-ifld*nfld)
      ipar = ipar+1
      ifld = ifld+ibfld
      jfld = jfld+ibfld
      if(bUseQ) then
         write(iw,1001) DIR(ifld),DIR(jfld),ipar,conv
      else
         iat = (ipar-1)/3
         icoord = (ipar-iat*3)
         iat = iat+1
         write(iw,1002) DIR(ifld),DIR(jfld),iat,DIR(icoord),conv
      end if
C     Format
 1001 FORMAT('  Field directions ',A1,A1,' mode ',I4,
     &       ' converged within ',  E13.6)
 1002 FORMAT('  Field directions ',A1,A1,' coord. ',I3,A1,
     &       ' converged within ',  E13.6)
      end
C
C*MODULE TDXITR  *DECK MakeFFfga
      subroutine MakeFFfga(FF,w,norb,nfld,ngeom)
      implicit none
C     Arguments
      integer norb,nfld,ngeom
      double precision w(2)
      double precision FF(norb*norb,nfld,nfld,ngeom)
C     Used common block
      DOUBLE PRECISION XX
      COMMON /FMCOM / XX(1)
C     Variables
      integer lDfg
      double precision dummy
C     Called subprogram
      logical GetGUXDEfg
C     Body of the subprogram
      call NEWMEM(lDfg,norb*norb*nfld*nfld)
      dummy = 0.0d+00
      if(.not.GetGUXDEfg(w(1),w(2),dummy,.false.,dummy,.false.,
     &                   dummy,.false.,XX(lDfg),.true.,dummy,.false.,
     &                   norb,nfld)) then
         write(*,*) 'Could not find Ufg(',w(1),',',w(2),')'
         call abrtx("MakeFFfga: not GetGUXDEfg")
      end if
      call GetDJKa(XX(lDfg),FF,nfld*nfld)
      call DELETEMEM(lDfg,norb*norb*nfld*nfld)
      end
C
C*MODULE TDXITR  *DECK MakedBetadXPerms
      subroutine MakedBetadXPerms(w,nfld,nterms,iflds,mults,ws,bgws)
      implicit none
C     Arguments
      double precision w(2)
      integer nfld
      integer nterms(2)
      integer iflds(4,*)
      double precision mults(*)
      double precision ws(2,2)
      logical bgws(2,2)
C     Variables and constants
      logical bw12,b12
      integer iterm,ifld,jfld,iifld,jjfld
      double precision mult,ONE,TWO
      parameter (ONE=1.0D+00,TWO=2.0D+00)
      integer i,iperm,nperm
      parameter (nperm=2)
C     Body of the subprogram
C     ------------------
C     Smart permutations
C     ------------------
      bw12=w(1).eq.w(2)
      ws(1,1) = w(1)
      ws(2,1) = w(2)
      ws(1,2) = w(2)
      ws(2,2) = w(1)
      iterm = 0
      do iperm=1,nperm
         do jjfld=1,nfld
            do 1011 iifld=1,nfld
               b12=(iifld.eq.jjfld).and.bw12
C
               if(iperm.eq.1) then
                  ifld = iifld
                  jfld = jjfld
                  if(b12) then
                     mult = TWO
                  else
                     mult = ONE
                  end if
               else if(iperm.eq.2) then
                  if(b12) go to 1011
                  ifld = jjfld
                  jfld = iifld
                  mult = ONE
               end if
               iterm=iterm+1
               iflds(1,iterm) = ifld
               iflds(2,iterm) = jfld
               iflds(3,iterm) = iifld
               iflds(4,iterm) = jjfld
               mults(iterm)   = mult
 1011       continue
         end do
         nterms(iperm) = iterm
      end do
      do i = 1,2
         bgws(i,1) = .true.
         bgws(i,2) = ws(i,1).ne.ws(i,2)
      end do
      end
C
C*MODULE TDXITR  *DECK MakeDDfga
      subroutine MakeDDfga(DD,w,C0,norb,nocc,nfld,ngeom)
      implicit none
C     Arguments
      integer norb,nocc,nfld,ngeom
      double precision w(2)
      double precision C0(norb*norb),DD(norb*norb,nfld,nfld,ngeom)
C     Used common block
      DOUBLE PRECISION XX
      COMMON /FMCOM / XX(1)
C     Variables
      integer lUf,lUga,lXga,lUa,lUfg,lXfg
C     Body of the subprogram
      call VCLR(DD,1,norb*norb*nfld*nfld*ngeom)
      call NEWMEM(lUa,norb*norb*ngeom)
      call NEWMEM(lUfg,norb*norb*nfld*nfld)
      call NEWMEM(lXfg,norb*norb*nfld*nfld)
      call CompDDfga1(DD,w,XX(lUa),XX(lUfg),XX(lXfg),
     &                norb,nocc,nfld,ngeom)
      call DELETEMEM(lXfg,norb*norb*nfld*nfld)
      call DELETEMEM(lUfg,norb*norb*nfld*nfld)
      call DELETEMEM(lUa,norb*norb*ngeom)
      call NEWMEM(lUf,norb*norb*nfld)
      call NEWMEM(lUga,norb*norb*nfld*ngeom)
      call NEWMEM(lXga,norb*norb*nfld*ngeom)
      call CompDDfga2(DD,w,XX(lUf),XX(lUga),XX(lXga),
     &                norb,nocc,nfld,ngeom)
      call DELETEMEM(lXga,norb*norb*nfld*ngeom)
      call DELETEMEM(lUga,norb*norb*nfld*ngeom)
      call DELETEMEM(lUf,norb*norb*nfld)
      call CompDDterm(DD,C0,DD,norb,nfld*nfld*ngeom)
      end
C
C*MODULE TDXITR  *DECK CompDDfga1
      subroutine CompDDfga1(DD,w,Ua,Ufg,Udfg,norb,nocc,nfld,ngeom)
      implicit none
C     Arguments
      integer norb,nocc,nfld,ngeom
      double precision w(2)
      double precision DD(norb*norb,nfld,nfld,ngeom)
      double precision Ua(norb*norb,ngeom)
      double precision Ufg(norb*norb,nfld,nfld)
      double precision Udfg(norb*norb,nfld,nfld)
C     Variables
      double precision dummy
      integer ifld,jfld,igeom
      double precision ONE,MINUSONE
      parameter (ONE=1.0D+00,MINUSONE=-1.0D+00)
C     Called subprogram
      logical GetGUIDEa,GetGUXDEfg
C     Body of the subprogram
      dummy = 0.0d+00
      if(.not.GetGUXDEfg(w(1),w(2),dummy,.false.,Ufg,.true.,
     &                   Udfg,.true.,dummy,.false.,dummy,.false.,
     &                   norb,nfld)) then
         write(*,*) 'Could not find Ufg(',w(1),',',w(2),')'
         call abrtx("MakeDDfga1: not GetGUXDEfg")
      end if
      call DAXPY(norb*norb*nfld*nfld,MINUSONE,Ufg,1,Udfg,1)
      if(.not.GetGUIDEa(dummy,.false.,Ua,.true.,dummy,.false.,
     &                  dummy,.false.,dummy,.false.,norb,ngeom)) then
         write(*,*) 'Could not find Ua'
         call abrtx("MakeDDfga1: not GetGUIDEa")
      end if
      do igeom=1,ngeom
         do jfld=1,nfld
            do ifld=1,nfld
               call CompUsum(DD(1,ifld,jfld,igeom),
     &                       Ufg(1,ifld,jfld),Ua(1,igeom),'a',
     &                       norb,nocc,ONE)
               call CompUsum(DD(1,ifld,jfld,igeom),
     &                       Ua(1,igeom),Udfg(1,ifld,jfld),'o',
     &                       norb,nocc,ONE)
            end do
         end do
      end do
      end
C
C*MODULE TDXITR  *DECK CompDDfga2
      subroutine CompDDfga2(DD,w,Uf,Uga,Udga,norb,nocc,nfld,ngeom)
      implicit none
C     Arguments
      integer norb,nocc,nfld,ngeom
      double precision w(2)
      double precision DD(norb*norb,nfld,nfld,ngeom)
      double precision Uf(norb*norb,nfld)
      double precision Uga(norb*norb,nfld,ngeom)
      double precision Udga(norb*norb,nfld,ngeom)
C     Variables
      double precision dummy
      integer iterm,ibterm,ieterm,ifld,jfld,iifld,jjfld,igeom
      double precision mult
      integer iperm,nperm
      parameter (nperm=2)
      integer nterms(nperm)
      integer iflds(4,nperm*9)
      double precision mults(nperm*9)
      double precision ws(2,nperm)
      logical bgws(2,nperm)
      double precision MINUSONE
      parameter (MINUSONE=-1.0D+00)
C     Called subprogram
      logical GetGUDEf,GetGUXDEfa
C     Body of the subprogram
      call MakedBetadXPerms(w,nfld,nterms,iflds,mults,ws,bgws)
      ieterm = 0
      dummy = 0.0d+00
      do iperm = 1,nperm
         if(bgws(1,iperm)) then
            if(.not.GetGUDEf(ws(1,iperm),
     &                       dummy,.false.,Uf,.true.,
     &                       dummy,.false.,dummy,.false.,
     &                       norb,nfld)) then
               write(*,*) 'Could not find Uf(',ws(1,iperm),')'
               call abrtx("MakeDDfga2: not GetGUDEf")
            end if
         end if
         if(bgws(2,iperm)) then
            if(.not.GetGUXDEfa(ws(2,iperm),
     &                         dummy,.false.,Uga,.true.,Udga,.true.,
     &                         dummy,.false.,dummy,.false.,
     &                         norb,nfld,ngeom)) then
               write(*,*) 'Could not find Ufa(',ws(2,iperm),')'
               call abrtx("MakeDDfga2: not GetGUXDEfa")
            end if
            call DAXPY(norb*norb*nfld*ngeom,MINUSONE,Uga,1,Udga,1)
         end if
         ibterm=ieterm+1
         ieterm=nterms(iperm)
         do iterm = ibterm,ieterm
            ifld = iflds(1,iterm)
            jfld = iflds(2,iterm)
            iifld = iflds(3,iterm)
            jjfld = iflds(4,iterm)
            mult = mults(iterm)
            do igeom=1,ngeom
               call CompUsum(DD(1,iifld,jjfld,igeom),
     &                       Uga(1,jfld,igeom),Uf(1,ifld),'f',
     &                       norb,nocc,mult)
               call CompUsum(DD(1,iifld,jjfld,igeom),
     &                       Uf(1,ifld),Udga(1,jfld,igeom),'o',
     &                       norb,nocc,mult)
            end do
         end do
      end do
      end
C
C*MODULE TDXITR  *DECK MakeQfga
      subroutine MakeQfga(Q,w,E0,norb,nfld,ngeom)
      implicit none
C     Arguments
      integer norb,nfld,ngeom
      double precision w(2)
      double precision E0(norb),Q(norb*norb,nfld,nfld,ngeom)
C     Used common block
      DOUBLE PRECISION XX
      COMMON /FMCOM / XX(1)
C     Variables
      integer lGfg,lUfg,lEfg,lGa,lUa,lIa,lEa
      integer lGf,lUf,lEf,lGga,lUga,lEga,lUg,lEg
C     Body of the subprogram
      call VCLR(Q,1,norb*norb*nfld*nfld*ngeom)
      call NEWMEM(lGa,norb*norb*ngeom)
      call NEWMEM(lUa,norb*norb*ngeom)
      call NEWMEM(lEa,norb*norb*ngeom)
      call NEWMEM(lGfg,norb*norb*nfld*nfld)
      call NEWMEM(lUfg,norb*norb*nfld*nfld)
      call NEWMEM(lEfg,norb*norb*nfld*nfld)
      call CompQfga1(Q,w,XX(lGa),XX(lUa),XX(lEa),
     &               XX(lGfg),XX(lUfg),XX(lEfg),
     &               norb,nfld,ngeom)
      call DELETEMEM(lEfg,norb*norb*nfld*nfld)
      call DELETEMEM(lUfg,norb*norb*nfld*nfld)
      call DELETEMEM(lGfg,norb*norb*nfld*nfld)
      call DELETEMEM(lEa,norb*norb*ngeom)
      call DELETEMEM(lUa,norb*norb*ngeom)
      call DELETEMEM(lGa,norb*norb*ngeom)
      call NEWMEM(lGf,norb*norb*nfld)
      call NEWMEM(lUf,norb*norb*nfld)
      call NEWMEM(lEf,norb*norb*nfld)
      call NEWMEM(lGga,norb*norb*nfld*ngeom)
      call NEWMEM(lUga,norb*norb*nfld*ngeom)
      call NEWMEM(lEga,norb*norb*nfld*ngeom)
      call CompQfga2(Q,w,XX(lGf),XX(lUf),XX(lEf),
     &               XX(lGga),XX(lUga),XX(lEga),
     &               norb,nfld,ngeom)
      call DELETEMEM(lEga,norb*norb*nfld*ngeom)
      call DELETEMEM(lUga,norb*norb*nfld*ngeom)
      call DELETEMEM(lGga,norb*norb*nfld*ngeom)
      call DELETEMEM(lEf,norb*norb*nfld)
      call DELETEMEM(lUf,norb*norb*nfld)
      call DELETEMEM(lGf,norb*norb*nfld)
      call NEWMEM(lIa,norb*norb*ngeom)
      call NEWMEM(lUfg,norb*norb*nfld*nfld)
      call NEWMEM(lUf,norb*norb*nfld)
      call NEWMEM(lEf,norb*norb*nfld)
      call NEWMEM(lUg,norb*norb*nfld)
      call NEWMEM(lEg,norb*norb*nfld)
      call CompQfga3(Q,w,E0,XX(lIa),
     &               XX(lUfg),XX(lUf),XX(lEf),XX(lUg),XX(lEg),
     &               norb,nfld,ngeom)
      call DELETEMEM(lEg,norb*norb*nfld)
      call DELETEMEM(lUg,norb*norb*nfld)
      call DELETEMEM(lEf,norb*norb*nfld)
      call DELETEMEM(lUf,norb*norb*nfld)
      call DELETEMEM(lUfg,norb*norb*nfld*nfld)
      call DELETEMEM(lIa,norb*norb*ngeom)
      end
C
C*MODULE TDXITR  *DECK CompQfga1
      subroutine CompQfga1(Q,w,Ga,Ua,Ea,Gfg,Ufg,Efg,
     &                     norb,nfld,ngeom)
      implicit none
C     Arguments
      integer norb,nfld,ngeom
      double precision w(2)
      double precision Q(norb*norb,nfld,nfld,ngeom)
      double precision Ga(norb*norb,ngeom),Gfg(norb*norb,nfld,nfld)
      double precision Ua(norb*norb,ngeom),Ufg(norb*norb,nfld,nfld)
      double precision Ea(norb*norb,ngeom),Efg(norb*norb,nfld,nfld)
C     Variables
      double precision dummy
      integer ifld,jfld,igeom
      double precision ONE
      parameter (ONE=1.0D+00)
C     Called subprogram
      logical GetGUIDEa,GetGUXDEfg
C     Body of the subprogram
      dummy = 0.0d+00
      if(.not.GetGUXDEfg(w(1),w(2),Gfg,.true.,Ufg,.true.,
     &                   dummy,.false.,dummy,.false.,Efg,.true.,
     &                   norb,nfld)) then
         write(*,*) 'Could not find Ufg(',w(1),',',w(2),')'
         call abrtx("CompQfga1: not GetGUXDEfg")
      end if
      if(.not.GetGUIDEa(Ga,.true.,Ua,.true.,dummy,.false.,
     &                  dummy,.false.,Ea,.true.,norb,ngeom)) then
         write(*,*) 'Could not find Ua)'
         call abrtx("CompQfga1: not GetGUIDEa")
      end if
      do igeom=1,ngeom
         do jfld=1,nfld
            do ifld=1,nfld
               call CompQterm(Q(1,ifld,jfld,igeom),
     &              Ga(1,igeom),Ua(1,igeom),'a',
     &              Ufg(1,ifld,jfld),Efg(1,ifld,jfld),norb,ONE)
               call CompQterm(Q(1,ifld,jfld,igeom),
     &              Gfg(1,ifld,jfld),Ufg(1,ifld,jfld),'o',
     &              Ua(1,igeom),Ea(1,igeom),norb,ONE)
            end do
         end do
      end do
      end
C
C*MODULE TDXITR  *DECK CompQfga2
      subroutine CompQfga2(Q,w,Gf,Uf,Ef,Gga,Uga,Ega,
     &                     norb,nfld,ngeom)
      implicit none
C     Arguments
      integer norb,nfld,ngeom
      double precision w(2)
      double precision Q(norb*norb,nfld,nfld,ngeom)
      double precision Gf(norb*norb,nfld),Gga(norb*norb,nfld,ngeom)
      double precision Uf(norb*norb,nfld),Uga(norb*norb,nfld,ngeom)
      double precision Ef(norb*norb,nfld),Ega(norb*norb,nfld,ngeom)
C     Variables
      double precision dummy
      integer iterm,ibterm,ieterm,ifld,jfld,iifld,jjfld,igeom
      double precision mult
      integer iperm,nperm
      parameter (nperm=2)
      integer nterms(nperm)
      integer iflds(4,nperm*9)
      double precision mults(nperm*9)
      double precision ws(2,nperm)
      logical bgws(2,nperm)
C     Called subprogram
      logical GetGUDEf,GetGUXDEfa
C     Body of the subprogram
      call MakedBetadXPerms(w,nfld,nterms,iflds,mults,ws,bgws)
      ieterm = 0
      dummy = 0.0d+00
      do iperm = 1,nperm
         if(bgws(1,iperm)) then
            if(.not.GetGUDEf(ws(1,iperm),
     &                       Gf,.true.,Uf,.true.,
     &                       dummy,.false.,Ef,.true.,
     &                       norb,nfld)) then
               write(*,*) 'Could not find Uf(',ws(1,iperm),')'
               call abrtx("CompQfga2: not GetGUDEf")
            end if
         end if
         if(bgws(2,iperm)) then
            if(.not.GetGUXDEfa(ws(2,iperm),
     &                         Gga,.true.,Uga,.true.,dummy,.false.,
     &                         dummy,.false.,Ega,.true.,
     &                         norb,nfld,ngeom)) then
               write(*,*) 'Could not find Ufa(',ws(2,iperm),')'
               call abrtx("CompQfga2: not GetGUXDEfa")
            end if
         end if
         ibterm=ieterm+1
         ieterm=nterms(iperm)
         do iterm = ibterm,ieterm
            ifld = iflds(1,iterm)
            jfld = iflds(2,iterm)
            iifld = iflds(3,iterm)
            jjfld = iflds(4,iterm)
            mult = mults(iterm)
            do igeom=1,ngeom
               call CompQterm(Q(1,iifld,jjfld,igeom),
     &              Gf(1,ifld),Uf(1,ifld),'o',
     &              Uga(1,jfld,igeom),Ega(1,jfld,igeom),norb,mult)
               call CompQterm(Q(1,iifld,jjfld,igeom),
     &              Gga(1,jfld,igeom),Uga(1,jfld,igeom),'o',
     &              Uf(1,ifld),Ef(1,ifld),norb,mult)
            end do
         end do
      end do
      end
C
C*MODULE TDXITR  *DECK CompQfga3
      subroutine CompQfga3(Q,w,E0,Ia,Ufg,Uf,Ef,Ug,Eg,
     &                     norb,nfld,ngeom)
      implicit none
C     Arguments
      integer norb,nfld,ngeom
      double precision w(2)
      double precision Q(norb*norb,nfld,nfld,ngeom),E0(norb)
      double precision Ia(norb*norb,ngeom),Ufg(norb*norb,nfld,nfld)
      double precision Uf(norb*norb,nfld),Ug(norb*norb,nfld)
      double precision Ef(norb*norb,nfld),Eg(norb*norb,nfld)
C     Variables
      double precision dummy,ws
      integer ifld,jfld,igeom
      double precision ONE,MINUSONE
      parameter (ONE=1.0D+00,MINUSONE=-1.0D+00)
      character*1 sN
      parameter (sN='N')
C     Called subprogram
      logical GetGUIDEa,GetGUDEf,GetGUXDEfg
C     Body of the subprogram
      ws=w(1)+w(2)
      dummy = 0.0d+00
      if(.not.GetGUXDEfg(w(1),w(2),dummy,.false.,Ufg,.true.,
     &                   dummy,.false.,dummy,.false.,dummy,.false.,
     &                   norb,nfld)) then
         write(*,*) 'Could not find Ufg(',w(1),',',w(2),')'
         call abrtx("CompQfga3: not GetGUXDEfg")
      end if
      if(.not.GetGUDEf(w(1),dummy,.false.,Uf,.true.,
     &                 dummy,.false.,Ef,.true.,norb,nfld)) then
         write(*,*) 'Could not find Uf(',w(1),')'
         call abrtx("CompQfga3: not GetGUDEf:1")
      end if
      if(.not.GetGUDEf(w(2),dummy,.false.,Ug,.true.,
     &                 dummy,.false.,Eg,.true.,norb,nfld)) then
         write(*,*) 'Could not find Uf(',w(2),')'
         call abrtx("CompQfga3: not GetGUDEf:2")
      end if
      if(.not.GetGUIDEa(dummy,.false.,dummy,.false.,Ia,.true.,
     &                  dummy,.false.,dummy,.false.,norb,ngeom)) then
         write(*,*) 'Could not find Ua'
         call abrtx("CompQfga3: not GetGUIDEa")
      end if
      do jfld=1,nfld
         do ifld=1,nfld
            call CompQterm2(Ufg(1,ifld,jfld),ws,E0,norb)
            call DGEMM(sN,sN,norb,norb,norb,
     &                 MINUSONE,Uf(1,ifld),norb,Eg(1,jfld),norb,
     &                 ONE,Ufg(1,ifld,jfld),norb)
            call DGEMM(sN,sN,norb,norb,norb,
     &                 MINUSONE,Ug(1,jfld),norb,Ef(1,ifld),norb,
     &                 ONE,Ufg(1,ifld,jfld),norb)
            do igeom=1,ngeom
               call DGEMM(sN,sN,norb,norb,norb,
     &                    ONE,Ia(1,igeom),norb,Ufg(1,ifld,jfld),norb,
     &                    ONE,Q(1,ifld,jfld,igeom),norb)
            end do
         end do
      end do
      end
C
C*MODULE TDXITR  *DECK MakeXfga
      subroutine MakeXfga(X,w,norb,nfld,ngeom)
      implicit none
C     Arguments
      integer norb,nfld,ngeom
      double precision w(2)
      double precision X(norb*norb,nfld,nfld,ngeom)
C     Used common block
      DOUBLE PRECISION XX
      COMMON /FMCOM / XX(1)
C     Variables
      integer lUfg,lUa,lUf,lUg,lUga
C     Body of the subprogram
      call VCLR(X,1,norb*norb*nfld*nfld*ngeom)
      call NEWMEM(lUa,norb*norb*ngeom)
      call NEWMEM(lUfg,norb*norb*nfld*nfld)
      call CompXfga1(X,w,XX(lUa),XX(lUfg),norb,nfld,ngeom)
      call DELETEMEM(lUfg,norb*norb*nfld*nfld)
      call NEWMEM(lUf,norb*norb*nfld)
      call NEWMEM(lUg,norb*norb*nfld)
      call CompXfga2(X,w,XX(lUa),XX(lUf),XX(lUg),norb,nfld,ngeom)
      call DELETEMEM(lUg,norb*norb*nfld)
      call DELETEMEM(lUf,norb*norb*nfld)
      call DELETEMEM(lUa,norb*norb*ngeom)
      call NEWMEM(lUf,norb*norb*nfld)
      call NEWMEM(lUga,norb*norb*nfld*ngeom)
      call CompXfga3(X,w,XX(lUf),XX(lUga),norb,nfld,ngeom)
      call DELETEMEM(lUga,norb*norb*nfld*ngeom)
      call DELETEMEM(lUf,norb*norb*nfld)
      end
C
C*MODULE TDXITR  *DECK CompXfga1
      subroutine CompXfga1(X,w,Ua,Ufg,norb,nfld,ngeom)
      implicit none
C     Arguments
      integer norb,nfld,ngeom
      double precision w(2)
      double precision X(norb*norb,nfld,nfld,ngeom)
      double precision Ua(norb*norb,ngeom),Ufg(norb*norb,nfld,nfld)
C     Variables
      double precision dummy
      integer ifld,jfld,igeom
      double precision ZERO,ONE,MINUSONE
      parameter (ZERO=0.0D+00,ONE=1.0D+00,MINUSONE=-1.0D+00)
      character*1 sN,sT
      parameter (sN='N',sT='T')
C     Called subprogram
      logical GetGUIDEa,GetGUXDEfg
C     Body of the subprogram
      dummy = 0.0d+00
      if(.not.GetGUXDEfg(w(1),w(2),dummy,.false.,Ufg,.true.,
     &                   dummy,.false.,dummy,.false.,dummy,.false.,
     &                   norb,nfld)) then
         write(*,*) 'Could not find Ufg(',w(1),',',w(2),')'
         call abrtx("CompXfga1: not GetGUXDEfg")
      end if
      if(.not.GetGUIDEa(dummy,.false.,Ua,.true.,dummy,.false.,
     &                  dummy,.false.,dummy,.false.,norb,ngeom)) then
         write(*,*) 'Could not find Ua'
         call abrtx("CompXfga1: not GetGUIDEa")
      end if
      do igeom=1,ngeom
         do jfld=1,nfld
            do ifld=1,nfld
               call DGEMM(sN,sN,norb,norb,norb,
     &              ONE,Ua(1,igeom),norb,Ufg(1,ifld,jfld),norb,
     &              ZERO,X(1,ifld,jfld,igeom),norb)
               call DGEMM(sN,sT,norb,norb,norb,
     &              MINUSONE,Ufg(1,ifld,jfld),norb,Ua(1,igeom),norb,
     &              ONE,X(1,ifld,jfld,igeom),norb)
            end do
         end do
      end do
      end
C
C*MODULE TDXITR  *DECK CompXfga2
      subroutine CompXfga2(X,w,Ua,Uf,Ug,norb,nfld,ngeom)
      implicit none
C     Arguments
      integer norb,nfld,ngeom
      double precision w(2)
      double precision X(norb*norb,nfld,nfld,ngeom)
      double precision Ua(norb*norb,ngeom)
      double precision Uf(norb*norb,nfld),Ug(norb*norb,nfld)
C     Used common block
      DOUBLE PRECISION XX
      COMMON /FMCOM / XX(1)
C     Variables
      integer lWRK1,lWRK2
      double precision dummy
      integer iterm,ibterm,ieterm,ifld,jfld,iifld,jjfld,igeom
      double precision mult
      integer iperm,nperm
      parameter (nperm=2)
      integer nterms(nperm)
      integer iflds(4,nperm*9)
      double precision mults(nperm*9)
      double precision ws(2,nperm)
      logical bgws(2,nperm)
      double precision ZERO,ONE,MINUSONE
      parameter (ZERO=0.0D+00,ONE=1.0D+00,MINUSONE=-1.0D+00)
      character*1 sN,sT
      parameter (sN='N',sT='T')
C     Called subprogram
      logical GetGUDEf
C     Body of the subprogram
      call MakedBetadXPerms(w,nfld,nterms,iflds,mults,ws,bgws)
C
      call NEWMEM(lWRK1,norb*norb)
      call NEWMEM(lWRK2,norb*norb)
C
      ieterm = 0
      dummy = 0.0d+00
      do iperm = 1,nperm
         if(bgws(1,iperm)) then
            if(.not.GetGUDEf(ws(1,iperm),
     &                       dummy,.false.,Uf,.true.,
     &                       dummy,.false.,dummy,.false.,
     &                       norb,nfld)) then
               write(*,*) 'Could not find Uf(',ws(1,iperm),')'
               call abrtx("CompXfga2: not GetGUDEf:1")
            end if
         end if
         if(bgws(2,iperm)) then
            if(.not.GetGUDEf(ws(2,iperm),
     &                       dummy,.false.,Ug,.true.,
     &                       dummy,.false.,dummy,.false.,
     &                       norb,nfld)) then
               write(*,*) 'Could not find Uf(',ws(2,iperm),')'
               call abrtx("CompXfga2: not GetGUDEf:2")
            end if
         end if
         ibterm=ieterm+1
         ieterm=nterms(iperm)
         do iterm = ibterm,ieterm
            ifld = iflds(1,iterm)
            jfld = iflds(2,iterm)
            iifld = iflds(3,iterm)
            jjfld = iflds(4,iterm)
            mult = mults(iterm)
            call DGEMM(sN,sN,norb,norb,norb,
     &           mult,Uf(1,ifld),norb,Ug(1,jfld),norb,
     &           ZERO,XX(lWRK1),norb)
            do igeom=1,ngeom
               call DGEMM(sN,sT,norb,norb,norb,
     &              ONE,XX(lWRK1),norb,Ua(1,igeom),norb,
     &              ONE,X(1,iifld,jjfld,igeom),norb)
               call DGEMM(sN,sN,norb,norb,norb,
     &              MINUSONE,Ua(1,igeom),norb,XX(lWRK1),norb,
     &              ONE,X(1,iifld,jjfld,igeom),norb)
               call DGEMM(sN,sN,norb,norb,norb,
     &              mult,Ua(1,igeom),norb,Ug(1,jfld),norb,
     &              ZERO,XX(lWRK2),norb)
               call DGEMM(sN,sN,norb,norb,norb,
     &              MINUSONE,Uf(1,ifld),norb,XX(lWRK2),norb,
     &              ONE,X(1,iifld,jjfld,igeom),norb)
            end do
         end do
      end do
C
      call DELETEMEM(lWRK2,norb*norb)
      call DELETEMEM(lWRK1,norb*norb)
C
      end
C
C*MODULE TDXITR  *DECK CompXfga3
      subroutine CompXfga3(X,w,Uf,Uga,norb,nfld,ngeom)
      implicit none
C     Arguments
      integer norb,nfld,ngeom
      double precision w(2)
      double precision X(norb*norb,nfld,nfld,ngeom)
      double precision Uf(norb*norb,nfld),Uga(norb*norb,nfld,ngeom)
C     Variables
      double precision dummy
      integer iterm,ibterm,ieterm,ifld,jfld,iifld,jjfld,igeom
      double precision mult
      integer iperm,nperm
      parameter (nperm=2)
      integer nterms(nperm)
      integer iflds(4,nperm*9)
      double precision mults(nperm*9)
      double precision ws(2,nperm)
      logical bgws(2,nperm)
      double precision ONE
      parameter (ONE=1.0D+00)
      character*1 sN
      parameter (sN='N')
C     Called subprogram
      logical GetGUDEf,GetGUXDEfa
C     Body of the subprogram
      call MakedBetadXPerms(w,nfld,nterms,iflds,mults,ws,bgws)
C
      ieterm = 0
      dummy = 0.0d+00
      do iperm = 1,nperm
         if(bgws(1,iperm)) then
            if(.not.GetGUDEf(ws(1,iperm),
     &                       dummy,.false.,Uf,.true.,
     &                       dummy,.false.,dummy,.false.,
     &                       norb,nfld)) then
               write(*,*) 'Could not find Uf(',ws(1,iperm),')'
               call abrtx("CompXfga3: not GetGUDEf")
            end if
         end if
         if(bgws(2,iperm)) then
            if(.not.GetGUXDEfa(ws(2,iperm),
     &                         dummy,.false.,Uga,.true.,dummy,.false.,
     &                         dummy,.false.,dummy,.false.,
     &                         norb,nfld,ngeom)) then
               write(*,*) 'Could not find Ufa(',ws(2,iperm),')'
               call abrtx("CompXfga3: not GetGUXDEfa")
            end if
         end if
         ibterm=ieterm+1
         ieterm=nterms(iperm)
         do iterm = ibterm,ieterm
            ifld = iflds(1,iterm)
            jfld = iflds(2,iterm)
            iifld = iflds(3,iterm)
            jjfld = iflds(4,iterm)
            mult = mults(iterm)
            do igeom=1,ngeom
               call DGEMM(sN,sN,norb,norb,norb,
     &              mult,Uf(1,ifld),norb,Uga(1,jfld,igeom),norb,
     &              ONE,X(1,iifld,jjfld,igeom),norb)
               call DGEMM(sN,sN,norb,norb,norb,
     &              mult,Uga(1,jfld,igeom),norb,Uf(1,ifld),norb,
     &              ONE,X(1,iifld,jjfld,igeom),norb)
            end do
         end do
      end do
C
      end
C*MODULE TDXITR  *DECK SLVLEQSVD
      subroutine SLVLEQSVD(A,X,N,LDA,IERR)
      implicit none
C     Arguments
      integer N,LDA,IERR
      double precision A(LDA,N), X(N)
C     Used common blocks
      DOUBLE PRECISION XX
      COMMON /FMCOM / XX(1)
C     Variables and Constants
      double precision dummy
      integer i,lW,lwork,lU,lVT,lS
      double precision ZERO,ONE,TOL
      PARAMETER (ZERO=0.0D+00,ONE=1.0D+00,TOL=1D-16)
C     body of the subroutine
      call NEWMEM(lU,n*n)
      call NEWMEM(lVT,n*n)
      call NEWMEM(lS,n)
C     Retrieve right working size
      call dgesvd('S','S',n,n,A,LDA,XX(lS),XX(lU),n,XX(lVt),n,dummy,-1,
     *            ierr)
      lwork=int(dummy)
      if(ierr.ne.0) return
      call NEWMEM(lw,lwork)
C     Factorize
      call dgesvd('S','S',n,n,A,LDA,XX(lS),XX(lU),n,XX(lVt),n,xx(lW),
     *            lwork,ierr)
      if(ierr.ne.0) return
      call DELETEMEM(lw,lwork)
C     Solve
      call NEWMEM(lw,n)
      call dgemv('T',n,n,ONE,xx(lU),n,X,1,ZERO,xx(lW),1)
C     Pseudo invert
      do i=0,n-1
         if(xx(lS+i).lt.tol) then
            xx(lW+i) = ZERO
         else
            xx(lW+i) = xx(lW+i)/xx(lS+i)
         end if
      end do
      call dgemv('T',n,n,ONE,XX(lVt),n,XX(lW),1,ZERO,X,1)
C     Free memory
      call DELETEMEM(lw,n)
      call DELETEMEM(lS,n)
      call DELETEMEM(lVT,n*n)
      call DELETEMEM(lU,n*n)
      end
