C 19 Sep 05 - MWS - add true nuclear charge array to INFOA common
C  5 Jul 05 - MWS - SELECT NEW ATOM,BASIS,EFP,PCM,DAF DIMENSIONS
C  5 Feb 05 - OQ  - Noniterative procedures for extended TDHF properties
C
C*MODULE TDXNI   *DECK CompNIUsum
      subroutine CompNIUsum(Usum,Ux,Tx,Uy,Ty,norb,nocc,mult,bClear)
      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 Tx,Ty
      logical bClear
C     ------------------------------------------------
C     mult * op(Ux) * op(Uy) = Uxy
C     op depends on the type of Ux and Uy (Tx and Ty).
C
C     Tx = 'a', Ux=Ua     : UaUdy         =  UaUdy
C     Tx = 'f', Ux=Uf(w)  : Uf(w)Udy      =  Uf(w)Udy
C     Tx = 'm', Ux=Uf(w)  : Uf(-w)Udy     = -Udf(w)Udy
C     Tx = 'o', Ux=Ux(w)  : Ux(w)Udy      =  Ux(w)Udy
C
C     Ty = 'a', Uy=Ua     : UxUda         =  UxUda
C     Ty = 'f', Uy=Uf(w)  : UxUdf(-w)     = -UxUf(w)
C     Ty = 'm', Uy=Uf(w)  : UxUdf(w)      =  UxUdf(w)
C     Ty = 'o', Uy=Udy(-w): UxUdy(-w)     =  UxUdy(-w)
C     ------------------------------------------------
C     Variables
      double precision dummy,clear
      double precision ZERO,ONE
      parameter (ZERO=0.0d+00,ONE=1.0d+00)
      character*1 sN,sT,sY,sX
      parameter (sN='N',sT='T')
C     Body of the subprogram
      dummy = mult
      if(Tx.eq.'a') then
         sX = sN
      else if(Tx.eq.'f') then
         sX = sN
      else if(Tx.eq.'m') then
         dummy = -dummy
         sX = sT
      else if(Tx.eq.'o') then
         sX = sN
      else
         write(*,*) 'Wrong type for Tx [',Tx,'] in CompNIUsum'
         call abrtx("CompNIUsum:1")
      end if
      if(Ty.eq.'a') then
         sY = sT
      else if(Ty.eq.'f') then
         dummy = -dummy
         sY = sN
      else if(Ty.eq.'m') then
         sY = sT
      else if(Ty.eq.'o') then
         sY = sN
      else
         write(*,*) 'Wrong type for Ty [',Ty,'] in CompNIUsum'
         call abrtx("CompNIUsum:2")
      end if
      if(bClear) clear=ZERO
      if(.not.bClear) clear=ONE
      call DGEMM(sX,sY,norb,norb,nocc,
     &           dummy,Ux,norb,Uy,norb,clear,Usum,norb)
C
      end
C
C*MODULE TDXNI   *DECK CompNIXsum
      subroutine CompNIXsum(Xsum,Ux,Tx,Uy,Ty,norb,nocc,mult,bClear)
      implicit none
C     Arguments
      integer norb,nocc
      double precision mult
      double precision Xsum(nocc*nocc)
      double precision Ux(norb*norb)
      double precision Uy(norb*norb)
      character*1 Tx,Ty
      logical bClear
C     -------------------------------------------------
C     -mult * op(Ux) * op(Uy) = Xxy
C     op depends on the types of Ux and Uy (Tx and Ty)
C
C     Tx = 'a', Ux=Ua     : IaUy + UdaUy  = -UaUy
C     Tx = 'f', Ux=Uf(w)  : Udf(-w)Uy     = -Uf(w)Uy
C     Tx = 'm', Ux=Uf(w)  : Udf(w)Uy      =  Udf(w)Uy
C     Tx = 'o', Ux=Udx(-w): Udx(-w)Uy     =  Udx(-w)Uy
C     Tx = 'J', Ux=Jab    : JdabUy        =  JdabUy
C
C     Ty = 'a', Uy=Ua     : UdxIa + UdxUa = -UdxUda
C     Ty = 'f', Uy=Uf(w)  : UdxUf(w)      =  UdxUf(w)
C     Ty = 'm', Uy=Uf(w)  : UdxUf(-w)     = -UdxUdf(w)
C     Ty = 'o', Uy=Uy(w)  : UdxUy(w)      =  UdxUy(w)
C     Ty = 'J', Uy=Jab    : UdxJab        =  UdxJab
C     -------------------------------------------------
C     Variables
      double precision dummy,clear
      double precision ZERO,ONE
      parameter (ZERO=0.0d+00,ONE=1.0d+00)
      character*1 sN,sT,sX,sY
      parameter (sN='N',sT='T')
C     Body of the subprogram
      dummy = -mult
      if(Tx.eq.'a') then
         dummy = -dummy
         sX = sN
      else if(Tx.eq.'f') then
         dummy = -dummy
         sX = sN
      else if(Tx.eq.'m') then
         sX = sT
      else if(Tx.eq.'o') then
         sX = sN
      else if(Tx.eq.'J') then
         sX = sT
      else
         write(*,*) 'Wrong type for Tx [',Tx,'] in CompNIXsum'
         call abrtx("CompNIXsum:1")
      end if
      if(Ty.eq.'a') then
         dummy = -dummy
         sY = sT
      else if(Ty.eq.'f') then
         sY = sN
      else if(Ty.eq.'m') then
         dummy = -dummy
         sY = sT
      else if(Ty.eq.'o') then
         sY = sN
      else if(Ty.eq.'J') then
         sY = sN
      else
         write(*,*) 'Wrong type for TY [',TY,'] in CompNIXsum'
         call abrtx("CompNIXsum:1")
      end if
      if(bClear) clear=ZERO
      if(.not.bClear) clear=ONE
      call DGEMM(sX,sY,nocc,nocc,norb,
     &           dummy,Ux,norb,Uy,norb,clear,Xsum,nocc)
C
      end
C
C*MODULE TDXNI   *DECK CompNIXsum2
      subroutine CompNIXsum2(Xsum,Ux,Tx,Ia,Uy,Ty,norb,nocc,mult,bClear)
      implicit none
C     Arguments
      integer norb,nocc
      double precision mult
      double precision Xsum(nocc*nocc)
      double precision Ux(norb*norb)
      double precision Uy(norb*norb)
      double precision Ia(norb*norb)
      character*1 Tx,Ty
      logical bClear
C     Used common block
      DOUBLE PRECISION XX
      COMMON /FMCOM / XX(1)
C     ------------------------------------------------
C     -mult * op(Ux) * I * op(Uy) = Xxy
C     op depends on the type of Ux and Uy (Tx and Ty).
C
C     Tx = 'f', Ux=Uf(w)  : Udf(-w)Uy    = -Uf(w)Uy
C     Tx = 'm', Ux=Uf(w)  : Udf(w)Uy     =  Udf(w)Uy
C     Tx = 'o', Ux=Udx(-w): Udx(-w)Uy    =  Udx(-w)Uy
C
C     Ty = 'f', Uy=Uf(w)  : UdxUf(w)     =  UdxUf(w)
C     Ty = 'm', Uy=Uf(w)  : UdxUf(-w)    = -UdxUdf(w)
C     Ty = 'o', Uy=Uy(w)  : UdxUy(w)     =  UdxUy(w)
C     ------------------------------------------------
C     Variables
      integer lWork
      double precision dummy,clear
      double precision ZERO,ONE
      parameter (ZERO=0.0d+00,ONE=1.0d+00)
      character*1 sN,sT,sX,sY
      parameter (sN='N',sT='T')
C     Body of the subprogram
      call NEWMEM(lWork,nocc*norb)
      dummy = -mult
      if(Tx.eq.'f') then
         dummy = -dummy
         sX = sN
      else if(Tx.eq.'m') then
         sX = sT
      else if(Tx.eq.'o') then
         sX = sN
      else
         write(*,*) 'Wrong type for Tx [',Tx,'] in CompNIXsum2'
         call abrtx("CompNIXsum2:1")
      end if
      if(Ty.eq.'f') then
         sY = sN
      else if(Ty.eq.'m') then
         dummy = -dummy
         sY = sT
      else if(Ty.eq.'o') then
         sY = sN
      else
         write(*,*) 'Wrong type for TY [',TY,'] in CompNIXsum2'
         call abrtx("CompNIXsum2:2")
      end if
      if(bClear) clear=ZERO
      if(.not.bClear) clear=ONE
      call DGEMM(sX,sN,nocc,norb,norb,
     &           ONE,Ux,norb,Ia,norb,ZERO,XX(lWork),nocc)
      call DGEMM(sN,sY,nocc,nocc,norb,
     &           dummy,XX(lWork),nocc,Uy,norb,clear,Xsum,nocc)
C
      call DELETEMEM(lWork,nocc*norb)
      end
C
C*MODULE TDXNI   *DECK CompNIUfaSum
      subroutine CompNIUfaSum(Usum,Uf,Tf,Ua,nocc,norb,mult,bclear)
      implicit none
C     Arguments
      character*1 Tf
      integer nocc,norb
      double precision mult
      logical bclear
      double precision Usum(norb*norb)
      double precision Uf(norb*norb),Ua(norb*norb)
C     Body of the subprogram
      call CompNIUsum(Usum,Ua,'a',Uf,Tf,norb,nocc,mult,bClear)
      call CompNIUsum(Usum,Uf,Tf,Ua,'a',norb,nocc,mult,.false.)
      end
C
C*MODULE TDXNI   *DECK CompNIXfaSum
      subroutine CompNIXfaSum(Xsum,Uf,Tf,Ua,nocc,norb,mult,bclear)
      implicit none
C     Arguments
      character*1 Tf
      integer nocc,norb
      double precision mult
      logical bclear
      double precision Xsum(nocc*nocc)
      double precision Uf(norb*norb),Ua(norb*norb)
C     Body of the subprogram
C     ------------------------------
C     -[Udf(-w)(Ia+Ua)+(Ia+Uda)Uf(w)]
C     = -Uf(w)Uda+UaUf(w)
C     ------------------------------
      call CompNIXsum(Xsum,Uf,Tf,Ua,'a',norb,nocc,mult,bclear)
      call CompNIXsum(Xsum,Ua,'a',Uf,Tf,norb,nocc,mult,.false.)
      end
C
C*MODULE TDXNI   *DECK CompNIUfgaSum1
      subroutine CompNIUfgaSum1(Usum,Ufg,Udfg,Ua,nocc,norb,mult,bclear)
      implicit none
C     Arguments
      integer nocc,norb
      double precision mult
      logical bclear
      double precision Usum(norb*norb)
      double precision Ufg(norb*norb),Udfg(norb*norb),Ua(norb*norb)
C     Body of the subprogram
      call CompNIUsum(Usum,Ua,'a',Udfg,'o',norb,nocc,mult,bclear)
      call CompNIUsum(Usum,Ufg,'o',Ua,'a',norb,nocc,mult,.false.)
      end
C
C*MODULE TDXNI   *DECK CompNIUfgaSum2
      subroutine CompNIUfgaSum2(Usum,Uf,Tf,Uga,Udga,
     &                          nocc,norb,mult,bclear)
      implicit none
C     Arguments
      character*1 Tf
      integer nocc,norb
      double precision mult
      logical bclear
      double precision Usum(norb*norb)
      double precision Uf(norb*norb),Uga(norb*norb),Udga(norb*norb)
C     Body of the subprogram
      call CompNIUsum(Usum,Uf,Tf,Udga,'o',norb,nocc,mult,bclear)
      call CompNIUsum(Usum,Uga,'o',Uf,Tf,norb,nocc,mult,.false.)
      end
C
C*MODULE TDXNI   *DECK CompNIXfgaSum1
      subroutine CompNIXfgaSum1(Xsum,Ufg,Udfg,Ua,nocc,norb,mult,bclear)
      implicit none
C     Arguments
      integer nocc,norb
      double precision mult
      logical bclear
      double precision Xsum(nocc*nocc)
      double precision Ufg(norb*norb),Udfg(norb*norb),Ua(norb*norb)
C     Body of the subprogram
      call CompNIXsum(Xsum,Ua,'a',Ufg,'o',norb,nocc,mult,bclear)
      call CompNIXsum(Xsum,Udfg,'o',Ua,'a',norb,nocc,mult,.false.)
      end
C
C*MODULE TDXNI   *DECK CompNIXfgaSum2
      subroutine CompNIXfgaSum2(Xsum,Uf,Tf,Ug,Tg,bSame,Ia,
     &                          nocc,norb,mult,bclear)
      implicit none
C     Arguments
      character*1 Tf,Tg
      integer nocc,norb
      double precision mult
      logical bclear,bsame
      double precision Xsum(nocc*nocc)
      double precision Uf(norb*norb),Ug(norb*norb),Ia(norb*norb)
C     Variables
      double precision dummy
C     Body of the subprogram
      dummy=mult
      if(bSame) dummy=dummy+dummy
      call CompNIXsum2(Xsum,Ug,Tg,Ia,Uf,Tf,norb,nocc,dummy,bclear)
      if(.not.bSame) then
         call CompNIXsum2(Xsum,Uf,Tf,Ia,Ug,Tg,norb,nocc,dummy,.false.)
      end if
      end
C
C*MODULE TDXNI   *DECK CompNIXfgaSum3
      subroutine CompNIXfgaSum3(Xsum,Uf,Tf,Uga,Udga,
     &                          nocc,norb,mult,bclear)
      implicit none
C     Arguments
      character*1 Tf
      integer nocc,norb
      double precision mult
      logical bclear
      double precision Xsum(nocc*nocc)
      double precision Uf(norb*norb),Uga(norb*norb),Udga(norb*norb)
C     Body of the subprogram
      call CompNIXsum(Xsum,Uf,Tf,Uga,'o',norb,nocc,mult,bclear)
      call CompNIXsum(Xsum,Udga,'o',Uf,Tf,norb,nocc,mult,.false.)
      end
C
C*MODULE TDXNI   *DECK CompNIUfgSum
      subroutine CompNIUfgSum(Usum,Uf,Tf,Ug,Tg,bsame,
     &                        nocc,norb,mult,bclear)
      implicit none
C     Arguments
      character*1 Tf,Tg
      integer nocc,norb
      double precision mult
      logical bclear,bsame
      double precision Usum(norb*norb)
      double precision Uf(norb*norb),Ug(norb*norb)
C     Variables
      double precision dummy
C     Body of the subprogram
      dummy=mult
      if(bSame) dummy=dummy+dummy
      call CompNIUsum(Usum,Uf,Tf,Ug,Tg,norb,nocc,dummy,bclear)
      if(.not.bSame) then
         call CompNIUsum(Usum,Ug,Tg,Uf,Tf,norb,nocc,dummy,.false.)
      end if
      end
C
C*MODULE TDXNI   *DECK CompNIXfgSum
      subroutine CompNIXfgSum(Xsum,Uf,Tf,Ug,Tg,bsame,
     &                        nocc,norb,mult,bclear)
      implicit none
C     Arguments
      character*1 Tf,Tg
      integer nocc,norb
      double precision mult
      logical bclear,bsame
      double precision Xsum(nocc*nocc)
      double precision Uf(norb*norb),Ug(norb*norb)
C     Variables
      double precision dummy
C     Body of the subprogram
      dummy=mult
      if(bSame) dummy=dummy+dummy
      call CompNIXsum(Xsum,Uf,Tf,Ug,Tg,norb,nocc,dummy,bclear)
      if(.not.bSame) then
         call CompNIXsum(Xsum,Ug,Tg,Uf,Tf,norb,nocc,dummy,.false.)
      end if
      end
C
C*MODULE TDXNI   *DECK CompNIUfghSum
      subroutine CompNIUfghSum(Usum,Uf,Ugh,Udgh,nocc,norb,mult,bclear)
      implicit none
C     Arguments
      integer nocc,norb
      double precision mult
      logical bclear
      double precision Usum(norb*norb)
      double precision Uf(norb*norb),Ugh(norb*norb),Udgh(norb*norb)
C     Body of the subprogram
      call CompNIUsum(Usum,Uf,'f',Udgh,'o',norb,nocc,mult,bclear)
      call CompNIUsum(Usum,Ugh,'o',Uf,'f',norb,nocc,mult,.false.)
      end
C
C*MODULE TDXNI   *DECK CompNIXfghSum
      subroutine CompNIXfghSum(Xsum,Uf,Ugh,Udgh,nocc,norb,mult,bclear)
      implicit none
C     Arguments
      integer nocc,norb
      double precision mult
      logical bclear
      double precision Xsum(nocc*nocc)
      double precision Uf(norb*norb),Ugh(norb*norb),Udgh(norb*norb)
C     Body of the subprogram
      call CompNIXsum(Xsum,Uf,'f',Ugh,'o',norb,nocc,mult,bclear)
      call CompNIXsum(Xsum,Udgh,'o',Uf,'f',norb,nocc,mult,.false.)
      end
C
C*MODULE TDXNI   *DECK CompNIUabSum
      subroutine CompNIUabSum(Usum,Ua,Ub,bsame,nocc,norb,mult,bclear)
      implicit none
C     Arguments
      integer nocc,norb
      double precision mult
      logical bsame,bclear
      double precision Usum(norb*norb)
      double precision Ua(norb*norb),Ub(norb*norb)
C     Variable
      double precision dummy
C     Body of the subprogram
      dummy = mult
      if(bSame) dummy=dummy+dummy
      call CompNIUsum(Usum,Ua,'a',Ub,'a',norb,nocc,dummy,bClear)
      if(.not.bSame) then
         call CompNIUsum(Usum,Ub,'a',Ua,'a',norb,nocc,mult,.false.)
      end if
      end
C
C*MODULE TDXNI   *DECK CompNIUfabSum
      subroutine CompNIUfabSum(Usum,Ufa,Udfa,Ub,Ufb,Udfb,Ua,bSame,
     &                         nocc,norb,mult,bclear)
      implicit none
C     Arguments
      integer nocc,norb
      double precision mult
      logical bclear,bSame
      double precision Usum(norb*norb)
      double precision Ufa(norb*norb),Udfa(norb*norb),Ub(norb*norb)
      double precision Ufb(norb*norb),Udfb(norb*norb),Ua(norb*norb)
C     Variable
      double precision dummy
C     Body of the subprogram
      dummy = mult
      if(bSame) dummy=dummy+dummy
C     -----------------------------------------
C     Uga(w)Udb+UbUdga(-w) + (a,b) permutations
C     -----------------------------------------
      call CompNIUsum(Usum,Ufa,'o',Ub,'a',norb,nocc,dummy,bclear)
      call CompNIUsum(Usum,Ub,'a',Udfa,'o',norb,nocc,dummy,.false.)
      if(.not.bSame) then
         call CompNIUsum(Usum,Ufb,'o',Ua,'a',norb,nocc,dummy,.false.)
         call CompNIUsum(Usum,Ua,'a',Udfb,'o',norb,nocc,dummy,.false.)
      end if
      end
C
C
C*MODULE TDXNI   *DECK CompNIXfabSum1
      subroutine CompNIXfabSum1(Xsum,Ufa,Udfa,Ub,Ufb,Udfb,Ua,
     &                          bSame,nocc,norb,mult,bclear)
      implicit none
C     Arguments
      integer nocc,norb
      double precision mult
      logical bclear,bSame
      double precision Xsum(nocc*nocc)
      double precision Ufa(norb*norb),Udfa(norb*norb),Ub(norb*norb)
      double precision Ufb(norb*norb),Udfb(norb*norb),Ua(norb*norb)
C     Variable
      double precision dummy
C     Body of the subprogram
C     -----------------------------------------------------------
C     -[Udfa(-w)Ib+Udfa(-w)Ub]=[Udfa(-w)Udb] + (a,b) permutations
C     -[IbUfa(w)+UdbUfa(w)]=[UbUfa(w)]       + (a,b) permutations
C     -----------------------------------------------------------
      dummy = mult
      if(bSame) dummy=dummy+dummy
      call CompNIXsum(Xsum,Udfa,'o',Ub,'a',norb,nocc,dummy,bclear)
      call CompNIXsum(Xsum,Ub,'a',Ufa,'o',norb,nocc,dummy,.false.)
      if(.not.bSame) then
         call CompNIXsum(Xsum,Udfb,'o',Ua,'a',norb,nocc,dummy,.false.)
         call CompNIXsum(Xsum,Ua,'a',Ufb,'o',norb,nocc,dummy,.false.)
      end if
      end
C
C*MODULE TDXNI   *DECK CompNIXfabSum2
      subroutine CompNIXfabSum2(Xsum,Jab,Uf,Tf,nocc,norb,mult,bclear)
      implicit none
C     Arguments
      character*1 Tf
      integer nocc,norb
      double precision mult
      logical bclear
      double precision Xsum(nocc*nocc)
      double precision Jab(norb*norb)
      double precision Uf(norb*norb)
C     Body of the subprogram
C     --------------------------------------
C     -[Udf(-w)(Iab+IaUb+IbUa)]=[Uf(w)Jab]
C     -[(Iab+UdbIa+UdaIb)Uf(w)]=-[JdabUf(w)]
C     --------------------------------------
      call CompNIXsum(Xsum,Uf,Tf,Jab,'J',norb,nocc,mult,bclear)
      call CompNIXsum(Xsum,Jab,'J',Uf,Tf,norb,nocc,mult,.false.)
      end
C
C*MODULE TDXNI   *DECK CompNIUfghaSum1
      subroutine CompNIUfghaSum1(Usum,Ua,Ufgh,Udfgh,
     &                           nocc,norb,mult,bclear)
      implicit none
C     Arguments
      integer nocc,norb
      double precision mult
      logical bclear
      double precision Usum(norb*norb)
      double precision Ua(norb*norb),Ufgh(norb*norb),Udfgh(norb*norb)
C     Body of the subprogram
      call CompNIUsum(Usum,Ua,'a',Udfgh,'o',norb,nocc,mult,bclear)
      call CompNIUsum(Usum,Ufgh,'o',Ua,'a',norb,nocc,mult,.false.)
      end
C
C*MODULE TDXNI   *DECK CompNIUfghaSum2
      subroutine CompNIUfghaSum2(Usum,Ufa,Udfa,Ugh,Udgh,
     &                           nocc,norb,mult,bclear)
      implicit none
C     Arguments
      integer nocc,norb
      double precision mult
      logical bclear
      double precision Usum(norb*norb)
      double precision Ufa(norb*norb),Udfa(norb*norb)
      double precision Ugh(norb*norb),Udgh(norb*norb)
C     Body of the subprogram
      call CompNIUsum(Usum,Ufa,'o',Udgh,'o',norb,nocc,mult,bclear)
      call CompNIUsum(Usum,Ugh,'o',Udfa,'o',norb,nocc,mult,.false.)
      end
C
C*MODULE TDXNI   *DECK CompNIXfghaSum1
      subroutine CompNIXfghaSum1(Xsum,Ua,Ufgh,Udfgh,
     &                           nocc,norb,mult,bclear)
      implicit none
C     Arguments
      integer nocc,norb
      double precision mult
      logical bclear
      double precision Xsum(nocc*nocc)
      double precision Ua(norb*norb),Ufgh(norb*norb),Udfgh(norb*norb)
C     Body of the subprogram
      call CompNIXsum(Xsum,Ua,'a',Ufgh,'o',norb,nocc,mult,bclear)
      call CompNIXsum(Xsum,Udfgh,'o',Ua,'a',norb,nocc,mult,.false.)
      end
C
C*MODULE TDXNI   *DECK CompNIXfghaSum2
      subroutine CompNIXfghaSum2(Xsum,Ia,Uf,Ufa,Udfa,Ugh,Udgh,
     &                           nocc,norb,mult,bclear)
      implicit none
C     Arguments
      integer nocc,norb
      double precision mult
      logical bclear
      double precision Xsum(nocc*nocc)
      double precision Ia(norb*norb),Uf(norb*norb)
      double precision Ufa(norb*norb),Udfa(norb*norb)
      double precision Ugh(norb*norb),Udgh(norb*norb)
C     Body of the subprogram
      call CompNIXsum(Xsum,Udfa,'o',Ugh,'o',norb,nocc,mult,bclear)
      call CompNIXsum2(Xsum,Uf,'f',Ia,Ugh,'o',norb,nocc,mult,.false.)
      call CompNIXsum(Xsum,Udgh,'o',Ufa,'o',norb,nocc,mult,.false.)
      call CompNIXsum2(Xsum,Udgh,'o',Ia,Uf,'f',norb,nocc,mult,.false.)
      end
C
C*MODULE TDXNI   *DECK CompNIUfgabSum
      subroutine CompNIUfgabSum(Usum,Ufa,Udfa,Ugb,Udgb,bSame,
     &                          nocc,norb,mult,bclear)
      implicit none
C     Arguments
      integer nocc,norb
      double precision mult
      logical bclear,bSame
      double precision Usum(norb*norb)
      double precision Ufa(norb*norb),Udfa(norb*norb)
      double precision Ugb(norb*norb),Udgb(norb*norb)
C     Variable
      double precision dummy
C     Body of the subprogram
      dummy = mult
      if(bSame) dummy=dummy+dummy
C     ----------------------------
C     UfaUdgb + (fa,gb) permutations
C     ----------------------------
      call CompNIUsum(Usum,Ufa,'o',Udgb,'o',norb,nocc,dummy,bclear)
      if(.not.bSame) then
         call CompNIUsum(Usum,Ugb,'o',Udfa,'o',norb,nocc,dummy,.false.)
      end if
      end
C*MODULE TDXNI   *DECK CompNIXfgabSum1
      subroutine CompNIXfgabSum1(Xsum,Ia,Uf,Ufa,Udfa,Ib,Ug,Ugb,Udgb,
     &                           bSame,nocc,norb,mult,bclear)
      implicit none
C     Arguments
      integer nocc,norb
      double precision mult
      logical bclear,bSame
      double precision Xsum(nocc*nocc)
      double precision Uf(norb*norb),Ia(norb*norb)
      double precision Ufa(norb*norb),Udfa(norb*norb)
      double precision Ug(norb*norb),Ib(norb*norb)
      double precision Ugb(norb*norb),Udgb(norb*norb)
C     Variable
      double precision dummy
C     Body of the subprogram
C     ------------------------------------------------------
C     -Udf(-w)IaUgb(w) - Udgb(-w)IaUf(w)
C     -Udfa(-w)[IbUg(w)+Ugb(w)] - [Udgb(-w)+Udg(-w)Ib]Ufa(w)
C     + (fa,gb) permutations
C     ------------------------------------------------------
      dummy = mult
      if(bSame) dummy=dummy+dummy
      call CompNIXsum(Xsum,Udfa,'o',Ugb,'o',norb,nocc,dummy,bclear)
      call CompNIXsum2(Xsum,Uf,'f',Ia,Ugb,'o',norb,nocc,dummy,.false.)
      call CompNIXsum2(Xsum,Udfa,'o',Ib,Ug,'f',norb,nocc,dummy,.false.)
      if(.not.bSame) then
         call CompNIXsum(Xsum,Udgb,'o',Ufa,'o',norb,nocc,dummy,.false.)
         call CompNIXsum2(Xsum,Ug,'f',Ib,Ufa,'o',norb,nocc,dummy,
     &                    .false.)
         call CompNIXsum2(Xsum,Udgb,'o',Ia,Uf,'f',norb,nocc,dummy,
     &                    .false.)
      end if
      end
C
C*MODULE TDXNI   *DECK CompNIXfgabSum2
      subroutine CompNIXfgabSum2(Xsum,Jab,Ufg,Udfg,
     &                           nocc,norb,mult,bclear)
      implicit none
C     Arguments
      integer nocc,norb
      double precision mult
      logical bclear
      double precision Xsum(nocc*nocc)
      double precision Jab(norb*norb)
      double precision Ufg(norb*norb),Udfg(norb*norb)
C     Body of the subprogram
C     -----------------------------------------
C     -[Udfg(-w)(Iab+IaUb+IbUa)]=-[Udfg(-w)Jab]
C     -[(Iab+UdbIa+UdaIb)Ufg(w)]=-[JdabUfg(w)]
C     -----------------------------------------
      call CompNIXsum(Xsum,Udfg,'o',Jab,'J',norb,nocc,mult,bclear)
      call CompNIXsum(Xsum,Jab,'J',Ufg,'o',norb,nocc,mult,.false.)
      end
C
C*MODULE TDXNI   *DECK CompNIXfgabSum3
      subroutine CompNIXfgabSum3(Xsum,Uf,Tf,Ug,Tg,bSame,Iab,
     &                           nocc,norb,mult,bclear)
      implicit none
C     Arguments
      character*1 Tf,Tg
      integer nocc,norb
      double precision mult
      logical bclear,bsame
      double precision Xsum(nocc*nocc)
      double precision Uf(norb*norb),Ug(norb*norb),Iab(norb*norb)
C     Variables
      double precision dummy
C     Body of the subprogram
      dummy=mult
      if(bSame) dummy=dummy+dummy
      call CompNIXsum2(Xsum,Ug,Tg,Iab,Uf,Tf,norb,nocc,dummy,bclear)
      if(.not.bSame) then
         call CompNIXsum2(Xsum,Uf,Tf,Iab,Ug,Tg,norb,nocc,dummy,.false.)
      end if
      end
C*MODULE TDXNI   *DECK MakeDipole
      subroutine MakeDipole
      implicit none
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 dipole(3)
      integer lC0,lHf,lD0
C     Called subprograms
      logical GetHf
C     Body of the subprogram
      call VCLR(dipole,1,3)
C
      call NEWMEM(lD0,norb*norb)
      call NEWMEM(lC0,norb*norb)
      call GetLCAO(XX(lC0),norb)
      call MakeD0(XX(lD0),XX(lC0),norb,nocc)
      call DELETEMEM(lC0,norb*norb)
C
      call NEWMEM(lHf,norb*norb*nfld)
      if(.not.GetHf(XX(lHF),.true.)) then
         call MakeHf()
         if(.not.GetHf(XX(lHF),.true.)) then
            call abrtx("MakeDipole: not GetHF")
         end if
      end if
C      call PrintMM(XX(lHf),norb,nfld,'OLIX:Hf')
C
      call CompDipole(dipole,xx(lD0),xx(lHf),norb,ibfld,nfld)
C
      call DELETEMEM(lHf,norb*norb*nfld)
      call DELETEMEM(lD0,norb*norb)
C
      call SetDipole(dipole)
      call PrintMuX(iw,dipole,'Dipole moment [in au]')
      CALL TIMIT(1)
      end
C
C*MODULE TDXNI   *DECK CompDipole
      subroutine CompDipole(mu,D0,Hf,norb,ibfld,nfld)
      implicit none
C     Arguments
      integer norb,ibfld,nfld
      double precision mu(3)
      double precision D0(norb*norb)
      double precision Hf(norb*norb,nfld)
C     Used common block
      integer IR, IW, IP, IS,IPK,IDAF,NAV,IODA
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      integer MXATM
      parameter (MXATM=2000)
      integer NAT,ICH,MUL,NUM,NQMT,NE,NA,NB,IAN
      double precision ZAN,C
      COMMON /INFOA / NAT,ICH,MUL,NUM,NQMT,NE,NA,NB,
     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
C     Variables
      integer ifld,iifld,iat
      double precision dummy
      double precision cm(3)
C     Called subprograms
      double precision dtrace
C     Body of the subprogram
      write(iw,3000)
      call CALCOM(cm(1),cm(2),cm(3))
      iifld = ibfld
      do ifld = 1,nfld
         dummy = -dtrace('S','S',D0,norb,Hf(1,ifld),norb,norb,norb)
C
         do iat = 1,nat
            dummy = dummy + ZAN(iat)*(C(iifld,iat)-cm(iifld))
         enddo
         call SetMuElem(mu,iifld,dummy)
         iifld = iifld+1
      end do
 3000 FORMAT('NonIterative procedure to compute Mu')
       end
C
C*MODULE TDXNI   *DECK MakeBetaNI
      subroutine MakeBetaNI(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)
C     Variables and constants
      double precision mu(3)
      double precision beta(27)
      double precision w(3)
      double precision dummy(1)
C     Called subprograms
      logical GotDip,GetGUDEf
C     Body of the subprogram
      call VCLR(beta,1,27)
C
      w(1)=-w1-w2
      w(2)=w1
      w(3)=w2
C
      dummy(1) = 0.0d+00
      call GetDipole(dummy,.false.,GotDip)
      if(.not.GotDip) call MakeDipole()
      if(.not.GetGUDEf(w(1),dummy,.false.,dummy,.false.,dummy,.false.,
     &     dummy,.false.,norb,nfld)) call MakeAlpha(w(1))
      if(.not.GetGUDEf(w(2),dummy,.false.,dummy,.false.,dummy,.false.,
     &     dummy,.false.,norb,nfld)) call MakeAlpha(w(2))
      if(.not.GetGUDEf(w(3),dummy,.false.,dummy,.false.,dummy,.false.,
     &     dummy,.false.,norb,nfld)) call MakeAlpha(w(3))
C
      call CompBetaNI(beta,w,norb,nocc,ibfld,nfld)
C
      call SetBeta(w1,w2,beta)
      call GetDipole(mu,.true.,GotDip)
      if(.not.GotDip) then
         call abrtx("MakeBetaNI: not GetDop")
      end if
      call PrintBetaX(iw,beta,mu,w,'Beta tensor [in au]')
      CALL TIMIT(1)
      end
C
C*MODULE TDXNI   *DECK CompBetaNI
      subroutine CompBetaNI(beta,w,norb,nocc,ibfld,nfld)
      implicit none
C     Arguments
      integer norb,nocc,ibfld,nfld
      double precision w(3)
      double precision beta(27)
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 lGf,lEf,lUg,lUh
C     Body of the subprogram
      write(iw,3000) w(1),w(2),w(3)
C     -----
C     Step1
C     -----
      call NEWMEM(lUh,norb*norb*nfld)
      call NEWMEM(lUg,norb*norb*nfld)
      call NEWMEM(lGf,norb*norb*nfld)
      call NEWMEM(lEf,norb*norb*nfld)
      call MakeBetaNI1(beta,w,XX(lGf),XX(lEf),XX(lUg),XX(lUh),
     &                 norb,nocc,ibfld,nfld)
      call DELETEMEM(lEf,norb*norb*nfld)
      call DELETEMEM(lGf,norb*norb*nfld)
      call DELETEMEM(lUg,norb*norb*nfld)
      call DELETEMEM(lUh,norb*norb*nfld)
C     Format
 3000 FORMAT('Non-Iterative procedure to compute Beta(',
     &       F13.6,';',F13.6,',',F13.6,')')
      end
C
C*MODULE TDXNI   *DECK MakeBetaNIPerms
      subroutine MakeBetaNIPerms(w,ibfld,nfld,nterms,iflds,mults,
     &                           ws,bgws)
      implicit none
C     Arguments
      double precision w(3)
      integer ibfld,nfld
      integer nterms(3)
      integer iflds(4,*)
      double precision mults(*)
      double precision ws(3,3)
      logical bgws(3,3)
C     Used common block
      integer NELEMB,IELEMB,INDEXB
      common/NLO_BETA/NELEMB,IELEMB(3,27),INDEXB(3,3,3)
C     Variables and constants
      logical bw13
      integer index,iterm,ifld,jfld,kfld,iifld,jjfld,kkfld
      double precision mult,ONE
      parameter (ONE=1.0d+00)
      integer i,iperm,nperm,permlist(3)
      parameter (nperm=3)
C     Body of the subprogram
C     ------------------
C     Smart permutations
C     ------------------
C--   bw12=w(1).eq.w(2)
      bw13=w(1).eq.w(3)
C--   bw23=w(2).eq.w(3)
      if(bw13) then
         permlist(1) = 1
         permlist(2) = 3
         permlist(3) = 2
      else
         permlist(1) = 1
         permlist(2) = 2
         permlist(3) = 3
      end if
      iterm = 0
      do i=1,nperm
         iperm = permlist(i)
         if(iperm.eq.1) then
            ws(1,i) = w(1)
            ws(2,i) = w(2)
            ws(3,i) = w(3)
         else if(iperm.eq.2) then
            ws(1,i) = w(2)
            ws(2,i) = w(1)
            ws(3,i) = w(3)
         else if(iperm.eq.3) then
            ws(1,i) = w(3)
            ws(2,i) = w(1)
            ws(3,i) = w(2)
         end if
         do 1011 index=1,NELEMB
            iifld = IELEMB(1,index)-ibfld+1
            jjfld = IELEMB(2,index)-ibfld+1
            kkfld = IELEMB(3,index)-ibfld+1
            if((iifld.lt.1).or.(iifld.gt.nfld)) go to 1011
            if((jjfld.lt.1).or.(jjfld.gt.nfld)) go to 1011
            if((kkfld.lt.1).or.(kkfld.gt.nfld)) go to 1011
            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
C               b12=(iifld.eq.jjfld).and.bw12
C               b13=(iifld.eq.kkfld).and.bw13
C               b23=(jjfld.eq.kkfld).and.bw23
            if(iperm.eq.1) then
C               if(b12.and.b13) then
C                  mult = THREE
C               else if(b12.or.b13) then
C                  mult = TWO
C               else
                  mult = ONE
C               end if
            else if(iperm.eq.2) then
C               if((b12.and.b13).or.b12) go to 1011
C               if(b23) then
C                  mult = TWO
C               else
                  mult = ONE
C               end if
            else if(iperm.eq.3) then
C               if((b12.and.b13).or.b13.or.b23) go to 1011
               mult = ONE
            end if
            iterm=iterm+1
            iflds(1,iterm) = ifld
            iflds(2,iterm) = jfld
            iflds(3,iterm) = kfld
            iflds(4,iterm) = index
            mults(iterm)   = mult
 1011    continue
         nterms(i) = iterm
      end do
      do i = 1,3
         bgws(i,1) = .true.
C         bgws(i,2) = ws(i,2).ne.ws(i,1)
C         bgws(i,3) = ws(i,3).ne.ws(i,2)
         bgws(i,2) = .true.
         bgws(i,3) = .true.
      end do
      end
C
C*MODULE TDXNI   *DECK MakeBetaNI1
      subroutine MakeBetaNI1(beta,w,Gf,Ef,Ug,Uh,
     &                       norb,nocc,ibfld,nfld)
      implicit none
C     Arguments
      double precision beta(27)
      double precision w(3)
      integer norb,nocc,ibfld,nfld
      double precision Gf(norb*norb,nfld)
      double precision Ef(norb*norb,nfld)
      double precision Ug(norb*norb,nfld)
      double precision Uh(norb*norb,nfld)
C     Used common block
      DOUBLE PRECISION XX
      COMMON /FMCOM / XX(1)
C     Variables
      integer iterm,ibterm,ieterm,ifld,jfld,kfld,index
      double precision mult,dummy
      integer iperm,nperm
      parameter (nperm=3)
      integer nterms(nperm)
      integer iflds(4,nperm*27)
      double precision mults(nperm*27)
      double precision ws(3,nperm)
      logical bgws(3,nperm)
      integer lWork1,lWork2
      double precision ONE,TWO
      parameter (ONE=1.0d+00,TWO=2.0d+00)
      logical bSame,bCompjk,bCompijk
      integer iold,jold,kold
C     Called subprograms
      logical GetGUDEf
      double precision DTrace
C     Body of the subprogram
      call MakeBetaNIPerms(w,ibfld,nfld,nterms,iflds,mults,
     &                     ws,bgws)
      call NEWMEM(lWork1,norb*norb)
      call NEWMEM(lWork2,nocc*nocc)
C
C      write(*,*) "BetaNI:iflds",iflds
C      write(*,*) "BetaNI:mutls",mults
C      write(*,*) "BetaNI:ws",ws
C      write(*,*) "BetaNI:bgws",bgws
C
      ieterm = 0
      dummy = 0.0d+00
      do iperm = 1,nperm
         if(bgws(1,iperm)) then
            iold = -1
            if(.not.GetGUDEf(ws(1,iperm),
     &                       Gf,.true.,dummy,.false.,
     &                       dummy,.false.,Ef,.true.,
     &                       norb,nfld)) then
               call abrtx("MakeBetaNI1: not GetGUDEf:1")
            end if
         end if
         if(bgws(2,iperm)) then
            jold = -1
            if(.not.GetGUDEf(ws(2,iperm),
     &                       dummy,.false.,Ug,.true.,
     &                       dummy,.false.,dummy,.false.,
     &                       norb,nfld)) then
               call abrtx("MakeBetaNI1: not GetGUDEf:2")
            end if
         end if
         if(bgws(3,iperm)) then
            kold = -1
            if(.not.GetGUDEf(ws(3,iperm),
     &                       dummy,.false.,Uh,.true.,
     &                       dummy,.false.,dummy,.false.,
     &                       norb,nfld)) then
               call abrtx("MakeBetaNI1: not GetGUDEf:3")
            end if
         end if
C
         ibterm=ieterm+1
         ieterm=nterms(iperm)
         do iterm = ibterm,ieterm
            ifld = iflds(1,iterm)
            jfld = iflds(2,iterm)
            kfld = iflds(3,iterm)
            index = iflds(4,iterm)
            mult = mults(iterm)
            bSame = (jfld.eq.kfld).and.(ws(2,iperm).eq.ws(3,iperm))
C
            bCompjk = (jfld.ne.jold).or.(kfld.ne.kold)
            bCompijk = (ifld.ne.iold).or.bCompjk
            iold = ifld
            jold = jfld
            kold = kfld
C
            if(bCompjk) then
               call CompNIUfgSum(XX(lWork1),
     &                           Ug(1,jfld),'f',Uh(1,kfld),'f',bsame,
     &                           nocc,norb,ONE,.true.)
               call CompNIXfgSum(XX(lWork2),
     &                           Ug(1,jfld),'f',Uh(1,kfld),'f',bsame,
     &                           nocc,norb,ONE,.true.)
            end if
            if(bCompijk) then
               dummy = DTrace('N','N',
     &                        XX(lWork1),norb,Gf(1,ifld),norb,norb,norb)
     &                +DTrace('N','N',
     &                        XX(lWork2),nocc,Ef(1,ifld),norb,nocc,nocc)
            end if
            mult = mult*TWO
            beta(index)=beta(index)-mult*dummy
         end do
      end do
C
      call DELETEMEM(lWork2,nocc*nocc)
      call DELETEMEM(lWork1,norb*norb)
C
      end
C*MODULE TDXNI   *DECK MakeGammaNI
      subroutine MakeGammaNI(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)
C     Variables and constants
      double precision gamma(21)
      double precision w(4)
      double precision dummy(1)
C     Called subprograms
      logical 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
C
      dummy(1) = 0.0d+00
      if(.not.GetGUDEf(w(1),dummy,.false.,dummy,.false.,dummy,.false.,
     &     dummy,.false.,norb,nfld)) call MakeAlpha(w(1))
      if(.not.GetGUDEf(w(2),dummy,.false.,dummy,.false.,dummy,.false.,
     &     dummy,.false.,norb,nfld)) call MakeAlpha(w(2))
      if(.not.GetGUDEf(w(3),dummy,.false.,dummy,.false.,dummy,.false.,
     &     dummy,.false.,norb,nfld)) call MakeAlpha(w(3))
      if(.not.GetGUDEf(w(4),dummy,.false.,dummy,.false.,dummy,.false.,
     &     dummy,.false.,norb,nfld)) call MakeAlpha(w(4))
      if(.not.GetGUXDEfg(w(2),w(3),dummy,.false.,dummy,.false.,
     &     dummy,.false.,dummy,.false.,dummy,.false.,norb,nfld))
     &     call MakeBeta(w(2),w(3))
      if(.not.GetGUXDEfg(w(2),w(4),dummy,.false.,dummy,.false.,
     &     dummy,.false.,dummy,.false.,dummy,.false.,norb,nfld))
     &     call MakeBeta(w(2),w(4))
      if(.not.GetGUXDEfg(w(3),w(4),dummy,.false.,dummy,.false.,
     &     dummy,.false.,dummy,.false.,dummy,.false.,norb,nfld))
     &     call MakeBeta(w(3),w(4))
C
      call CompGammaNI(gamma,w,norb,nocc,ibfld,nfld)
C
      call SetGamma(w1,w2,w3,gamma)
      call PrintGammaX(iw,gamma,w,'Gamma tensor [in au]')
      CALL TIMIT(1)
      end
C
C*MODULE TDXNI   *DECK CompGamma
      subroutine CompGammaNI(gamma,w,norb,nocc,ibfld,nfld)
      implicit none
C     Arguments
      integer norb,nocc,ibfld,nfld
      double precision w(4)
      double precision gamma(21)
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 lGf,LUf,lEf,lGg,lUg,lEg,lGhi,lUhi,lXhi,lEhi
C     Body of the subprogram
      write(iw,3000) w(1),w(2),w(3),w(4)
C     -----
C     Step1
C     -----
      call NEWMEM(lGf,norb*norb*nfld)
      call NEWMEM(lUf,norb*norb*nfld)
      call NEWMEM(lEf,norb*norb*nfld)
      call NEWMEM(lGg,norb*norb*nfld)
      call NEWMEM(lUg,norb*norb*nfld)
      call NEWMEM(lEg,norb*norb*nfld)
      call NEWMEM(lGhi,norb*norb*nfld*nfld)
      call NEWMEM(lUhi,norb*norb*nfld*nfld)
      call NEWMEM(lXhi,norb*norb*nfld*nfld)
      call NEWMEM(lEhi,norb*norb*nfld*nfld)
      call MakeGammaNI1(gamma,w,
     &                  XX(lUf),XX(lGf),XX(lEf),
     &                  XX(lUg),XX(lGg),XX(lEg),
     &                  XX(lUhi),XX(lXhi),XX(lGhi),XX(lEhi),
     &                  norb,nocc,ibfld,nfld)
      call DELETEMEM(lEhi,norb*norb*nfld*nfld)
      call DELETEMEM(lXhi,norb*norb*nfld*nfld)
      call DELETEMEM(lUhi,norb*norb*nfld*nfld)
      call DELETEMEM(lGhi,norb*norb*nfld*nfld)
      call DELETEMEM(lEg,norb*norb*nfld)
      call DELETEMEM(lUg,norb*norb*nfld)
      call DELETEMEM(lGg,norb*norb*nfld)
      call DELETEMEM(lEf,norb*norb*nfld)
      call DELETEMEM(lUf,norb*norb*nfld)
      call DELETEMEM(lGf,norb*norb*nfld)
C     Format
 3000 FORMAT('Non-Iterative procedure to compute Gamma(',
     &       F13.6,';',F13.6,',',F13.6,',',F13.6,')')
      end
C
C*MODULE TDXNI   *DECK MakeGammaaNIPerms
      subroutine MakeGammaNIPerms(w,ibfld,nfld,nterms,iflds,mults,
     &                            ws,bgws)
      implicit none
C     Arguments
      double precision w(4)
      integer ibfld,nfld
      integer nterms(3)
      integer iflds(5,*)
      double precision mults(*)
      double precision ws(4,3)
      logical bgws(4,3)
C     Used common block
      integer NELEMG,IELEMG,INDEXG
      common/NLO_GAMMA/NELEMG,IELEMG(4,21),INDEXG(3,3,3,3)
C     Variables and constants
      integer iterm,ifld,jfld,kfld,lfld,index,iifld,jjfld,kkfld,llfld
      double precision mult,ONE
      parameter (ONE=1.0d+00)
      integer i,iperm,nperm,permlist(3)
      parameter (nperm=3)
C     Body of the subprogram
C     ------------------
C     Smart permutations
C     ------------------
C--      bw23=w(2).eq.w(3)
C--      bw24=w(2).eq.w(4)
C--      bw34=w(3).eq.w(4)
      permlist(1) = 1
      permlist(2) = 2
      permlist(3) = 3
      iterm = 0
      do i=1,nperm
         iperm = permlist(i)
         if(iperm.eq.1) then
            ws(1,i) = w(1)
            ws(2,i) = w(2)
            ws(3,i) = w(3)
            ws(4,i) = w(4)
         else if(iperm.eq.2) then
            ws(1,i) = w(1)
            ws(2,i) = w(3)
            ws(3,i) = w(2)
            ws(4,i) = w(4)
         else if(iperm.eq.3) then
            ws(1,i) = w(1)
            ws(2,i) = w(4)
            ws(3,i) = w(2)
            ws(4,i) = w(3)
         end if
         do 1011 index=1,NELEMG
            iifld = IELEMG(1,index)-ibfld+1
            jjfld = IELEMG(2,index)-ibfld+1
            kkfld = IELEMG(3,index)-ibfld+1
            llfld = IELEMG(4,index)-ibfld+1
            if((iifld.lt.1).or.(iifld.gt.nfld)) go to 1011
            if((jjfld.lt.1).or.(jjfld.gt.nfld)) go to 1011
            if((kkfld.lt.1).or.(kkfld.gt.nfld)) go to 1011
            if((llfld.lt.1).or.(llfld.gt.nfld)) go to 1011
            if(iperm.eq.1) then
               ifld = iifld
               jfld = jjfld
               kfld = kkfld
               lfld = llfld
            else if(iperm.eq.2) then
               ifld = iifld
               jfld = kkfld
               kfld = jjfld
               lfld = llfld
            else if(iperm.eq.3) then
               ifld = iifld
               jfld = llfld
               kfld = jjfld
               lfld = kkfld
            end if
Cc$$$ logical bw23,bw24,bw34,b23,b24,b34
Cc$$$            b23=(jjfld.eq.kkfld).and.bw23
Cc$$$            b24=(jjfld.eq.llfld).and.bw24
Cc$$$            b34=(kkfld.eq.llfld).and.bw34
Cc$$$            if(iperm.eq.1) then
Cc$$$               if(b23.and.b24) then
Cc$$$                  mult = THREE
Cc$$$               else if(b23.or.b24) then
Cc$$$                  mult = TWO
Cc$$$               else
Cc$$$                  mult = ONE
Cc$$$               end if
Cc$$$            else if(iperm.eq.2) then
Cc$$$               if((b23.and.b24).or.b23) go to 1011
Cc$$$               if(b24) then
Cc$$$                  mult = TWO
Cc$$$               else
Cc$$$                  mult = ONE
Cc$$$               end if
Cc$$$            else if(iperm.eq.3) then
Cc$$$               if((b23.and.b24).or.b23.or.b24) go to 1011
Cc$$$               mult = ONE
Cc$$$            end if
            mult = ONE
            iterm=iterm+1
            iflds(1,iterm) = ifld
            iflds(2,iterm) = jfld
            iflds(3,iterm) = kfld
            iflds(4,iterm) = lfld
            iflds(5,iterm) = index
            mults(iterm)   = mult
 1011    continue
         nterms(i) = iterm
      end do
      do i = 1,4
         bgws(i,1) = .true.
         bgws(i,2) = ws(i,2).ne.ws(i,1)
         bgws(i,3) = ws(i,3).ne.ws(i,2)
      end do
      end
C
C*MODULE TDXNI   *DECK MakeGammaNI1
      subroutine MakeGammaNI1(gamma,w,Uf,Gf,Ef,Ug,Gg,Eg,
     &                        Uhi,Udhi,Ghi,Ehi,
     &                        norb,nocc,ibfld,nfld)
      implicit none
C     Arguments
      integer norb,nocc,ibfld,nfld
      double precision gamma(21)
      double precision w(4)
      double precision Uf(norb*norb,nfld)
      double precision Gf(norb*norb,nfld)
      double precision Ef(norb*norb,nfld)
      double precision Ug(norb*norb,nfld)
      double precision Gg(norb*norb,nfld)
      double precision Eg(norb*norb,nfld)
      double precision Uhi(norb*norb,nfld,nfld)
      double precision Udhi(norb*norb,nfld,nfld)
      double precision Ghi(norb*norb,nfld,nfld)
      double precision Ehi(norb*norb,nfld,nfld)
C     Used common block
      DOUBLE PRECISION XX
      COMMON /FMCOM / XX(1)
C     Variables and contants
      integer iterm,ibterm,ieterm
      integer ifld,jfld,kfld,lfld,index
      double precision mult,dummy
      integer iperm,nperm
      parameter (nperm=3)
      integer nterms(nperm)
      integer iflds(5,nperm*21)
      double precision mults(nperm*21)
      double precision ws(4,nperm)
      logical bgws(4,nperm)
      integer lWork
      double precision ONE,TWO,MINUSONE
      parameter (ONE=1.0d+00,TWO=2.0d+00,MINUSONE=-1.0d+00)
      logical bSameij
C     Called subprograms
      logical GetGUDEf,GetGUXDEfg
      double precision DTrace
C     Body of the subprogram
      call MakeGammaNIPerms(w,ibfld,nfld,nterms,iflds,mults,
     &                      ws,bgws)
C
      call NEWMEM(lWork,norb*norb)
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("MakeGammaNI1: not GetGUDEf:1")
            end if
         end if
         if(bgws(2,iperm)) then
            if(.not.GetGUDEf(ws(2,iperm),
     &                       Gg,.true.,Ug,.true.,
     &                       dummy,.false.,Eg,.true.,
     &                       norb,nfld)) then
               call abrtx("MakeGammaNI1: not GetGUDEf:2")
            end if
         end if
         if(bgws(3,iperm).or.bgws(4,iperm)) then
            if(.not.GetGUXDEfg(ws(3,iperm),ws(4,iperm),
     &                         Ghi,.true.,Uhi,.true.,Udhi,.true.,
     &                         dummy,.false.,Ehi,.true.,
     &                         norb,nfld)) then
               call abrtx("MakeGammaNI1: not GetGUXDEfg")
            end if
C           form Udagger(-w2,-w3)=X(w2,w3)-U(w2,w3)
            call DAXPY(norb*norb*nfld*nfld,MINUSONE,Uhi,1,Udhi,1)
         end if
C
         ibterm=ieterm+1
         ieterm=nterms(iperm)
         do iterm = ibterm,ieterm
            ifld = iflds(1,iterm)
            jfld = iflds(2,iterm)
            kfld = iflds(3,iterm)
            lfld = iflds(4,iterm)
            index = iflds(5,iterm)
            mult = mults(iterm)
C
            bSameij = (ifld.eq.jfld).and.(ws(1,iperm).eq.ws(2,iperm))
C            bSameij = .false.
C
C
            call CompNIUfghSum(XX(lWork),Ug(1,jfld),
     &                         Uhi(1,kfld,lfld),Udhi(1,kfld,lfld),
     &                         nocc,norb,ONE,.true.)
            dummy = DTrace('N','N',XX(lWork),norb,Gf(1,ifld),norb,
     &                     norb,norb)
C
            call CompNIXfghSum(XX(lWork),Ug(1,jfld),
     &                         Uhi(1,kfld,lfld),Udhi(1,kfld,lfld),
     &                         nocc,norb,ONE,.true.)
            dummy = dummy
     &             +DTrace('N','N',XX(lWork),nocc,Ef(1,ifld),norb,
     &                     nocc,nocc)
C
            call CompNIUfghSum(XX(lWork),Uf(1,ifld),
     &                         Uhi(1,kfld,lfld),Udhi(1,kfld,lfld),
     &                         nocc,norb,ONE,.true.)
            dummy = dummy
     &             +DTrace('N','N',XX(lWork),norb,Gg(1,jfld),norb,
     &                     norb,norb)
C
            call CompNIXfghSum(XX(lWork),Uf(1,ifld),
     &                         Uhi(1,kfld,lfld),Udhi(1,kfld,lfld),
     &                         nocc,norb,ONE,.true.)
            dummy = dummy
     &             +DTrace('N','N',XX(lWork),nocc,Eg(1,jfld),norb,
     &                     nocc,nocc)
C
            call CompNIUfgSum(XX(lWork),Uf(1,ifld),'f',Ug(1,jfld),'f',
     &                        bSameij,nocc,norb,ONE,.true.)
            dummy = dummy
     &            +DTrace('N','N',XX(lWork),norb,Ghi(1,kfld,lfld),norb,
     &                    norb,norb)
C
            call CompNIXfgSum(XX(lWork),Uf(1,ifld),'f',Ug(1,jfld),'f',
     &                        bSameij,nocc,norb,ONE,.true.)
            dummy = dummy
     &            +DTrace('N','N',XX(lWork),nocc,Ehi(1,kfld,lfld),norb,
     &                    nocc,nocc)
C
            mult=mult*TWO
            gamma(index) = gamma(index)-mult*dummy
         end do
      end do
C
      call DELETEMEM(lWork,norb*norb)
C
      end
C
C*MODULE TDXNI   *DECK MakedMudXNI
      subroutine MakedMudXNI()
      implicit none
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 ldMdX
      double precision dummy(1)
C     Called subprograms
      logical GetGUIDEa,GetHf,GetHfa
C     Body of the subprogram
      dummy(1) = 0.0d+00
      if(.not.GetGUIDEa(dummy,.false.,dummy,.false.,dummy,.false.,
     &                  dummy,.false.,dummy,.false.,norb,ngeom)) then
         call abrtx("MakedMudXNI: not GetGUIDEa")
      end if
      if(.not.GetHf(dummy,.false.)) call MakeHf()
      if(.not.GetHfa(dummy,.false.)) call MakeHfa()
C
      call NEWMEM(ldMdX,ngeom*3)
      call VCLR(XX(ldMdX),1,ngeom*3)
C
      call CompdMudXNI(XX(ldMdX),norb,nocc,ibfld,nfld,ngeom)
      call PostdMudX(XX(ldMdX))
C
      call DELETEMEM(ldMdX,ngeom*3)
      CALL TIMIT(1)
C
      end
C
C*MODULE TDXNI   *DECK PostdMudX
      subroutine PostdMudX(dmdx)
      implicit none
C     Variables
      double precision dmdx(3,*)
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
      integer ldMdQ,lQ,lW
      logical bUseQ
C     Called subprograms
      logical GetFNM
      integer IGetOLI
C     Body of the subprogram
      bUseQ = IGetOLI(4).ne.0
C
      if(bUseQ) then
         call NEWMEM(lW,nmodes)
         if(.not.GetFNM(XX(lW),.true.,XX(lW),.false.,nxyz,nmodes)) then
            call abrtx("PostdMudX: not GetFNM:1")
         end if
         call PrintdMudQ(iw,XX(lW),dMdX,nmodes,
     &                   'Dipole derivatives [in au]')
         call DELETEMEM(lW,nmodes)
         call SetdMudQ(dMdX,nmodes)
      else
         call PrintdMudX(iw,dMdX,ngeom,
     &                   'Dipole derivatives [in au]')
         call NEWMEM(ldMdQ,nmodes*3)
         call NEWMEM(lQ,nxyz*nmodes)
         call NEWMEM(lW,nmodes)
         if(.not.GetFNM(XX(lW),.true.,XX(lQ),.true.,nxyz,nmodes)) then
            call abrtx("PostdMudX: not GetFNM:2")
         end if
         call Project2Q(XX(lQ),dMdX,XX(ldMdQ),nxyz,nmodes,3,1)
         call PrintdMudQ(iw,XX(lW),XX(ldMdQ),nmodes,
     &                   'Dipole derivatives [in au]')
         call DELETEMEM(lW,nmodes)
         call DELETEMEM(lQ,nxyz*nmodes)
         call SetdMudQ(XX(ldMdQ),nmodes)
         call DELETEMEM(ldMdQ,nmodes*3)
      end if
C
      end
C
C*MODULE TDXNI   *DECK CompdMudXNI
      subroutine CompdMudXNI(dmdx,norb,nocc,ibfld,nfld,ngeom)
      implicit none
C     Arguments
      integer norb,nocc,ibfld,nfld,ngeom
      double precision dmdx(3,ngeom)
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 lHf,lDa,lHfa,lD0
C     Body of the subprogram
      write(iw,3000)
C     --------------
C     Step 1 : Hf*Da
C     --------------
      call NEWMEM(lHf,norb*norb*nfld)
      call NEWMEM(lDa,norb*norb*ngeom)
      call MakedMuDx1(dmdx,XX(lHf),XX(lDa),
     &                norb,ibfld,nfld,ngeom)
      call DELETEMEM(lDa,norb*norb*ngeom)
      call DELETEMEM(lHf,norb*norb*nfld)
C     ---------------
C     Step 2 : Hfa*D0
C     ---------------
      call NEWMEM(lHfa,norb*norb*ngeom*nfld)
      call NEWMEM(lD0,norb*norb)
      call MakedMuDx2(dmdx,XX(lHfa),XX(lD0),
     &                norb,nocc,ibfld,nfld,ngeom)
      call DELETEMEM(lD0,norb*norb)
      call DELETEMEM(lHfa,norb*norb*ngeom*nfld)
C     ---------------------------------
C     Step 3 : Add nuclear contribution
C     ---------------------------------
      call MakedMuDx3(dmdx)
C     Format
 3000 FORMAT('Non-Iterative procedure to compute dMudX')
      end
C
C*MODULE TDXNI   *DECK MakedMuDx1
      subroutine MakedMuDx1(dmdx,Hf,Da,
     &                      norb,ibfld,nfld,ngeom)
      implicit none
C     Arguments
      integer norb,ibfld,nfld,ngeom
      double precision dmdx(3,ngeom)
      double precision Hf(norb*norb,nfld)
      double precision Da(norb*norb,ngeom)
C     Variables and constants
      integer igeom,ifld
C     Called subprograms
      double precision DTrace
      logical GetHf,GetGUIDEa
C     Body of the subprogram
      if(.not.GetHf(Hf,.true.)) then
         call abrtx("MakedMuDx1: not GetHf")
      end if
      if(.not.GetGUIDEa(Da,.false.,Da,.false.,Da,.false.,
     &                  Da,.true.,Da,.false.,norb,ngeom)) then
         call abrtx("MakedMuDx1: not GetGUIDEa")
      end if
C
C      call PrintMM(Hf,norb,nfld,'Hf')
C      call PrintMM(Da,norb,ngeom,'Da')
C
      do igeom=1,ngeom
         do ifld=1,nfld
            dmdx(ifld+ibfld-1,igeom)=
     &           -DTrace('S','S',Hf(1,ifld),norb,Da(1,igeom),norb,
     &                   norb,norb)
         end do
      end do
      end
C
C*MODULE TDXNI   *DECK MakedMuDx2
      subroutine MakedMuDx2(dmdx,Hfa,D0,
     &                      norb,nocc,ibfld,nfld,ngeom)
      implicit none
C     Arguments
      integer norb,nocc,ibfld,nfld,ngeom
      double precision dmdx(3,ngeom)
      double precision Hfa(norb*norb,nfld,ngeom)
      double precision D0(norb*norb)
C     Variables and constants
      integer igeom,ifld
C     Called subprograms
      double precision DTrace
      logical GetHfa
C     Body of the subprogram
      call GetLCAO(Hfa,norb)
      call MakeD0(D0,Hfa,norb,nocc)
      if(.not.GetHfa(Hfa,.true.)) then
         call abrtx("MakedMuDx2: not GetHfa")
      end if
C
C      call PrintMM(Hfa,norb,ngeom*nfld,'Hfa')
C      call PrintM(D0,norb,'D0')
C
      do igeom=1,ngeom
         do ifld=1,nfld
            dmdx(ifld+ibfld-1,igeom)=dmdx(ifld+ibfld-1,igeom)
     &        -DTrace('S','S',Hfa(1,ifld,igeom),norb,D0,norb,norb,norb)
         end do
      end do
      end
C
C*MODULE TDXNI   *DECK MakedMuDx3
      subroutine MakedMuDx3(dmdx)
      implicit none
C     Arguments
      double precision dmdx(3,*)
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 MXATM
      parameter (MXATM=2000)
      integer NAT,ICH,MUL,NUM,NQMT,NE,NA,NB,IAN
      double precision ZAN,C
      COMMON /INFOA / NAT,ICH,MUL,NUM,NQMT,NE,NA,NB,
     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
      DOUBLE PRECISION XX
      COMMON /FMCOM / XX(1)
C     Variables and constants
      integer lQ
      integer iat,ifld,ixyz,imode
      double precision Z,dummy
      double precision ZERO
      parameter (ZERO=0.0d+00)
      logical bUseQ
C     Called subprograms
      integer IGetOLI
      logical GetFNM
C     Body of the subprogram
      bUseQ = IGetOLI(4).ne.0
      if(bUseQ) then
         call NEWMEM(lQ,nxyz*nmodes)
         if(.not.GetFNM(XX(lQ),.false.,XX(lQ),.true.,nxyz,nmodes)) then
            call abrtx("MakedMuDx3: not GetFNM")
         end if
         do imode=1,nmodes
            do ifld = ibfld,iefld
               dummy = ZERO
               do iat=1,nat
                  Z = zan(iat)
                  ixyz=(iat-1)*3+ifld
                  dummy=dummy+Z*XX(lQ+(ixyz-1)+(imode-1)*nxyz)
               end do
               dMdX(ifld,imode)=dMdX(ifld,imode)+dummy
            end do
         end do
         call DELETEMEM(lQ,nxyz*nmodes)
      else
         do iat = 1,nat
            Z = zan(iat)
            do ifld = ibfld,iefld
               ixyz=(iat-1)*3+ifld
               dMdX(ifld,ixyz)=dMdX(ifld,ixyz)+Z
            end do
         end do
      end if
C
      end
C
C
C*MODULE TDXNI   *DECK MakedAlphadXNI
      subroutine MakedAlphadXNI(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
      DOUBLE PRECISION XX
      COMMON /FMCOM / XX(1)
C     Variables and constants
      integer ldAdX
      double precision dummy(1)
      double precision w(2)
C     Called subprograms
      logical GetGUIDEa,GetGUDEf,GetHf,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("MakedAlphadXNI: 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 CompdAlphadXNI(XX(ldAdX),w1,
     &                    norb,nocc,ibfld,nfld,ngeom)
      call PostdAlphadX(w,XX(ldAdX))
C
      call DELETEMEM(ldAdX,9*ngeom)
      CALL TIMIT(1)
C
      end
C
C*MODULE TDXNI   *DECK PostdAlphadX
      subroutine PostdAlphadX(w,dAdX)
      implicit none
C     Argument
      double precision w(2)
      double precision dAdX(9,*)
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
      integer ldAdQ,lQ,lW
      logical bUseQ
C     Called subprograms
      logical GetFNM
      integer IGetOLI
C     Body of the subprogram
      bUseQ = IGetOLI(4).ne.0
C
      if(bUseQ) then
         call NEWMEM(lW,nmodes)
         if(.not.GetFNM(XX(lW),.true.,XX(lW),.false.,nxyz,nmodes)) then
            call abrtx("PostdAlphadXNI: not GetFNM:1")
         end if
         call PrintdAlphadQ(iw,XX(lW),dAdX,w,nmodes,
     &                      'Alpha derivatives[in au]')
         call DELETEMEM(lW,nmodes)
         call SetdAlphadQ(w(2),dAdX,nmodes)
      else
         call PrintdAlphadX(iw,dAdX,w,ngeom,
     &                      'Alpha derivatives [in au]')
         call NEWMEM(ldAdQ,9*nmodes)
         call NEWMEM(lQ,nxyz*nmodes)
         call NEWMEM(lW,nmodes)
         if(.not.GetFNM(XX(lW),.true.,XX(lQ),.true.,nxyz,nmodes)) then
            call abrtx("PostdAlphadXNI: not GetFNM:2")
         end if
         call Project2Q(XX(lQ),dAdX,XX(ldAdQ),nxyz,nmodes,9,1)
         call PrintdAlphadQ(iw,XX(lW),XX(ldAdQ),w,nmodes,
     &                      'Alpha derivatives [in au]')
         call DELETEMEM(lW,nmodes)
         call DELETEMEM(lQ,nxyz*nmodes)
         call SetdAlphadQ(w(2),XX(ldAdQ),nmodes)
         call DELETEMEM(ldAdQ,9*nmodes)
      end if
C
      end
C
C*MODULE TDXNI   *DECK CompdAlphadXNI
      subroutine CompdAlphadXNI(dAdX,w,
     &                          norb,nocc,ibfld,nfld,ngeom)
      implicit none
C     Arguments
      double precision w
      integer norb,nocc,ibfld,nfld,ngeom
      double precision dAdX(9,ngeom)
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 lHfa,lDg,lUa,lGa,lIa,lEa,lE0,lUf,lGf,lEf
C     Body of the subprogram
      write(iw,3000) -w,w
C      call MakedAlphaDxPerms(iperms,nterms,ibfld,iefld,nfld)
C     ------------------------
C     Step 1 : Hfa*Dg + Fga*Df
C     ------------------------
      call NEWMEM(lHfa,norb*norb*ngeom*nfld)
      call NEWMEM(lDg,norb*norb*nfld)
      call MakedAlphaDxNI1(dadx,w,XX(lHfa),XX(lDg),
     &                     norb,ibfld,nfld,ngeom)
      call DELETEMEM(lDg,norb*norb*nfld)
      call DELETEMEM(lHfa,norb*norb*ngeom*nfld)
C     --------------------------------------------------
C     Step 2 : w*[Uf+*Ug-Ug+*Uf]Ia - [Uf+*Ug-Ug+*Uf]IaE0
C     --------------------------------------------------
      call NEWMEM(lE0,norb)
      call NEWMEM(lIa,norb*norb*ngeom)
      call NEWMEM(lUf,norb*norb*nfld)
      call MakedAlphaDxNI2(dadx,w,XX(lE0),XX(lIa),XX(lUf),
     &                     norb,nocc,ibfld,nfld,ngeom)
      call DELETEMEM(lUf,norb*norb*nfld)
      call DELETEMEM(lIa,norb*norb*ngeom)
      call DELETEMEM(lE0,norb)
C     --------------------------------------
C     Step 3 : [UxUy+UyUx]Gz + [UxUy+UyUx]Ez
C     --------------------------------------
      call NEWMEM(lUf,norb*norb*nfld)
      call NEWMEM(lGf,norb*norb*nfld)
      call NEWMEM(lEf,norb*norb*nfld)
      call NEWMEM(lUa,norb*norb*ngeom)
      call NEWMEM(lGa,norb*norb*ngeom)
      call NEWMEM(lEa,norb*norb*ngeom)
      call MakedAlphaDxNI3(dadx,w,XX(lUf),XX(lGf),XX(lEf),
     &                     XX(lUa),XX(lGa),XX(lEa),
     &                     norb,nocc,ibfld,nfld,ngeom)
      call DELETEMEM(lEa,norb*norb*ngeom)
      call DELETEMEM(lGa,norb*norb*ngeom)
      call DELETEMEM(lUa,norb*norb*ngeom)
      call DELETEMEM(lEf,norb*norb*nfld)
      call DELETEMEM(lGf,norb*norb*nfld)
      call DELETEMEM(lUf,norb*norb*nfld)
C     Format
 3000 FORMAT('Non-Iterative procedure to compute dAlphadX(',
     &       F13.6,';',F13.6,')')
      end
C
C*MODULE TDXNI   *DECK MakedAlphaDxNI1
      subroutine MakedAlphaDxNI1(dadx,w,Hfa,Dg,
     &                           norb,ibfld,nfld,ngeom)
      implicit none
C     Arguments
      double precision w
      integer norb,ibfld,nfld,ngeom
      double precision dadx(9,ngeom)
      double precision Hfa(norb*norb,nfld,ngeom)
      double precision Dg(norb*norb,nfld)
C     Used common block
      integer NELEMA,IELEMA,INDEXA
      common/NLO_ALPHA/NELEMA,IELEMA(2,9),INDEXA(3,3)
C     Variables and constants
      integer igeom,ifld,jfld,index
      double precision dummy
C     Called subprograms
      logical GetGUDEf,GetHfa,GetFFfa
      double precision DTrace
C     Body of the subprogram
      dummy = 0.0d+00
      if(.not.GetGUDEf(w,dummy,.false.,dummy,.false.,Dg,.true.,
     &                 dummy,.false.,norb,nfld)) then
         call abrtx("MakedAlphaDxNI1: not GetGUDEf")
      end if
      if(.not.GetHfa(Hfa,.true.)) then
         call abrtx("MakedAlphaDxNI1: not GetHfa")
      end if
C     ---------
C     Dg(w)*Hfa
C     ---------
      do igeom=1,ngeom
         do jfld=1,nfld
            do ifld=1,nfld
               dummy=DTrace('N','S',Dg(1,jfld),norb,
     &                      Hfa(1,ifld,igeom),norb,norb,norb)
               index = INDEXA(ifld+ibfld-1,jfld+ibfld-1)
               dadx(index,igeom)=-dummy
            end do
         end do
      end do
      if(.not.GetFFfa(w,Hfa,.true.)) then
         call abrtx("MakedAlphaDxNI1: not GetFFfa")
      end if
C     ---------------------------
C     Df(-w)FFga(w)=Dfd(w)FFga(w)
C     ---------------------------
      do igeom=1,ngeom
         do jfld=1,nfld
            do ifld=1,nfld
               dummy=DTrace('T','S',Dg(1,ifld),norb,
     &                      Hfa(1,jfld,igeom),norb,norb,norb)
               index = INDEXA(ifld+ibfld-1,jfld+ibfld-1)
               dadx(index,igeom)=dadx(index,igeom)-dummy
            end do
         end do
      end do
      end
C
C*MODULE TDXNI   *DECK MakedAlphaDxNI2
      subroutine MakedAlphaDxNI2(dadx,w,E0,Ia,Ug,
     &                           norb,nocc,ibfld,nfld,ngeom)
      implicit none
C     Arguments
      double precision w
      integer norb,nocc,ibfld,nfld,ngeom
      double precision dadx(9,ngeom)
      double precision E0(norb),Ia(norb*norb,ngeom),Ug(norb*norb,nfld)
C     Used common block
      DOUBLE PRECISION XX
      COMMON /FMCOM / XX(1)
      integer NELEMA,IELEMA,INDEXA
      common/NLO_ALPHA/NELEMA,IELEMA(2,9),INDEXA(3,3)
C     Variables and constants
      integer igeom,ifld,jfld,index
      integer lTMP,LVTR
      double precision dummy
      double precision ZERO,ONE,TWO
      parameter (ZERO=0.0d+00,ONE=1.0d+00,TWO=2.0d+00)
      character*1 sN
      parameter (sN='N')
C     Called subprograms
      logical GetGUDEf,GetGUIDEa
      double precision DDOT,DSUM
C     Body of the subprogram
      call GetEnergy(E0,norb)
      dummy = 0.0d+00
      if(.not.GetGUDEf(w,dummy,.false.,Ug,.true.,dummy,.false.,
     &                 dummy,.false.,norb,nfld)) then
         call abrtx("MakedAlphaDxNI2: not GetGUDEf")
      end if
      if(.not.GetGUIDEa(dummy,.false.,dummy,.false.,Ia,.true.,
     &                  dummy,.false.,dummy,.false.,norb,ngeom)) then
         call abrtx("MakedAlphaDxNI2: not GetGUIDEa")
      end if
C
      call NEWMEM(lTMP,norb*nocc)
      call NEWMEM(lVTR,nocc)
C
C     ------------------------------------------
C     -[Udf(w)IaUg(w)+Udg(-w)IaUf(-w)]E0 =
C           -[Udf(w)(IaUg(w))+(Ug(w)Ia)Udf(w)]E0
C          +w[Udf(w)IaUg(w)-Udg(-w)IaUf(-w)]
C          +w[Udf(w)(IaUg(w))-(Ug(w)Ia)Udf(w)]
C     ------------------------------------------
C
      do igeom=1,ngeom
         do jfld=1,nfld
C     -------------
C     [IaUg(w)]=TMP
C     -------------
            call DGEMM(sN,sN,norb,nocc,norb,
     &                 ONE,Ia(1,igeom),norb,Ug(1,jfld),norb,
     &                 ZERO,XX(lTMP),norb)
            do ifld=1,nfld
C     ----------------
C     Udf(w) TMP = VTR
C     ----------------
               call VTRACE(XX(lVTR),'T','N',
     &                     Ug(1,ifld),norb,XX(lTMP),norb,nocc,norb)
C     --------------------
C     -VTR E0 + w sum(VTR)
C     --------------------
               dummy=-DDOT(nocc,XX(lVTR),1,E0,1)
     &               +DSUM(nocc,XX(lVTR),1)*w
               index = INDEXA(ifld+ibfld-1,jfld+ibfld-1)
               dadx(index,igeom)=dadx(index,igeom)-TWO*dummy
            end do
C     -------------
C     [Ug(w)Ia]=TMP
C     -------------
            call DGEMM(sN,sN,nocc,norb,norb,
     &                 ONE,Ug(1,jfld),norb,Ia(1,igeom),norb,
     &                 ZERO,XX(lTMP),nocc)
            do ifld=1,nfld
C     ----------------
C     TMP Udf(w) = VTR
C     ----------------
               call VTRACE(XX(lVTR),'N','T',
     &                     XX(lTMP),nocc,Ug(1,ifld),norb,nocc,norb)
C     --------------------
C     -VTR E0 - w sum(VTR)
C     --------------------
               dummy=-DDOT(nocc,XX(lVTR),1,E0,1)
     &               -DSUM(nocc,XX(lVTR),1)*w
               index = INDEXA(ifld+ibfld-1,jfld+ibfld-1)
               dadx(index,igeom)=dadx(index,igeom)-TWO*dummy
            end do
         end do
      end do
C
      call DELETEMEM(lVTR,nocc)
      call DELETEMEM(lTMP,norb*nocc)
      end
C
C*MODULE TDXNI   *DECK MakedAlphaDxNI3
      subroutine MakedAlphaDxNI3(dadx,w,Ug,Gg,Eg,Ua,Ga,Ea,
     &                           norb,nocc,ibfld,nfld,ngeom)
      implicit none
C     Arguments
      double precision w
      integer norb,nocc,ibfld,nfld,ngeom
      double precision dadx(9,ngeom)
      double precision Ug(norb*norb,nfld),Ua(norb*norb,ngeom)
      double precision Gg(norb*norb,nfld),Ga(norb*norb,ngeom)
      double precision Eg(norb*norb,nfld),Ea(norb*norb,ngeom)
C     Used common block
      DOUBLE PRECISION XX
      COMMON /FMCOM / XX(1)
      integer NELEMA,IELEMA,INDEXA
      common/NLO_ALPHA/NELEMA,IELEMA(2,9),INDEXA(3,3)
C     Variables and constants
      integer lUsum
      integer igeom,ifld,jfld,index
      double precision dummy
      double precision ONE,TWO
      parameter (ONE=1.0d+00,TWO=2.0d+00)
C     Called subprograms
      logical GetGUDEf,GetGUIDEa
      double precision DTrace
C     Body of the subprogram
      dummy = 0.0d+00
      if(.not.GetGUDEf(w,Gg,.true.,Ug,.true.,dummy,.false.,
     &                 Eg,.true.,norb,nfld)) then
         call abrtx("MakedAlphaDxNI3: not GetGUDEf")
      end if
      if(.not.GetGUIDEa(Ga,.true.,Ua,.true.,dummy,.false.,
     &                  dummy,.false.,Ea,.true.,norb,ngeom)) then
         call abrtx("MakedAlphaDxNI3: not GetGUIDEa")
      end if
C
      call NEWMEM(lUsum,norb*norb)
C
C     -----------------------------------------------------------
C     [Ugd(-w)*Ua+Uad*Ug(w)]Gf(-w) + [Ugd(-w)*Uad+Ua*Ug(w)]Ef(-w)
C     -----------------------------------------------------------
      do igeom=1,ngeom
         do jfld=1,nfld
            call CompNIUfaSum(XX(lUsum),Ug(1,jfld),'f',Ua(1,igeom),
     &                        nocc,norb,ONE,.true.)
            do ifld=1,nfld
               dummy=DTrace('N','T',XX(lUsum),norb,Gg(1,ifld),norb,
     &                      norb,norb)
               index = INDEXA(ifld+ibfld-1,jfld+ibfld-1)
               dadx(index,igeom)=dadx(index,igeom)-TWO*dummy
            end do
            call CompNIXfaSum(XX(lUsum),Ug(1,jfld),'f',Ua(1,igeom),
     &                        nocc,norb,ONE,.true.)
            do ifld=1,nfld
               dummy=DTrace('N','T',XX(lUsum),nocc,Eg(1,ifld),norb,
     &                      nocc,nocc)
               index = INDEXA(ifld+ibfld-1,jfld+ibfld-1)
               dadx(index,igeom)=dadx(index,igeom)-TWO*dummy
            end do
         end do
      end do
C     ---------------------------------------------------------
C     [Ufd(w)*Ua+Uad*Uf(-w)]Gg(w) + [Ufd(w)*Uad+Ua*Uf(-w)]Eg(w)
C     ---------------------------------------------------------
      do igeom=1,ngeom
         do ifld=1,nfld
            call CompNIUfaSum(XX(lUsum),Ug(1,ifld),'m',Ua(1,igeom),
     &                        nocc,norb,ONE,.true.)
            do jfld=1,nfld
               dummy=DTrace('N','N',XX(lUsum),norb,Gg(1,jfld),norb,
     &                       norb,norb)
               index = INDEXA(ifld+ibfld-1,jfld+ibfld-1)
               dadx(index,igeom)=dadx(index,igeom)-TWO*dummy
            end do
            call CompNIXfaSum(XX(lUsum),Ug(1,ifld),'m',Ua(1,igeom),
     &                        nocc,norb,ONE,.true.)
            do jfld=1,nfld
               dummy=DTrace('N','N',XX(lUsum),nocc,Eg(1,jfld),norb,
     &                       nocc,nocc)
               index = INDEXA(ifld+ibfld-1,jfld+ibfld-1)
               dadx(index,igeom)=dadx(index,igeom)-TWO*dummy
            end do
         end do
      end do
C     -----------------------------------------------------------------
C     [Ufd(w)*Ug(w)+Ugd(-w)*Uf(-w)]Ga + [Ufd(w)*Ug(w)+Ugd(-w)*Uf(-w)]Ea
C     -----------------------------------------------------------------
      do jfld=1,nfld
         do ifld=1,nfld
            call CompNIUfgSum(XX(lUsum),Ug(1,ifld),'m',Ug(1,jfld),'f',
     &                        .false.,nocc,norb,ONE,.true.)
            do igeom=1,ngeom
               dummy=DTrace('N','S',XX(lUsum),norb,Ga(1,igeom),norb,
     &                      norb,norb)
               index = INDEXA(ifld+ibfld-1,jfld+ibfld-1)
               dadx(index,igeom)=dadx(index,igeom)-TWO*dummy
            end do
            call CompNIXfgSum(XX(lUsum),Ug(1,ifld),'m',Ug(1,jfld),'f',
     &                        .false.,nocc,norb,ONE,.true.)
            do igeom=1,ngeom
               dummy=DTrace('N','S',XX(lUsum),nocc,Ea(1,igeom),norb,
     &                      nocc,nocc)
               index = INDEXA(ifld+ibfld-1,jfld+ibfld-1)
               dadx(index,igeom)=dadx(index,igeom)-TWO*dummy
            end do
         end do
      end do
C
      call DELETEMEM(lUsum,norb*norb)
      end
C*MODULE TDXNI   *DECK MakedBetadXNI
      subroutine MakedBetadXNI(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 ldBdX
      double precision w(3)
      double precision dummy(1)
C     Called subprograms
      logical GetHf,GetGUIDEa,GetGUDEf,GetGUXDEfg,GetHfa,GetFFfa
C     Body of the subprogram
      w(1)=-w1-w2
      w(2)=w1
      w(3)=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("MakedBetadXNI: not GetGUIDEa")
      end if
      if(.not.GetGUDEf(w(1),dummy,.false.,dummy,.false.,dummy,.false.,
     &                 dummy,.false.,norb,nfld)) call MakeAlpha(w(1))
      if(.not.GetGUDEf(w(2),dummy,.false.,dummy,.false.,dummy,.false.,
     &                 dummy,.false.,norb,nfld)) call MakeAlpha(w(2))
      if(.not.GetGUDEf(w(3),dummy,.false.,dummy,.false.,dummy,.false.,
     &                 dummy,.false.,norb,nfld)) call MakeAlpha(w(3))
      if(.not.GetGUXDEfg(w(1),w(2),dummy,.false.,dummy,.false.,dummy,
     &                   .false.,dummy,.false.,dummy,.false.,norb,nfld))
     &     call MakeBeta(w(1),w(2))
      if(.not.GetGUXDEfg(w(1),w(3),dummy,.false.,dummy,.false.,dummy,
     &                   .false.,dummy,.false.,dummy,.false.,norb,nfld))
     &     call MakeBeta(w(1),w(3))
      if(.not.GetGUXDEfg(w(2),w(3),dummy,.false.,dummy,.false.,dummy,
     &                   .false.,dummy,.false.,dummy,.false.,norb,nfld))
     &     call MakeBeta(w(2),w(3))
      if(.not.GetHf(dummy,.false.)) call MakeHf()
      if(.not.GetHfa(dummy,.false.)) call MakeHfa()
      if(.not.GetFFfa(w(1),dummy,.false.)) call MakeFFfa(w(1))
      if(.not.GetFFfa(w(2),dummy,.false.)) call MakeFFfa(w(2))
      if(.not.GetFFfa(w(3),dummy,.false.)) call MakeFFfa(w(3))
C
      call NEWMEM(ldBdX,27*ngeom)
      call VCLR(XX(ldBdX),1,27*ngeom)
C
      call CompdBetadXNI(XX(ldBdX),w,
     &                   norb,nocc,ibfld,nfld,ngeom)
C
      call PostdBetadX(w,XX(ldBdX))
C
      call DELETEMEM(ldBdX,27*ngeom)
      CALL TIMIT(1)
C
      end
C
C*MODULE TDXNI   *DECK PostdBetadX
      subroutine PostdBetadX(w,dBdX)
      implicit none
C     Argument
      double precision w(3)
      double precision dBdX(27,*)
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
      integer ldBdQ,lQ,lW
      logical bUseQ
C     Called subprograms
      logical GetFNM
      integer IGetOLI
C     Body of the subprogram
      bUseQ = IGetOLI(4).ne.0
C
      if(bUseQ) then
         call NEWMEM(lW,nmodes)
         if(.not.GetFNM(XX(lW),.true.,XX(lW),.false.,nxyz,nmodes)) then
            call abrtx("PostdBetadX: not GetFNM:1")
         end if
         call PrintdBetadQ(iw,XX(lW),dBdX,w,nmodes,
     &                      'Beta derivatives[in au]')
         call DELETEMEM(lW,nmodes)
         call SetdBetadQ(w(2),w(3),dBdX,nmodes)
      else
         call PrintdBetadX(iw,dBdX,w,ngeom,
     &                      'Beta derivatives [in au]')
         call NEWMEM(ldBdQ,27*nmodes)
         call NEWMEM(lQ,nxyz*nmodes)
         call NEWMEM(lW,nmodes)
         if(.not.GetFNM(XX(lW),.true.,XX(lQ),.true.,nxyz,nmodes)) then
            call abrtx("PostdBetadX: not GetFNM:2")
         end if
         call Project2Q(XX(lQ),dBdX,XX(ldBdQ),nxyz,nmodes,27,1)
         call PrintdBetadQ(iw,XX(lW),XX(ldBdQ),w,nmodes,
     &                     'Beta derivatives [in au]')
         call DELETEMEM(lW,nmodes)
         call DELETEMEM(lQ,nxyz*nmodes)
         call SetdBetadQ(w(2),w(3),XX(ldBdQ),nmodes)
         call DELETEMEM(ldBdQ,27*nmodes)
      end if
C
      end
C
C*MODULE TDXNI   *DECK CompdBetadXNI
      subroutine CompdBetadXNI(dBdX,w,
     &                         norb,nocc,ibfld,nfld,ngeom)
      implicit none
C     Arguments
      integer norb,nocc,ibfld,nfld,ngeom
      double precision w(3)
      double precision dBdX(27,ngeom)
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 lFfa,lDgh,lGa,lUa,lIa,lEa,lGf,lUf,lEf,lE0
      integer lGgh,lEgh,lUgh,lXgh,lUg,lUh
C     Body of the subprogram
      write(iw,3000) w(1),w(2),w(3)
C     ------------------------------
C     Step 1 : FxaDyz + permutations
C     ------------------------------
      call NEWMEM(lFfa,norb*norb*nfld*ngeom)
      call NEWMEM(lDgh,norb*norb*nfld*nfld)
      call MakedBetaDxNI1(dbdx,w,XX(lFfa),XX(lDgh),
     &                    norb,ibfld,nfld,ngeom)
      call DELETEMEM(lDgh,norb*norb*nfld*nfld)
      call DELETEMEM(lFfa,norb*norb*nfld*ngeom)
C     ---------------------------------------
C     Step 2 : UxaGyz + XfaEyz + permutations
C     ---------------------------------------
      call NEWMEM(lUa,norb*norb*ngeom)
      call NEWMEM(lUf,norb*norb*nfld)
      call NEWMEM(lGgh,norb*norb*nfld*nfld)
      call NEWMEM(lEgh,norb*norb*nfld*nfld)
      call MakedBetaDxNI2(dbdx,w,
     &                    XX(lUa),XX(lUf),XX(lGgh),XX(lEgh),
     &                    norb,nocc,ibfld,nfld,ngeom)
      call DELETEMEM(lEgh,norb*norb*nfld*nfld)
      call DELETEMEM(lGgh,norb*norb*nfld*nfld)
      call DELETEMEM(lUf,norb*norb*nfld)
      call DELETEMEM(lUa,norb*norb*ngeom)
C     ---------------------------------------
C     Step 3 : UyzaGx + XyzaEx + permutations
C     ---------------------------------------
      call NEWMEM(lUa,norb*norb*ngeom)
      call NEWMEM(lIa,norb*norb*ngeom)
      call NEWMEM(lGf,norb*norb*nfld)
      call NEWMEM(lEf,norb*norb*nfld)
      call NEWMEM(lUgh,norb*norb*nfld*nfld)
      call NEWMEM(lXgh,norb*norb*nfld*nfld)
      call NEWMEM(lUg,norb*norb*nfld)
      call NEWMEM(lUh,norb*norb*nfld)
      call MakedBetaDxNI3(dbdx,w,
     &                    XX(lUa),XX(lIa),
     &                    XX(lGf),XX(lEf),XX(lUgh),XX(lXgh),
     &                    XX(lUg),XX(lUh),
     &                    norb,nocc,ibfld,nfld,ngeom)
      call DELETEMEM(lUh,norb*norb*nfld)
      call DELETEMEM(lUg,norb*norb*nfld)
      call DELETEMEM(lXgh,norb*norb*nfld*nfld)
      call DELETEMEM(lUgh,norb*norb*nfld*nfld)
      call DELETEMEM(lEf,norb*norb*nfld)
      call DELETEMEM(lGf,norb*norb*nfld)
      call DELETEMEM(lIa,norb*norb*ngeom)
      call DELETEMEM(lUa,norb*norb*ngeom)
C     ------------------------
C     Step 4 : UxyzGa + XxyzEa
C     ------------------------
      call NEWMEM(lGa,norb*norb*ngeom)
      call NEWMEM(lEa,norb*norb*ngeom)
      call NEWMEM(lUf,norb*norb*nfld)
      call NEWMEM(lUgh,norb*norb*nfld*nfld)
      call NEWMEM(lXgh,norb*norb*nfld*nfld)
      call MakedBetaDxNI4(dbdx,w,
     &                    XX(lGa),XX(lEa),XX(lUf),XX(lUgh),XX(lXgh),
     &                    norb,nocc,ibfld,nfld,ngeom)
      call DELETEMEM(lXgh,norb*norb*nfld*nfld)
      call DELETEMEM(lUgh,norb*norb*nfld*nfld)
      call DELETEMEM(lUf,norb*norb*nfld)
      call DELETEMEM(lEa,norb*norb*ngeom)
      call DELETEMEM(lGa,norb*norb*ngeom)
C     -------------------------------------------
C     Step 5 : XxyzaE0 + (w Sxyza + permutations)
C     -------------------------------------------
      call NEWMEM(lE0,norb)
      call NEWMEM(lIa,norb*norb*ngeom)
      call NEWMEM(lUf,norb*norb*nfld)
      call NEWMEM(lUgh,norb*norb*nfld*nfld)
      call NEWMEM(lXgh,norb*norb*nfld*nfld)
      call MakedBetaDxNI5(dbdx,w,
     &                    XX(lE0),XX(lIa),XX(lUf),XX(lUgh),XX(lXgh),
     &                    norb,nocc,ibfld,nfld,ngeom)
      call DELETEMEM(lXgh,norb*norb*nfld*nfld)
      call DELETEMEM(lUgh,norb*norb*nfld*nfld)
      call DELETEMEM(lUf,norb*norb*nfld)
      call DELETEMEM(lIa,norb*norb*ngeom)
      call DELETEMEM(lE0,norb)
C     Format
 3000 FORMAT('Non-Iterative procedure to compute dBetadX(',
     &       F13.6,';',F13.6,',',F13.6,')')
      end
C
C*MODULE TDXNI   *DECK MakedBetadXNIPerms
      subroutine MakedBetadXNIPerms(w,ibfld,nfld,nterms,iflds,
     &                              mults,ws,bgws)
      implicit none
C     Arguments
      double precision w(3)
      integer ibfld,nfld
      integer nterms(3)
      integer iflds(4,*)
      double precision mults(*)
      double precision ws(3,3)
      logical bgws(3,3)
C     Body of the subprogram
      call MakeBetaNIPerms(w,ibfld,nfld,nterms,iflds,mults,
     &                     ws,bgws)
      end
C
C*MODULE TDXNI   *DECK MakedBetaDxNI1
      subroutine MakedBetaDxNI1(dbdx,w,Ffa,Dgh,
     &                          norb,ibfld,nfld,ngeom)
      implicit none
C     Arguments
      integer norb,ibfld,nfld,ngeom
      double precision w(3)
      double precision dBdX(27,ngeom)
      double precision Ffa(norb*norb,nfld,ngeom)
      double precision Dgh(norb*norb,nfld,nfld)
C     Variables and constants
      double precision dummy,mult
      integer igeom,iterm,ibterm,ieterm,ifld,jfld,kfld,index
      integer iperm,nperm
      parameter (nperm=3)
      integer nterms(nperm)
      integer iflds(4,nperm*27)
      double precision mults(nperm*27)
      double precision ws(3,nperm)
      logical bgws(3,nperm)
C     Called subprograms
      logical GetFFfa,GetGUXDEfg
      double precision dtrace
C     Body of the subprogram
      call MakedBetadXNIPerms(w,ibfld,nfld,nterms,iflds,mults,
     &                        ws,bgws)
C
C     ---------------------------------
C     FFfa(wf)Dgh(wg,wh) + permutations
C     ---------------------------------
      ieterm = 0
      dummy = 0.0d+00
      do iperm = 1,nperm
         if(bgws(1,iperm)) then
            if(.not.GetFFfa(ws(1,iperm),Ffa,.true.)) then
               call abrtx("MakedBetaDxNI1: not GetFFfa")
            end if
         end if
         if(bgws(2,iperm).or.bgws(3,iperm)) then
            if(.not.GetGUXDEfg(ws(2,iperm),ws(3,iperm),
     &                         dummy,.false.,dummy,.false.,
     &                         dummy,.false.,Dgh,.true.,dummy,.false.,
     &                         norb,nfld)) then
               call abrtx("MakedBetaDxNI1: not GetGUXDEfg")
            end if
         end if
C
         ibterm=ieterm+1
         ieterm=nterms(iperm)
         do iterm = ibterm,ieterm
            ifld = iflds(1,iterm)
            jfld = iflds(2,iterm)
            kfld = iflds(3,iterm)
            index = iflds(4,iterm)
            mult = mults(iterm)
            do igeom=1,ngeom
               dummy=dtrace('N','N',Ffa(1,ifld,igeom),norb,
     &                      Dgh(1,jfld,kfld),norb,norb,norb)
               dbdx(index,igeom)=dbdx(index,igeom)-mult*dummy
            end do
         end do
      end do
      end
C
C*MODULE TDXNI   *DECK MakedBetaDxNI2
      subroutine MakedBetaDxNI2(dbdx,w,Ua,Uf,Ggh,Egh,
     &                          norb,nocc,ibfld,nfld,ngeom)
      implicit none
C     Arguments
      integer norb,nocc,ibfld,nfld,ngeom
      double precision w(3)
      double precision dBdX(27,ngeom)
      double precision Ua(norb*norb,ngeom)
      double precision Uf(norb*norb,nfld)
      double precision Ggh(norb*norb,nfld,nfld),Egh(norb*norb,nfld,nfld)
C     Used common block
      DOUBLE PRECISION XX
      COMMON /FMCOM / XX(1)
C     Variables and constants
      integer lWork1,lWork2
      double precision dummy,mult
      integer igeom,iterm,ibterm,ieterm,ifld,jfld,kfld,index
      integer iperm,nperm
      parameter (nperm=3)
      integer nterms(nperm)
      integer iflds(4,nperm*27)
      double precision mults(nperm*27)
      double precision ws(3,nperm)
      logical bgws(3,nperm)
      double precision ONE,TWO
      parameter (ONE=1.0d+00,TWO=2.0d+00)
C     Called subprograms
      logical GetGUIDEa,GetGUDEf,GetGUXDEfg
      double precision dtrace
C     Body of the subprogram
      call MakedBetadXNIPerms(w,ibfld,nfld,nterms,iflds,mults,
     &                        ws,bgws)
      dummy = 0.0d+00
      if(.not.GetGUIDEa(dummy,.false.,Ua,.true.,dummy,.false.,
     &                  dummy,.false.,dummy,.false.,norb,ngeom)) then
         call abrtx("MakedBetaDxNI2: not GetGUIDEa")
      end if
C
      call NEWMEM(lWork1,norb*norb)
      call NEWMEM(lWork2,nocc*nocc)
C
C     ----------------------------------------------------
C     Ufa(wf)Ggh(wg,wh) + Xfa(wf)Egh(wg,wh) + permutations
C     ----------------------------------------------------
C
      ieterm = 0
      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("MakedBetaDxNI2: not GetGUDEf")
            end if
         end if
         if(bgws(2,iperm).or.bgws(3,iperm)) then
            if(.not.GetGUXDEfg(ws(2,iperm),ws(3,iperm),
     &                         Ggh,.true.,dummy,.false.,dummy,.false.,
     &                         dummy,.false.,Egh,.true.,
     &                         norb,nfld)) then
               call abrtx("MakedBetaDxNI2: not GetGUXDEfg")
            end if
         end if
C
         ibterm=ieterm+1
         ieterm=nterms(iperm)
         do iterm = ibterm,ieterm
            ifld = iflds(1,iterm)
            jfld = iflds(2,iterm)
            kfld = iflds(3,iterm)
            index = iflds(4,iterm)
            mult = mults(iterm)
            mult = mult*TWO
            do igeom=1,ngeom
               call CompNIUfaSum(XX(lWork1),Uf(1,ifld),'f',Ua(1,igeom),
     &                           nocc,norb,ONE,.true.)
               call CompNIXfaSum(XX(lWork2),Uf(1,ifld),'f',Ua(1,igeom),
     &                           nocc,norb,ONE,.true.)
               dummy= dtrace('N','N',
     &                       XX(lWork1),norb,Ggh(1,jfld,kfld),norb,
     &                       norb,norb)
     &               +dtrace('N','N',
     &                       XX(lWork2),nocc,Egh(1,jfld,kfld),norb,
     &                       nocc,nocc)
               dbdx(index,igeom)=dbdx(index,igeom)-mult*dummy
            end do
         end do
      end do
C
      call DELETEMEM(lWork2,nocc*nocc)
      call DELETEMEM(lWork1,norb*norb)
C
      end
C
C*MODULE TDXNI   *DECK MakedBetaDxNI3
      subroutine MakedBetaDxNI3(dbdx,w,Ua,Ia,Gf,Ef,Ugh,Udgh,Ug,Uh,
     &                          norb,nocc,ibfld,nfld,ngeom)
      implicit none
C     Arguments
      integer norb,nocc,ibfld,nfld,ngeom
      double precision w(3)
      double precision dBdX(27,ngeom)
      double precision Ua(norb*norb,ngeom),Ia(norb*norb,ngeom)
      double precision Gf(norb*norb,nfld),Ef(norb*norb,nfld)
      double precision Ug(norb*norb,nfld),Uh(norb*norb,nfld)
      double precision Ugh(norb*norb,nfld,nfld)
      double precision Udgh(norb*norb,nfld,nfld)
C     Used common block
      DOUBLE PRECISION XX
      COMMON /FMCOM / XX(1)
C     Variables and constants
      integer lWork1,lWork2
      double precision dummy,mult
      integer igeom,iterm,ibterm,ieterm,ifld,jfld,kfld,index
      integer iperm,nperm
      parameter (nperm=3)
      integer nterms(nperm)
      integer iflds(4,nperm*27)
      double precision mults(nperm*27)
      double precision ws(3,nperm)
      logical bgws(3,nperm)
      double precision ONE,TWO,MINUSONE
      parameter (ONE=1.0d+00,TWO=2.0d+00,MINUSONE=-1.0d+00)
      logical bSamejk
C     Called subprograms
      logical GetGUIDEa,GetGUDEf,GetGUXDEfg
      double precision dtrace
C     Body of the subprogram
      call MakedBetadXNIPerms(w,ibfld,nfld,nterms,iflds,mults,
     &                        ws,bgws)
      dummy = 0.0d+00
      if(.not.GetGUIDEa(dummy,.false.,Ua,.true.,Ia,.true.,
     &                  dummy,.false.,dummy,.false.,norb,ngeom)) then
         call abrtx("MakedBetaDxNI3: not GetGUIDEa")
      end if
C
      call NEWMEM(lWork1,norb*norb)
      call NEWMEM(lWork2,nocc*nocc)
C
C     ----------------------------------------------------
C     Ugha(wg,wh)Gf(wf) + Xgha(wg,wh)Ef(wf) + permutations
C     ----------------------------------------------------
C
      ieterm = 0
      do iperm = 1,nperm
         if(bgws(1,iperm)) then
            if(.not.GetGUDEf(ws(1,iperm),
     &                       Gf,.true.,dummy,.false.,
     &                       dummy,.false.,Ef,.true.,
     &                       norb,nfld)) then
               call abrtx("MakedBetaDxNI3: 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
               call abrtx("MakedBetaDxNI3: not GetGUDEf:2")
            end if
         end if
         if(bgws(3,iperm)) then
            if(.not.GetGUDEf(ws(3,iperm),
     &                       dummy,.false.,Uh,.true.,
     &                       dummy,.false.,dummy,.false.,
     &                       norb,nfld)) then
               call abrtx("MakedBetaDxNI3: not GetGUDEf:3")
            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("MakedBetaDxNI3: not GetGUXDEfg")
            end if
            call DAXPY(norb*norb*nfld*nfld,MINUSONE,Ugh,1,Udgh,1)
         end if
C
         ibterm=ieterm+1
         ieterm=nterms(iperm)
         do iterm = ibterm,ieterm
            ifld = iflds(1,iterm)
            jfld = iflds(2,iterm)
            kfld = iflds(3,iterm)
            index = iflds(4,iterm)
            mult = mults(iterm)
            mult = mult*TWO
            bSameJK = (jfld.eq.kfld).and.(ws(2,iperm).eq.ws(3,iperm))
            do igeom=1,ngeom
               call CompNIUfgaSum1(XX(lWork1),
     &                             Ugh(1,jfld,kfld),Udgh(1,jfld,kfld),
     &                             Ua(1,igeom),nocc,norb,ONE,.true.)
               call CompNIXfgaSum1(XX(lWork2),
     &                             Ugh(1,jfld,kfld),Udgh(1,jfld,kfld),
     &                             Ua(1,igeom),nocc,norb,ONE,.true.)
               call CompNIXfgaSum2(XX(lWork2),
     &                             Ug(1,jfld),'f',Uh(1,kfld),'f',
     &                             bSameJK,Ia(1,igeom),
     &                             nocc,norb,ONE,.false.)
C
               dummy= dtrace('N','N',
     &                       XX(lWork1),norb,Gf(1,ifld),norb,norb,norb)
     &               +dtrace('N','N',
     &                       XX(lWork2),nocc,Ef(1,ifld),norb,nocc,nocc)
               dbdx(index,igeom)=dbdx(index,igeom)-mult*dummy
            end do
         end do
      end do
C
      call DELETEMEM(lWork2,nocc*nocc)
      call DELETEMEM(lWork1,norb*norb)
C
      end
C
C*MODULE TDXNI   *DECK MakedBetaDxNI4
      subroutine MakedBetaDxNI4(dbdx,w,Ga,Ea,Uf,Ugh,Udgh,
     &                          norb,nocc,ibfld,nfld,ngeom)
      implicit none
C     Arguments
      integer norb,nocc,ibfld,nfld,ngeom
      double precision w(3)
      double precision dBdX(27,ngeom)
      double precision Ga(norb*norb,ngeom),Ea(norb*norb,ngeom)
      double precision Uf(norb*norb,nfld)
      double precision Ugh(norb*norb,nfld,nfld)
      double precision Udgh(norb*norb,nfld,nfld)
C     Used common block
      DOUBLE PRECISION XX
      COMMON /FMCOM / XX(1)
C     Variables and constants
      integer lWork1,lWork2
      double precision dummy,mult
      integer igeom,iterm,ibterm,ieterm,ifld,jfld,kfld,index
      integer iperm,nperm
      parameter (nperm=3)
      integer nterms(nperm)
      integer iflds(4,nperm*27)
      double precision mults(nperm*27)
      double precision ws(3,nperm)
      logical bgws(3,nperm)
      double precision ONE,TWO,MINUSONE
      parameter (ONE=1.0d+00,TWO=2.0d+00,MINUSONE=-1.0d+00)
C     Called subprograms
      logical GetGUIDEa,GetGUDEf,GetGUXDEfg
      double precision dtrace
C     Body of the subprogram
      call MakedBetadXNIPerms(w,ibfld,nfld,nterms,iflds,mults,
     &                        ws,bgws)
      dummy = 0.0d+00
      if(.not.GetGUIDEa(Ga,.true.,dummy,.false.,dummy,.false.,
     &                  dummy,.false.,Ea,.true.,norb,ngeom)) then
         call abrtx("MakedBetaDxNI4: not GetGUIDEa")
      end if
C
      call NEWMEM(lWork1,norb*norb)
      call NEWMEM(lWork2,nocc*nocc)
C
C     ------------------------------------
C     Ufgh(wf,wg,wh)Ga + Xfgha(wf,wg,wh)Ea
C     ------------------------------------
C
      ieterm = 0
      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("MakedBetaDxNI4: not GetGUDEf")
            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("MakedBetaDxNI4: not GetGUXDEfg")
            end if
            call DAXPY(norb*norb*nfld*nfld,MINUSONE,Ugh,1,Udgh,1)
         end if
C
         ibterm=ieterm+1
         ieterm=nterms(iperm)
         do iterm = ibterm,ieterm
            ifld = iflds(1,iterm)
            jfld = iflds(2,iterm)
            kfld = iflds(3,iterm)
            index = iflds(4,iterm)
            mult = mults(iterm)
            mult = mult*TWO
C
            call CompNIUfghSum(XX(lWork1),Uf(1,ifld),
     &                         Ugh(1,jfld,kfld),Udgh(1,jfld,kfld),
     &                         nocc,norb,ONE,.true.)
            call CompNIXfghSum(XX(lWork2),Uf(1,ifld),
     &                         Ugh(1,jfld,kfld),Udgh(1,jfld,kfld),
     &                         nocc,norb,ONE,.true.)
C
            do igeom=1,ngeom
               dummy= dtrace('N','S',
     &                       XX(lWork1),norb,Ga(1,igeom),norb,norb,norb)
     &               +dtrace('N','S',
     &                       XX(lWork2),nocc,Ea(1,igeom),norb,nocc,nocc)
               dbdx(index,igeom)=dbdx(index,igeom)-mult*dummy
            end do
         end do
      end do
C
      call DELETEMEM(lWork2,nocc*nocc)
      call DELETEMEM(lWork1,norb*norb)
C
      end
C
C*MODULE TDXNI   *DECK MakedBetaDxNI5
      subroutine MakedBetaDxNI5(dbdx,w,E0,Ia,Uf,Ugh,Udgh,
     &                          norb,nocc,ibfld,nfld,ngeom)
      implicit none
C     Arguments
      integer norb,nocc,ibfld,nfld,ngeom
      double precision w(3)
      double precision dBdX(27,ngeom)
      double precision E0(norb)
      double precision Ia(norb*norb,ngeom)
      double precision Uf(norb*norb,nfld)
      double precision Ugh(norb*norb,nfld,nfld)
      double precision Udgh(norb*norb,nfld,nfld)
C     Used common block
      DOUBLE PRECISION XX
      COMMON /FMCOM / XX(1)
C     Variables and constants
      integer lVTR,lTMP
      double precision dummy,mult
      integer igeom,iterm,ibterm,ieterm,ifld,jfld,kfld,index
      integer iperm,nperm
      parameter (nperm=3)
      integer nterms(nperm)
      integer iflds(4,nperm*27)
      double precision mults(nperm*27)
      double precision ws(3,nperm)
      logical bgws(3,nperm)
      double precision ZERO,ONE,TWO,MINUSONE
      parameter (ZERO=0.0d+00,ONE=1.0d+00,TWO=2.0d+00,MINUSONE=-1.0d+00)
      character*1 sN
      parameter (sN='N')
C     Called subprograms
      logical GetGUDEf,GetGUXDEfg,GetGUIDEa
      double precision DDOT,DSUM
C     Body of the subprogram
      call MakedBetadXNIPerms(w,ibfld,nfld,nterms,iflds,mults,
     &                        ws,bgws)
      call GetEnergy(E0,norb)
      dummy = 0.0d+00
      if(.not.GetGUIDEa(dummy,.false.,dummy,.false.,Ia,.true.,
     &                  dummy,.false.,dummy,.false.,norb,ngeom)) then
         call abrtx("MakedBetaDxNI5: not GetGUIDEa")
      end if
C
      call NEWMEM(lTMP,norb*nocc)
      call NEWMEM(lVTR,nocc)
C
C     ------------------------------------------------------------
C     -[Ufd(-wf)(IaUgh(wg,wh))+(Udgh(-wg,-wh)Ia)Uf(wf) + perms]E0
C     +wf[-Ufd(-wf)(IaUgh(wg,wh))+(Udgh(-wg,-wh)Ia)Uf(wf)] + perms
C     ------------------------------------------------------------
      ieterm = 0
      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("MakedBetaDxNI5: not GetGUDEf")
            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("MakedBetaDxNI5: not GetGUXDEfg")
            end if
            call DAXPY(norb*norb*nfld*nfld,MINUSONE,Ugh,1,Udgh,1)
         end if
C
         ibterm=ieterm+1
         ieterm=nterms(iperm)
         do iterm = ibterm,ieterm
            ifld = iflds(1,iterm)
            jfld = iflds(2,iterm)
            kfld = iflds(3,iterm)
            index = iflds(4,iterm)
            mult = mults(iterm)
            mult = mult*TWO
C
            do igeom = 1,ngeom
C     ------------------
C     [IaUgh(wg,wh)]=TMP
C     ------------------
               call DGEMM(sN,sN,norb,nocc,norb,
     &                    ONE,Ia(1,igeom),norb,Ugh(1,jfld,kfld),norb,
     &                    ZERO,XX(lTMP),norb)
C     ---------------------------------
C     Ufd(-wf) TMP = -Uf(wf) TMP = -VTR
C     ---------------------------------
               call VTRACE(XX(lVTR),'N','N',
     &                     Uf(1,ifld),norb,XX(lTMP),norb,nocc,norb)
C     ---------------------
C     -VTR E0 + wf sum(VTR)
C     ---------------------
               dummy= DDOT(nocc,XX(lVTR),1,E0,1)
     &               +DSUM(nocc,XX(lVTR),1)*ws(1,iperm)
               dbdx(index,igeom)=dbdx(index,igeom)-mult*dummy
C     ---------------------
C     [Udgh(-wg,-wh)Ia]=TMP
C     ---------------------
               call DGEMM(sN,sN,nocc,norb,norb,
     &                    ONE,Udgh(1,jfld,kfld),norb,Ia(1,igeom),norb,
     &                    ZERO,XX(lTMP),nocc)
C     ---------------
C     TMP Uf(w) = VTR
C     ---------------
               call VTRACE(XX(lVTR),'N','N',
     &                     XX(lTMP),nocc,Uf(1,ifld),norb,nocc,norb)
C     ---------------------
C     -VTR E0 + wf sum(VTR)
C     ---------------------
               dummy=-DDOT(nocc,XX(lVTR),1,E0,1)
     &               +DSUM(nocc,XX(lVTR),1)*ws(1,iperm)
C
               dbdx(index,igeom)=dbdx(index,igeom)-mult*dummy
            end do
         end do
      end do
C
      call DELETEMEM(lVTR,nocc)
      call DELETEMEM(lTMP,norb*nocc)
C
      end
C
C*MODULE TDXNI   *DECK Maked2MudX2NI
      subroutine Maked2MudX2NI()
      implicit none
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 ld2MdX2
      integer ngeom2
      double precision dummy(1)
      double precision ZERO
      parameter (ZERO=0.0d+00)
C     Called subprograms
      logical GetGUIDEa,GetGUDEf,GetHf,GetHfa,GetHfab,GetFFab,GetXab
C     Body of the subprogram
C
      ngeom2=(ngeom*ngeom+ngeom)/2
C
      dummy(1) = 0.0d+00
      if(.not.GetGUIDEa(dummy,.false.,dummy,.false.,dummy,.false.,
     &                  dummy,.false.,dummy,.false.,norb,ngeom)) then
         call abrtx("Maked2MudX2NI : not GetGUIDEa")
      end if
      if(.not.GetHf(dummy,.false.)) call MakeHf()
      if(.not.GetHfa(dummy,.false.)) call MakeHfa()
      if(.not.GetHfab(dummy,.false.)) call MakeHfab()
      if(.not.GetFFab(dummy,.false.)) call MakeFFab()
      if(.not.GetXab(dummy,.false.)) call MakeXab()
      if(.not.GetGUDEf(ZERO,dummy,.false.,dummy,.false.,dummy,.false.,
     &                 dummy,.false.,norb,nfld)) call MakeAlpha(ZERO)
C
      call NEWMEM(ld2MdX2,ngeom2*3)
      call VCLR(XX(ld2MdX2),1,ngeom2*3)
C
      call Compd2MudX2NI(XX(ld2MdX2),norb,nocc,
     &                   ibfld,nfld,ngeom,ngeom2)
      call Postd2MudX2(XX(ld2MdX2))
C
      call DELETEMEM(ld2MdX2,ngeom2*3)
      CALL TIMIT(1)
C
      end
C
C*MODULE TDXNI   *DECK Postd2MudX2
      subroutine Postd2MudX2(d2mdx2)
      implicit none
C     Variables
      double precision d2mdx2(3,*)
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
      integer ld2MdQ2,lQ,lW
      logical bUseQ
      integer nmode2
C     Called subprograms
      logical GetFNM
      integer IGetOLI
C     Body of the subprogram
      bUseQ = IGetOLI(4).ne.0
      nmode2=(nmodes*nmodes+nmodes)/2
C
      if(bUseQ) then
         call NEWMEM(lW,nmodes)
         if(.not.GetFNM(XX(lW),.true.,XX(lW),.false.,nxyz,nmodes)) then
            call abrtx("Postd2MudX2: not GetFNM:1")
         end if
         call Printd2MudQ2(iw,XX(lW),d2MdX2,nmodes,
     &                     'Dipole second order derivatives [in au]')
         call Setd2MudQ2(d2MdX2,nmode2)
         call DELETEMEM(lW,nmodes)
      else
         call Printd2MudX2(iw,d2MdX2,ngeom,
     &                     'Dipole second order derivatives [in au]')
         call NEWMEM(ld2MdQ2,nmode2*3)
         call NEWMEM(lQ,nxyz*nmodes)
         call NEWMEM(lW,nmodes)
         if(.not.GetFNM(XX(lW),.true.,XX(lQ),.true.,nxyz,nmodes)) then
            call abrtx("Postd2MudX2: not GetFNM:2")
         end if
         call Project2QQ(XX(lQ),d2MdX2,XX(ld2MdQ2),nxyz,nmodes,3,1)
         call Printd2MudQ2(iw,XX(lW),XX(ld2MdQ2),nmodes,
     &                     'Dipole second order derivatives [in au]')
         call DELETEMEM(lW,nmodes)
         call DELETEMEM(lQ,nxyz*nmodes)
         call Setd2MudQ2(d2MdX2,nmode2)
         call DELETEMEM(ld2MdQ2,nmode2*3)
      end if
C
      end
C
C*MODULE TDXNI   *DECK Compd2MudX2NI
      subroutine Compd2MudX2NI(d2mdx2,norb,nocc,
     &                         ibfld,nfld,ngeom,ngeom2)
      implicit none
C     Arguments
      integer norb,nocc,ibfld,nfld,ngeom,ngeom2
      double precision d2mdx2(3,ngeom2)
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 lD0,lHfab,lDa,lHfb,lDf,lFFab,lJab,lE0,lXab
      integer lUa,lGa,lEa,lUf,lGf,lEf
C     Body of the subprogram
      write(iw,3000)
C     ----------------
C     Step 1 : D0*Hfab
C     ----------------
      call NEWMEM(lD0,norb*norb)
      call NEWMEM(lHfab,norb*norb*ngeom2*nfld)
      call Maked2MuDx2A(d2mdx2,XX(lD0),XX(lHfab),
     &                  norb,nocc,ibfld,nfld,ngeom2)
      call DELETEMEM(lHfab,norb*norb*ngeom2*nfld)
      call DELETEMEM(lD0,norb*norb)
C     ----------------------
C     Step 2 : Da*Hfb+Db*Hfa
C     ----------------------
      call NEWMEM(lDa,norb*norb*ngeom)
      call NEWMEM(lHfb,norb*norb*ngeom*nfld)
      call Maked2MuDx2B(d2mdx2,XX(lDa),XX(lHfb),
     &                  norb,ibfld,nfld,ngeom)
      call DELETEMEM(lHfb,norb*norb*ngeom*nfld)
      call DELETEMEM(lDa,norb*norb*ngeom)
C     ----------------
C     Step 3 : Df*FFab
C     ----------------
      call NEWMEM(lDf,norb*norb*nfld)
      call NEWMEM(lFFab,norb*norb*ngeom2)
      call Maked2MuDx2C(d2mdx2,XX(lDf),XX(lFFab),
     &                  norb,ibfld,nfld,ngeom2)
      call DELETEMEM(lFFab,norb*norb*ngeom2)
      call DELETEMEM(lDf,norb*norb*nfld)
C     ---------------------------------------------
C     Step 4 : Ufa*Gb + Xfa*Eb + (a,b) permutations
C     ---------------------------------------------
      call NEWMEM(lUf,norb*norb*nfld)
      call NEWMEM(lUa,norb*norb*ngeom)
      call NEWMEM(lGa,norb*norb*ngeom)
      call NEWMEM(lEa,norb*norb*ngeom)
      call Maked2MuDx2D(d2mdx2,XX(lUf),XX(lUa),XX(lGa),XX(lEa),
     &                  norb,nocc,ibfld,nfld,ngeom)
      call DELETEMEM(lEa,norb*norb*ngeom)
      call DELETEMEM(lGa,norb*norb*ngeom)
      call DELETEMEM(lUa,norb*norb*ngeom)
      call DELETEMEM(lUf,norb*norb*nfld)
C     ------------------------
C     Step 5 : Uab*Gf + Xab*Ef
C     ------------------------
      call NEWMEM(lGf,norb*norb*nfld)
      call NEWMEM(lEf,norb*norb*nfld)
      call NEWMEM(lUa,norb*norb*ngeom)
      call NEWMEM(lXab,norb*norb*ngeom2)
      call Maked2MuDx2E(d2mdx2,XX(lGf),XX(lEf),XX(lUa),XX(lXab),
     &                  norb,nocc,ibfld,nfld,ngeom,ngeom2)
      call DELETEMEM(lXab,norb*norb*ngeom2)
      call DELETEMEM(lUa,norb*norb*ngeom)
      call DELETEMEM(lEf,norb*norb*nfld)
      call DELETEMEM(lGf,norb*norb*nfld)
C     ------------------------
C     Step 6 : Xabf*E0
C     ------------------------
      call NEWMEM(lE0,norb)
      call NEWMEM(lUf,norb*norb*nfld)
      call NEWMEM(lJab,norb*norb*ngeom2)
      call Maked2MuDx2F(d2mdx2,XX(lE0),XX(lUf),XX(lJab),
     &                  norb,nocc,ibfld,nfld,ngeom2)
      call DELETEMEM(lJab,norb*norb*ngeom2)
      call DELETEMEM(lUf,norb*norb*nfld)
      call DELETEMEM(lE0,norb)
C     Format
 3000 FORMAT('Non-Iterative procedure to compute d2MudX2')
      end
C
C*MODULE TDXNI   *DECK Maked2MuDx2A
      subroutine Maked2MuDx2A(d2mdx2,D0,Hfab,
     &                        norb,nocc,ibfld,nfld,ngeom2)
      implicit none
C     Arguments
      integer norb,nocc,ibfld,nfld,ngeom2
      double precision d2mdx2(3,ngeom2)
      double precision D0(norb*norb)
      double precision Hfab(norb*norb,nfld,ngeom2)
C     Variables and constants
      integer ijgeom,ifld
C     Called subprograms
      double precision dtrace
      logical GetHfab
C     Body of the subprogram
C
      call GetLCAO(Hfab,norb)
      call MakeD0(D0,Hfab,norb,nocc)
      if(.not.GetHfab(Hfab,.true.)) then
         call abrtx("Maked2MuDx2A: not GetHfab")
      end if
C
C     ------
C     D0Hfab
C     ------
      do ijgeom=1,ngeom2
         do ifld=1,nfld
            d2mdx2(ifld+ibfld-1,ijgeom)=
     &           -dtrace('S','S',D0,norb,Hfab(1,ifld,ijgeom),norb,
     &                   norb,norb)
         end do
      end do
      end
C
C*MODULE TDXNI   *DECK MakedMuDx2B
      subroutine Maked2MuDx2B(d2mdx2,Da,Hfb,
     &                        norb,ibfld,nfld,ngeom)
      implicit none
C     Arguments
      integer norb,ibfld,nfld,ngeom
      double precision d2mdx2(3,*)
      double precision Da(norb*norb,ngeom)
      double precision Hfb(norb*norb,nfld,ngeom)
C     Variables and constants
      integer igeom,jgeom,ijgeom,ifld
      double precision dummy
C     Called subprograms
      double precision dtrace
      logical GetGUIDEa,GetHfa
C     Body of the subprogram
C
      dummy = 0.0d+00
      if(.not.GetGUIDEa(dummy,.false.,dummy,.false.,dummy,.false.,
     &                  Da,.true.,dummy,.false.,norb,ngeom)) then
         call abrtx("Maked2MuDx2B: GetGUIDEa")
      end if
      if(.not.GetHfa(Hfb,.true.)) then
         call abrtx("Maked2MuDx2B: GetHfa")
      end if
C
C     -----------
C     DaHfb+DbHfa
C     -----------
      ijgeom=0
      do jgeom=1,ngeom
         do igeom=1,jgeom
            ijgeom=ijgeom+1
            do ifld=1,nfld
               dummy=dtrace('S','S',
     &                      Da(1,igeom),norb,Hfb(1,ifld,jgeom),
     &                      norb,norb,norb)
C
               if(igeom.eq.jgeom) then
                  dummy=dummy+dummy
               else
                  dummy=dummy+
     &                  dtrace('S','S',
     &                         Da(1,jgeom),norb,Hfb(1,ifld,igeom),
     &                         norb,norb,norb)
               end if
C
               d2mdx2(ifld+ibfld-1,ijgeom)=d2mdx2(ifld+ibfld-1,ijgeom)
     &              -dummy
            end do
         end do
      end do
      end
C
C
C*MODULE TDXNI   *DECK MakedMuDx2C
      subroutine Maked2MuDx2C(d2mdx2,Df,FFab,
     &                        norb,ibfld,nfld,ngeom2)
      implicit none
C     Arguments
      integer norb,ibfld,nfld,ngeom2
      double precision d2mdx2(3,*)
      double precision Df(norb*norb,nfld)
      double precision FFab(norb*norb,ngeom2)
C     Variables and constants
      integer ijgeom,ifld
      double precision dummy
      double precision ZERO
      parameter (ZERO=0.0d+00)
C     Called subprograms
      double precision dtrace
      logical GetGUDEf,GetFFab
C     Body of the subprogram
C
      dummy = 0.0d+00
      if(.not.GetGUDEf(ZERO,dummy,.false.,dummy,.false.,
     &                 Df,.true.,dummy,.false.,norb,nfld)) then
         call abrtx("Maked2MuDx2C: not GetGUDEf")
      end if
      if(.not.GetFFab(FFab,.true.)) then
         call abrtx("Maked2MuDx2C: not GetFFab")
      end if
C
C     ---------
C     Df(0)FFab
C     ---------
      do ijgeom=1,ngeom2
         do ifld=1,nfld
            dummy= dtrace('S','S',Df(1,ifld),norb,FFab(1,ijgeom),
     &                    norb,norb,norb)
            d2mdx2(ifld+ibfld-1,ijgeom)=d2mdx2(ifld+ibfld-1,ijgeom)
     &           -dummy
         end do
      end do
      end
C
C*MODULE TDXNI   *DECK MakedMuDx2D
      subroutine Maked2MuDx2D(d2mdx2,Uf,Ua,Ga,Ea,
     &                        norb,nocc,ibfld,nfld,ngeom)
      implicit none
C     Arguments
      integer norb,nocc,ibfld,nfld,ngeom
      double precision d2mdx2(3,*)
      double precision Uf(norb*norb,nfld)
      double precision Ua(norb*norb,ngeom)
      double precision Ga(norb*norb,ngeom)
      double precision Ea(norb*norb,ngeom)
C     Used common block
      DOUBLE PRECISION XX
      COMMON /FMCOM / XX(1)
C     Variables and constants
      integer lWork1,lWork2
      integer ifld,igeom,jgeom,ijgeom
      double precision dummy
      double precision ZERO,ONE,TWO
      parameter (ZERO=0.0d+00,ONE=1.0d+00,TWO=2.0d+00)
C     Called subprograms
      logical GetGUIDEa,GetGUDEf
      double precision dtrace
C     Body of the subprogram
C
      dummy = 0.0d+00
      if(.not.GetGUIDEa(Ga,.true.,Ua,.true.,dummy,.false.,
     &                  dummy,.false.,Ea,.true.,norb,ngeom)) then
         call abrtx("Maked2MuDx2D: not GetGUIDEa")
      end if
      if(.not.GetGUDEf(ZERO,dummy,.false.,Uf,.true.,dummy,
     &                 .false.,dummy,.false.,norb,nfld)) then
         call abrtx("Maked2MuDx2D: not GetGUDEf")
      end if
C
      call NEWMEM(lWork1,norb*norb)
      call NEWMEM(lWork2,nocc*nocc)
C
C     -------------------------
C     UfbGa+XfbEa + (a,b) perms
C     -------------------------
C
      do ifld=1,nfld
         do jgeom=1,ngeom
            call CompNIUfaSum(XX(lWork1),Uf(1,ifld),'f',Ua(1,jgeom),
     &                        nocc,norb,ONE,.true.)
            call CompNIXfaSum(XX(lWork2),Uf(1,ifld),'f',Ua(1,jgeom),
     &                        nocc,norb,ONE,.true.)
            do igeom=1,ngeom
               dummy= dtrace('N','S',
     &                       XX(lWork1),norb,Ga(1,igeom),norb,norb,norb)
     &               +dtrace('N','S',
     &                       XX(lWork2),nocc,Ea(1,igeom),norb,nocc,nocc)
               dummy=dummy+dummy
C
               if(igeom.eq.jgeom) then
                  ijgeom = (igeom*igeom+igeom)/2
                  dummy=dummy*TWO
               else if(igeom.gt.jgeom) then
                  ijgeom = (igeom*igeom-igeom)/2+jgeom
               else
                  ijgeom = (jgeom*jgeom-jgeom)/2+igeom
               end if
C
               d2mdx2(ifld+ibfld-1,ijgeom)=d2mdx2(ifld+ibfld-1,ijgeom)
     &              - dummy
            end do
         end do
      end do
C
      call DELETEMEM(lWork2,nocc*nocc)
      call DELETEMEM(lWork1,norb*norb)
C
      end
C
C
C*MODULE TDXNI   *DECK MakedMuDx2E
      subroutine Maked2MuDx2E(d2mdx2,Gf,Ef,Ua,Xab,
     &                        norb,nocc,ibfld,nfld,ngeom,ngeom2)
      implicit none
C     Arguments
      integer norb,nocc,ibfld,nfld,ngeom,ngeom2
      double precision d2mdx2(3,*)
      double precision Gf(norb*norb,nfld)
      double precision Ef(norb*norb,nfld)
      double precision Ua(norb*norb,ngeom)
      double precision Xab(norb*norb,ngeom2)
C     Used common block
      DOUBLE PRECISION XX
      COMMON /FMCOM / XX(1)
C     Variables and constants
      logical bSame
      integer lWork1
      integer ifld,igeom,jgeom,ijgeom
      double precision dummy
      double precision ZERO,ONE
      parameter (ZERO=0.0d+00,ONE=1.0d+00)
C     Called subprograms
      logical GetGUIDEa,GetGUDEf,GetXab
      double precision dtrace
C     Body of the subprogram
C
      dummy = 0.0d+00
      if(.not.GetGUIDEa(dummy,.false.,Ua,.true.,dummy,.false.,
     &                  dummy,.false.,dummy,.false.,norb,ngeom)) then
         call abrtx("Maked2MuDx2E: not GetGUIDEa")
      end if
      if(.not.GetGUDEf(ZERO,Gf,.true.,dummy,.false.,
     &                 dummy,.false.,Ef,.true.,norb,nfld)) then
         call abrtx("Maked2MuDx2E: not GetGUDEf")
      end if
      if(.not.GetXab(Xab,.true.)) then
         call abrtx("Maked2MuDx2E: not GetXab")
      end if
C
      call NEWMEM(lWork1,norb*norb)
C
C     -----------------------
C     [UaUb+UbUa]Gf + [Xab]Ef
C     -----------------------
C
      ijgeom=0
      do jgeom=1,ngeom
         do igeom=1,jgeom
            ijgeom=ijgeom+1
C
            bsame = igeom.eq.jgeom
            call CompNIUabSum(XX(lWork1),Ua(1,igeom),Ua(1,jgeom),bsame,
     &                        nocc,norb,ONE,.true.)
            do ifld=1,nfld
               dummy= dtrace('S','N',
     &                       XX(lWork1),norb,Gf(1,ifld),norb,norb,norb)
     &               +dtrace('S','N',
     &                     Xab(1,ijgeom),norb,Ef(1,ifld),norb,nocc,nocc)
C
               dummy=dummy+dummy
               d2mdx2(ifld+ibfld-1,ijgeom)=d2mdx2(ifld+ibfld-1,ijgeom)
     &              - dummy
            end do
         end do
      end do
C
      call DELETEMEM(lWork1,norb*norb)
C
      end
C
C
C*MODULE TDXNI   *DECK MakedMuDx2F
      subroutine Maked2MuDx2F(d2mdx2,E0,Uf,Jab,
     &                        norb,nocc,ibfld,nfld,ngeom2)
      implicit none
C     Arguments
      integer norb,nocc,ibfld,nfld,ngeom2
      double precision d2mdx2(3,*)
      double precision E0(norb)
      double precision Uf(norb*norb,nfld)
      double precision Jab(norb*norb,ngeom2)
C     Used common block
      DOUBLE PRECISION XX
      COMMON /FMCOM / XX(1)
C     Variables and constants
      integer LVTR
      integer ifld,ijgeom
      double precision dummy
      double precision ZERO
      parameter (ZERO=0.0d+00)
C     Called subprograms
      logical GetGUDEf,GetJab
      double precision DDOT
C     Body of the subprogram
C
      call GetEnergy(E0,norb)
      dummy = 0.0d+00
      if(.not.GetGUDEf(ZERO,dummy,.false.,Uf,.true.,
     &                 dummy,.false.,dummy,.false.,norb,nfld)) then
         call abrtx("Maked2MuDx2F: not GetGUDEf")
      end if
      if(.not.GetJab(Jab,.true.)) then
         call abrtx("Maked2MuDx2F: not GetJab")
      end if
C
      call NEWMEM(lVTR,nocc)
C
C     ------------------------------------------------------------------
C     -[Udf(0)Iab+Udf(0)IaUb+Udf(0)IbUa+IabUf(0)+UdbIaUf(0)UdaIbUf(0)]E0
C     =-[Udf(0)][Iab+IaUb+IbUa=Jab]E0-[Iab+UdbIa+UdaIb=Jdab][Uf(0)]E0
C     ------------------------------------------------------------------
      do ifld=1,nfld
         do ijgeom=1,ngeom2
            call VTRACE(XX(lVTR),'T','N',
     &                  Uf(1,ifld),norb,Jab(1,ijgeom),norb,
     &                  nocc,norb)
            dummy = -DDOT(nocc,XX(lVTR),1,E0,1)
            call VTRACE(XX(lVTR),'T','N',
     &                  Jab(1,ijgeom),norb,Uf(1,ifld),norb,
     &                  nocc,norb)
            dummy = dummy-DDOT(nocc,XX(lVTR),1,E0,1)
C
            dummy=dummy+dummy
            d2mdx2(ifld+ibfld-1,ijgeom)=d2mdx2(ifld+ibfld-1,ijgeom)
     &           - dummy
         end do
      end do
C
      call DELETEMEM(lVTR,nocc)
C
      end
C
C*MODULE TDXNI   *DECK Maked2AlphadX2NI
      subroutine Maked2AlphadX2NI(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
      DOUBLE PRECISION XX
      COMMON /FMCOM / XX(1)
C     Variables and constants
      integer ld2AdX2
      integer ngeom2
      double precision dummy(1)
      double precision w(2)
C     Called subprograms
      logical GetGUIDEa,GetGUDEf,GetGUXDEfg,GetHf,GetHfa
      logical GetHfab,GetFFab,GetXab,GetFFfa,GetFFfab,GetGUXDEfa
C     Body of the subprogram
C
      ngeom2=(ngeom*ngeom+ngeom)/2
C
      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("Maked2AlphadX2NI: not GetGUIDEa")
      end if
      if(.not.GetHf(dummy,.false.)) call MakeHf()
      if(.not.GetHfa(dummy,.false.)) call MakeHfa()
      if(.not.GetHfab(dummy,.false.)) call MakeHfab()
      if(.not.GetFFab(dummy,.false.)) call MakeFFab()
      if(.not.GetXab(dummy,.false.)) call MakeXab()
      if(.not.GetGUDEf(w1,dummy,.false.,dummy,.false.,dummy,.false.,
     &     dummy,.false.,norb,nfld))
     &     call MakeAlpha(w1)
      if(.not.GetGUXDEfg(w(1),w(2),dummy,.false.,dummy,.false.,
     &     dummy,.false.,dummy,.false.,dummy,.false.,norb,nfld))
     &     call MakeBeta(w(1),w(2))
      if(.not.GetFFfa(w1,dummy,.false.)) call MakeFFfa(w1)
      if(.not.GetGUXDEfa(w1,dummy,.false.,dummy,.false.,dummy,.false.,
     &     dummy,.false.,dummy,.false.,norb,nfld,ngeom))
     &     call MakedAlphadX(w1)
      if(.not.GetFFfab(w1,dummy,.false.)) call MakeFFfab(w1)
C
      call NEWMEM(ld2AdX2,ngeom2*9)
      call VCLR(XX(ld2AdX2),1,ngeom2*9)
C
      call Compd2AlphadX2NI(XX(ld2AdX2),w1,norb,nocc,
     &                      ibfld,nfld,ngeom,ngeom2)
      call Postd2AlphadX2(w,XX(ld2AdX2))
C
      call DELETEMEM(ld2AdX2,ngeom2*9)
      CALL TIMIT(1)
C
      end
C
C*MODULE TDXNI   *DECK Postd2AlphadX2
      subroutine Postd2AlphadX2(w,d2AdX2)
      implicit none
C     Argument
      double precision w(2)
      double precision d2AdX2(9,*)
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
      integer ld2AdQ2,lQ,lW
      integer nmode2
      logical bUseQ
C     Called subprograms
      logical GetFNM
      integer IGetOLI
C     Body of the subprogram
      bUseQ = IGetOLI(4).ne.0
C
      nmode2=(nmodes*nmodes+nmodes)/2
C
      if(bUseQ) then
         call NEWMEM(lW,nmodes)
         if(.not.GetFNM(XX(lW),.true.,XX(lW),.false.,nxyz,nmodes)) then
            call abrtx("Postd2AlphadX2: not GetFNM:1")
         end if
         call Printd2AlphadQ2(iw,XX(lW),d2AdX2,w,nmodes,
     &        'Polarizability second order derivatives [in au]')
         call DELETEMEM(lW,nmodes)
         call Setd2AlphadQ2(w(2),d2AdX2,nmode2)
      else
         call Printd2AlphadX2(iw,d2AdX2,w,ngeom,
     &        'Polarizability second order derivatives [in au]')
         call NEWMEM(ld2AdQ2,nmode2*9)
         call NEWMEM(lQ,nxyz*nmodes)
         call NEWMEM(lW,nmodes)
         if(.not.GetFNM(XX(lW),.true.,XX(lQ),.true.,nxyz,nmodes)) then
            call abrtx("Postd2AlphadX2: not GetFNM:2")
         end if
         call Project2QQ(XX(lQ),d2AdX2,XX(ld2AdQ2),nxyz,nmodes,9,1)
         call Printd2AlphadQ2(iw,XX(lW),XX(ld2AdQ2),w,nmodes,
     &        'Polarizability second order derivatives [in au]')
         call DELETEMEM(lW,nmodes)
         call DELETEMEM(lQ,nxyz*nmodes)
         call Setd2AlphadQ2(w(2),XX(ld2AdQ2),nmode2)
         call DELETEMEM(ld2AdQ2,nmode2*9)
      end if
      end
C
C*MODULE TDXNI   *DECK Compd2AlphadX2NI
      subroutine Compd2AlphadX2NI(d2AdX2,w,norb,nocc,
     &                            ibfld,nfld,ngeom,ngeom2)
      implicit none
C     Arguments
      integer norb,nocc,ibfld,nfld,ngeom,ngeom2
      double precision w
      double precision d2adx2(9,ngeom2)
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 lE0
      integer lDg,lHfab,lDga,lHfb,lFFab
      integer lGa,lUa,lIa,lEa
      integer lXab
      integer lGfg,lUfg,lXfg,lDfg,lEfg
      integer lGf,lUf,lEf
      integer lGga,lUga,lXga,lEga
      integer lIab,lJab
C     Body of the subprogram
      write(iw,3000) -w,w
C     -------------------------------------
C     Step 1 : Dg(w)*Hfab + Df(-w)*FFgab(w)
C     -------------------------------------
      call NEWMEM(lDg,norb*norb*nfld)
      call NEWMEM(lHfab,norb*norb*nfld*ngeom2)
      call Maked2AlphadX2NI1(d2adx2,w,XX(lDg),XX(lHfab),
     &                       norb,ibfld,nfld,ngeom2)
      call DELETEMEM(lHfab,norb*norb*nfld*ngeom2)
      call DELETEMEM(lDg,norb*norb*nfld)
C     ----------------------------------------
C     Step 2 : Dga(w)*Hfb + (a,b) permutations
C     ----------------------------------------
      call NEWMEM(lDga,norb*norb*nfld*ngeom)
      call NEWMEM(lHfb,norb*norb*nfld*ngeom)
      call Maked2AlphadX2NI2(d2adx2,w,XX(lDga),XX(lHfb),
     &                       norb,ibfld,nfld,ngeom,ngeom2)
      call DELETEMEM(lHfb,norb*norb*nfld*ngeom)
      call DELETEMEM(lDga,norb*norb*nfld*ngeom)
C     ----------------------
C     Step 3 : Dfg(-w,w)FFab
C     ----------------------
      call NEWMEM(lDfg,norb*norb*nfld*nfld)
      call NEWMEM(lFFab,norb*norb*ngeom2)
      call Maked2AlphadX2NI3(d2adx2,w,XX(lDfg),XX(lFFab),
     &                       norb,ibfld,nfld,ngeom2)
      call DELETEMEM(lFFab,norb*norb*ngeom2)
      call DELETEMEM(lDfg,norb*norb*nfld*nfld)
C     ------------------------------------
C     Step 4 : UabGfg(-w,w) + XabEfg(-w,w)
C     ------------------------------------
      call NEWMEM(lUa,norb*norb*ngeom)
      call NEWMEM(lXab,norb*norb*ngeom2)
      call NEWMEM(lGfg,norb*norb*nfld*nfld)
      call NEWMEM(lEfg,norb*norb*nfld*nfld)
      call Maked2AlphadX2NI4(d2adx2,w,XX(lUa),XX(lXab),
     &                       XX(lGfg),XX(lEfg),
     &                       norb,nocc,ibfld,nfld,ngeom,ngeom2)
      call DELETEMEM(lEfg,norb*norb*nfld*nfld)
      call DELETEMEM(lGfg,norb*norb*nfld*nfld)
      call DELETEMEM(lXab,norb*norb*ngeom2)
      call DELETEMEM(lUa,norb*norb*ngeom)
C     ------------------------------------------------------------
C     Step 5 : Ufb(-w)Gga(w) + Xfb(-w)*Ega(w) + (a,b) permutations
C     ------------------------------------------------------------
      call NEWMEM(lUa,norb*norb*ngeom)
      call NEWMEM(lUf,norb*norb*nfld)
      call NEWMEM(lGga,norb*norb*nfld*ngeom)
      call NEWMEM(lEga,norb*norb*nfld*ngeom)
      call Maked2AlphadX2NI5(d2adx2,w,XX(lUa),XX(lUf),XX(lGga),XX(lEga),
     &                       norb,nocc,ibfld,nfld,ngeom,ngeom2)
      call DELETEMEM(lEga,norb*norb*nfld*ngeom)
      call DELETEMEM(lGga,norb*norb*nfld*ngeom)
      call DELETEMEM(lUf,norb*norb*nfld)
      call DELETEMEM(lUa,norb*norb*ngeom)
C     --------------------------------------
C     Step 6 : Ugab(w)Gf(-w) + Xgab(w)Ef(-w)
C     --------------------------------------
      call NEWMEM(lUa,norb*norb*ngeom)
      call NEWMEM(lUga,norb*norb*nfld*ngeom)
      call NEWMEM(lXga,norb*norb*nfld*ngeom)
      call NEWMEM(lGf,norb*norb*nfld)
      call NEWMEM(lUf,norb*norb*nfld)
      call NEWMEM(lEf,norb*norb*nfld)
      call NEWMEM(lJab,norb*norb*ngeom2)
      call Maked2AlphadX2NI6(d2adx2,w,XX(lUa),XX(lUga),XX(lXga),
     &                       XX(lGf),XX(lUf),XX(lEf),XX(lJab),
     &                       norb,nocc,ibfld,nfld,ngeom,ngeom2)
      call DELETEMEM(lJab,norb*norb*ngeom2)
      call DELETEMEM(lEf,norb*norb*nfld)
      call DELETEMEM(lUf,norb*norb*nfld)
      call DELETEMEM(lGf,norb*norb*nfld)
      call DELETEMEM(lXga,norb*norb*nfld*ngeom)
      call DELETEMEM(lUga,norb*norb*nfld*ngeom)
      call DELETEMEM(lUa,norb*norb*ngeom)
C     ----------------------
C     Step 7 : Xfab(-w)Eg(w)
C     ----------------------
      call NEWMEM(lUf,norb*norb*nfld)
      call NEWMEM(lEf,norb*norb*nfld)
      call NEWMEM(lJab,norb*norb*ngeom2)
      call Maked2AlphadX2NI7(d2adx2,w,XX(lUf),XX(lEf),XX(lJab),
     &                       norb,nocc,ibfld,nfld,ngeom2)
      call DELETEMEM(lJab,norb*norb*ngeom2)
      call DELETEMEM(lEf,norb*norb*nfld)
      call DELETEMEM(lUf,norb*norb*nfld)
C     ---------------------------------------------------------
C     Step 8 : Ufgb(-w,w)Ga + Xfgb(-w,w)Ea + (a,b) permutations
C     ---------------------------------------------------------
      call NEWMEM(lGa,norb*norb*ngeom)
      call NEWMEM(lUa,norb*norb*ngeom)
      call NEWMEM(lIa,norb*norb*ngeom)
      call NEWMEM(lEa,norb*norb*ngeom)
      call NEWMEM(lUf,norb*norb*nfld)
      call NEWMEM(lUfg,norb*norb*nfld*nfld)
      call NEWMEM(lXfg,norb*norb*nfld*nfld)
      call NEWMEM(lUga,norb*norb*nfld*ngeom)
      call NEWMEM(lXga,norb*norb*nfld*ngeom)
      call Maked2AlphadX2NI8(d2adx2,w,XX(lGa),XX(lUa),XX(lIa),XX(lEa),
     &                       XX(lUf),XX(lUfg),XX(lXfg),
     &                       XX(lUga),XX(lXga),
     &                       norb,nocc,ibfld,nfld,ngeom,ngeom2)
      call DELETEMEM(lXga,norb*norb*nfld*ngeom)
      call DELETEMEM(lUga,norb*norb*nfld*ngeom)
      call DELETEMEM(lXfg,norb*norb*nfld*nfld)
      call DELETEMEM(lUfg,norb*norb*nfld*nfld)
      call DELETEMEM(lUf,norb*norb*nfld)
      call DELETEMEM(lEa,norb*norb*ngeom)
      call DELETEMEM(lIa,norb*norb*ngeom)
      call DELETEMEM(lUa,norb*norb*ngeom)
      call DELETEMEM(lGa,norb*norb*ngeom)
C     --------------------------------
C     Step 9 : Xfgab(-w,w)E0 + w Sfgab
C     --------------------------------
      call NEWMEM(lE0,norb)
      call NEWMEM(lJab,norb*norb*ngeom2)
      call NEWMEM(lUfg,norb*norb*nfld*nfld)
      call NEWMEM(lXfg,norb*norb*nfld*nfld)
      call Maked2AlphadX2NI9a(d2adx2,w,XX(lE0),XX(lJab),
     &                        XX(lUfg),XX(lXfg),
     &                        norb,nocc,ibfld,nfld,ngeom2)
      call DELETEMEM(lXfg,norb*norb*nfld*nfld)
      call DELETEMEM(lUfg,norb*norb*nfld*nfld)
      call DELETEMEM(lJab,norb*norb*ngeom2)
      call DELETEMEM(lE0,norb)
C
      call NEWMEM(lE0,norb)
      call NEWMEM(lIa,norb*norb*ngeom)
      call NEWMEM(lIab,norb*norb*ngeom2)
      call NEWMEM(lUf,norb*norb*nfld)
      call NEWMEM(lUga,norb*norb*nfld*ngeom)
      call NEWMEM(lXga,norb*norb*nfld*ngeom)
      call Maked2AlphadX2NI9b(d2adx2,w,XX(lE0),XX(lIa),XX(lIab),XX(lUf),
     &                        XX(lUga),XX(lXga),
     &                        norb,nocc,ibfld,nfld,ngeom,ngeom2)
      call DELETEMEM(lXga,norb*norb*nfld*ngeom)
      call DELETEMEM(lUga,norb*norb*nfld*ngeom)
      call DELETEMEM(lUf,norb*norb*nfld)
      call DELETEMEM(lIab,norb*norb*ngeom2)
      call DELETEMEM(lIa,norb*norb*ngeom)
      call DELETEMEM(lE0,norb)
C
C     Format
 3000 FORMAT('Non-Iterative procedure to compute d2AlphadX2(',
     &       F13.6,';',F13.6,')')
      end
C
C*MODULE TDXNI   *DECK Maked2AlphadX2NI1
      subroutine Maked2AlphadX2NI1(d2adx2,w,Dg,Hfab,
     &                             norb,ibfld,nfld,ngeom2)
      implicit none
C     Arguments
      integer norb,ibfld,nfld,ngeom2
      double precision w
      double precision d2adx2(9,ngeom2)
      double precision Dg(norb*norb,nfld)
      double precision Hfab(norb*norb,nfld,ngeom2)
C     Used common block
      integer NELEMA,IELEMA,INDEXA
      common/NLO_ALPHA/NELEMA,IELEMA(2,9),INDEXA(3,3)
C     Variables and constants
      integer ijgeom,ifld,jfld,index
      double precision dummy
C     Called subprograms
      double precision Dtrace
      logical GetGUDEf,GetHfab,GetFFfab
C     Body of the subprogram
C
      dummy = 0.0d+00
      if(.not.GetGUDEf(w,dummy,.false.,dummy,.false.,Dg,.true.,
     &     dummy,.false.,norb,nfld)) then
         call abrtx("Maked2AlphadX2NI1: not GetGUDEf")
      end if
C     ----------
C     Dg(w)*Hfab
C     ----------
      if(.not.GetHfab(Hfab,.true.)) then
         call abrtx("Maked2AlphadX2NI1: not GetHfab")
      end if
C
      do ijgeom=1,ngeom2
         do jfld=1,nfld
            do ifld=1,nfld
               dummy=Dtrace('N','S',
     &                      Dg(1,jfld),norb,
     &                      Hfab(1,ifld,ijgeom),norb,
     &                       norb,norb)
               index = INDEXA(ifld+ibfld-1,jfld+ibfld-1)
               d2adx2(index,ijgeom)=d2adx2(index,ijgeom)-dummy
            end do
         end do
      end do
C     ---------------------------------
C     Df(-w)*FFgab(w) = Dfd(w)*FFgab(w)
C     ---------------------------------
      if(.not.GetFFfab(w,Hfab,.true.)) then
         call abrtx("Maked2AlphadX2NI1: not GetFFfab")
      end if
C
      do ijgeom=1,ngeom2
         do jfld=1,nfld
            do ifld=1,nfld
               dummy=Dtrace('T','N',
     &                      Dg(1,ifld),norb,
     &                      Hfab(1,jfld,ijgeom),norb,
     &                      norb,norb)
               index = INDEXA(ifld+ibfld-1,jfld+ibfld-1)
               d2adx2(index,ijgeom)=d2adx2(index,ijgeom)-dummy
            end do
         end do
      end do
C
      end
C
C*MODULE TDXNI   *DECK Maked2AlphadX2NI2
      subroutine Maked2AlphadX2NI2(d2adx2,w,Dga,Hfb,
     &                             norb,ibfld,nfld,ngeom,ngeom2)
      implicit none
C     Arguments
      integer norb,ibfld,nfld,ngeom,ngeom2
      double precision w
      double precision d2adx2(9,ngeom2)
      double precision Dga(norb*norb,nfld,ngeom)
      double precision Hfb(norb*norb,nfld,ngeom)
C     Used common block
      integer NELEMA,IELEMA,INDEXA
      common/NLO_ALPHA/NELEMA,IELEMA(2,9),INDEXA(3,3)
C     Variables and constants
      integer ijgeom,igeom,jgeom,ifld,jfld,index
      double precision dummy
C     Called subprograms
      double precision Dtrace
      logical GetGUXDEfa,GetHfa
C     Body of the subprogram
      dummy = 0.0d+00
      if(.not.GetGUXDEfa(w,dummy,.false.,dummy,.false.,dummy,.false.,
     &     Dga,.true.,dummy,.false.,norb,nfld,ngeom)) then
         call abrtx("Maked2AlphadX2NI2: not GetGUXDEfa")
      end if
      if(.not.GetHfa(Hfb,.true.)) then
         call abrtx("Maked2AlphadX2NI2: not GetHfa")
      end if
C     -------------------------------
C     Dga(w)*Hfb + (a,b) permutations
C     -------------------------------
      ijgeom=0
      do jgeom=1,ngeom
         do igeom=1,jgeom
            ijgeom=ijgeom+1
            do jfld=1,nfld
               do ifld=1,nfld
                  dummy=Dtrace('N','S',
     &                         Dga(1,jfld,igeom),norb,
     &                         Hfb(1,ifld,jgeom),norb,
     &                          norb,norb)
                  if(igeom.eq.jgeom) then
                     dummy=dummy+dummy
                  else
                     dummy=dummy+
     &                    Dtrace('N','S',
     &                           Dga(1,jfld,jgeom),norb,
     &                           Hfb(1,ifld,igeom),norb,
     &                           norb,norb)
                  end if
C
                  index = INDEXA(ifld+ibfld-1,jfld+ibfld-1)
                  d2adx2(index,ijgeom)=d2adx2(index,ijgeom)-dummy
               end do
            end do
         end do
      end do
C
      end
C
C*MODULE TDXNI   *DECK Maked2AlphadX2NI3
      subroutine Maked2AlphadX2NI3(d2adx2,w,Dfg,FFab,
     &                             norb,ibfld,nfld,ngeom2)
      implicit none
C     Arguments
      integer norb,ibfld,nfld,ngeom2
      double precision w
      double precision d2adx2(9,ngeom2)
      double precision Dfg(norb*norb,nfld,nfld)
      double precision FFab(norb*norb,ngeom2)
C     Used common block
      integer NELEMA,IELEMA,INDEXA
      common/NLO_ALPHA/NELEMA,IELEMA(2,9),INDEXA(3,3)
C     Variables and constants
      integer ijgeom,ifld,jfld,index
      double precision dummy
C     Called subprograms
      double precision Dtrace
      logical GetGUXDEfg,GetFFab
C     Body of the subprogram
      dummy = 0.0d+00
      if(.not.GetGUXDEfg(-w,w,dummy,.false.,dummy,.false.,dummy,.false.,
     &     Dfg,.true.,dummy,.false.,norb,nfld)) then
         call abrtx("Maked2AlphadX2NI3: not GetGUXDEfg")
      end if
      if(.not.GetFFab(FFab,.true.)) then
         call abrtx("Maked2AlphadX2NI3: not GetFFab")
      end if
C     -------------
C     Dfg(-w,w)FFab
C     -------------
      do ijgeom=1,ngeom2
         do jfld=1,nfld
            do ifld=1,nfld
               dummy= Dtrace('N','S',
     &                       Dfg(1,ifld,jfld),norb,
     &                       FFab(1,ijgeom),norb,
     &                       norb,norb)
               index = INDEXA(ifld+ibfld-1,jfld+ibfld-1)
               d2adx2(index,ijgeom)=d2adx2(index,ijgeom)-dummy
            end do
         end do
      end do
C
      end
C
C*MODULE TDXNI   *DECK Maked2AlphadX2NI4
      subroutine Maked2AlphadX2NI4(d2adx2,w,Ua,Xab,Gfg,Efg,
     &                             norb,nocc,ibfld,nfld,ngeom,ngeom2)
      implicit none
C     Arguments
      integer norb,nocc,ibfld,nfld,ngeom,ngeom2
      double precision w
      double precision d2adx2(9,ngeom2)
      double precision Ua(norb*norb,ngeom)
      double precision Xab(norb*norb,ngeom2)
      double precision Gfg(norb*norb,nfld,nfld)
      double precision Efg(norb*norb,nfld,nfld)
C     Used common blocks
      integer NELEMA,IELEMA,INDEXA
      common/NLO_ALPHA/NELEMA,IELEMA(2,9),INDEXA(3,3)
      DOUBLE PRECISION XX
      COMMON /FMCOM / XX(1)
C     Variables and constants
      integer lWork1
      integer ijgeom,igeom,jgeom,ifld,jfld,index
      double precision dummy
      logical bSame
      double precision ONE
      parameter (ONE=1.0d+00)
C     Called subprograms
      double precision Dtrace
      logical GetGUIDEa,GetXab,GetGUXDEfg
C     Body of the subprogram
      dummy = 0.0d+00
      if(.not.GetGUIDEa(dummy,.false.,Ua,.true.,dummy,.false.,
     &                  dummy,.false.,dummy,.false.,norb,ngeom)) then
         call abrtx("Maked2AlphadX2NI4: not GetGUIDEa")
      end if
      if(.not.GetXab(Xab,.true.)) then
         call abrtx("Maked2AlphadX2NI4: not GetXab")
      end if
      if(.not.GetGUXDEfg(-w,w,Gfg,.true.,dummy,.false.,
     &     dummy,.false.,dummy,.false.,Efg,.true.,norb,nfld)) then
         call abrtx("Maked2AlphadX2NI4: not GetGUXDEfg")
      end if
C
      call NEWMEM(lWork1,norb*norb)
C
C     ---------------------------
C     UabGfg(-w,w) + XabEfg(-w,w)
C     ---------------------------
C
      ijgeom=0
      do jgeom=1,ngeom
         do igeom=1,jgeom
            ijgeom=ijgeom+1
C
            bsame = igeom.eq.jgeom
            call CompNIUabSum(XX(lWork1),Ua(1,igeom),Ua(1,jgeom),bsame,
     &                        nocc,norb,ONE,.true.)
C
            do jfld=1,nfld
               do ifld=1,nfld
                  dummy= Dtrace('S','N',
     &                          XX(lWork1),norb,
     &                          Gfg(1,ifld,jfld),norb,
     &                          norb,norb)
     &                  +Dtrace('S','N',
     &                          Xab(1,ijgeom),norb,
     &                          Efg(1,ifld,jfld),norb,
     &                          nocc,nocc)
C
                  dummy=dummy+dummy
                  index = INDEXA(ifld+ibfld-1,jfld+ibfld-1)
                  d2adx2(index,ijgeom)=d2adx2(index,ijgeom)-dummy
               end do
            end do
         end do
      end do
C
      call DELETEMEM(lWork1,norb*norb)
C
      end
C
C*MODULE TDXNI   *DECK Maked2AlphadX2NI5
      subroutine Maked2AlphadX2NI5(d2adx2,w,Ub,Ug,Gga,Ega,
     &                             norb,nocc,ibfld,nfld,ngeom,ngeom2)
      implicit none
C     Arguments
      integer norb,nocc,ibfld,nfld,ngeom,ngeom2
      double precision w
      double precision d2adx2(9,ngeom2)
      double precision Ub(norb*norb,ngeom)
      double precision Ug(norb*norb,nfld)
      double precision Gga(norb*norb,nfld,ngeom)
      double precision Ega(norb*norb,nfld,ngeom)
C     Used common blocks
      integer NELEMA,IELEMA,INDEXA
      common/NLO_ALPHA/NELEMA,IELEMA(2,9),INDEXA(3,3)
      DOUBLE PRECISION XX
      COMMON /FMCOM / XX(1)
C     Variables and constants
      integer lWork1,lWork2
      integer ijgeom,igeom,jgeom,ifld,jfld,index
      double precision dummy
      double precision ONE,TWO
      parameter (ONE=1.0d+00,TWO=2.0d+00)
C     Called subprograms
      double precision Dtrace
      logical GetGUIDEa,GetGUDEf,GetGUXDEfa
C     Body of the subprogram
      dummy = 0.0d+00
      if(.not.GetGUIDEa(dummy,.false.,Ub,.true.,dummy,.false.,
     &                  dummy,.false.,dummy,.false.,norb,ngeom)) then
         call abrtx("Maked2AlphadX2NI5: not GetGUIDEa")
      end if
      if(.not.GetGUDEf(w,dummy,.false.,Ug,.true.,dummy,.false.,
     &     dummy,.false.,norb,nfld)) then
         call abrtx("Maked2AlphadX2NI5: not GetGUDEf")
      end if
      if(.not.GetGUXDEfa(w,Gga,.true.,dummy,.false.,dummy,.false.,
     &     dummy,.false.,Ega,.true.,norb,nfld,ngeom)) then
         call abrtx("Maked2AlphadX2NI5: not GetGUXDEfa")
      end if
C
C
      call NEWMEM(lWork1,norb*norb)
      call NEWMEM(lWork2,nocc*nocc)
C
C     ---------------------------------------------------
C     Ufb(-w)Gga(w) + Xfb(-w)*Ega(w) + (a,b) permutations
C     ---------------------------------------------------
      ijgeom=0
      do ifld=1,nfld
         do jgeom=1,ngeom
            call CompNIUfaSum(XX(lWork1),Ug(1,ifld),'m',Ub(1,jgeom),
     &                        nocc,norb,ONE,.true.)
            call CompNIXfaSum(XX(lWork2),Ug(1,ifld),'m',Ub(1,jgeom),
     &                        nocc,norb,ONE,.true.)
C
            do jfld=1,nfld
               do igeom=1,ngeom
                  dummy= Dtrace('N','N',
     &                          XX(lWork1),norb,
     &                          Gga(1,jfld,igeom),norb,
     &                          norb,norb)
     &                  +Dtrace('N','N',
     &                          XX(lWork2),nocc,
     &                          Ega(1,jfld,igeom),norb,
     &                          nocc,nocc)
                  dummy=dummy+dummy
C
                  if(igeom.eq.jgeom) then
                     ijgeom = (igeom*igeom+igeom)/2
                     dummy=dummy*TWO
                  else if(igeom.gt.jgeom) then
                     ijgeom = (igeom*igeom-igeom)/2+jgeom
                  else
                     ijgeom = (jgeom*jgeom-jgeom)/2+igeom
                  end if
C
                  index = INDEXA(ifld+ibfld-1,jfld+ibfld-1)
                  d2adx2(index,ijgeom)=d2adx2(index,ijgeom)-dummy
               end do
            end do
         end do
      end do
C
      call DELETEMEM(lWork2,nocc*nocc)
      call DELETEMEM(lWork1,norb*norb)
C
      end
C
C*MODULE TDXNI   *DECK Maked2AlphadX2NI6
      subroutine Maked2AlphadX2NI6(d2adx2,w,Ua,Uga,Udga,Gg,Ug,Eg,Jab,
     &                             norb,nocc,ibfld,nfld,ngeom,ngeom2)
      implicit none
C     Arguments
      integer norb,nocc,ibfld,nfld,ngeom,ngeom2
      double precision w
      double precision d2adx2(9,ngeom2)
      double precision Ua(norb*norb,ngeom)
      double precision Uga(norb*norb,nfld,ngeom)
      double precision Udga(norb*norb,nfld,ngeom)
      double precision Gg(norb*norb,nfld)
      double precision Ug(norb*norb,nfld)
      double precision Eg(norb*norb,nfld)
      double precision Jab(norb*norb,ngeom2)
C     Used common blocks
      integer NELEMA,IELEMA,INDEXA
      common/NLO_ALPHA/NELEMA,IELEMA(2,9),INDEXA(3,3)
      DOUBLE PRECISION XX
      COMMON /FMCOM / XX(1)
C     Variables and constants
      integer lWork1,lWork2
      integer ijgeom,igeom,jgeom,ifld,jfld,index
      double precision dummy
      logical bSame
      double precision ONE,MINUSONE
      parameter (ONE=1.0d+00,MINUSONE=-1.0d+00)
C     Called subprograms
      double precision Dtrace
      logical GetGUIDEa,GetGUDEf,GetGUXDEfa,GetJab
C     Body of the subprogram
      dummy = 0.0d+00
      if(.not.GetGUIDEa(dummy,.false.,Ua,.true.,dummy,.false.,
     &                  dummy,.false.,dummy,.false.,norb,ngeom)) then
         call abrtx("Maked2AlphadX2NI6: not GetGUIDEa")
      end if
      if(.not.GetGUDEf(w,Gg,.true.,Ug,.true.,dummy,.false.,
     &     Eg,.true.,norb,nfld)) then
         call abrtx("Maked2AlphadX2NI6: not GetGUDEf")
      end if
      if(.not.GetGUXDEfa(w,dummy,.false.,Uga,.true.,Udga,.true.,
     &     dummy,.false.,dummy,.false.,norb,nfld,ngeom)) then
         call abrtx("Maked2AlphadX2NI6: not GetGUXDEfa")
      end if
      call DAXPY(norb*norb*nfld*ngeom,MINUSONE,Uga,1,Udga,1)
      if(.not.GetJab(Jab,.true.)) then
         call abrtx("Maked2AlphadX2NI6: not GetJab")
      end if
C
      call NEWMEM(lWork1,norb*norb)
      call NEWMEM(lWork2,nocc*nocc)
C
C     -------------------------------------------------------------
C     Ugab(w)Gf(-w) + Xgab(w)Ef(-w) = Ugab(w)Gdf(w) + Xgab(w)Edf(w)
C     -------------------------------------------------------------
      do jfld=1,nfld
         ijgeom = 0
         do jgeom=1,ngeom
            do igeom=1,jgeom
               ijgeom=ijgeom+1
C
               bSame = igeom.eq.jgeom
               call CompNIUfabSum(XX(lWork1),
     &                            Uga(1,jfld,igeom),Udga(1,jfld,igeom),
     &                            Ua(1,jgeom),
     &                            Uga(1,jfld,jgeom),Udga(1,jfld,jgeom),
     &                            Ua(1,igeom),
     &                            bSame,nocc,norb,ONE,.true.)
               call CompNIXfabSum1(XX(lWork2),
     &                             Uga(1,jfld,igeom),Udga(1,jfld,igeom),
     &                             Ua(1,jgeom),
     &                             Uga(1,jfld,jgeom),Udga(1,jfld,jgeom),
     &                             Ua(1,igeom),
     &                             bSame,nocc,norb,ONE,.true.)
               call CompNIXfabSum2(XX(lWork2),
     &                             Jab(1,ijgeom),Ug(1,jfld),'f',
     &                             nocc,norb,ONE,.false.)
C
               do ifld=1,nfld
                  dummy= Dtrace('N','T',
     &                          XX(lWork1),norb,
     &                          Gg(1,ifld),norb,
     &                          norb,norb)
     &                  +Dtrace('N','T',
     &                          XX(lWork2),nocc,
     &                          Eg(1,ifld),norb,
     &                          nocc,nocc)
                  dummy=dummy+dummy
C
                  index = INDEXA(ifld+ibfld-1,jfld+ibfld-1)
                  d2adx2(index,ijgeom)=d2adx2(index,ijgeom)-dummy
               end do
            end do
         end do
      end do
C
      call DELETEMEM(lWork2,nocc*nocc)
      call DELETEMEM(lWork1,norb*norb)
C
      end
C
C*MODULE TDXNI   *DECK Maked2AlphadX2NI7
      subroutine Maked2AlphadX2NI7(d2adx2,w,Ug,Eg,Jab,
     &                             norb,nocc,ibfld,nfld,ngeom2)
      implicit none
C     Arguments
      integer norb,nocc,ibfld,nfld,ngeom2
      double precision w
      double precision d2adx2(9,ngeom2)
      double precision Ug(norb*norb,nfld)
      double precision EG(norb*norb,nfld)
      double precision Jab(norb*norb,ngeom2)
C     Used common blocks
      integer NELEMA,IELEMA,INDEXA
      common/NLO_ALPHA/NELEMA,IELEMA(2,9),INDEXA(3,3)
      DOUBLE PRECISION XX
      COMMON /FMCOM / XX(1)
C     Variables and constants
      integer lWork2
      integer ijgeom,ifld,jfld,index
      double precision dummy
      double precision ONE
      parameter (ONE=1.0d+00)
C     Called subprograms
      double precision Dtrace
      logical GetGUDEf,GetJab
C     Body of the subprogram
      dummy = 0.0d+00
      if(.not.GetGUDEf(w,dummy,.false.,Ug,.true.,dummy,.false.,
     &     Eg,.true.,norb,nfld)) then
         call abrtx("Maked2AlphadX2NI7: not GetGUDEf")
      end if
      if(.not.GetJab(Jab,.true.)) then
         call abrtx("Maked2AlphadX2NI7: not GetJab")
      end if
C
      call NEWMEM(lWork2,nocc*nocc)
C
C     ---------------------------------------
C     -[Udf(w)(Iab+IaUb+IbUa)]=-[Udf(w)Jab]
C     -[(Iab+UdbIa+UdaIb)Uf(-w)]=[JdabUdf(w)]
C     ---------------------------------------
C
      do ifld=1,nfld
         do ijgeom=1,ngeom2
            call CompNIXfabSum2(XX(lWork2),
     &                          Jab(1,ijgeom),Ug(1,ifld),'m',
     &                          nocc,norb,ONE,.true.)
C
            do jfld=1,nfld
               dummy= Dtrace('N','N',
     &                       XX(lWork2),nocc,
     &                       Eg(1,jfld),norb,
     &                       nocc,nocc)
               dummy=dummy+dummy
C
               index = INDEXA(ifld+ibfld-1,jfld+ibfld-1)
               d2adx2(index,ijgeom)=d2adx2(index,ijgeom)-dummy
            end do
         end do
      end do
C
      call DELETEMEM(lWork2,nocc*nocc)
C
      end
C
C*MODULE TDXNI   *DECK Maked2AlphadX2NI8
      subroutine Maked2AlphadX2NI8(d2adx2,w,Ga,Ua,Ea,Ia,Ug,
     &                             Ufg,Udfg,Uga,Udga,
     &                             norb,nocc,ibfld,nfld,ngeom,ngeom2)
      implicit none
C     Arguments
      integer norb,nocc,ibfld,nfld,ngeom,ngeom2
      double precision w
      double precision d2adx2(9,ngeom2)
      double precision Ga(norb*norb,ngeom)
      double precision Ua(norb*norb,ngeom)
      double precision Ia(norb*norb,ngeom)
      double precision Ea(norb*norb,ngeom)
      double precision Uga(norb*norb,nfld,ngeom)
      double precision Udga(norb*norb,nfld,ngeom)
      double precision Ug(norb*norb,nfld)
      double precision Ufg(norb*norb,nfld,nfld)
      double precision Udfg(norb*norb,nfld,nfld)
C     Used common blocks
      integer NELEMA,IELEMA,INDEXA
      common/NLO_ALPHA/NELEMA,IELEMA(2,9),INDEXA(3,3)
      DOUBLE PRECISION XX
      COMMON /FMCOM / XX(1)
C     Variables and constants
      integer lWork1,lWork2
      integer ijgeom,igeom,jgeom,ifld,jfld,index
      double precision dummy
      double precision ONE,TWO,MINUSONE
      parameter (ONE=1.0d+00,TWO=2.0d+00,MINUSONE=-1.0d+00)
C     Called subprograms
      double precision Dtrace
      logical GetGUIDEa,GetGUDEf,GetGUXDEfg,GetGUXDEfa
C     Body of the subprogram
      dummy = 0.0d+00
      if(.not.GetGUIDEa(Ga,.true.,Ua,.true.,Ia,.true.,
     &                  dummy,.false.,Ea,.true.,norb,ngeom)) then
         call abrtx("Maked2AlphadX2NI8: not GetGUIDEa")
      end if
      if(.not.GetGUDEf(w,dummy,.false.,Ug,.true.,dummy,.false.,
     &     dummy,.false.,norb,nfld)) then
         call abrtx("Maked2AlphadX2NI8: not GetGUDEf")
      end if
      if(.not.GetGUXDEfg(-w,w,dummy,.false.,Ufg,.true.,Udfg,.true.,
     &     dummy,.false.,dummy,.false.,norb,nfld)) then
         call abrtx("Maked2AlphadX2NI8: not GetGUXDEfg")
      end if
      call DAXPY(norb*norb*nfld*nfld,MINUSONE,Ufg,1,Udfg,1)
      if(.not.GetGUXDEfa(w,dummy,.false.,Uga,.true.,Udga,.true.,
     &     dummy,.false.,dummy,.false.,norb,nfld,ngeom)) then
         call abrtx("Maked2AlphadX2NI8: not GetGUXDEfa")
      end if
      call DAXPY(norb*norb*nfld*ngeom,MINUSONE,Uga,1,Udga,1)
C
      call NEWMEM(lWork1,norb*norb)
      call NEWMEM(lWork2,nocc*nocc)
C
C     ------------------------------------------------
C     Ufgb(-w,w)Ga + Xfgb(-w,w)Ea + (a,b) permutations
C     ------------------------------------------------
      do ifld=1,nfld
         do jfld=1,nfld
            ijgeom = 0
            do jgeom=1,ngeom
C
               call CompNIUfgaSum1(XX(lWork1),
     &                             Ufg(1,ifld,jfld),Udfg(1,ifld,jfld),
     &                             Ua(1,jgeom),nocc,norb,ONE,.true.)
               call CompNIUfgaSum2(XX(lWork1),
     &                             Ug(1,ifld),'m',Uga(1,jfld,jgeom),
     &                             Udga(1,jfld,jgeom),
     &                             nocc,norb,ONE,.false.)
C
               call CompNIXfgaSum1(XX(lWork2),
     &                             Ufg(1,ifld,jfld),Udfg(1,ifld,jfld),
     &                             Ua(1,jgeom),nocc,norb,ONE,.true.)
               call CompNIXfgaSum2(XX(lWork2),
     &                             Ug(1,ifld),'m',Ug(1,jfld),'f',
     &                             .false.,Ia(1,jgeom),
     &                             nocc,norb,ONE,.false.)
               call CompNIXfgaSum3(XX(lWork2),
     &                             Ug(1,ifld),'m',Uga(1,jfld,jgeom),
     &                             Udga(1,jfld,jgeom),
     &                             nocc,norb,ONE,.false.)
C
               do igeom=1,ngeom
C
                  dummy= Dtrace('N','N',
     &                          XX(lWork1),norb,
     &                          Ga(1,igeom),norb,
     &                          norb,norb)
     &                  +Dtrace('N','N',
     &                          XX(lWork2),nocc,
     &                          Ea(1,igeom),norb,
     &                          nocc,nocc)
                  dummy=dummy+dummy
C
                  if(igeom.eq.jgeom) then
                     ijgeom = (igeom*igeom+igeom)/2
                     dummy=dummy*TWO
                  else if(igeom.gt.jgeom) then
                     ijgeom = (igeom*igeom-igeom)/2+jgeom
                  else
                     ijgeom = (jgeom*jgeom-jgeom)/2+igeom
                  end if
C
                  index = INDEXA(ifld+ibfld-1,jfld+ibfld-1)
                  d2adx2(index,ijgeom)=d2adx2(index,ijgeom)-dummy
C
               end do
            end do
         end do
      end do
C
      call DELETEMEM(lWork2,nocc*nocc)
      call DELETEMEM(lWork1,norb*norb)
C
      end
C
C*MODULE TDXNI   *DECK Maked2AlphadX2NI9a
      subroutine Maked2AlphadX2NI9a(d2adx2,w,E0,Jab,Ufg,Udfg,
     &                              norb,nocc,ibfld,nfld,ngeom2)
      implicit none
C     Arguments
      integer norb,nocc,ibfld,nfld,ngeom2
      double precision w
      double precision d2adx2(9,ngeom2)
      double precision E0(norb)
      double precision Jab(norb*norb,ngeom2)
      double precision Ufg(norb*norb,nfld,nfld)
      double precision Udfg(norb*norb,nfld,nfld)
C     Used common blocks
      integer NELEMA,IELEMA,INDEXA
      common/NLO_ALPHA/NELEMA,IELEMA(2,9),INDEXA(3,3)
      DOUBLE PRECISION XX
      COMMON /FMCOM / XX(1)
C     Variables and constants
      integer lVTR
      integer ijgeom,ifld,jfld,index
      double precision dummy
      double precision MINUSONE
      parameter (MINUSONE=-1.0d+00)
C     Called subprograms
      double precision DDOT
      logical GetGUXDEfg,GetJab
C     Body of the subprogram
      call GetEnergy(E0,norb)
      dummy = 0.0d+00
      if(.not.GetGUXDEfg(-w,w,dummy,.false.,Ufg,.true.,Udfg,.true.,
     &     dummy,.false.,dummy,.false.,norb,nfld)) then
         call abrtx("Maked2AlphadX2NI9a: not GetGUXDEfg")
      end if
      call DAXPY(norb*norb*nfld*nfld,MINUSONE,Ufg,1,Udfg,1)
      if(.not.GetJab(Jab,.true.)) then
         call abrtx("Maked2AlphadX2NI9a: not GetJab")
      end if
C
      call NEWMEM(lVTR,nocc)
C
C     ----------------------------------
C     -[Udfg(w,-w) Jab+Jdab Ufg(-w,w)]E0
C     ----------------------------------
      do ifld = 1,nfld
         do jfld = 1,nfld
            do ijgeom = 1,ngeom2
               call VTRACE(XX(lVTR),'N','N',
     &                     Udfg(1,ifld,jfld),norb,Jab(1,ijgeom),norb,
     &                     nocc,norb)
               dummy = -DDOT(nocc,XX(lVTR),1,E0,1)
               call VTRACE(XX(lVTR),'T','N',
     &                     Jab(1,ijgeom),norb,Ufg(1,ifld,jfld),norb,
     &                     nocc,norb)
               dummy = dummy-DDOT(nocc,XX(lVTR),1,E0,1)
               dummy=dummy+dummy
C
               index = INDEXA(ifld+ibfld-1,jfld+ibfld-1)
               d2adx2(index,ijgeom)=d2adx2(index,ijgeom)-dummy
            end do
         end do
      end do
C
      call DELETEMEM(lVTR,nocc)
C
      end
C
C*MODULE TDXNI   *DECK Maked2AlphadX2NI9b
      subroutine Maked2AlphadX2NI9b(d2adx2,w,E0,Ia,Iab,Ug,Uga,Udga,
     &                              norb,nocc,ibfld,nfld,ngeom,ngeom2)
      implicit none
C     Arguments
      integer norb,nocc,ibfld,nfld,ngeom,ngeom2
      double precision w
      double precision d2adx2(9,ngeom2)
      double precision E0(norb)
      double precision Ia(norb*norb,ngeom)
      double precision Iab(norb*norb,ngeom2)
      double precision Ug(norb*norb,nfld)
      double precision Uga(norb*norb,nfld,ngeom)
      double precision Udga(norb*norb,nfld,ngeom)
C     Used common blocks
      integer NELEMA,IELEMA,INDEXA
      common/NLO_ALPHA/NELEMA,IELEMA(2,9),INDEXA(3,3)
      DOUBLE PRECISION XX
      COMMON /FMCOM / XX(1)
C     Variables and constants
      integer lTMP,lVTR
      integer ijgeom,igeom,jgeom,ifld,jfld,index
      double precision dummy
      double precision ZERO,ONE,MINUSONE
      parameter (ZERO=0.0d+00,ONE=1.0d+00,MINUSONE=-1.0d+00)
      character*1 sN
      parameter (sN='N')
C     Called subprograms
      double precision DDOT,DSUM
      logical GetGUIDEa,GetGUDEf,GetGUXDEfa,GetIab
C     Body of the subprogram
      call GetEnergy(E0,norb)
      dummy = 0.0d+00
      if(.not.GetGUIDEa(dummy,.false.,dummy,.false.,Ia,.true.,
     &                  dummy,.false.,dummy,.false.,norb,ngeom)) then
         call abrtx("Maked2AlphadX2NI9b: not GetGUIDEa")
      end if
      if(.not.GetGUDEf(w,dummy,.false.,Ug,.true.,dummy,.false.,
     &     dummy,.false.,norb,nfld)) then
         call abrtx("Maked2AlphadX2NI9b: not GetGUDEf")
      end if
      if(.not.GetGUXDEfa(w,dummy,.false.,Uga,.true.,Udga,.true.,
     &     dummy,.false.,dummy,.false.,norb,nfld,ngeom)) then
         call abrtx("Maked2AlphadX2NI9b: not GetGUXDEfa")
      end if
      call DAXPY(norb*norb*nfld*ngeom,MINUSONE,Uga,1,Udga,1)
      if(.not.GetIab(Iab,.true.)) then
         call abrtx("Maked2AlphadX2NI9b: not GetIab")
      end if
C
      call NEWMEM(lTMP,norb*nocc)
      call NEWMEM(lVTR,nocc)
C
      do jfld = 1,nfld
         ijgeom=0
         do jgeom = 1,ngeom
            do igeom = 1,jgeom
C
               ijgeom=ijgeom+1
C     --------------------------------------
C     [IabUg(w) + IaUgb(w) + IbUga(w)] = TMP
C     --------------------------------------
               call DGEMM(sN,sN,norb,nocc,norb,
     &                    ONE,Iab(1,ijgeom),norb,Ug(1,jfld),norb,
     &                    ZERO,XX(lTMP),norb)
               call DGEMM(sN,sN,norb,nocc,norb,
     &                    ONE,Ia(1,igeom),norb,Uga(1,jfld,jgeom),norb,
     &                    ONE,XX(lTMP),norb)
               call DGEMM(sN,sN,norb,nocc,norb,
     &                    ONE,Ia(1,jgeom),norb,Uga(1,jfld,igeom),norb,
     &                    ONE,XX(lTMP),norb)
               do ifld = 1,nfld
C     ----------------
C     [Udf(w) TMP]=VTR
C     ----------------
                  call VTRACE(XX(lVTR),'T','N',
     &                        Ug(1,ifld),norb,XX(lTMP),norb,nocc,norb)
C     --------------------
C     -VTR E0 + w sum(VTR)
C     --------------------
                  dummy=-DDOT(nocc,XX(lVTR),1,E0,1)
     &                  +DSUM(nocc,XX(lVTR),1)*w
                  dummy=dummy+dummy
                  index = INDEXA(ifld+ibfld-1,jfld+ibfld-1)
                  d2adx2(index,ijgeom)=d2adx2(index,ijgeom)-dummy
               end do
C     ---------------------------------------------
C     [Udg(-w)Iab + Udga(-w)Ib + Udgb(-w)Ia]
C     = [-Ug(w)Iab + Udga(-w)Ib + Udgb(-w)Ia] = TMP
C     ---------------------------------------------
               call DGEMM(sN,sN,nocc,norb,norb,
     &                    MINUSONE,Ug(1,jfld),norb,Iab(1,ijgeom),norb,
     &                    ZERO,XX(lTMP),nocc)
               call DGEMM(sN,sN,nocc,norb,norb,
     &                    ONE,Udga(1,jfld,jgeom),norb,Ia(1,igeom),norb,
     &                    ONE,XX(lTMP),nocc)
               call DGEMM(sN,sN,nocc,norb,norb,
     &                    ONE,Udga(1,jfld,igeom),norb,Ia(1,jgeom),norb,
     &                    ONE,XX(lTMP),nocc)
               do ifld = 1,nfld
C     ---------------------------------
C     [TMP Uf(-w)]=[TMP (-Udf(w))]=-VTR
C     ---------------------------------
                  call VTRACE(XX(lVTR),'N','T',
     &                        XX(lTMP),nocc,Ug(1,ifld),norb,nocc,norb)
C     -------------------
C     VTR E0 - w sum(VTR)
C     -------------------
                  dummy= DDOT(nocc,XX(lVTR),1,E0,1)
     &                  +DSUM(nocc,XX(lVTR),1)*w
                  dummy=dummy+dummy
                  index = INDEXA(ifld+ibfld-1,jfld+ibfld-1)
                  d2adx2(index,ijgeom)=d2adx2(index,ijgeom)-dummy
               end do
C
            end do
         end do
      end do
C
      call DELETEMEM(lVTR,nocc)
      call DELETEMEM(lTMP,norb*nocc)
C
      end
C*MODULE TDXNI   *DECK Maked2BetadX2NI
      subroutine Maked2BetadX2NI(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
      DOUBLE PRECISION XX
      COMMON /FMCOM / XX(1)
C     Variables and constants
      integer ld2BdX2
      integer ngeom2
      double precision dummy(1)
      double precision w(3)
C     Called subprograms
      logical GetGUIDEa,GetGUDEf,GetGUXDEfg,GetGUXDEfgh,GetHf,GetHfa
      logical GetHfab,GetFFab,GetXab,GetFFfa,GetFFfab,GetGUXDEfa
C     Body of the subprogram
      ngeom2=(ngeom*ngeom+ngeom)/2
C
      w(1)=-w1-w2
      w(2)=w1
      w(3)=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("Maked2BetadX2NI: not GetGUIDEa")
      end if
      if(.not.GetHf(dummy,.false.)) call MakeHf()
      if(.not.GetHfa(dummy,.false.)) call MakeHfa()
      if(.not.GetHfab(dummy,.false.)) call MakeHfab()
      if(.not.GetFFab(dummy,.false.)) call MakeFFab()
      if(.not.GetXab(dummy,.false.)) call MakeXab()
      if(.not.GetGUDEf(w(1),dummy,.false.,dummy,.false.,dummy,.false.,
     &     dummy,.false.,norb,nfld))
     &     call MakeAlpha(w(1))
      if(.not.GetGUDEf(w(2),dummy,.false.,dummy,.false.,dummy,.false.,
     &     dummy,.false.,norb,nfld))
     &     call MakeAlpha(w(2))
      if(.not.GetGUDEf(w(3),dummy,.false.,dummy,.false.,dummy,.false.,
     &     dummy,.false.,norb,nfld))
     &     call MakeAlpha(w(3))
      if(.not.GetGUXDEfg(w(1),w(2),dummy,.false.,dummy,.false.,
     &     dummy,.false.,dummy,.false.,dummy,.false.,norb,nfld))
     &     call MakeBeta(w(1),w(2))
      if(.not.GetGUXDEfg(w(1),w(3),dummy,.false.,dummy,.false.,
     &     dummy,.false.,dummy,.false.,dummy,.false.,norb,nfld))
     &     call MakeBeta(w(1),w(3))
      if(.not.GetGUXDEfg(w(2),w(3),dummy,.false.,dummy,.false.,
     &     dummy,.false.,dummy,.false.,dummy,.false.,norb,nfld))
     &     call MakeBeta(w(2),w(3))
      if(.not.GetGUXDEfgh(w(1),w(2),w(3),dummy,.false.,dummy,.false.,
     &     dummy,.false.,dummy,.false.,dummy,.false.,norb,nfld))
     &     call MakeGamma(w(1),w(2),w(3))
      if(.not.GetFFfa(w(1),dummy,.false.)) call MakeFFfa(w(1))
      if(.not.GetFFfa(w(2),dummy,.false.)) call MakeFFfa(w(2))
      if(.not.GetFFfa(w(3),dummy,.false.)) call MakeFFfa(w(3))
      if(.not.GetGUXDEfa(w(1),dummy,.false.,dummy,.false.,dummy,.false.,
     &     dummy,.false.,dummy,.false.,norb,nfld,ngeom))
     &     call MakedAlphadX(w(1))
      if(.not.GetGUXDEfa(w(2),dummy,.false.,dummy,.false.,dummy,.false.,
     &     dummy,.false.,dummy,.false.,norb,nfld,ngeom))
     &     call MakedAlphadX(w(2))
      if(.not.GetGUXDEfa(w(3),dummy,.false.,dummy,.false.,dummy,.false.,
     &     dummy,.false.,dummy,.false.,norb,nfld,ngeom))
     &     call MakedAlphadX(w(3))
      if(.not.GetFFfab(w(1),dummy,.false.)) call MakeFFfab(w(1))
      if(.not.GetFFfab(w(2),dummy,.false.)) call MakeFFfab(w(2))
      if(.not.GetFFfab(w(3),dummy,.false.)) call MakeFFfab(w(3))
C
      call NEWMEM(ld2BdX2,ngeom2*27)
      call VCLR(XX(ld2BdX2),1,ngeom2*27)
C
      call Compd2BetadX2NI(XX(ld2BdX2),w,norb,nocc,nvir,
     &                     ibfld,iefld,nfld,ngeom,ngeom2)
      call Postd2BetadX2NI(w,XX(ld2BdX2))
C
      call DELETEMEM(ld2BdX2,ngeom2*27)
      CALL TIMIT(1)
C
      end
C
C*MODULE TDXNI   *DECK Postd2BetadX2NI
      subroutine Postd2BetadX2NI(w,d2bdx2)
      implicit none
C     Arguments
      double precision w(3)
      double precision d2bdx2(27,*)
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
      integer ld2BdQ2,lQ,lW
      integer nmode2
      logical bUseQ
C     Called subprograms
      logical GetFNM
      integer IGetOLI
C     Body of the subprogram
      bUseQ = IGetOLI(4).ne.0
C
      nmode2=(nmodes*nmodes+nmodes)/2
C
      if(bUseQ) then
         call NEWMEM(lW,nmodes)
         if(.not.GetFNM(XX(lW),.true.,XX(lW),.false.,nxyz,nmodes)) then
            call abrtx("Postd2BetadX2NI: not GetFNM:1")
         end if
         call Printd2BetadQ2(iw,XX(lW),d2BdX2,w,nmodes,
     &     'First Hyperpolarizability second order derivatives [in au]')
         call DELETEMEM(lW,nmodes)
         call Setd2BetadQ2(w(2),w(3),d2bdx2,nmode2)
      else
         call Printd2BetadX2(iw,d2BdX2,w,ngeom,
     &     'First Hyperpolarizability second order derivatives [in au]')
         call NEWMEM(ld2BdQ2,nmode2*27)
         call NEWMEM(lQ,nxyz*nmodes)
         call NEWMEM(lW,nmodes)
         if(.not.GetFNM(XX(lW),.true.,XX(lQ),.true.,nxyz,nmodes)) then
            call abrtx("Postd2BetadX2NI: not GetFNM:2")
         end if
         call Project2QQ(XX(lQ),d2BdX2,XX(ld2BdQ2),nxyz,nmodes,27,1)
         call Printd2BetadQ2(iw,XX(lW),XX(ld2BdQ2),w,nmodes,
     &     'First Hyperpolarizability second order derivatives [in au]')
         call DELETEMEM(lW,nmodes)
         call DELETEMEM(lQ,nxyz*nmodes)
         call Setd2BetadQ2(w(2),w(3),XX(ld2BdQ2),nmode2)
         call DELETEMEM(ld2BdQ2,nmode2*27)
      end if
C
      end
C
C*MODULE TDXNI   *DECK Compd2BetadX2NI
      subroutine Compd2BetadX2NI(d2BdX2,w,norb,nocc,nvir,
     &                           ibfld,iefld,nfld,ngeom,ngeom2)
      implicit none
C     Arguments
      integer norb,nocc,nvir,ibfld,iefld,nfld,ngeom,ngeom2
      double precision w(3)
      double precision d2bdx2(27,ngeom2)
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 lE0
      integer lGa,lUa,lIa,lEa
      integer lFFab,lXab,lIab,lJab
      integer lFFfab
      integer lGf,lUf,lEf,lUg,lUh
      integer lGgh,lUgh,lXgh,lDgh,lEgh
      integer lGfgh,lUfgh,lXfgh,lDfgh,lEfgh
      integer lGfa,lUfa,lXfa,lEfa,lUga,lXga,lUha,lXha
C     Body of the subprogram
      write(iw,3000) w(1),w(2),w(3)
      write(*,*) norb,nocc,nvir,ibfld,iefld,nfld,ngeom,ngeom2
C     ----------------------------
C     Step 1 : Dfgh(w0,w1,w2)*FFab
C     ----------------------------
      write(*,*) "Step 1"
      call NEWMEM(lDfgh,norb*norb*nfld*nfld*nfld)
      call NEWMEM(lFFab,norb*norb*ngeom2)
C     Verified
      if(.true.)
     &call Maked2BetaDx2NI1(d2bdx2,w,XX(lFFab),XX(lDfgh),
     &                      norb,ibfld,nfld,ngeom,ngeom2)
      call DELETEMEM(lFFab,norb*norb*ngeom2)
      call DELETEMEM(lDfgh,norb*norb*nfld*nfld*nfld)
C     ----------------------------------------------------
C     Step 2 : Dgh(w1,w2)*FFfab(w0) + (f,g,h) permutations
C     ----------------------------------------------------
      write(*,*) "Step 2"
      call NEWMEM(lDgh,norb*norb*nfld*nfld)
      call NEWMEM(lFFfab,norb*norb*nfld*ngeom2)
C     Verified
      if(.true.)
     &call Maked2BetaDx2NI2(d2bdx2,w,XX(lFFfab),XX(lDgh),
     &                      norb,ibfld,nfld,ngeom2)
      call DELETEMEM(lFFfab,norb*norb*nfld*ngeom2)
      call DELETEMEM(lDgh,norb*norb*nfld*nfld)
C     ----------------------------------------------
C     Step 3 : UabGfgh(w0,w1,w2) + XabEfgh(w0,w1,w2)
C     ----------------------------------------------
      write(*,*) "Step 3"
      call NEWMEM(lUa,norb*norb*ngeom)
      call NEWMEM(lXab,norb*norb*ngeom2)
      call NEWMEM(lGfgh,norb*norb*nfld*nfld*nfld)
      call NEWMEM(lEfgh,norb*norb*nfld*nfld*nfld)
C     Verified
      if(.true.)
     &call Maked2BetaDx2NI3(d2bdx2,w,XX(lUa),XX(lXab),
     &                      XX(lGfgh),XX(lEfgh),
     &                      norb,nocc,ibfld,
     &                      nfld,ngeom,ngeom2)
      call DELETEMEM(lEfgh,norb*norb*nfld*nfld*nfld)
      call DELETEMEM(lGfgh,norb*norb*nfld*nfld*nfld)
      call DELETEMEM(lXab,norb*norb*ngeom2)
      call DELETEMEM(lUa,norb*norb*ngeom)
C     ---------------------------------------------------------------
C     Step 4 : Ughb(w1,w2)Gfa(w0) + Xghb(w1,w2)Efa(w0) + permutations
C     ---------------------------------------------------------------
      write(*,*) "Step 4"
      call NEWMEM(lGfa,norb*norb*nfld*ngeom)
      call NEWMEM(lEfa,norb*norb*nfld*ngeom)
      call NEWMEM(lUgh,norb*norb*nfld*nfld)
      call NEWMEM(lXgh,norb*norb*nfld*nfld)
      call NEWMEM(lUa,norb*norb*ngeom)
      call NEWMEM(lIa,norb*norb*ngeom)
      call NEWMEM(lUg,norb*norb*nfld)
      call NEWMEM(lUh,norb*norb*nfld)
      call NEWMEM(lUga,norb*norb*nfld*ngeom)
      call NEWMEM(lXga,norb*norb*nfld*ngeom)
      call NEWMEM(lUha,norb*norb*nfld*ngeom)
      call NEWMEM(lXha,norb*norb*nfld*ngeom)
C     Verified
      if(.true.)
     &call Maked2BetaDx2NI4(d2bdx2,w,XX(lUa),XX(lIa),XX(lUg),XX(lUh),
     &                      XX(lUgh),XX(lXgh),XX(lUga),XX(lXga),
     &                      XX(lUha),XX(lXha),XX(lGfa),XX(lEfa),
     &                      norb,nocc,ibfld,nfld,ngeom,ngeom2)
      call DELETEMEM(lXha,norb*norb*nfld*ngeom)
      call DELETEMEM(lUha,norb*norb*nfld*ngeom)
      call DELETEMEM(lXga,norb*norb*nfld*ngeom)
      call DELETEMEM(lUga,norb*norb*nfld*ngeom)
      call DELETEMEM(lUh,norb*norb*nfld)
      call DELETEMEM(lUg,norb*norb*nfld)
      call DELETEMEM(lIa,norb*norb*ngeom)
      call DELETEMEM(lUa,norb*norb*ngeom)
      call DELETEMEM(lXgh,norb*norb*nfld*nfld)
      call DELETEMEM(lUgh,norb*norb*nfld*nfld)
      call DELETEMEM(lEfa,norb*norb*nfld*ngeom)
      call DELETEMEM(lGfa,norb*norb*nfld*ngeom)
C     ---------------------------------------------------------------
C     Step 5 : Ufab(w0)Ggh(w1,w2) + Xfab(w0)Egh(w1,w2) + permutations
C     ---------------------------------------------------------------
      write(*,*) "Step 5"
      call NEWMEM(lGgh,norb*norb*nfld*nfld)
      call NEWMEM(lEgh,norb*norb*nfld*nfld)
      call NEWMEM(lUa,norb*norb*ngeom)
      call NEWMEM(lUf,norb*norb*nfld)
      call NEWMEM(lUfa,norb*norb*nfld*ngeom)
      call NEWMEM(lXfa,norb*norb*nfld*ngeom)
      call NEWMEM(lJab,norb*norb*ngeom2)
C     Verified
      if(.true.)
     &call Maked2BetaDx2NI5(d2bdx2,w,XX(lUa),XX(lJab),XX(lUf),
     &                      XX(lUfa),XX(lXfa),XX(lGgh),XX(lEgh),
     &                      norb,nocc,ibfld,nfld,ngeom,ngeom2)
      call DELETEMEM(lJab,norb*norb*ngeom2)
      call DELETEMEM(lXfa,norb*norb*nfld*ngeom)
      call DELETEMEM(lUfa,norb*norb*nfld*ngeom)
      call DELETEMEM(lUf,norb*norb*nfld)
      call DELETEMEM(lUa,norb*norb*ngeom)
      call DELETEMEM(lEgh,norb*norb*nfld*nfld)
      call DELETEMEM(lGgh,norb*norb*nfld*nfld)
C     ----------------------------------------------
C     Step 6 : Ufghb(w0,w1,w2)Ga + Xfghb(w0,w1,w2)Ea
C                    + (a,b) permutations
C     ----------------------------------------------
      write(*,*) "Step 6a"
      call NEWMEM(lGa,norb*norb*ngeom)
      call NEWMEM(lEa,norb*norb*ngeom)
      call NEWMEM(lUa,norb*norb*ngeom)
      call NEWMEM(lUfgh,norb*norb*nfld*nfld*nfld)
      call NEWMEM(lXfgh,norb*norb*nfld*nfld*nfld)
C     Verified
      if(.true.)
     &call Maked2BetaDx2NI6a(d2bdx2,w,XX(lUa),XX(lUfgh),XX(lXfgh),
     &                       XX(lGa),XX(lEa),
     &                       norb,nocc,ibfld,nfld,ngeom,ngeom2)
      call DELETEMEM(lXfgh,norb*norb*nfld*nfld*nfld)
      call DELETEMEM(lUfgh,norb*norb*nfld*nfld*nfld)
      call DELETEMEM(lUa,norb*norb*ngeom)
      call DELETEMEM(lEa,norb*norb*ngeom)
      call DELETEMEM(lGa,norb*norb*ngeom)
      write(*,*) "Step 6b"
      call NEWMEM(lGa,norb*norb*ngeom)
      call NEWMEM(lEa,norb*norb*ngeom)
      call NEWMEM(lUa,norb*norb*ngeom)
      call NEWMEM(lUf,norb*norb*nfld)
      call NEWMEM(lUfa,norb*norb*nfld*ngeom)
      call NEWMEM(lXfa,norb*norb*nfld*ngeom)
      call NEWMEM(lUgh,norb*norb*nfld*nfld)
      call NEWMEM(lXgh,norb*norb*nfld*nfld)
C     Verified
      if(.true.)
     &call Maked2BetaDx2NI6b(d2bdx2,w,XX(lUa),XX(lUf),XX(lUfa),XX(lXfa),
     &                       XX(lUgh),XX(lXgh),
     &                       XX(lGa),XX(lEa),
     &                       norb,nocc,ibfld,nfld,ngeom,ngeom2)
      call DELETEMEM(lXgh,norb*norb*nfld*nfld)
      call DELETEMEM(lUgh,norb*norb*nfld*nfld)
      call DELETEMEM(lXfa,norb*norb*nfld*ngeom)
      call DELETEMEM(lUfa,norb*norb*nfld*ngeom)
      call DELETEMEM(lUf,norb*norb*nfld)
      call DELETEMEM(lUa,norb*norb*ngeom)
      call DELETEMEM(lEa,norb*norb*ngeom)
      call DELETEMEM(lGa,norb*norb*ngeom)
C     ---------------------------------------------------------------
C     Step 7 : Ughab(w1,w2)Gf(w0) + Xghab(w1,w2)Ef(w0) + permutations
C     ---------------------------------------------------------------
      write(*,*) "Step 7a"
      call NEWMEM(lGf,norb*norb*nfld)
      call NEWMEM(lEf,norb*norb*nfld)
      call NEWMEM(lUga,norb*norb*nfld*ngeom)
      call NEWMEM(lXga,norb*norb*nfld*ngeom)
      call NEWMEM(lUha,norb*norb*nfld*ngeom)
      call NEWMEM(lXha,norb*norb*nfld*ngeom)
      call NEWMEM(lIa,norb*norb*ngeom)
      call NEWMEM(lUg,norb*norb*nfld)
      call NEWMEM(lUh,norb*norb*nfld)
      !Verified
      if(.true.)
     &call Maked2BetaDx2NI7a(d2bdx2,w,XX(lIa),XX(lUg),XX(lUh),
     &                       XX(lUga),XX(lXga),XX(lUha),XX(lXha),
     &                       XX(lGf),XX(lEf),
     &                       norb,nocc,ibfld,nfld,ngeom,ngeom2)
      call DELETEMEM(lUh,norb*norb*nfld)
      call DELETEMEM(lUg,norb*norb*nfld)
      call DELETEMEM(lIa,norb*norb*ngeom)
      call DELETEMEM(lXha,norb*norb*nfld*ngeom)
      call DELETEMEM(lUha,norb*norb*nfld*ngeom)
      call DELETEMEM(lXga,norb*norb*nfld*ngeom)
      call DELETEMEM(lUga,norb*norb*nfld*ngeom)
      call DELETEMEM(lEf,norb*norb*nfld)
      call DELETEMEM(lGf,norb*norb*nfld)
      write(*,*) "Step 7b"
      call NEWMEM(lGf,norb*norb*nfld)
      call NEWMEM(lEf,norb*norb*nfld)
      call NEWMEM(lUgh,norb*norb*nfld*nfld)
      call NEWMEM(lXgh,norb*norb*nfld*nfld)
      call NEWMEM(lJab,norb*norb*ngeom2)
      !Verified
      if(.true.)
     &call Maked2BetaDx2NI7b(d2bdx2,w,XX(lJab),XX(lUgh),XX(lXgh),
     &                       XX(lGf),XX(lEf),
     &                       norb,nocc,ibfld,nfld,ngeom2)
      call DELETEMEM(lJab,norb*norb*ngeom2)
      call DELETEMEM(lXgh,norb*norb*nfld*nfld)
      call DELETEMEM(lUgh,norb*norb*nfld*nfld)
      call DELETEMEM(lEf,norb*norb*nfld)
      call DELETEMEM(lGf,norb*norb*nfld)
      write(*,*) "Step 7c"
      call NEWMEM(lGf,norb*norb*nfld)
      call NEWMEM(lEf,norb*norb*nfld)
      call NEWMEM(lUg,norb*norb*nfld)
      call NEWMEM(lUh,norb*norb*nfld)
      call NEWMEM(lIab,norb*norb*ngeom2)
      !Verified
      if(.true.)
     &call Maked2BetaDx2NI7c(d2bdx2,w,XX(lIab),XX(lUg),XX(lUh),
     &                       XX(lGf),XX(lEf),
     &                       norb,nocc,ibfld,nfld,ngeom2)
      call DELETEMEM(lIab,norb*norb*ngeom2)
      call DELETEMEM(lUh,norb*norb*nfld)
      call DELETEMEM(lUg,norb*norb*nfld)
      call DELETEMEM(lEf,norb*norb*nfld)
      call DELETEMEM(lGf,norb*norb*nfld)
C     ---------------------------------------
C     Step 8 : Xfghab(w0,w1,w2)E0 + wi Sfghab
C     ---------------------------------------
      write(*,*) "Step 8a"
      call NEWMEM(lE0,norb)
      call NEWMEM(lJab,norb*norb*ngeom2)
      call NEWMEM(lUfgh,norb*norb*nfld*nfld*nfld)
      call NEWMEM(lXfgh,norb*norb*nfld*nfld*nfld)
C     Verified
      if(.true.)
     &call Maked2BetadX2NI8a(d2bdx2,w,XX(lE0),XX(lJab),
     &                       XX(lUfgh),XX(lXfgh),
     &                       norb,nocc,ibfld,nfld,ngeom2)
      call DELETEMEM(lXfgh,norb*norb*nfld*nfld*nfld)
      call DELETEMEM(lUfgh,norb*norb*nfld*nfld*nfld)
      call DELETEMEM(lJab,norb*norb*ngeom2)
      call DELETEMEM(lE0,norb)
      write(*,*) "Step 8b"
      call NEWMEM(lE0,norb)
      call NEWMEM(lIa,norb*norb*ngeom)
      call NEWMEM(lIab,norb*norb*ngeom2)
      call NEWMEM(lUf,norb*norb*nfld)
      call NEWMEM(lUfa,norb*norb*nfld*ngeom)
      call NEWMEM(lXfa,norb*norb*nfld*ngeom)
      call NEWMEM(lUgh,norb*norb*nfld*nfld)
      call NEWMEM(lXgh,norb*norb*nfld*nfld)
C     Verified
      if(.true.)
     &call Maked2BetadX2NI8b(d2bdx2,w,XX(lE0),XX(lIa),XX(lIab),XX(lUf),
     &                       XX(lUfa),XX(lXfa),XX(lUgh),XX(lXgh),
     &                       norb,nocc,ibfld,
     &                       nfld,ngeom,ngeom2)
      call DELETEMEM(lXgh,norb*norb*nfld*nfld)
      call DELETEMEM(lUgh,norb*norb*nfld*nfld)
      call DELETEMEM(lXfa,norb*norb*nfld*ngeom)
      call DELETEMEM(lUfa,norb*norb*nfld*ngeom)
      call DELETEMEM(lUf,norb*norb*nfld)
      call DELETEMEM(lIab,norb*norb*ngeom2)
      call DELETEMEM(lIa,norb*norb*ngeom)
      call DELETEMEM(lE0,norb)
C     Format
 3000 FORMAT('Non-Iterative procedure to compute d2BetadX2(',
     &       F13.6,';',F13.6,',',F13.6,')')
      end
C
C*MODULE TDXNI   *DECK Maked2BetadX2NIPerms
      subroutine Maked2BetadX2NIPerms(w,ibfld,nfld,nterms,iflds,
     &                                mults,ws,bgws)
      implicit none
C     Arguments
      double precision w(3)
      integer ibfld,nfld
      integer nterms(3)
      integer iflds(4,*)
      double precision mults(*)
      double precision ws(3,3)
      logical bgws(3,3)
C     Body of the subprogram
      call MakeBetaNIPerms(w,ibfld,nfld,nterms,iflds,mults,
     &                     ws,bgws)
      end
C
C*MODULE TDXNI   *DECK Maked2BetaDx2NI1
      subroutine Maked2BetaDx2NI1(d2bdx2,w,FFab,Dfgh,
     &                            norb,ibfld,nfld,ngeom,ngeom2)
      implicit none
C     Arguments
      integer norb,ibfld,nfld,ngeom,ngeom2
      double precision w(3)
      double precision d2BdX2(27,ngeom)
      double precision FFab(norb*norb,ngeom2)
      double precision Dfgh(norb*norb,nfld,nfld,nfld)
C     Used common block
      integer NELEMB,IELEMB,INDEXB
      common/NLO_BETA/NELEMB,IELEMB(3,27),INDEXB(3,3,3)
C     Variables and constants
      double precision dummy
      integer ijgeom,ifld,jfld,kfld,index
C     Called subprograms
      logical GetFFab,GetGUXDEfgh
      double precision dtrace
C     Body of the subprogram
      if(.not.GetFFab(FFab,.true.)) then
         call abrtx("Maked2BetaDx2NI1: not GetFFab")
      end if
      dummy = 0.0d+00
      if(.not.GetGUXDEfgh(w(1),w(2),w(3),dummy,.false.,dummy,.false.,
     &     dummy,.false.,Dfgh,.true.,dummy,.false.,norb,nfld)) then
         call abrtx("Maked2BetaDx2NI1: not GetGUXDEfgh")
      end if
C     ------------------
C     FFabDfgh(wf,wg,wh)
C     ------------------
      do ijgeom = 1,ngeom2
         do kfld=1,nfld
            do jfld=1,nfld
               do ifld=1,nfld
                  dummy = dtrace('S','N',FFab(1,ijgeom),norb,
     &                           Dfgh(1,ifld,jfld,kfld),norb,norb,norb)
C
                  index = INDEXB(ifld+ibfld-1,jfld+ibfld-1,kfld+ibfld-1)
                  d2bdx2(index,ijgeom)=d2bdx2(index,ijgeom)-dummy
               end do
            end do
         end do
      end do
      end
C
C*MODULE TDXNI   *DECK MakedBetaDxNI2
      subroutine Maked2BetaDx2NI2(d2bdx2,w,FFfab,Dgh,
     &                            norb,ibfld,nfld,ngeom2)
      implicit none
C     Arguments
      integer norb,ibfld,nfld,ngeom2
      double precision w(3)
      double precision d2BdX2(27,ngeom2)
      double precision FFfab(norb*norb,nfld,ngeom2)
      double precision Dgh(norb*norb,nfld,nfld)
C     Variables and constants
      double precision dummy,mult
      integer ijgeom,iterm,ibterm,ieterm,ifld,jfld,kfld,index
      integer iperm,nperm
      parameter (nperm=3)
      integer nterms(nperm)
      integer iflds(4,nperm*27)
      double precision mults(nperm*27)
      double precision ws(3,nperm)
      logical bgws(3,nperm)
C     Called subprograms
      logical GetFFfab,GetGUXDEfg
      double precision dtrace
C     Body of the subprogram
      call Maked2BetadX2NIPerms(w,ibfld,nfld,nterms,iflds,mults,
     &                          ws,bgws)
C
C     ----------------------------------
C     FFfab(wf)Dgh(wg,wh) + permutations
C     ----------------------------------
      ieterm = 0
      dummy = 0.0d+00
      do iperm = 1,nperm
         if(bgws(1,iperm)) then
            if(.not.GetFFfab(ws(1,iperm),FFfab,.true.)) then
               call abrtx("Maked2BetaDx2NI2: not GetFFfab")
            end if
         end if
         if(bgws(2,iperm).or.bgws(3,iperm)) then
            if(.not.GetGUXDEfg(ws(2,iperm),ws(3,iperm),
     &                         dummy,.false.,dummy,.false.,
     &                         dummy,.false.,Dgh,.true.,dummy,.false.,
     &                         norb,nfld)) then
               call abrtx("Maked2BetaDx2NI2: not GetGUXDEfg")
            end if
         end if
C
         ibterm=ieterm+1
         ieterm=nterms(iperm)
         do iterm = ibterm,ieterm
            ifld = iflds(1,iterm)
            jfld = iflds(2,iterm)
            kfld = iflds(3,iterm)
            index = iflds(4,iterm)
            mult = mults(iterm)
            do ijgeom=1,ngeom2
               dummy=dtrace('N','N',FFfab(1,ifld,ijgeom),norb,
     &                      Dgh(1,jfld,kfld),norb,norb,norb)
               d2bdx2(index,ijgeom)=d2bdx2(index,ijgeom)-mult*dummy
            end do
         end do
      end do
      end
C
C*MODULE TDXNI   *DECK Maked2BetaDxNI3
      subroutine Maked2BetaDx2NI3(d2bdx2,w,Ua,Xab,Gfgh,Efgh,
     &                            norb,nocc,ibfld,nfld,ngeom,ngeom2)
      implicit none
C     Arguments
      integer norb,nocc,ibfld,nfld,ngeom,ngeom2
      double precision w(3)
      double precision d2BdX2(27,ngeom2)
      double precision Ua(norb*norb,ngeom)
      double precision Xab(norb*norb,ngeom2)
      double precision Gfgh(norb*norb,nfld,nfld,nfld)
      double precision Efgh(norb*norb,nfld,nfld,nfld)
C     Used common blocks
      integer NELEMB,IELEMB,INDEXB
      common/NLO_BETA/NELEMB,IELEMB(3,27),INDEXB(3,3,3)
      DOUBLE PRECISION XX
      COMMON /FMCOM / XX(1)
C     Variables and constants
      integer lWork1
      double precision dummy
      integer igeom,jgeom,ijgeom,ifld,jfld,kfld,index
      logical bsame
      double precision ONE
      parameter (ONE=1.0d+00)
C     Called subprograms
      logical GetXab,GetGUXDEfgh,GetGUIDEa
      double precision dtrace
C     Body of the subprogram
      dummy = 0.0d+00
      if(.not.GetGUIDEa(dummy,.false.,Ua,.true.,dummy,.false.,
     &                  dummy,.false.,dummy,.false.,norb,ngeom)) then
         call abrtx("Maked2BetaDx2NI3: not GetGUIDEa")
      end if
      if(.not.GetXab(Xab,.true.)) then
         call abrtx("Maked2BetaDx2NI3: not GetXab")
      end if
      if(.not.GetGUXDEfgh(w(1),w(2),w(3),Gfgh,.true.,dummy,.false.,
     &     dummy,.false.,dummy,.false.,Efgh,.true.,norb,nfld)) then
         call abrtx("Maked2BetaDx2NI3: not GetGUXDEfgh")
      end if
C
      call NEWMEM(lWork1,norb*norb)
C
C     -------------------------------------
C     UabGfgh(wf,wg,wh) + XabEfgh(wf,wg,wh)
C     -------------------------------------
      ijgeom = 0
      do jgeom = 1,ngeom
         do igeom = 1,jgeom
            ijgeom=ijgeom+1
C
            bsame = igeom.eq.jgeom
            call CompNIUabSum(XX(lWork1),Ua(1,igeom),Ua(1,jgeom),bsame,
     &                        nocc,norb,ONE,.true.)
C
            do kfld=1,nfld
               do jfld=1,nfld
                  do ifld=1,nfld
                     dummy = dtrace('S','N',XX(lWork1),norb,
     &                            Gfgh(1,ifld,jfld,kfld),norb,norb,norb)
     &                      +dtrace('S','N',Xab(1,ijgeom),norb,
     &                            Efgh(1,ifld,jfld,kfld),norb,nocc,nocc)
                     dummy=dummy+dummy
                     index = INDEXB(ifld+ibfld-1,jfld+ibfld-1,
     &                              kfld+ibfld-1)
                     d2bdx2(index,ijgeom)=d2bdx2(index,ijgeom)-dummy
                  end do
               end do
            end do
         end do
      end do
C
      call DELETEMEM(lWork1,norb*norb)
C
      end
C
C*MODULE TDXNI   *DECK Maked2BetaDxNI4
      subroutine Maked2BetaDx2NI4(d2bdx2,w,Ua,Ia,Ug,Uh,Ugh,Udgh,
     &                            Uga,Udga,Uha,Udha,Gfa,Efa,
     &                            norb,nocc,ibfld,nfld,ngeom,ngeom2)
      implicit none
C     Arguments
      integer norb,nocc,ibfld,nfld,ngeom,ngeom2
      double precision w(3)
      double precision d2BdX2(27,ngeom2)
      double precision Ua(norb*norb,ngeom)
      double precision Ia(norb*norb,ngeom)
      double precision Ug(norb*norb,nfld)
      double precision Uh(norb*norb,nfld)
      double precision Ugh(norb*norb,nfld,nfld)
      double precision Udgh(norb*norb,nfld,nfld)
      double precision Uga(norb*norb,nfld,ngeom)
      double precision Udga(norb*norb,nfld,ngeom)
      double precision Uha(norb*norb,nfld,ngeom)
      double precision Udha(norb*norb,nfld,ngeom)
      double precision Gfa(norb*norb,nfld,ngeom)
      double precision Efa(norb*norb,nfld,ngeom)
C     Used common blocks
      DOUBLE PRECISION XX
      COMMON /FMCOM / XX(1)
C     Variables and constants
      integer lWork1,lWork2
      double precision dummy,mult
      integer ijgeom,igeom,jgeom
      integer iterm,ibterm,ieterm,ifld,jfld,kfld,index
      integer iperm,nperm
      parameter (nperm=3)
      integer nterms(nperm)
      integer iflds(4,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)
C     Called subprograms
      double precision Dtrace
      logical GetGUIDEa,GetGUDEf,GetGUXDEfa,GetGUXDEfg
C     Body of the subprogram
      call Maked2BetadX2NIPerms(w,ibfld,nfld,nterms,iflds,mults,
     &                          ws,bgws)
C
      dummy = 0.0d+00
      if(.not.GetGUIDEa(dummy,.false.,Ua,.true.,Ia,.true.,
     &                  dummy,.false.,dummy,.false.,norb,ngeom)) then
         call abrtx("Maked2BetaDx2NI4: not GetGUIDEa")
      end if
C
C     ------------------------------------------------------
C     Ughb(w1,w2)Gfa(w0) + Xghb(w1,w2)Efa(w0) + permutations
C     ------------------------------------------------------
C
      call NEWMEM(lWork1,norb*norb)
      call NEWMEM(lWork2,nocc*nocc)
C
      ieterm = 0
      do iperm = 1,nperm
         if(bgws(1,iperm)) then
            if(.not.GetGUXDEfa(ws(1,iperm),Gfa,.true.,dummy,.false.,
     &                         dummy,.false.,dummy,.false.,Efa,.true.,
     &                         norb,nfld,ngeom)) then
               call abrtx("Maked2BetaDx2NI4: not GetGUXDEfa")
            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
               call abrtx("Maked2BetaDx2NI4: not GetGUDEf:1")
            end if
            if(.not.GetGUXDEfa(ws(2,iperm),dummy,.false.,Uga,.true.,
     &                         Udga,.true.,dummy,.false.,dummy,.false.,
     &                         norb,nfld,ngeom)) then
               call abrtx("Maked2BetaDx2NI4: not GetGUXDEfa:1")
            end if
            call DAXPY(norb*norb*nfld*ngeom,MINUSONE,Uga,1,Udga,1)
         end if
         if(bgws(3,iperm)) then
            if(.not.GetGUDEf(ws(3,iperm),dummy,.false.,Uh,.true.,
     &                       dummy,.false.,dummy,.false.,
     &                       norb,nfld)) then
               call abrtx("Maked2BetaDx2NI4: not GetGUDEf:2")
            end if
            if(.not.GetGUXDEfa(ws(3,iperm),dummy,.false.,Uha,.true.,
     &                         Udha,.true.,dummy,.false.,dummy,.false.,
     &                         norb,nfld,ngeom)) then
               call abrtx("Maked2BetaDx2NI4: not GetGUXDEfa:2")
            end if
            call DAXPY(norb*norb*nfld*ngeom,MINUSONE,Uha,1,Udha,1)
         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("Maked2BetaDx2NI4: not GetGUXDEfg")
            end if
            call DAXPY(norb*norb*nfld*nfld,MINUSONE,Ugh,1,Udgh,1)
         end if
C
         ibterm=ieterm+1
         ieterm=nterms(iperm)
         do iterm = ibterm,ieterm
            ifld = iflds(1,iterm)
            jfld = iflds(2,iterm)
            kfld = iflds(3,iterm)
            index = iflds(4,iterm)
            mult = mults(iterm)
C
            do jgeom=1,ngeom
               call CompNIUfgaSum1(XX(lWork1),
     &                             Ugh(1,jfld,kfld),Udgh(1,jfld,kfld),
     &                             Ua(1,jgeom),nocc,norb,ONE,.true.)
               call CompNIUfgaSum2(XX(lWork1),
     &                             Ug(1,jfld),'f',Uha(1,kfld,jgeom),
     &                             Udha(1,kfld,jgeom),
     &                             nocc,norb,ONE,.false.)
               call CompNIUfgaSum2(XX(lWork1),
     &                             Uh(1,kfld),'f',Uga(1,jfld,jgeom),
     &                             Udga(1,jfld,jgeom),
     &                             nocc,norb,ONE,.false.)
C
               call CompNIXfgaSum1(XX(lWork2),
     &                             Ugh(1,jfld,kfld),Udgh(1,jfld,kfld),
     &                             Ua(1,jgeom),nocc,norb,ONE,.true.)
               call CompNIXfgaSum2(XX(lWork2),
     &                             Ug(1,jfld),'f',Uh(1,kfld),'f',
     &                             .false.,Ia(1,jgeom),
     &                             nocc,norb,ONE,.false.)
               call CompNIXfgaSum3(XX(lWork2),
     &                             Ug(1,jfld),'f',Uha(1,kfld,jgeom),
     &                             Udha(1,kfld,jgeom),
     &                             nocc,norb,ONE,.false.)
               call CompNIXfgaSum3(XX(lWork2),
     &                             Uh(1,kfld),'f',Uga(1,jfld,jgeom),
     &                             Udga(1,jfld,jgeom),
     &                             nocc,norb,ONE,.false.)
C
               do igeom=1,ngeom
                  dummy=dtrace('N','N',XX(lWork1),norb,
     &                         Gfa(1,ifld,igeom),norb,norb,norb)
     &                 +dtrace('N','N',XX(lWork2),nocc,
     &                         Efa(1,ifld,igeom),norb,nocc,nocc)
C
                  if(igeom.eq.jgeom) then
                     ijgeom = (igeom*igeom+igeom)/2
                     dummy=dummy+dummy
                  else if(igeom.gt.jgeom) then
                     ijgeom = (igeom*igeom-igeom)/2+jgeom
                  else
                     ijgeom = (jgeom*jgeom-jgeom)/2+igeom
                  end if
C
                  dummy=dummy+dummy
                  d2bdx2(index,ijgeom)=d2bdx2(index,ijgeom)-mult*dummy
               end do
            end do
         end do
      end do
C
      call DELETEMEM(lWork2,nocc*nocc)
      call DELETEMEM(lWork1,norb*norb)
C
      end
C
C*MODULE TDXNI   *DECK Maked2BetaDxNI5
      subroutine Maked2BetaDx2NI5(d2bdx2,w,Ua,Jab,Uf,Ufa,Udfa,Ggh,Egh,
     &                            norb,nocc,ibfld,nfld,ngeom,ngeom2)
      implicit none
C     Arguments
      integer norb,nocc,ibfld,nfld,ngeom,ngeom2
      double precision w(3)
      double precision d2BdX2(27,ngeom2)
      double precision Ua(norb*norb,ngeom)
      double precision Jab(norb*norb,ngeom2)
      double precision Uf(norb*norb,nfld)
      double precision Ufa(norb*norb,nfld,ngeom)
      double precision Udfa(norb*norb,nfld,ngeom)
      double precision Ggh(norb*norb,nfld,nfld)
      double precision Egh(norb*norb,nfld,nfld)
C     Used common blocks
      DOUBLE PRECISION XX
      COMMON /FMCOM / XX(1)
C     Variables and constants
      integer lWork1,lWork2
      double precision dummy,mult
      integer ijgeom,igeom,jgeom
      integer iterm,ibterm,ieterm,ifld,jfld,kfld,index
      logical bSame
      integer iperm,nperm
      parameter (nperm=3)
      integer nterms(nperm)
      integer iflds(4,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)
C     Called subprograms
      double precision Dtrace
      logical GetGUIDEa,GetJab,GetGUDEf,GetGUXDEfa,GetGUXDEfg
C     Body of the subprogram
      call Maked2BetadX2NIPerms(w,ibfld,nfld,nterms,iflds,mults,
     &                          ws,bgws)
C
      dummy = 0.0d+00
      if(.not.GetGUIDEa(dummy,.false.,Ua,.true.,dummy,.false.,
     &                  dummy,.false.,dummy,.false.,norb,ngeom)) then
         call abrtx("Maked2BetaDx2NI5: not GetGUIDEa")
      end if
      if(.not.GetJab(Jab,.true.)) then
         call abrtx("Maked2BetaDx2NI5: not GetJab")
      end if
C
C     ------------------------------------------------------
C     Ufab(w0)Ggh(w1,w2) + Xfab(w0)Egh(w1,w2) + permutations
C     ------------------------------------------------------
C
      call NEWMEM(lWork1,norb*norb)
      call NEWMEM(lWork2,nocc*nocc)
C
      ieterm = 0
      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("Maked2BetaDx2NI5: not GetGUDEf")
            end if
            if(.not.GetGUXDEfa(ws(1,iperm),dummy,.false.,Ufa,.true.,
     &                         Udfa,.true.,dummy,.false.,dummy,.false.,
     &                         norb,nfld,ngeom)) then
               call abrtx("Maked2BetaDx2NI5: not GetGUXDEfa")
            end if
            call DAXPY(norb*norb*nfld*ngeom,MINUSONE,Ufa,1,Udfa,1)
         end if
         if(bgws(2,iperm).or.bgws(3,iperm)) then
            if(.not.GetGUXDEfg(ws(2,iperm),ws(3,iperm),
     &                         Ggh,.true.,dummy,.false.,dummy,.false.,
     &                         dummy,.false.,Egh,.true.,
     &                         norb,nfld)) then
               call abrtx("Maked2BetaDx2NI5: not GetGUXDEfg")
            end if
         end if
C
         ibterm=ieterm+1
         ieterm=nterms(iperm)
         do iterm = ibterm,ieterm
            ifld = iflds(1,iterm)
            jfld = iflds(2,iterm)
            kfld = iflds(3,iterm)
            index = iflds(4,iterm)
            mult = mults(iterm)
C
            ijgeom = 0
C
            do jgeom = 1,ngeom
               do igeom = 1,jgeom
                  ijgeom=ijgeom+1
                  bSame=igeom.eq.jgeom
C
                  call CompNIUfabSum(XX(lWork1),Ufa(1,ifld,igeom),
     &                               Udfa(1,ifld,igeom),Ua(1,jgeom),
     &                               Ufa(1,ifld,jgeom),
     &                               Udfa(1,ifld,jgeom),Ua(1,igeom),
     &                               bSame,nocc,norb,ONE,.true.)
                  call CompNIXfabSum1(XX(lWork2),Ufa(1,ifld,igeom),
     &                                Udfa(1,ifld,igeom),Ua(1,jgeom),
     &                                Ufa(1,ifld,jgeom),
     &                                Udfa(1,ifld,jgeom),Ua(1,igeom),
     &                                bSame,nocc,norb,ONE,.true.)
                  call CompNIXfabSum2(XX(lWork2),Jab(1,ijgeom),
     &                                Uf(1,ifld),'f',
     &                                nocc,norb,ONE,.false.)
C
                  dummy=dtrace('N','N',XX(lWork1),norb,
     &                         Ggh(1,jfld,kfld),norb,norb,norb)
     &                 +dtrace('N','N',XX(lWork2),nocc,
     &                         Egh(1,jfld,kfld),norb,nocc,nocc)
C
                  dummy=dummy+dummy
                  d2bdx2(index,ijgeom)=d2bdx2(index,ijgeom)-mult*dummy
               end do
            end do
         end do
      end do
C
      call DELETEMEM(lWork2,nocc*nocc)
      call DELETEMEM(lWork1,norb*norb)
C
      end
C
C*MODULE TDXNI   *DECK Maked2BetaDxNI6a
      subroutine Maked2BetaDx2NI6a(d2bdx2,w,Ua,Ufgh,Udfgh,Ga,Ea,
     &                             norb,nocc,ibfld,nfld,ngeom,ngeom2)
      implicit none
C     Arguments
      integer norb,nocc,ibfld,nfld,ngeom,ngeom2
      double precision w(3)
      double precision d2BdX2(27,ngeom2)
      double precision Ua(norb*norb,ngeom)
      double precision Ufgh(norb*norb,nfld,nfld,nfld)
      double precision Udfgh(norb*norb,nfld,nfld,nfld)
      double precision Ga(norb*norb,ngeom)
      double precision Ea(norb*norb,ngeom)
C     Used common blocks
      integer NELEMB,IELEMB,INDEXB
      common/NLO_BETA/NELEMB,IELEMB(3,27),INDEXB(3,3,3)
      DOUBLE PRECISION XX
      COMMON /FMCOM / XX(1)
C     Variables and constants
      integer lWork1,lWork2
      double precision dummy
      integer igeom,jgeom,ijgeom,ifld,jfld,kfld,index
      double precision ONE,MINUSONE
      parameter (ONE=1.0d+00,MINUSONE=-1.0d+00)
C     Called subprograms
      logical GetGUIDEa,GetGUXDEfgh
      double precision dtrace
C     Body of the subprogram
      dummy = 0.0d+00
      if(.not.GetGUIDEa(Ga,.true.,Ua,.true.,dummy,.false.,
     &                  dummy,.false.,Ea,.true.,norb,ngeom)) then
         call abrtx("Maked2BetaDx2NI6a: not GetGUIDEa")
      end if
      if(.not.GetGUXDEfgh(w(1),w(2),w(3),dummy,.false.,Ufgh,.true.,
     &     Udfgh,.true.,dummy,.false.,dummy,.false.,norb,nfld)) then
         call abrtx("Maked2BetaDx2NI6a: not GetGUXDEfgh")
      end if
      call DAXPY(norb*norb*nfld*nfld*nfld,MINUSONE,Ufgh,1,Udfgh,1)
C
      call NEWMEM(lWork1,norb*norb)
      call NEWMEM(lWork2,nocc*nocc)
C
C     ----------------------------------------------------------
C     Ufghb(w0,w1,w2)Ga + Xfghb(w0,w1,w2)Ea + (a,b) permutations
C     ----------------------------------------------------------
      do kfld=1,nfld
         do jfld=1,nfld
            do ifld=1,nfld
               do jgeom = 1,ngeom
C
                  call CompNIUfghaSum1(XX(lWork1),Ua(1,jgeom),
     &                                 Ufgh(1,ifld,jfld,kfld),
     &                                 Udfgh(1,ifld,jfld,kfld),
     &                                 nocc,norb,ONE,.true.)
                  call CompNIXfghaSum1(XX(lWork2),Ua(1,jgeom),
     &                                 Ufgh(1,ifld,jfld,kfld),
     &                                 Udfgh(1,ifld,jfld,kfld),
     &                                 nocc,norb,ONE,.true.)
C
                  do igeom = 1,ngeom
C
                     dummy = dtrace('N','S',XX(lWork1),norb,
     &                              Ga(1,igeom),norb,norb,norb)
     &                      +dtrace('N','S',XX(lWork2),nocc,
     &                              Ea(1,igeom),norb,nocc,nocc)
C
                     if(igeom.eq.jgeom) then
                        ijgeom = (igeom*igeom+igeom)/2
                        dummy=dummy+dummy
                     else if(igeom.gt.jgeom) then
                        ijgeom = (igeom*igeom-igeom)/2+jgeom
                     else
                        ijgeom = (jgeom*jgeom-jgeom)/2+igeom
                     end if
C
                     dummy=dummy+dummy
                     index=INDEXB(ifld+ibfld-1,jfld+ibfld-1,
     &                            kfld+ibfld-1)
                     d2bdx2(index,ijgeom)=d2bdx2(index,ijgeom)-dummy
                  end do
               end do
            end do
         end do
      end do
C
      call DELETEMEM(lWork2,nocc*nocc)
      call DELETEMEM(lWork1,norb*norb)
C
      end
C
C*MODULE TDXNI   *DECK Maked2BetaDxNI6b
      subroutine Maked2BetaDx2NI6b(d2bdx2,w,Ia,Uf,Ufa,Udfa,Ugh,Udgh,
     &                             Ga,Ea,
     &                             norb,nocc,
     &                             ibfld,nfld,ngeom,ngeom2)
      implicit none
C     Arguments
      integer norb,nocc,ibfld,nfld,ngeom,ngeom2
      double precision w(3)
      double precision d2BdX2(27,ngeom2)
      double precision Ia(norb*norb,ngeom)
      double precision Uf(norb*norb,nfld)
      double precision Ufa(norb*norb,nfld,ngeom)
      double precision Udfa(norb*norb,nfld,ngeom)
      double precision Ugh(norb*norb,nfld,nfld)
      double precision Udgh(norb*norb,nfld,nfld)
      double precision Ga(norb*norb,ngeom)
      double precision Ea(norb*norb,ngeom)
C     Used common blocks
      DOUBLE PRECISION XX
      COMMON /FMCOM / XX(1)
C     Variables and constants
      integer lWork1,lWork2
      double precision dummy,mult
      integer ijgeom,igeom,jgeom
      integer iterm,ibterm,ieterm,ifld,jfld,kfld,index
      integer iperm,nperm
      parameter (nperm=3)
      integer nterms(nperm)
      integer iflds(4,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)
C     Called subprograms
      double precision Dtrace
      logical GetGUIDEa,GetGUDEf,GetGUXDEfa,GetGUXDEfg
C     Body of the subprogram
      call Maked2BetadX2NIPerms(w,ibfld,nfld,nterms,iflds,mults,
     &                          ws,bgws)
C
      dummy = 0.0d+00
      if(.not.GetGUIDEa(Ga,.true.,dummy,.false.,Ia,.true.,
     &                  dummy,.false.,Ea,.true.,norb,ngeom)) then
         call abrtx("Maked2BetaDx2NI6b: not GetGUIDEa")
      end if
C
C     ----------------------------------------------------------
C     Ufghb(w0,w1,w2)Ga + Xfghb(w0,w1,w2)Ea + (a,b) permutations
C     ----------------------------------------------------------
C
      call NEWMEM(lWork1,norb*norb)
      call NEWMEM(lWork2,nocc*nocc)
C
      ieterm = 0
      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("Maked2BetaDx2NI6b: not GetGUDEf")
            end if
            if(.not.GetGUXDEfa(ws(1,iperm),dummy,.false.,Ufa,.true.,
     &                         Udfa,.true.,dummy,.false.,dummy,.false.,
     &                         norb,nfld,ngeom)) then
               call abrtx("Maked2BetaDx2NI6b: not GetGUXDEfa")
            end if
            call DAXPY(norb*norb*nfld*ngeom,MINUSONE,Ufa,1,Udfa,1)
         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("Maked2BetaDx2NI6b: not GetGUXDEfg")
            end if
            call DAXPY(norb*norb*nfld*nfld,MINUSONE,Ugh,1,Udgh,1)
         end if
C
         ibterm=ieterm+1
         ieterm=nterms(iperm)
         do iterm = ibterm,ieterm
            ifld = iflds(1,iterm)
            jfld = iflds(2,iterm)
            kfld = iflds(3,iterm)
            index = iflds(4,iterm)
            mult = mults(iterm)
C
            do jgeom = 1,ngeom
               call CompNIUfghaSum2(XX(lWork1),Ufa(1,ifld,jgeom),
     &                              Udfa(1,ifld,jgeom),
     &                              Ugh(1,jfld,kfld),Udgh(1,jfld,kfld),
     &                              nocc,norb,ONE,.true.)
               call CompNIXfghaSum2(XX(lWork2),Ia(1,jgeom),Uf(1,ifld),
     &                              Ufa(1,ifld,jgeom),
     &                              Udfa(1,ifld,jgeom),
     &                              Ugh(1,jfld,kfld),Udgh(1,jfld,kfld),
     &                              nocc,norb,ONE,.true.)
C
               do igeom = 1,ngeom
                  dummy = dtrace('N','S',XX(lWork1),norb,
     &                           Ga(1,igeom),norb,norb,norb)
     &                   +dtrace('N','S',XX(lWork2),nocc,
     &                           Ea(1,igeom),norb,nocc,nocc)
C
                  if(igeom.eq.jgeom) then
                     ijgeom = (igeom*igeom+igeom)/2
                     dummy=dummy+dummy
                  else if(igeom.gt.jgeom) then
                     ijgeom = (igeom*igeom-igeom)/2+jgeom
                  else
                     ijgeom = (jgeom*jgeom-jgeom)/2+igeom
                  end if
C
                  dummy=dummy+dummy
                  d2bdx2(index,ijgeom)=d2bdx2(index,ijgeom)-dummy*mult
               end do
            end do
         end do
      end do
C
      call DELETEMEM(lWork2,nocc*nocc)
      call DELETEMEM(lWork1,norb*norb)
C
      end
C
C*MODULE TDXNI   *DECK Maked2BetaDxNI7a
      subroutine Maked2BetaDx2NI7a(d2bdx2,w,Ia,Ug,Uh,Uga,Udga,Uha,Udha,
     &                             Gf,Ef,
     &                             norb,nocc,ibfld,nfld,ngeom,ngeom2)
      implicit none
C     Arguments
      integer norb,nocc,ibfld,nfld,ngeom,ngeom2
      double precision w(3)
      double precision d2BdX2(27,ngeom2)
      double precision Ia(norb*norb,ngeom)
      double precision Ug(norb*norb,nfld)
      double precision Uh(norb*norb,nfld)
      double precision Uga(norb*norb,nfld,ngeom)
      double precision Udga(norb*norb,nfld,ngeom)
      double precision Uha(norb*norb,nfld,ngeom)
      double precision Udha(norb*norb,nfld,ngeom)
      double precision Gf(norb*norb,nfld)
      double precision Ef(norb*norb,nfld)
C     Used common blocks
      DOUBLE PRECISION XX
      COMMON /FMCOM / XX(1)
C     Variables and constants
      integer lWork1,lWork2
      double precision dummy,mult
      integer ijgeom,igeom,jgeom
      integer iterm,ibterm,ieterm,ifld,jfld,kfld,index
      integer iperm,nperm
      parameter (nperm=3)
      integer nterms(nperm)
      integer iflds(4,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)
C     Called subprograms
      double precision Dtrace
      logical GetGUIDEa,GetGUDEf,GetGUXDEfa
C     Body of the subprogram
      call Maked2BetadX2NIPerms(w,ibfld,nfld,nterms,iflds,mults,
     &                          ws,bgws)
C
      dummy = 0.0d+00
      if(.not.GetGUIDEa(dummy,.false.,dummy,.false.,Ia,.true.,
     &                  dummy,.false.,dummy,.false.,norb,ngeom)) then
         call abrtx("Maked2BetaDx2NI7a: not GetGUIDEa")
      end if
C
C     ------------------------------------------------------
C     Ughab(w1,w2)Gf(w0) + Xghab(w1,w2)Ef(w0) + permutations
C     ------------------------------------------------------
C
      call NEWMEM(lWork1,norb*norb)
      call NEWMEM(lWork2,nocc*nocc)
C
      ieterm = 0
      do iperm = 1,nperm
         if(bgws(1,iperm)) then
            if(.not.GetGUDEf(ws(1,iperm),Gf,.true.,dummy,.false.,
     &                       dummy,.false.,Ef,.true.,
     &                       norb,nfld)) then
               call abrtx("Maked2BetaDx2NI7a: 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
               call abrtx("Maked2BetaDx2NI7a: not GetGUDEf:2")
            end if
            if(.not.GetGUXDEfa(ws(2,iperm),dummy,.false.,Uga,.true.,
     &                         Udga,.true.,dummy,.false.,dummy,.false.,
     &                         norb,nfld,ngeom)) then
               call abrtx("Maked2BetaDx2NI7a: not GetGUXDEfa:2")
            end if
            call DAXPY(norb*norb*nfld*ngeom,MINUSONE,Uga,1,Udga,1)
         end if
         if(bgws(3,iperm)) then
            if(.not.GetGUDEf(ws(3,iperm),dummy,.false.,Uh,.true.,
     &                       dummy,.false.,dummy,.false.,
     &                       norb,nfld)) then
               call abrtx("Maked2BetaDx2NI7a: not GetGUDEf:3")
            end if
            if(.not.GetGUXDEfa(ws(3,iperm),dummy,.false.,Uha,.true.,
     &                         Udha,.true.,dummy,.false.,dummy,.false.,
     &                         norb,nfld,ngeom)) then
               call abrtx("Maked2BetaDx2NI7a: not GetGUXDEfa:3")
            end if
            call DAXPY(norb*norb*nfld*ngeom,MINUSONE,Uha,1,Udha,1)
         end if
C
         ibterm=ieterm+1
         ieterm=nterms(iperm)
         do iterm = ibterm,ieterm
            ifld = iflds(1,iterm)
            jfld = iflds(2,iterm)
            kfld = iflds(3,iterm)
            index = iflds(4,iterm)
            mult = mults(iterm)
C
            do jgeom = 1,ngeom
               do igeom = 1,ngeom
C
                  call CompNIUfgabSum(XX(lWork1),
     &                                Uga(1,jfld,igeom),
     &                                Udga(1,jfld,igeom),
     &                                Uha(1,kfld,jgeom),
     &                                Udha(1,kfld,jgeom),.false.,
     &                                nocc,norb,ONE,.true.)
                  call CompNIXfgabSum1(XX(lWork2),
     &                                 Ia(1,igeom),Ug(1,jfld),
     &                                 Uga(1,jfld,igeom),
     &                                 Udga(1,jfld,igeom),
     &                                 Ia(1,jgeom),Uh(1,kfld),
     &                                 Uha(1,kfld,jgeom),
     &                                 Udha(1,kfld,jgeom),.false.,
     &                                 nocc,norb,ONE,.true.)
C
                  dummy = dtrace('N','N',XX(lWork1),norb,
     &                           Gf(1,ifld),norb,norb,norb)
     &                   +dtrace('N','N',XX(lWork2),nocc,
     &                           Ef(1,ifld),norb,nocc,nocc)
C
                  if(igeom.eq.jgeom) then
                     ijgeom = (igeom*igeom+igeom)/2
                     dummy=dummy+dummy
                  else if(igeom.gt.jgeom) then
                     ijgeom = (igeom*igeom-igeom)/2+jgeom
                  else
                     ijgeom = (jgeom*jgeom-jgeom)/2+igeom
                  end if
C
                  dummy=dummy+dummy
                  d2bdx2(index,ijgeom)=d2bdx2(index,ijgeom)-dummy*mult
               end do
            end do
         end do
      end do
C
      call DELETEMEM(lWork2,nocc*nocc)
      call DELETEMEM(lWork1,norb*norb)
C
      end
C
C*MODULE TDXNI   *DECK Maked2BetaDxNI7b
      subroutine Maked2BetaDx2NI7b(d2bdx2,w,Jab,Ugh,Udgh,Gf,Ef,
     &                             norb,nocc,ibfld,nfld,ngeom2)
      implicit none
C     Arguments
      integer norb,nocc,ibfld,nfld,ngeom2
      double precision w(3)
      double precision d2BdX2(27,ngeom2)
      double precision Jab(norb*norb,ngeom2)
      double precision Ugh(norb*norb,nfld,nfld)
      double precision Udgh(norb*norb,nfld,nfld)
      double precision Gf(norb*norb,nfld)
      double precision Ef(norb*norb,nfld)
C     Used common blocks
      DOUBLE PRECISION XX
      COMMON /FMCOM / XX(1)
C     Variables and constants
      integer lWork2
      double precision dummy,mult
      integer ijgeom
      integer iterm,ibterm,ieterm,ifld,jfld,kfld,index
      integer iperm,nperm
      parameter (nperm=3)
      integer nterms(nperm)
      integer iflds(4,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)
C     Called subprograms
      double precision Dtrace
      logical GetJab,GetGUDEf,GetGUXDEfg
C     Body of the subprogram
      call Maked2BetadX2NIPerms(w,ibfld,nfld,nterms,iflds,mults,
     &                          ws,bgws)
C
      if(.not.GetJab(Jab,.true.)) then
         call abrtx("Maked2BetaDx2NI7b: not GetJab")
      end if
C
C     ------------------------------------------------------
C     Ughab(w1,w2)Gf(w0) + Xghab(w1,w2)Ef(w0) + permutations
C     ------------------------------------------------------
C
      call NEWMEM(lWork2,nocc*nocc)
C
      ieterm = 0
      dummy = 0.0d+00
      do iperm = 1,nperm
         if(bgws(1,iperm)) then
            if(.not.GetGUDEf(ws(1,iperm),Gf,.true.,dummy,.false.,
     &                       dummy,.false.,Ef,.true.,
     &                       norb,nfld)) then
               call abrtx("Maked2BetaDx2NI7b: not GetGUDEf")
            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("Maked2BetaDx2NI7b: not GetGUXDEfg")
            end if
            call DAXPY(norb*norb*nfld*nfld,MINUSONE,Ugh,1,Udgh,1)
         end if
C
         ibterm=ieterm+1
         ieterm=nterms(iperm)
         do iterm = ibterm,ieterm
            ifld = iflds(1,iterm)
            jfld = iflds(2,iterm)
            kfld = iflds(3,iterm)
            index = iflds(4,iterm)
            mult = mults(iterm)
C
            do ijgeom = 1,ngeom2
               call CompNIXfgabSum2(XX(lWork2),Jab(1,ijgeom),
     &                              Ugh(1,jfld,kfld),Udgh(1,jfld,kfld),
     &                              nocc,norb,ONE,.true.)
C
               dummy =dtrace('N','N',XX(lWork2),nocc,
     &                       Ef(1,ifld),norb,nocc,nocc)
C
               dummy=dummy+dummy
               d2bdx2(index,ijgeom)=d2bdx2(index,ijgeom)-dummy*mult
            end do
         end do
      end do
C
      call DELETEMEM(lWork2,nocc*nocc)
C
      end
C
C*MODULE TDXNI   *DECK Maked2BetaDxNI7c
      subroutine Maked2BetaDx2NI7c(d2bdx2,w,Iab,Ug,Uh,Gf,Ef,
     &                             norb,nocc,ibfld,nfld,ngeom2)
      implicit none
C     Arguments
      integer norb,nocc,ibfld,nfld,ngeom2
      double precision w(3)
      double precision d2BdX2(27,ngeom2)
      double precision Iab(norb*norb,ngeom2)
      double precision Ug(norb*norb,nfld)
      double precision Uh(norb*norb,nfld)
      double precision Gf(norb*norb,nfld)
      double precision Ef(norb*norb,nfld)
C     Used common blocks
      DOUBLE PRECISION XX
      COMMON /FMCOM / XX(1)
C     Variables and constants
      integer lWork2
      double precision dummy,mult
      integer ijgeom
      integer iterm,ibterm,ieterm,ifld,jfld,kfld,index
      integer iperm,nperm
      parameter (nperm=3)
      integer nterms(nperm)
      integer iflds(4,nperm*27)
      double precision mults(nperm*27)
      double precision ws(3,nperm)
      logical bgws(3,nperm)
      double precision ONE
      parameter (ONE=1.0d+00)
C     Called subprograms
      double precision Dtrace
      logical GetIab,GetGUDEf
C     Body of the subprogram
      call Maked2BetadX2NIPerms(w,ibfld,nfld,nterms,iflds,mults,
     &                          ws,bgws)
C
      if(.not.GetIab(Iab,.true.)) then
         call abrtx("Maked2BetaDx2NI7c: not GetIab")
      end if
C
C     ------------------------------------------------------
C     Ughab(w1,w2)Gf(w0) + Xghab(w1,w2)Ef(w0) + permutations
C     ------------------------------------------------------
C
      call NEWMEM(lWork2,nocc*nocc)
C
      ieterm = 0
      dummy = 0.0d+00
      do iperm = 1,nperm
         if(bgws(1,iperm)) then
            if(.not.GetGUDEf(ws(1,iperm),Gf,.true.,dummy,.false.,
     &                       dummy,.false.,Ef,.true.,
     &                       norb,nfld)) then
               call abrtx("Maked2BetaDx2NI7c: 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
               call abrtx("Maked2BetaDx2NI7c: not GetGUDEf:2")
            end if
         end if
         if(bgws(3,iperm)) then
            if(.not.GetGUDEf(ws(3,iperm),dummy,.false.,Uh,.true.,
     &                       dummy,.false.,dummy,.false.,
     &                       norb,nfld)) then
               call abrtx("Maked2BetaDx2NI7c: not GetGUDEf:3")
            end if
         end if
C
         ibterm=ieterm+1
         ieterm=nterms(iperm)
         do iterm = ibterm,ieterm
            ifld = iflds(1,iterm)
            jfld = iflds(2,iterm)
            kfld = iflds(3,iterm)
            index = iflds(4,iterm)
            mult = mults(iterm)
C
            do ijgeom = 1,ngeom2
               call CompNIXfgabSum3(XX(lWork2),Ug(1,jfld),'f',
     &                              Uh(1,kfld),'f',.false.,
     &                              Iab(1,ijgeom),
     &                              nocc,norb,ONE,.true.)
C
               dummy =dtrace('N','N',XX(lWork2),nocc,
     &                       Ef(1,ifld),norb,nocc,nocc)
C
               dummy=dummy+dummy
               d2bdx2(index,ijgeom)=d2bdx2(index,ijgeom)-dummy*mult
            end do
         end do
      end do
C
      call DELETEMEM(lWork2,nocc*nocc)
C
      end
C
C*MODULE TDXNI   *DECK Maked2BetadX2NI8a
      subroutine Maked2BetadX2NI8a(d2bdx2,w,E0,Jab,Ufgh,Udfgh,
     &                             norb,nocc,ibfld,nfld,ngeom2)
      implicit none
C     Arguments
      integer norb,nocc,ibfld,nfld,ngeom2
      double precision w(3)
      double precision d2bdx2(27,ngeom2)
      double precision E0(norb)
      double precision Jab(norb*norb,ngeom2)
      double precision Ufgh(norb*norb,nfld,nfld,nfld)
      double precision Udfgh(norb*norb,nfld,nfld,nfld)
C     Used common blocks
      integer NELEMB,IELEMB,INDEXB
      common/NLO_BETA/NELEMB,IELEMB(3,27),INDEXB(3,3,3)
      DOUBLE PRECISION XX
      COMMON /FMCOM / XX(1)
C     Variables and constants
      integer lVTR
      integer ijgeom,ifld,jfld,kfld,index
      double precision dummy
      double precision MINUSONE
      parameter (MINUSONE=-1.0d+00)
C     Called subprograms
      double precision DDOT
      logical GetGUXDEfgh,GetJab
C     Body of the subprogram
      call GetEnergy(E0,norb)
      dummy = 0.0d+00
      if(.not.GetGUXDEfgh(w(1),w(2),w(3),
     &                    dummy,.false.,Ufgh,.true.,Udfgh,.true.,
     &                    dummy,.false.,dummy,.false.,norb,nfld)) then
         call abrtx("Maked2BetadX2NI8a: not GetGUXDEfgh")
      end if
      call DAXPY(norb*norb*nfld*nfld*nfld,MINUSONE,Ufgh,1,Udfgh,1)
      if(.not.GetJab(Jab,.true.)) then
         call abrtx("Maked2BetadX2NI8a: not GetJab")
      end if
C
      call NEWMEM(lVTR,nocc)
C
C     -------------------------------
C     -[Udfgh(-w) Jab+Jdab Ufgh(w)]E0
C     -------------------------------
      do ifld = 1,nfld
         do jfld = 1,nfld
            do kfld = 1,nfld
               do ijgeom = 1,ngeom2
                  call VTRACE(XX(lVTR),'N','N',
     &                        Udfgh(1,ifld,jfld,kfld),norb,
     &                        Jab(1,ijgeom),norb,
     &                        nocc,norb)
                  dummy = -DDOT(nocc,XX(lVTR),1,E0,1)
                  call VTRACE(XX(lVTR),'T','N',
     &                        Jab(1,ijgeom),norb,
     &                        Ufgh(1,ifld,jfld,kfld),norb,
     &                        nocc,norb)
                  dummy = dummy-DDOT(nocc,XX(lVTR),1,E0,1)
                  dummy=dummy+dummy
C
                  index = INDEXB(ifld+ibfld-1,jfld+ibfld-1,kfld+ibfld-1)
                  d2bdx2(index,ijgeom)=d2bdx2(index,ijgeom)-dummy
               end do
            end do
         end do
      end do
C
      call DELETEMEM(lVTR,nocc)
      return
      end
C
C*MODULE TDXNI   *DECK Maked2BetadX2NI8b
      subroutine Maked2BetadX2NI8b(d2bdx2,w,E0,Ia,Iab,Uf,Ufa,Udfa,
     &                             Ugh,Udgh,
     &                             norb,nocc,
     &                             ibfld,nfld,ngeom,ngeom2)
      implicit none
C     Arguments
      integer norb,nocc,ibfld,nfld,ngeom,ngeom2
      double precision w(3)
      double precision d2BdX2(27,ngeom2)
      double precision E0(norb)
      double precision Ia(norb*norb,ngeom)
      double precision Iab(norb*norb,ngeom2)
      double precision Uf(norb*norb,nfld)
      double precision Ufa(norb*norb,nfld,ngeom)
      double precision Udfa(norb*norb,nfld,ngeom)
      double precision Ugh(norb*norb,nfld,nfld)
      double precision Udgh(norb*norb,nfld,nfld)
C     Used common blocks
      DOUBLE PRECISION XX
      COMMON /FMCOM / XX(1)
C     Variables and constants
      integer lTMP,lVTR
      double precision dummy,mult
      integer ijgeom,igeom,jgeom
      integer iterm,ibterm,ieterm,ifld,jfld,kfld,index
      integer iperm,nperm
      parameter (nperm=3)
      integer nterms(nperm)
      integer iflds(4,nperm*27)
      double precision mults(nperm*27)
      double precision ws(3,nperm)
      logical bgws(3,nperm)
      double precision ZERO,ONE,MINUSONE
      parameter (ZERO=0.0d+00,ONE=1.0d+00,MINUSONE=-1.0d+00)
      character*1 sN
      parameter (sN='N')
C     Called subprograms
      double precision DSUM,DDOT
      logical GetGUIDEa,GetIab,GetGUDEf,GetGUXDEfa,GetGUXDEfg
C     Body of the subprogram
      call Maked2BetadX2NIPerms(w,ibfld,nfld,nterms,iflds,mults,
     &                          ws,bgws)
C
      call GetEnergy(E0,norb)
      dummy = 0.0d+00
      if(.not.GetGUIDEa(dummy,.false.,dummy,.false.,Ia,.true.,
     &                  dummy,.false.,dummy,.false.,norb,ngeom)) then
         call abrtx("Maked2BetadX2NI8b: not GetGUIDEa")
      end if
      if(.not.GetIab(Iab,.true.)) then
         call abrtx("Maked2BetadX2NI8b: not GetIab")
      end if
C
      call NEWMEM(lTMP,norb*nocc)
      call NEWMEM(lVTR,nocc)
C
      ieterm = 0
      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("Maked2BetadX2NI8b: not GetGUDEf")
            end if
            if(.not.GetGUXDEfa(ws(1,iperm),dummy,.false.,Ufa,.true.,
     &                         Udfa,.true.,dummy,.false.,dummy,.false.,
     &                         norb,nfld,ngeom)) then
               call abrtx("Maked2BetadX2NI8b: not GetGUXDEfa")
            end if
            call DAXPY(norb*norb*nfld*ngeom,MINUSONE,Ufa,1,Udfa,1)
         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("Maked2BetadX2NI8b: not GetGUXDEfg")
            end if
            call DAXPY(norb*norb*nfld*nfld,MINUSONE,Ugh,1,Udgh,1)
         end if
C
         ibterm=ieterm+1
         ieterm=nterms(iperm)
         do iterm = ibterm,ieterm
            ifld = iflds(1,iterm)
            jfld = iflds(2,iterm)
            kfld = iflds(3,iterm)
            index = iflds(4,iterm)
            mult = mults(iterm)
C
            ijgeom = 0
            do jgeom = 1,ngeom
               do igeom = 1,jgeom
C
                  ijgeom=ijgeom+1
C     --------------------------------------
C     [IabUf(w) + IaUfb(w) + IbUfa(w)] = TMP
C     --------------------------------------
                  call DGEMM(sN,sN,norb,nocc,norb,
     &                       ONE,Iab(1,ijgeom),norb,
     &                       Uf(1,ifld),norb,
     &                       ZERO,XX(lTMP),norb)
                  call DGEMM(sN,sN,norb,nocc,norb,
     &                       ONE,Ia(1,igeom),norb,
     &                       Ufa(1,ifld,jgeom),norb,
     &                       ONE,XX(lTMP),norb)
                  call DGEMM(sN,sN,norb,nocc,norb,
     &                       ONE,Ia(1,jgeom),norb,
     &                       Ufa(1,ifld,igeom),norb,
     &                       ONE,XX(lTMP),norb)
C     ------------------
C     [Udgh(-w) TMP]=VTR
C     ------------------
                  call VTRACE(XX(lVTR),'N','N',
     &                        Udgh(1,jfld,kfld),norb,XX(lTMP),norb,
     &                        nocc,norb)
C     --------------------
C     -VTR E0 + w sum(VTR)
C     --------------------
                  dummy=-DDOT(nocc,XX(lVTR),1,E0,1)
     &                  +DSUM(nocc,XX(lVTR),1)*ws(1,iperm)
                  dummy=dummy+dummy
                  d2bdx2(index,ijgeom)=d2bdx2(index,ijgeom)-dummy*mult
C     ---------------------------------------------
C     [Udf(-w)Iab + Udfa(-w)Ib + Udfb(-w)Ia]
C     = [-Uf(w)Iab + Udfa(-w)Ib + Udfb(-w)Ia] = TMP
C     ---------------------------------------------
                  call DGEMM(sN,sN,nocc,norb,norb,
     &                       MINUSONE,Uf(1,ifld),norb,
     &                       Iab(1,ijgeom),norb,
     &                       ZERO,XX(lTMP),nocc)
                  call DGEMM(sN,sN,nocc,norb,norb,
     &                       ONE,Udfa(1,ifld,jgeom),norb,
     &                       Ia(1,igeom),norb,
     &                       ONE,XX(lTMP),nocc)
                  call DGEMM(sN,sN,nocc,norb,norb,
     &                       ONE,Udfa(1,ifld,igeom),norb,
     &                       Ia(1,jgeom),norb,
     &                       ONE,XX(lTMP),nocc)
C     ----------------
C     [TMP Ugh(w)]=VTR
C     ----------------
                  call VTRACE(XX(lVTR),'N','N',
     &                        XX(lTMP),nocc,Ugh(1,jfld,kfld),norb,
     &                        nocc,norb)
C     --------------------
C     -VTR E0 - w sum(VTR)
C     --------------------
                  dummy=-DDOT(nocc,XX(lVTR),1,E0,1)
     &                  -DSUM(nocc,XX(lVTR),1)*ws(1,iperm)
                  dummy=dummy+dummy
                  d2bdx2(index,ijgeom)=d2bdx2(index,ijgeom)-dummy*mult
               end do
            end do
         end do
      end do
C
      call DELETEMEM(lVTR,nocc)
      call DELETEMEM(lTMP,norb*nocc)
C
      end
