C 21 APR 08 - SHY - INITIALIZATION FOR TD-DFT/EFP1, U-TDDT ON LEBEDEV
C 11 APR 08 - HPTI - GRIDS AND PARALLELIZATION
C  4 MAR 08 - MWS - TDDFTCALC: FIX DAWRIT LENGTH CHANGING PROBLEM
C  4 MAR 08 - SHY - CHANGES FOR OPEN SHELL TDDFT, AND TO EFP1 SOLVATION
C 28 AUG 07 - SHY - TDDFTCALC: ADD CALL FOR EFP SOLVATION CORRECTIONS
C 20 AUG 07 - MC  - TDDFTCALC,TDPRNT,INFOTD: ADD ARGUMENTS FOR FMO-TDDFT
C 27 DEC 06 - DGF - TDPRNT: SYMMETRY PRINTING CHANGES FOR FMO
C 19 NOV 06 - MC  - RPARESVEC: FIXES FOR THE SOLVER
C  6 NOV 06 - MC  - NEW MODULE FOR TDDFT COMPUTATION
C
C*MODULE TDDFT  *DECK MEMGRD
      SUBROUTINE MEMGRD
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      COMMON /FMCOM / X(1)
      COMMON /INFGRD/ RHOMIN,ILENG,MAXGRD
      COMMON /MACHIN/ NWDVAR,MAXFM,MAXSM,LIMFM,LIMSM
C
C     ----- GRDDFT
C
C      NOTE: THIS ROUTINE IS FAST, SO CHECK RUNS CAN JUST DO IT ALL.
C
      CALL VALFM(LOADFM)
      CALL GRDDFT(L2,NGRDMEM)
      IGRDDFT = LOADFM  + 1
      LAST    = IGRDDFT + NGRDMEM
      NEED = LAST-LOADFM-1
      CALL GETFM(NEED)
C
C     ----- DFTSET
C     SETTING FOR GRID DFT
C
      CALL DFTSET(X(IGRDDFT),1,.FALSE.)
C
C     ----- MAXGRD: NUMBER OF GRID POITS
C
      CALL MGRDMNG(X(IGRDDFT),DUM,DUM,DUM,MAXGRD,0)
C
C     ----- GRID INFORMATION
C
      CALL VALFM(LOADFM)
      M10  = LOADFM + 1
      M20  = M10    + MAXGRD*3
      M30  = M20    + MAXGRD
      LAST = M30    + MAXGRD/NWDVAR + 1
      NEED2 = LAST-LOADFM-1
      CALL GETFM(NEED2)
      CALL MGRDMNG(X(IGRDDFT),X(M10),X(M20),X(M30),MAXGRD,1)
      CALL GRDWT(22,X(M10),X(M20),X(M30),MAXGRD)
C
      CALL RETFM(NEED+NEED2)
      RETURN
      END
C*MODULE TDDFT  *DECK TDDINP
      SUBROUTINE TDDINP
C
C     READ INFORMATION FROM $TDDFT
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      PARAMETER (NUMLEB=32)
      DIMENSION LEBSOK(NUMLEB)
C
      PARAMETER (MXRT=100, MXATM=2000)
C
      LOGICAL TRIPLET,GOPARR,DSKWRK,MASWRK,SG1T,OK
C
      COMMON /INFOA / NAT,ICH,MUL,NUM,NQMT,NE,NA,NB,
     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /INFGRD/ RHOMIN,ILENG,MAXGRD
      COMMON /INFOTD/ TDM(3,MXRT),CNVTOL,TRIPLET,MOCC(MXRT),MVIR(MXRT),
     *                NRADT,NTHET,NPHIT,NSTAT,NTRIAL,MAXVEC,NTHST,
     *                IRECTD,NLEBT,SG1T
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /WFNOPT/ SCFTYP,VBTYP,DFTYPE,TDDFTYP,CITYP,CCTYP,
     *                MPLEVL,MPCTYP
C
      PARAMETER (NNAM=14)
      DIMENSION QNAM(NNAM),KQNAM(NNAM)
      DATA TDDFT/8HTDDFT   /
      DATA QNAM /8HNSTATE  ,8HIROOT   ,8HMULT    ,8HMAXVEC  ,
     *           8HNTRIAL  ,8HCNVTOL  ,8HMAXGRD  ,8HRHOMIN  ,
     *           8HILENG   ,8HNRAD    ,8HNTHE    ,8HNPHI    ,
     *           8HNLEB    ,8HSG1     /
      DATA KQNAM/1,1,1,1,  1,3,1,3,  1,1,1,1,  1,0/
C
C          74, 230, 266 ARE FORBIDDEN DUE TO NEGATIVE QUADRATURE WEIGHTS
C          THE FIRST THREE ARE JUST RIDICULOUSLY SMALL.
      DATA LEBSOK/-6,-14,-26,38,50,-74,86,110,146,170,194,-230,-266,302,
     *            350,434,590,770,974,1202,1454,1730,2030,2354,2702,
     *            3074,3470,3890,4334,4802,5294,5810/
C
      DATA RHF,UHF,ROHF/8HRHF     ,8HUHF     ,8HROHF    /
C
C     READ THE $TDDFT INPUT GROUP
C        THERE ARE TWO UNDOCUMENTED KEYWORDS:
C     RHOMIN = THRESHOLD FOR DENSITY CUTOFF IN GRID CALCULATIONS
C     ILENG  = SIZE OF SMALL BATCH IN GRID CALCULATIONS
C        OBSOLETE INPUT VALUES:
C     MAXGRD = NOW COMPUTED BY THE RUN ITSELF
C
C          -- GRID INFO
      RHOMIN= 1.0D-12
      ILENG = 200
      MAXGRD= 1000000
C
C          -- TDDFT INFO
      NSTAT  = 1         ! CARRIES KEYWORD -NSTATE-
      MAXVEC =50
      NTHST  = 1         ! CARRIES KEYWORD -IROOT-
      NTRIAL = 5
      MULT   = MUL
      TRIPLET=.FALSE.  ! CARRIES KEYWORD -MULT-, AFTER COPYING BELOW
      CNVTOL =1.0D-07
C                 GRID POINT SELECTION,
C                 IF NO ANGULAR GRID IS INPUT, LEBEDEV WILL BE CHOSEN.
      NRADT  = 48
      NTHET  =  0
      NPHIT  =  0
      NLEBT  =  0
      SG1T   = .FALSE.
C
      JRET     = 0
      CALL NAMEIO (IR,JRET,TDDFT,NNAM,QNAM,KQNAM,
     *             NSTAT,NTHST,MULT,MAXVEC,NTRIAL,CNVTOL,
     *             MAXGRD,RHOMIN,ILENG,NRADT,NTHET,NPHIT,NLEBT,SGIT,
     *     0,0,0,0,0,    0,0,0,0,0,
     *     0,0,0,0,0,    0,0,0,0,0,    0,0,0,0,0,   0,0,0,0,0,
     *     0,0,0,0,0,    0,0,0,0,0,    0,0,0,0,0,   0,0,0,0,0)
      IF (JRET.GT.1) THEN
         WRITE(IW,*) 'ERROR IN $TDDFT GROUP'
         CALL ABRT
         STOP
      END IF
C
      IF (SCFTYP.EQ.RHF) THEN
         IF(MULT.NE.1  .AND.  MULT.NE.3) THEN
            IF(MASWRK) WRITE(IW,9020) MULT
            CALL ABRT
         END IF
         IF(MULT.EQ.3) TRIPLET = .TRUE.
      END IF
      IF (SCFTYP.EQ.ROHF .OR. SCFTYP.EQ.UHF) THEN
         IF (MULT.NE.MUL) THEN
            IF(MASWRK) WRITE(IW,9021) MUL,MULT
            CALL ABRT
         END IF
      END IF
C
C     ECHO THE INPUT
C
      IF(MASWRK) THEN
         IF(SG1T) THEN
            WRITE(IW,9018) NSTAT,NTHST,MULT,
     *                     MAXVEC,NTRIAL,CNVTOL
         ELSE IF (NLEBT.NE.0) THEN
            WRITE(IW,9015) NSTAT,NTHST,MULT,
     *                     MAXVEC,NTRIAL,CNVTOL,
     *                     NRADT,NLEBT
         ELSE
            WRITE(IW,9010) NSTAT,NTHST,MULT,
     *                     MAXVEC,NTRIAL,CNVTOL,
     *                     NRADT,NTHET,NPHIT
        ENDIF
      ENDIF
C
C        SEE IF A SENSIBLE GRID WAS REQUESTED
C
      NERR=0
      IF(SG1T) THEN
         NTHET=0
         NPHIT=0
         NLEBT=0
         IF(SCFTYP.NE.RHF) THEN
            IF(MASWRK) WRITE(IW,9080) SCFTYP
            NERR=NERR+1
         ENDIF
      ELSE
C              NEXT LINE ESTABLISHES THE DEFAULT ANGULAR GRID
         IF(NLEBT.EQ.0  .AND.  (NTHET+NPHIT).EQ.0) NLEBT=110
C----------------DEBUG------------
C---     IF(SCFTYP.EQ.RHF) THEN
C---        IF(NLEBT.EQ.0  .AND.  (NTHET+NPHIT).EQ.0) NLEBT=110
C---     ELSE
C---        IF(NLEBT.EQ.0  .AND.  (NTHET+NPHIT).EQ.0) THEN
C---           NTHET= 8
C---           NPHIT=16
C---        END IF
C---     END IF
C----------------DEBUG------------
C
         IF(NLEBT.GT.0  .AND.  (NTHET+NPHIT).GT.0) THEN
            IF(MASWRK) WRITE(IW,9100) NLEBT,NTHET,NPHIT
            NERR=NERR+1
         END IF
         IF(NLEBT.GT.0) THEN
            OK=.FALSE.
            DO ILEB=1,NUMLEB
               IF(NLEBT.EQ.LEBSOK(ILEB)) OK=.TRUE.
            ENDDO
            IF(.NOT.OK) THEN
               NERR=NERR+1
               IF(MASWRK) WRITE(IW,9110) NLEBT,LEBSOK
            END IF
            IF(NRADT.LE.10.OR.NLEBT.LE.26) THEN
               IF(MASWRK) WRITE(IW,9075) NRADT,NLEBT
               NERR=NERR+1
            ENDIF
         ELSE
            IF(NRADT.LE.10.OR.NTHET.LE.4.OR.NPHIT.LE.8) THEN
               IF(MASWRK) WRITE(IW,9070) NRADT,NTHET,NPHIT
               NERR=NERR+1
            ENDIF
         END IF
      END IF
C
      IF(NERR.GT.0) CALL ABRT
C
C     -- SET XC INDEX
C
      CALL XCINDX
      IRECTD=472
      RETURN
C
 9010 FORMAT(/10X,22(1H-)/10X,'TDDFT INPUT PARAMETERS'/10X,22(1H-)/
     *        1X,'NSTATE=',I8,'  IROOT=',I8,'   MULT=',I8/
     *        1X,'MAXVEC=',I8,' NTRIAL=',I8,' CNVTOL=',1P,E10.2/
     *        1X,'  NRAD=',I8,'   NTHE=',I8,'   NPHI=',I8)
 9015 FORMAT(/10X,22(1H-)/10X,'TDDFT INPUT PARAMETERS'/10X,22(1H-)/
     *        1X,'NSTATE=',I8,'  IROOT=',I8,'   MULT=',I8/
     *        1X,'MAXVEC=',I8,' NTRIAL=',I8,' CNVTOL=',1P,E10.2/
     *        1X,'  NRAD=',I8,'   NLEB=',I8)
 9018 FORMAT(/10X,22(1H-)/10X,'TDDFT INPUT PARAMETERS'/10X,22(1H-)/
     *        1X,'NSTATE=',I8,'  IROOT=',I8,'   MULT=',I8/
     *        1X,'MAXVEC=',I8,' NTRIAL=',I8,' CNVTOL=',1P,E10.2/
     *        1X,'  WITH THE SG-1 PRUNED LEBEDEV GRID')
 9020 FORMAT(/1X,'*** ERROR IN $TDDFT INPUT GROUP: MULT=',I4/
     *       1X,'ONLY MULT=1 OR 3 CAN BE REACHED FROM A CLOSED',
     *          ' SHELL REFERENCE.'/)
 9021 FORMAT(/1X,'*** ERROR IN $TDDFT INPUT GROUP:'/
     *       1X,'OPEN SHELL TD-DFT EXCITATIONS CAN ONLY REACH THE'/
     *       1X,'SAME MULTIPLICITY AS THE GROUND STATE.'/
     *       1X,'YOUR $CONTRL HAS MULT=',I4,' BUT $TDDFT HAS MULT=',I4)
 9070 FORMAT(/1X,'***ERROR: $TDDFT''S NRAD,NTHE,NPHI=',3I5,
     *           ' ARE INSENSIBLY SMALL.')
 9075 FORMAT(/1X,'***ERROR: $TDDFT''S NRAD,NLEB=',2I5,
     *           ' ARE INSENSIBLY SMALL.')
 9080 FORMAT(1X,'***ERROR, SCFTYP=',A8,
     *          ' REQUIRES POLAR COORDINATE GRID, ONLY')
 9100 FORMAT(/1X,'*** ERROR: $TDDFT GRID CONFUSION, NLEB=',I5,
     *           ' -AND- NTHE/NPHI=',2I5,'?'/
     *       1X,'LEBEDEV AND POLAR COORDINATE ANGULAR GRIDS ARE',
     *          ' MUTUALLY INCOMPATIBLE')
 9110 FORMAT(/1X,'*** ERROR: INVALID $TDDFT NLEB=',I5,
     *          ', THE LEGAL VALUES ARE'/
     *       1X,12I6/1X,12I6/1X,8I6/
     *       1X,'PICK ONLY FROM THE POSITIVE VALUES!')
      END
C*MODULE TDDFT  *DECK TDDFTCALC
      SUBROUTINE TDDFTCALC
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      PARAMETER (MXATM=2000, MXRT=100)
      PARAMETER (MXFRG=1050, MXDFG=5, MXDPPT=MXFRG*MXDFG*12)
      PARAMETER (MXGRID=10)
C
      PARAMETER (ZERO=0.0D+00,ONE=1.0D+00,TWO=2.0D+00)
C
      LOGICAL DBUG,TRIPLET,NONEQR,IPCFP,
     *        GOPARR,DSKWRK,MASWRK,PACK2E,SG1,SG1T
C
      COMMON /DFGRID/ DFTTHR,DFTGTHR,SWOFF,SW0,BSLRD(137),NDFTFG,
     *                NRAD,NTHE,NPHI,NRAD0,NTHE0,NPHI0,
     *                NANGPT(MXGRID),NANGPT0(MXGRID),SG1
      COMMON /DFTPAR/ DFTTYP(20),EXENA,EXENB,EXENC,IDFT34,NAUXFUN,
     *                                                    NAUXSHL
      COMMON /FMCOM / X(1)
      COMMON /FMOINF/ NFG,NLAYER,NATFMO,NBDFG,NAOTYP,NBODY
      COMMON /FRGINF/ NMPTS(MXFRG),NMTTPT,IEFC,IEFD,IEFQ,IEFO,
     *                NPPTS(MXFRG),NPTTPT,IEFP,
     *                NRPTS(MXFRG),NRTTPT,IREP,ICHGP,NFRG,
     *                NDPPTS(MXDPPT),NDPTTPT,IEFDP,LSTMPTS(MXFRG)
      COMMON /INFOA / NAT,ICH,MUL,NUM,NQMT,NE,NA,NB,
     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
      COMMON /INFGRD/ RHOMIN,ILENG,MAXGRD
      COMMON /INTFIL/ NINTMX,NHEX,NTUPL,PACK2E,INTTYP,IGRDTYP
      COMMON /INFOTD/ TDM(3,MXRT),CNVTOL,TRIPLET,MOCC(MXRT),MVIR(MXRT),
     *                NRADT,NTHET,NPHIT,NSTAT,NTRIAL,MAXVEC,NTHST,
     *                IRECTD,NLEBT,SG1T
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /PCKLAB/ LABSIZ
      COMMON /PCMDIM/ MXSP,MXTS,MEMPCM1,MEMPCM2,NTS
      COMMON /PCMITR/ RCUT(2),THRES,IPCMIT,IMUL,MXDIIS,NREG,MXITR1,
     *                MXITR2,MODPAP
      COMMON /PCMOPT/ RABI,RASC,REFPOL,THRSLS,DENSLS,WB,WA,ETA2,GD,EVAC,
     *                RHOW,PM,AREATL,AREAKP,BONDRY,OMEGA,RET,FRO,EPSINF,
     *                EPS,DR,RSOLV,VMOL,TCE,STEN,DSTEN,CMF,TABS,IDIRCT,
     *                IPCDER,IDP,ICOMP,IFIELD,ICAV,IDISP,IPRINT,IRETCAV,
     *                ICENT,IFAST,NEVAL,IEFPOL,KEEPSM,IMGABI,IMGASC,NADD
      COMMON /PCMPAR/ IPCM,NFT26,NFT27,IRPPCM,IEF,IP_F,ITDPCM
      COMMON /PCMPNT/ LXYZRE,LSSFE,LLIST,LALPCM,LRINPCM,LINAPCM,LINFPCM,
     *                LINIPCM,LMEPCM,LAXYZCT,LXYZCT2,LNVERT,LQSN,LQSE,
     *                LQFS,LQIND,LISPHE,LVAD,LQOR,LVECMUL,LAIPRJ,
     *                LFIPRJ,LPEL
      COMMON /RUNOPT/ RUNTYP,EXETYP,NEVALS,NGLEVL,NHLEVL
      COMMON /WFNOPT/ SCFTYP,VBTYP,DFTYPE,TDDFTYP,CITYP,CCTYP,
     *                MPLEVL,MPCTYP
C
      DATA CHECK/8HCHECK   /
      DATA RHF/8HRHF     /
      DATA MXITER/100/
C
C     ----- DRIVER FOR CLOSED SHELL + TDDFT EXCITATION ENERGIES -----
C
C        OF COURSE, MANY OF THE EXCITED STATES ARE NOT TOTALLY SYMMETRIC
C
      CALL SYMOFF
C
C        DO INITIALIZATIONS FOR POSSIBLE TDDFT/EFP1
C
      IF (NFRG.GT.0 .AND. .NOT.TRIPLET) THEN
         CALL TDEFPINIT
      END IF
C
      CALL VCLR(TDM,1,3*MXRT)
      DBUG=.FALSE.
      NONEQR=.TRUE.
C     --- DEFINE NOCC & NVIR
      LX=NQMT
      NOCC=NA
      NVIR=LX-NOCC
      L0=NINTMX
      L1=NUM
      L2=L1*(L1+1)/2
      L3=L1*L1
C        L7 USED FOR COMPUTATIONS AND STORAGE, L7MAX USED FOR DA-WRITE,
C        AS THE NUMBER OF VIRTUAL ORBITALS CAN FLUCTUATE DURING A RUN.
      L7=NOCC*NVIR
      L7MAX = NOCC*(NUM-NOCC)
      CSCAL=ONE
      IF(TRIPLET) CSCAL=ZERO
      IPF=1
      IMF=0
      IF(DFTTYP(3).NE.ZERO) IMF=1
      IF(DBUG) WRITE(IW,*)'NOCC,NVIR,L7=',NOCC,NVIR,L7
C     --- DEFINE MXVEC & NDSR & NVEC:NUMBER OF INITIAL VECORS
      NDSR=NSTAT
      NVEC=NTRIAL
      IF(NDSR.GT.L7) NDSR=L7
      MXVEC=MAXVEC*NDSR
      IF(MXVEC.GT.L7) MXVEC=L7
C     --LAST CHECK
      IF(NVEC.GT.MXVEC) NVEC=MXVEC
      IF(NVEC.LT.NDSR) NVEC=NDSR
      IF(DBUG) THEN
      WRITE(IW,*)'ILENG=',ILENG
      WRITE(IW,*)'NDSR=',NDSR
      WRITE(IW,*)'L7=',L7
      WRITE(IW,*)'NOCC=',NOCC
      WRITE(IW,*)'NVIR=',NVIR
      WRITE(IW,*)'NTHST=',NTHST
      WRITE(IW,*)'MXVEC=',MXVEC
      WRITE(IW,*)'MAXGRD=',MAXGRD
      ENDIF
C
      NMAX=MAX(2*NDSR,NVEC)
      MXVEC2=MXVEC*(MXVEC+1)/2
      IST=1
      IEND=NVEC
C
C     TDPCM
      IF(IPCM.EQ.1) THEN
       EPS=EPSINF
C
       ITDPCM=1
      ENDIF
C     TDPCM
      CALL GOTFM(NGOTMX)
      CALL VALFM(LOADFM)
C     -- MEM1
      I010=LOADFM+1
C     I010:IPR(2,2,L3)
      I020=I010+4*L3
C     I020:INXOV(2,L7)
      I10=I020+L7*2
C     I10:E(L1)
      I20=I10+L1
C     I20:XME(L7):ABD(L7)
      I30=I20+L7
C     I30:IM(L7) OR INDX(NVEC)
      I40=I30+L7
C     I40:BVEC(L7,MXVEC)
      K10=I40+L7*MXVEC
C     K10:BVEC_AO(L1,L1,NMAX) OR SCRATCH
      K20=K10+L1*L1*NMAX
C     K20:B1PROA(L2,NMAX)
      K30=K20+L2*NMAX
C     K30:B2PROA(L1,L1,NMAX) OR SCRACTH
      K40=K30+L3*NMAX
C     K40:ABXC(L2) OR SCRATCH
      K50=K40+L3
C     K50:XX(L0)
      K60=K50+L0
C     K60:IX(L0*LABSIZ)
      K70=K60+L0*LABSIZ
C     K70:MO(L1,L1)
      K80=K70+L3
C     K80:SCR(NOCC,L1) OR SCRACTH
      I50=K80+L1*L1
C     I50:B1PRO_MO(L7,MXVEC)
      I60=I50+L7*MXVEC
C     I60:B2PRO_MO(L7,MXVEC)
      I70=I60+MXVEC*L7
C     I70:APB(MXVEC,MXVEC)
      I80=I70+MXVEC*MXVEC
C     I80:AMB(MXVEC,MXVEC)
      I90=I80+MXVEC*MXVEC
C     I90:SCR(MXVEC*8)
      I100=I90+MXVEC*8
C     I100:VL(MXVEC,MXVEC)
      I110=I100+MXVEC*MXVEC
C     I110:VR(MXVEC,MXVEC)
      I120=I110+MXVEC*MXVEC
C     I120:EE(MXVEC)
      I130=I120+MXVEC
C     I130:IPVT(MXVEC*MXVEC)
      I140=I130+MXVEC*MXVEC
C     I140:VLO(L7,NDSR)
      I150=I140+L7*NDSR
C     I150:VRO(L7,NDSR)
      I160=I150+L7*NDSR
C     I160:SCR(MXVEC2)
      I170=I160+MXVEC2
C     I170:ABP_R(L7,NDSR) OR W_L(L7,NDSR)
      I180=I170+L7*NDSR
C     I180:ABM_L(L7,NDSR) OR W_R(L7,NDSR)
      I190=I180+L7*NDSR
C     I190:Q(L7,NDSR,2)
      I200=I190+L7*NDSR*2
      I220=I200
C     -- MEM2: FOR CONVERGED
C     I220:TEMP(L1*L1)
      I230=I220+L1*L1
C     I230:---
      I240=I230
C     I240:---
      I250=I240
C     I250:---
      I260=I250
C     I260:DIX(L2)
      I270=I260+L2
C     I270:DIY(L2)
      I280=I270+L2
C     I280:DIZ(L2)/T(L2)
      I290=I280+L2
C     I290:--
      I300=I290
C     I300:TXYZ(3,NDSR)
      I310=I300+3*NDSR
C     I310:OS(NDSR)
      I315=I310+NDSR
C     I315:STATE SYMMETRIES
      I320=I315+NDSR
C     I320:AOV(L1,L1)
      I330=I320+L1*L1
C     I330:--
      LASTI=I330
C
C     ----- MEMORY FOR DFT GRID
C     -- MEM3
      IF (NDFTFG.EQ.1) THEN
      J10=I220+1
C     J10:XYZGRD(MAXGRD*3)
      J20=J10+MAXGRD*3
C     J20:XYZWGT(MAXGRD)
      J30=J20+MAXGRD
C     J30:XYZDCH
      J40=J30+4*L1*ILENG
C     J40:D(L2)
      J50=J40+L2
C     J50:DCC(NAT*NAT)
      J60=J50+(NAT*NAT)
C     J60:ACC(NAT*NAT)
      J70=J60+(NAT*NAT)
C     J70:PCC(NAT*NAT)
      J80=J70+(NAT*NAT)
C     J80:PC(NAT*NAT)
      J90=J80+(NAT*NAT    )
C     J90:--
      J100=J90
C     J100:COEF(3,ILENG) OR COEF(3,2,ILENG)
      J110=J100+6*ILENG
C     J110:DTRAI_(4,ILENG)
      J120=J110+4*ILENG
C     J120:DRHOI_(MAXGRD,4)
      J130=J120+4*MAXGRD
C     J130:EX0(ILENG)
      J140=J130+ILENG
C     J140:EX(ILENG,10)
      J150=J140+10*ILENG
C     J150:EC0(ILENG)
      J160=J150+ILENG
C     J160:EC(ILENG,20)
      J170=J160+20*ILENG
      J180=J170
C     J180:AOMAX(L1)
      J190=J180+L1
C     J190:IAO(L1)
      LASTJ=J190+L1
      LAST=MAX(LASTI,LASTJ)
      ELSE
      LAST=LASTI
      ENDIF
C
C     -- MEMORY FOR THE PCM SOLVATION METHOD
C     MEMORY ALLOCATIONS UPDATED NOVEMBER 2006 BY MWS,
C     BUT WITHOUT ANY ATTEMPT TO VERIFY THE CORRECTNESS.
C
C     THIS IS ALMOST CERTAINLY NON-FUNCTIONAL CODE, AND SHOULD
C     PROBABLY BE REMOVED ENTIRELY.  OR PERHAPS FIXED TO WORK.
C
      ISD=1  ! NOT SURE WHAT THIS SHOULD BE...
      IPCFP=.FALSE.
      CALL PCMMEM3(IPCFP,L2,LAST,NTSPAR,LSOL,LCHG,LDMATM,LQPOT,LQFLD,
     *             LVPOT,LSE,LDE,LQ0,LQ1,LQ2,LQ3,LD0,LQA,LDIMAT,LQREP,
     *             LTMP0,LTMP1,LTMP2,LIPVT,LPOTTMP,LRMUL,LCHG2,
     *             LASCCRD,LASCCHG,LASCDIP,LASCQAD,LDISV,LDIS1,LDIS2,
     *             LCQEF,LQEFF,LELD,LXDINT,LYDINT,LZDINT,LEFLD,MADD,
     *             LABFLD,LVECMP,LVEC_2,LVEC_S,LAPROJ,LFPROJ)
      IF(IPCM.EQ.1) THEN
         LWRK1 = LAST
         LWRK2 = LWRK1 + L3
         LWRK3 = LWRK2 + L3
         LSCR  = LWRK3 + L3
         LAST  = LSCR  + 8*L1
      END IF
C
C         STORAGE FOR EFP SOLVATION MODEL:
C
      IF (NFRG.GT.0) THEN
         IF(SCFTYP.EQ.RHF) THEN
C                THE POLARIZATION TERM DROPS OUT OF THE EQUATION
C                IF THE EXCITATION IS TRIPLET COUPLED,
C                BUT SINGLET COUPLING DOES NEED EFP1 STORAGE.
            IF(TRIPLET) THEN
               CONTINUE
            ELSE
               LXEFI = LAST
               LYEFI = LXEFI + L2
               LZEFI = LYEFI + L2
               LAST  = LZEFI + L2
            END IF
         ELSE
C                OPEN SHELL HAS NO EFP SOLVENT CAPABILITY,
C                BUT THIS IS ALREADY TRAPPED WHEN THE JOB STARTS.
C                IT WOULD BE REPETITIVE AND REDUNDANT TO TRAP AGAIN.
            CONTINUE
         END IF
      END IF
C
      NEED=LAST-LOADFM-1
      IF(MASWRK) WRITE(IW,9010) NEED
      CALL GETFM(NEED)
C
C---  MEM1=I220-LOADFM-1
C---  MEM2=LASTI-I220
C---  IF(NDFTFG.EQ.1) MEMXC=LASTJ-I220
C---  IF(NDFTFG.EQ.1) THEN
C---     IF(MASWRK) WRITE(IW,*) 'MEM1/MEM2/MEMXC=',MEM1,MEM2,MEMXC
C---  ELSE
C---     IF(MASWRK) WRITE(IW,*) 'MEM1/MEM2=',MEM1,MEM2
C---  ENDIF
C
      IF(EXETYP.EQ.CHECK) GO TO 800
C
C        READ INVERSE D MATRIX FOR PCM (SKIP FOR ITERATIVE SOLVER)
C
      IF(IPCMIT.EQ.0  .AND.  IPCM.EQ.1) THEN
         NTS3 = NTS*NTS
         CALL SEQREW(NFT26)
         CALL DERCHK(NDER)
         IF(NDER.EQ.1.AND.IPCDER.EQ.1) THEN
            CALL SEQADV(NFT26)
            CALL SEQADV(NFT26)
            CALL SEQADV(NFT26)
            CALL SEQADV(NFT26)
         END IF
         IF(IEF.GT.0.AND.IEF.LT.3) THEN
            CALL SEQADV(NFT26)
            CALL SEQADV(NFT26)
            CALL SQREAD(NFT26,X(LSE),NTS3)
            CALL SQREAD(NFT26,X(LDE),NTS3)
         END IF
         CALL SQREAD(NFT26,X(LDMATM),NTS3)
         CALL SEQREW(NFT26)
      END IF
C
C     -- READ MO AND DENSITY
C
      CALL DAREAD(IDAF,IODA,X(K70),L1*LX,15,0)
      IF(NDFTFG.EQ.1) CALL DAREAD(IDAF,IODA,X(J40),L2,16,0)
C
C     --CREATE DRHOI_ & XC DERIVATIVES
C
      IF (NDFTFG.EQ.1) THEN
         IF(DBUG) WRITE(IW,*)'SET TDDFT INFORMATION '
         NPTGRD=MAXGRD
         CALL TDDFTSET(RHOSET,X(J10),X(J20),X(J30),X(K70),
     *                 X(J120),X(J180),X(I010),ILENG,NPTGRD,L1)
      ENDIF
C
C     --- INITIAL TRIAL VECTORS BVEC(L7,NVEC)
C
      CALL DAREAD(IDAF,IODA,X(I10),LX,17,0)
      CALL TDINIVEC(X(I10),X(I020),X(I20),X(I30),X(I40),LX,L7,
     *              NOCC,NVEC,MXVEC)
      IF(MASWRK) WRITE(IW,*)'INITIAL VECTORS / MXVEC ',NVEC,'/',MXVEC
C     --NEW IEND
      IEND=NVEC
C
C     ---- TDDFT LOOP START ----
C
      ITER=0
   10 CONTINUE
      ITER=ITER+1
      IF(MASWRK) WRITE(IW,*)'ITERATION # ',ITER
      IF(ITER.GT.MXITER) THEN
       IF(MASWRK) THEN
       WRITE(IW,*)'  .. OVER MAXIMUM ITERATION NUMBER ..'
       WRITE(IW,*)'ITER,MXITER=',ITER,MXITER
       ENDIF
       ICK=1
       GOTO 20
      ENDIF
C ********************************************************************
C
C     --- MATRIX-VECTOR PRODUCTS
C         (A+B)_B & (A-B)_B
C
      NV=IEND-IST+1
      IF(MASWRK) WRITE(IW,9020) IST,IEND,NV
C
C     -- AO BASIS
C
      IV=0
      DO IVEC=IST,IEND
      IV=IV+1
       CALL IATOGENTD(X(I40),X(K40),LX,L7,NOCC,IVEC,MXVEC)
C      SUBROUTINE IATOGENTD(PV,AV,LX,L7,NOCC,IVEC,MXVEC)
       CALL DGEMM('N','N',L1,LX,LX,ONE,X(K70),L1,X(K40),LX,ZERO
     * ,X(K30),L1)
       CALL DGEMM('N','T',L1,L1,LX,ONE,X(K30),L1,X(K70),L1,ZERO
     * ,X(K10+(IV-1)*L1*L1),L1)
C     S UBROUTINE DGEMM(FORMA,FORMB,L,N,M,ONE,A,LDA,B,LDB,ZERO,C,LDC)
      ENDDO
C
C     -- 2E : A+B AND A-B
C
      CALL VCLR(X(K20),1,L2*NMAX)
      CALL VCLR(X(K30),1,L1*L1*NMAX)
      CALL TD2E(X(K10),X(K20),X(K30),X(K50),X(K60),L1,CSCAL,IMF,IPF,NV)
C
C     --- EFP1 ---
C
      IF (NFRG.GT.0 .AND. .NOT.TRIPLET) THEN
         CALL TDEFP (X(K10),X(K20),X(LXEFI),X(LYEFI),X(LZEFI),
     *        L2,L1,NV)
      END IF
C
C     ****** VECTOR LOOP ENDS STARTS HERE *******
C
      IV=0
      DO IVEC=IST,IEND
      IV=IV+1
C
C     -- PCM EFFECT
C     ARGUMENTS TO PCMFLD UPDATED NOVEMBER 2006 BY MWS,
C     BUT WITHOUT ANY ATTEMPT TO VERIFY THE CORRECTNESS.
C
       IF(IPCM.EQ.1.AND.NONEQR) THEN
          CALL TDPCMDEN(X(K40),X(K10+(IV-1)*L1*L1),L1)
          CALL PCMFLD(X(K20+(IV-1)*L2),DUMMY,X(K40),X(LSOL),
     *               X(LCHG),X(LWRK1),X(LDISV),
     *               X(LDIS1),X(LDIS2),X(LWRK3),X(LSCR),X(LDMATM),
     *               X(LCQEF),X(LQEFF),X(LELD),X(LQPOT),X(LQFLD),
     *               X(LVPOT),X(LSE),X(LDE),
     *               X(LQ0),X(LQ1),X(LQ2),X(LQ3),
     *               X(LD0),X(LQA),X(LDIMAT),X(LQREP),
     *               X(LTMP0),X(LTMP1),X(LTMP2),X(LIPVT),
     *               X(LPOTTMP),X(LRMUL),
     *               TCH,QET,QETN,QESC,ISD,NFT27,IPCFP,L1,L2,NTSPAR,
     *               X(LAXYZCT),X(LAXYZCT+MXTS),X(LAXYZCT+MXTS*2),
     *               X(LXYZCT2),X(LXYZCT2+MXTS),X(LXYZCT2+MXTS*2),
     *               X(LAXYZCT+MXTS*3),X(LQSN),X(LQSE),X(LQFS),
     *               X(LVECMUL),X(LCHG2),X(LQIND+MXTS),X(LPEL),
     *               X(LASCCRD),X(LASCCHG),X(LASCDIP),X(LASCQAD),
     *               X(LXYZRE),X(LXYZRE+MXSP),X(LXYZRE+MXSP*2),
     *               X(LXYZRE+MXSP*3),X(LISPHE),X(LLIST))
       END IF
C
C     -- DFT FUNCTIONAL PART
C      -- HALF-TRANSFORMING
C
       CALL IATOGENTD(X(I40),X(K40),LX,L7,NOCC,IVEC,MXVEC)
       CALL DGEMM('N','T',LX,L1,LX,ONE,X(K40),LX,X(K70),
     *            L1,ZERO,X(K80),LX)
       IF(NDFTFG.EQ.1) THEN
        CALL VCLR(X(K40),1,L2)
        CALL TDFXCP(X(K40),RHO,X(J10),X(J20),X(J30),X(K80),X(K70),
     *              X(J120),X(J110),X(J100),X(J140),X(J160),X(J130),
     *              X(J150),X(J180),X(I010),X(J190),ILENG,NPTGRD,
     *              L1,L2,LX,2,1)
        IF(GOPARR.AND.IBTYP.EQ.1) THEN
           CALL DDI_GSUMF(2330,X(K40),L2)
           CALL DDI_GSUMF(2331,RHO,1)
        ENDIF
        CALL DAXPY(L2,TWO,X(K40),1,X(K20+(IV-1)*L2),1)
       ENDIF
C
C     -- PRODUCT A+B *X
C
       CALL CANTOGEN(X(K20+(IV-1)*L2),X(K40),L1)
       CALL MNTOIA(X(K40),X(I50),X(K70),X(K80),L1,LX,L7,NOCC,IVEC)
C      SUBROUTINE MNTOIA(PAO,PMO,V,SCR,L1,LX,L7,NOCC,IVEC)
       CALL ESUM(X(I10),X(I50),X(I40),LX,L7,NOCC,IVEC)
C      SUBROUTINE ESUM(E,PMO,Z,LX,L7,NOCC,IVEC)
C
C     -- PRODUCT A-B *X
C
       CALL MNTOIA(X(K30+(IV-1)*L1*L1),X(I60),X(K70),X(K80),
     * L1,LX,L7,NOCC,IVEC)
       CALL ESUM(X(I10),X(I60),X(I40),LX,L7,NOCC,IVEC)
      ENDDO
C
C     ****** VECTOR LOOP ENDS JUST ABOVE *******
C
C     --- SMALL REDUCED MATRICES A+B,A-B.
C
      CALL RPAREDMS(X(I40),X(I50),X(I60),X(I70),X(I80),NVEC,L7)
C
C     --- FORM REDUCED RPA MATRIX AND DIAGONALIZE IT
C
      CALL RPAEIG(X(I120),X(I100),X(I110),X(I70),X(I80),
     *            X(I90),X(I160),X(I30),NVEC)
C
C     --- IMPOSING BIORTHOGONALITY CONDITION ON V_L & V_R
C
      CALL RPAVNORM(X(I110),X(I100),NVEC)
C     SUBROUTINE RPAVNORM(VR,VL,NVEC)
C
C     -- PRINT CURRENT EXCITATION ENERGIES
C
      CALL RPAECHK(X(I120),X(I30),NVEC,NDSR,IMAX)
C     SUBROUTINE RPAECHK(EE,INDX,NVEC,NDSR,IMAX)
C
C     -- VRO(L7,*), VLO(L7,*), (A+B)VRO, (A-B)VLO
C
      CALL RPAEXPNDV(X(I110),X(I100),X(I150),X(I140),X(I40),X(I40),
     *               NDSR,NVEC,L7,X(I30))
      CALL RPAEXPNDV(X(I110),X(I100),X(I170),X(I180),X(I50),X(I60),
     *               NDSR,NVEC,L7,X(I30))
C
C     -- RESIDUAL VECTORS W AND CHECK CONVERGENCE
C        AND GET PERTURBED VECTORS Q IF REQUIRED
C
      CALL RPARESVEC(X(I190),X(I170),X(I180),X(I140),X(I150),
     *               X(I120),X(I20),L7,NDSR,ICK,ERROR,CNVTOL,IMAX)
C
C     -- FOR CONVERGED CASE
C
      IF (ICK.EQ.0) GOTO 20
      IF (ICK.EQ.1.AND.NVEC.EQ.MXVEC) THEN
      IF(MASWRK) WRITE(IW,*)'  ..SOMETHING IS WRONG.. NVEC=MXVEC'
      GOTO 20
      ENDIF
C
C     -- ORTHONORMALIZE Q(L7,NDSR*2) AND APPEND NEW B VECTORS
C
      CALL RPANEWB(L7,NDSR,MXVEC,X(I40),X(I190),NOVEC,NVEC,ICK)
C          RPANEWB(NDIM,NDSR,MXVEC,BVEC,Q,NOVEC,NVEC,ICK)
C     ICK=1 NVEC OVER MXVEC:NOT CONVERGED CASE
      IF(ICK.EQ.1) THEN
        IF(MASWRK) THEN
        WRITE(IW,*)'  ..SOMETHING IS WRONG.. NVEC.GT.MXVEC'
        WRITE(IW,*)'NVEC/MXVEC=',NVEC,'/',MXVEC
        ENDIF
      GOTO 20
      ENDIF
C
C     --- UPDATE IST AND IVEC
C
      IST=NOVEC+1
      IEND=NVEC
      GOTO 10
C
   20 CONTINUE
      IF(ICK.EQ.0) THEN
       IF(MASWRK) THEN
       WRITE(IW,*)' '
       WRITE(IW,*)'  ---------------------------------  '
       WRITE(IW,*)'    R-TDDFT CALCULATION CONVERGED'
       WRITE(IW,*)'  ---------------------------------  '
       ENDIF
      ELSE
       IF(MASWRK) THEN
       WRITE(IW,*)' '
       WRITE(IW,*)'  XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX  '
       WRITE(IW,*)'    R-TDDFT CALCULATION NOT CONVERGED'
       WRITE(IW,*)'  XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX  '
       ENDIF
      ENDIF
C
C     ---- TDDFT CALCULATION FINISHED ----
C     RECOVER POINT GROUP SYMMETRY SO THAT WE CAN LABEL THE STATES
C
      CALL SYMON
C
      IF(MASWRK) THEN
         WRITE(IW,9977) NVEC ,MXVEC
         WRITE(IW,9978) ITER, MXITER
         WRITE(IW,9979) L7
         WRITE(IW,9980) NUM
         IF(NDFTFG.EQ.1) WRITE(IW,2111) RHOSET
      ENDIF
C
C     --- OSCILLATOR STRENGTH
C
      IF(.NOT.TRIPLET) THEN
       CALL TDOSCALC(X(I310),X(I300),X(I320),X(K30),X(I150),X(I140),
     *               X(I120),X(K70),X(I260),X(I270),X(I280),
     *               L1,LX,L2,L7,NDSR,X(I020))
      ELSE
       CALL VCLR(X(I310),1,NDSR)
       CALL VCLR(X(I300),1,NDSR*3)
      ENDIF
C
C     --- STORE THE MOST DOMINANT TRANSITION
C
      NSTMAX=NDSR
      DO IST=1,NSTMAX
       CALL MOCVIR(X(I150),IST,L7,MOCC(IST),MVIR(IST),X(I020))
      ENDDO
C
      IF(ICK.EQ.1) THEN
         IF(MASWRK) WRITE(IW,*) 'FAILED IN TDDFT STEP. STOP'
      ELSE
         CALL TDPRNT(X(I310),X(I315),X(I300),X(I020),NDSR,
     *               X(I120),X(I150),L7,IMAX)
      ENDIF
      IF(NFG.NE.0) CALL DAWRIT(IDAF,IODA,X(I150),L7MAX*NDSR,471,0)
C
C     --- PREPARATION FOR GRADIENT CALCULATION
C
C     --- X+Y AND X-Y
      IPOS=L7*(NTHST-1)
      FNOR =DDOT(L7,X(I150+IPOS),1,X(I140+IPOS),1)
      SCAL=SQRT(0.5D+00/FNOR)
      CALL DSCAL(L7,SCAL,X(I150+IPOS),1)
      CALL DSCAL(L7,SCAL,X(I140+IPOS),1)
C          MAYBE WRITE MORE DATA THAN WE HAVE COMPUTED, L7MAX VS. L7
      CALL DAWRIT(IDAF,IODA,X(I150+IPOS),L7MAX,IRECTD,0)
      CALL DAWRIT(IDAF,IODA,X(I140+IPOS),L7MAX,IRECTD+1,0)
C     --- X+Y AND X-Y IN AO-BASIS @ 267 AND 268
C     X+Y
      CALL IATOGENTD(X(I150+IPOS),X(I220),LX,L7,NOCC,1,1)
      CALL DGEMM('N','N',L1,LX,LX,ONE,X(K70),L1,X(I220),LX,ZERO
     * ,X(I320),L1)
      CALL DGEMM('N','T',L1,L1,LX,ONE,X(I320),L1,X(K70),L1,ZERO
     * ,X(I220),L1)
      CALL DAWRIT(IDAF,IODA,X(I220),L1*L1,IRECTD+2,0)
C     X-Y
      CALL IATOGENTD(X(I140+IPOS),X(I220),LX,L7,NOCC,1,1)
      CALL DGEMM('N','N',L1,LX,LX,ONE,X(K70),L1,X(I220),LX,ZERO
     * ,X(I320),L1)
      CALL DGEMM('N','T',L1,L1,LX,ONE,X(I320),L1,X(K70),L1,ZERO
     * ,X(I220),L1)
      CALL DAWRIT(IDAF,IODA,X(I220),L1*L1,IRECTD+3,0)
C
  800 CONTINUE
      CALL RETFM(NEED)
      IF(MASWRK) WRITE(IW,9030)
      CALL TIMIT(1)
      RETURN
C
 2111 FORMAT(' NUMBER OF TOTAL ELECTRONS    =',F20.7)
 9977 FORMAT(' NUMBER OF USED VECTORS / MAX =',I4,1X,'/',I4)
 9978 FORMAT(' NUMBER OF ITERATIONS   / MAX =',I4,1X,'/',I4)
 9979 FORMAT(' NUMBER OF SINGLE EXCITATIONS =',I6)
 9980 FORMAT(' NUMBER OF ATOMIC ORBITALS    =',I4)
 9010 FORMAT(1X,'MEMORY REQUIRED FOR TD-DFT ITERATIONS IS',I10,
     *          ' WORDS.'/)
 9020 FORMAT(1X,'   CURRENT VECTORS ',I4,' TO',I4,' (',I4,' VECTORS)')
 9030 FORMAT(1X,'..... DONE WITH TD-DFT EXCITATION ENERGIES .....')
      END
C
C*MODULE TDDFT   *DECK GENTOCAN
      SUBROUTINE GENTOCAN(GEN,CAN,L1)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      COMMON /IJPAIR/ IA(8192)
      DIMENSION GEN(L1,*)
        DIMENSION CAN(*)
      DO I=1,L1
      DO J=1,I
      IJ=IA(I)+J
      CAN(IJ)=GEN(I,J)
      ENDDO
      ENDDO
      RETURN
        END
C*MODULE TDDFT  *DECK TDINIVEC
      SUBROUTINE TDINIVEC(E,INXOV,XM,IM,BVEC,LX,L7,NOCC,NVEC,MXVEC)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      LOGICAL DBUG
      DIMENSION E(LX)
      DIMENSION XM(L7)
      DIMENSION IM(L7)
      DIMENSION BVEC(L7,MXVEC)
      DIMENSION INXOV(2,*)
      DATA MXINIVEC/50/
C
      DBUG=.FALSE.
C     -- MXINIVEC
      IF(MXINIVEC.GT.MXVEC) MXINIVEC=MXVEC
C     -- SET INXOV(2,L7) & XM(L7)
      DO I=1,NOCC
      DO J=NOCC+1,LX
      JJ=J-NOCC
      IJ=(JJ-1)*NOCC+I
      INXOV(1,IJ)=I
      INXOV(2,IJ)=J
      XM(IJ)=E(J)-E(I)
      ENDDO
      ENDDO
C     -- SAVE PREVIOUS ORDERING
      DO I=1,L7
       IM(I)=I
      ENDDO
C     -- ORDERING XM(L7):XM(SMALL) =< XM(LARGE)
C     ** GET SMALLER DIAGONAL VALUES
      DO 20 J=1,L7-1
      DO 20 I=J+1,L7
       IF(XM(J).LE.XM(I)) GO TO 20
       XMJ=XM(J)
       IMJ=IM(J)
       XM(J)=XM(I)
       IM(J)=IM(I)
       XM(I)=XMJ
       IM(I)=IMJ
   20 CONTINUE
C     -- DBUG
      IF (DBUG) THEN
      WRITE(6,*)'PRINT XM(L7) ORDERING'
      DO I=1,L7
      WRITE(6,*)'I,XM(I),IM(I)=',I,XM(I),IM(I)
      ENDDO
      ENDIF
C
C     ----  GET INITIAL VECTORS: BVEC(L7,NVEC)
C
      CALL VCLR(BVEC,1,L7*NVEC)
      DO K=1,NVEC
         BVEC(IM(K),K)=1.0D+00
      ENDDO
C
C     -- SET XM(L7) AGAIN
C
      DO I=1,NOCC
      DO J=NOCC+1,LX
      JJ=J-NOCC
      IJ=(JJ-1)*NOCC+I
      XM(IJ)=E(J)-E(I)
      ENDDO
      ENDDO
C
      RETURN
      END
C*MODULE TDDFT  *DECK RPAREDMS
      SUBROUTINE RPAREDMS(B,AP_B,AM_B,XM_P,XM_M,NVEC,NDIM)
C
C     CONSTRUCT SMALL REDUCED (A+B) AND (A-B) MATRICES
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DIMENSION B(NDIM,NVEC)
      DIMENSION AP_B(NDIM,NVEC)
      DIMENSION AM_B(NDIM,NVEC)
      DIMENSION XM_P(NVEC,NVEC)
      DIMENSION XM_M(NVEC,NVEC)
      DO I=1,NVEC
      DO J=1,NVEC
       XM_P(I,J)=DDOT(NDIM,B(1,I),1,AP_B(1,J),1)
       XM_M(I,J)=DDOT(NDIM,B(1,I),1,AM_B(1,J),1)
      ENDDO
      ENDDO
      RETURN
      END
C*MODULE TDDFT  *DECK RPAEIG
      SUBROUTINE RPAEIG(EE,VL,VR,APB,AMB,SCR2,SCR,IWRK,NDIM)
C
C     DIAGONALIZE SMALL REDUCED RPA MATRIX
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DIMENSION APB(NDIM,NDIM),AMB(NDIM,NDIM)
      DIMENSION VL(NDIM,NDIM), VR(NDIM,NDIM)
      DIMENSION EE(NDIM)
      DIMENSION SCR(*)
      DIMENSION SCR2(NDIM*8)
      DIMENSION IWRK(*)
C
      DATA ONE /1.0D+00/
      DATA ZERO/0.0D+00/
C
C     1. SQRT(A-B)
C
      IF(NDIM.EQ.1) THEN
        AMB(1,1) = SQRT(AMB(1,1))
      ELSE
        CALL GENTOCAN(AMB,SCR,NDIM)
        CALL GLDIAG(NDIM,NDIM,NDIM,SCR,SCR2,EE,VR,IERR,IWRK)
        DO I = 1, NDIM
           EE(I) = SIGN(SQRT(ABS(EE(I))),EE(I))
        ENDDO
        DO J = 1, NDIM
        DO I = 1, NDIM
         VL(I,J) = VR(J,I) * EE(I)
        ENDDO
        ENDDO
        CALL DGEMM('N','N',NDIM,NDIM,NDIM,ONE,VR,NDIM,
     *            VL,NDIM,ZERO,AMB,NDIM)
      ENDIF
C
C     2. FORM SQRT(A-B)*(A+B)*SQRT(A-B)
C
C     -- (A+B)*SQRT(A-B) : VL
      CALL DGEMM('N','N',NDIM,NDIM,NDIM,ONE,APB,NDIM,
     *            AMB,NDIM,ZERO,VL,NDIM)
C     -- SQRT(A-B)*(A+B)*SQRT(A-B) : VR
      CALL DGEMM('N','N',NDIM,NDIM,NDIM,ONE,AMB,NDIM,
     *            VL,NDIM,ZERO,VR,NDIM)
C
C     3. DIAGONAILZE SQRT(A-B)*(A+B)*SQRT(A-B) : VR
C
      IF(NDIM.EQ.1) THEN
        EE(1) = VR(1,1)
        VL(1,1) = ONE
      ELSE
        CALL GENTOCAN(VR,SCR,NDIM)
        CALL GLDIAG(NDIM,NDIM,NDIM,SCR,SCR2,EE,VL,IERR,IWRK)
      ENDIF
C
C     * CURRENT VECTOR VL IS SQRT(1/(A-B)))|X+Y>.
C
C     4. VL INTO RIGHT EIGENVECTOR  VR = |X+Y>
C
      CALL DGEMM('N','N',NDIM,NDIM,NDIM,ONE,AMB,NDIM,
     *            VL,NDIM,ZERO,VR,NDIM)
C
C     5. LEFT EIGENVECTOR   VL = = 1/E (A+B)|X+Y>
C
      CALL DGEMM('N','N',NDIM,NDIM,NDIM,ONE,APB,NDIM,
     *            VR,NDIM,ZERO,VL,NDIM)
      DO J=1,NDIM
       DUMJ=SIGN(SQRT(ABS(EE(J))),EE(J))
       DO I=1,NDIM
       VL(I,J)=VL(I,J)/DUMJ
       ENDDO
      ENDDO
C
      RETURN
      END
C*MODULE TDDFT  *DECK RPAECHK
      SUBROUTINE RPAECHK(EE,INDX,NVEC,NDSR,IMAX)
C
C     REMOVE MINUS EIGENVALUES AND
C     PRINT CURRENT EXCITATION ENERGIES
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      LOGICAL GOPARR,DSKWRK,MASWRK
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      DIMENSION EE(NVEC)
      DIMENSION INDX(NVEC)
      DATA ZERO/0.0D+00/
C
C     --- DEFINE INDX(NDIM)
C
      DO I=1,NVEC
      INDX(I)=I
      ENDDO
C
C     --- NUMBER OF MINUS EIGENVALUES : IMAX
C
      IMAX=0
      DO I=1,NVEC
         IF(EE(I).LT.ZERO) IMAX=I
      ENDDO
      IF(IMAX.NE.0) THEN
         IF(MASWRK) WRITE(IW,*)'NUMBER OF MINUS EIGENVALUES=',IMAX
      END IF
C
C     -- PRINT CURRENT EXCITATION ENERGIES
C
      DO ISTAT=1,NDSR
         EE(ISTAT)=SQRT(ABS(EE(ISTAT)))
         YDUM=EE(ISTAT)*27.212D+00
         IF(MASWRK) WRITE(IW,8888) ISTAT, YDUM
      ENDDO
C
      RETURN
 8888 FORMAT(1X,'   STATE NUMBER',I4,2X,'ENERGY =',F12.6,1X,'EV')
        END
C*MODULE TDDFT  *DECK RPAVNORM
      SUBROUTINE RPAVNORM(VR,VL,NVEC)
C
C     NORMALIZE V1 & V2 BY BIORTHOGONALITY CONDITION
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
        DIMENSION VR(NVEC,NVEC)
        DIMENSION VL(NVEC,NVEC)
      DO IVEC=1,NVEC
      VRL=DDOT(NVEC,VR(1,IVEC),1,VL(1,IVEC),1)
      SCAL=SQRT(1.0D+00/VRL)
      CALL DSCAL(NVEC,SCAL,VR(1,IVEC),1)
      CALL DSCAL(NVEC,SCAL,VL(1,IVEC),1)
      ENDDO
      RETURN
        END
C*MODULE TDDFT  *DECK RPAEXPNDV
      SUBROUTINE RPAEXPNDV(VR,VL,VRO,VLO,BR,BL,NDSR,NVEC,NDIM,INDX)
C
C     EXPAND REDUCED VECTORS TO REAL SIZE SPACE
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      LOGICAL DBUG
      DIMENSION VRO(NDIM,*)
      DIMENSION VLO(NDIM,*)
      DIMENSION VR(NVEC,*)
      DIMENSION VL(NVEC,*)
      DIMENSION INDX(*)
      DIMENSION BR(NDIM,*)
      DIMENSION BL(NDIM,*)
      DBUG=.FALSE.
C
      DO ISTAT=1,NDSR
      II=INDX(ISTAT)
       DO I=1,NDIM
       DUMJ_R=0.0D+00
       DUMJ_L=0.0D+00
        DO JVEC=1,NVEC
        DUMJ_R=DUMJ_R+BR(I,JVEC)*VR(JVEC,II)
        DUMJ_L=DUMJ_L+BL(I,JVEC)*VL(JVEC,II)
        ENDDO
       VRO(I,ISTAT)=DUMJ_R
       VLO(I,ISTAT)=DUMJ_L
       ENDDO
      ENDDO
C
      IF(DBUG) THEN
      DO ISTAT=1,NDSR
       DO I=1,NDIM
       DUMR=VRO(I,ISTAT)
       DUML=VLO(I,ISTAT)
       VRL=MAX(DUMR,DUML)
       IF(ABS(VRL).GT.5.0D+00) THEN
       WRITE(6,*)'ISTAT,I,VRO(I,ISTAT)=',ISTAT,I,DUMR
       WRITE(6,*)'ISTAT,I,VLO(I,ISTAT)=',ISTAT,I,DUML
       ENDIF
       ENDDO
      ENDDO
      ENDIF
C
      RETURN
      END
C*MODULE TDDFT  *DECK RPARESVEC
      SUBROUTINE RPARESVEC(Q,W_L,W_R,V_L,V_R,EE,ABD,NDIM,
     * NDSR,ICNVG,ERROR,TOL,IMAX)
C
C      CONSTRUCT RESIDUAL VECTORS AND CHECK CONVERGENCE
C
C     *INPUT W_L: ABP_R
C     *INPUT W_R: ABM_L
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      LOGICAL GOPARR,DSKWRK,MASWRK
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      DIMENSION Q(NDIM,NDSR,2)
      DIMENSION W_L(NDIM,NDSR)
      DIMENSION W_R(NDIM,NDSR)
      DIMENSION V_L(NDIM,NDSR)
      DIMENSION V_R(NDIM,NDSR)
      DIMENSION EE(*)
      DIMENSION ABD(*)
C
C     --- RESIDUAL VECTOR W_L & W_R
C
      DO IVEC=1,NDSR
      CALL DAXPY(NDIM,-EE(IVEC),V_L(1,IVEC),1,W_L(1,IVEC),1)
      CALL DAXPY(NDIM,-EE(IVEC),V_R(1,IVEC),1,W_R(1,IVEC),1)
      ENDDO
C
C     --NORMS OF W
C
      ERROR_MAX=0.0D+00
      DO 10 IVEC=1,NDSR
      IF(IVEC.LE.IMAX) GOTO 10
       ER_L=DDOT(NDIM,W_L(1,IVEC),1,W_L(1,IVEC),1)
       ER_R=DDOT(NDIM,W_R(1,IVEC),1,W_R(1,IVEC),1)
       ERR=MAX(ER_L,ER_R)
        IF(ERR.GT.1.0D+000) THEN
        WRITE(IW,*)'LARGE ERROR DETECTED'
        WRITE(IW,*)'STATE#=',IVEC
        WRITE(IW,*)'ERROR LEFT/RIGHT=',ER_L,'/',ER_R
        ERR=0.0D+00
        ELSE
        IF(ERR.GT.ERROR_MAX) ERROR_MAX=ERR
        ENDIF
C
C      --- GET NEW VECTORS Q
C
       IF(ERR.GT.TOL) THEN
        DO I=1,NDIM
        DUM=1.0D+00/(EE(IVEC)-ABD(I))
        Q(I,IVEC,1)=DUM*W_L(I,IVEC)
        Q(I,IVEC,2)=DUM*W_R(I,IVEC)
        ENDDO
       ELSE
        CALL VCLR(Q(1,IVEC,1),1,NDIM)
        CALL VCLR(Q(1,IVEC,2),1,NDIM)
       ENDIF
C
   10 CONTINUE
C
C     --CHECK CONVERGENCE OF TDDFT CALCULATION
C     ICNVG=0:CONVERGED
C
      IF(MASWRK) THEN
         WRITE(IW,1111) ERROR_MAX,TOL
         CALL FLSHBF(IW)
      ENDIF
C
      ERROR=ERROR_MAX
      IF (ERROR.LE.TOL) THEN
       ICNVG=0
       ELSE
       ICNVG=1
       CONTINUE
      ENDIF
C
      RETURN
 1111 FORMAT(1X,'   MAX ERROR =',1X,1P,E10.3,1X,'/',1P,E10.3)
      END
C*MODULE TDDFT  *DECK RPANEWB
      SUBROUTINE RPANEWB(NDIM,NDSR,MXVEC,BVEC,Q,NOVEC,NVEC,ICK)
C
C      ORTHONORMALIZE Q(NDIM,NDSR*2) AND APPEND TO BVEC
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      LOGICAL DBUG
      LOGICAL GOPARR,DSKWRK,MASWRK
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      DIMENSION Q(NDIM,NDSR*2)
      DIMENSION BVEC(NDIM,*)
      DATA THREB /1.0D-09/
      DBUG=.FALSE.
C
C     SAVE NVEC AS NOVEC
C
      NOVEC=NVEC
      ICK=0
C
C     SCHMIDT ORTHONORMALIZATION
C
      DO 100 K=1,NDSR*2
C
       DO 50 ISTAT=1,NVEC
       BQ=DDOT(NDIM,BVEC(1,ISTAT),1,Q(1,K),1)
   50  CALL DAXPY(NDIM,-BQ,BVEC(1,ISTAT),1,Q(1,K),1)
C
       FNOR =DDOT(NDIM,Q(1,K),1,Q(1,K),1)
       FNOR=SQRT(FNOR)
       IF(FNOR.LT.THREB) GO TO 100
C
C     --- APPEND NEW B VECTOR
C
       NVEC=NVEC+1
       IF(NVEC.GT.MXVEC) THEN
        ICK=1
        RETURN
       ELSE
        CALL DSCAL(NDIM,1.0D+00/FNOR,Q(1,K),1)
        CALL DCOPY(NDIM,Q(1,K),1,BVEC(1,NVEC),1)
       END IF
C
  100 CONTINUE
C
C     ERROR TERMINATION
C
      MS=NVEC-NOVEC
      IF(MS.EQ.0) THEN
       IF(MASWRK) WRITE(IW,*)
     *     'SOMETHING WRONG...THERE ARE NO VECTORS TO BE ADDED'
       CALL ABRT
       STOP
      ENDIF
C
      IF(DBUG) THEN
      WRITE(IW,*)'NUMBER OF VECTORS ',NVEC,'=',NOVEC,'(OLD)+',MS,'(NEW)'
      ENDIF
      RETURN
      END
C*MODULE TDDFT  *DECK IATOGENTD
      SUBROUTINE IATOGENTD(PV,AV,LX,L7,NOCC,IVEC,MXVEC)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DIMENSION PV(L7,MXVEC),AV(LX,LX)
C
      CALL VCLR(AV,1,LX*LX)
      DO I=1,NOCC
        DO J=NOCC+1,LX
          JJ=J-NOCC
          IJ=(JJ-1)*NOCC+I
          AV(I,J)=PV(IJ,IVEC)
        ENDDO
      ENDDO
      RETURN
      END
C*MODULE TDDFT  *DECK MNTOIA
      SUBROUTINE MNTOIA(PAO,PMO,V,SCR,L1,LX,L7,NOCC,IVEC)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      PARAMETER (ZERO=0.0D+00,ONE=1.0D+00)
      DIMENSION PAO(L1,*)
      DIMENSION PMO(L7,*)
      DIMENSION V(L1,LX)
      DIMENSION SCR(NOCC,*)
C
C     -- (M,N) TO (I,N)
C
      CALL DGEMM('T','N',NOCC,L1,L1,ONE,V,L1,PAO,L1,ZERO
     * ,SCR,NOCC)
C
C     -- (I,N) TO (IJ)
C
      DO I=1,NOCC
      DO J=NOCC+1,LX
       JJ=J-NOCC
       IJ=(JJ-1)*NOCC+I
       DUMN=ZERO
       DO N=1,L1
       DUMN=DUMN+SCR(I,N)*V(N,J)
       ENDDO
       PMO(IJ,IVEC)=DUMN
      ENDDO
      ENDDO
C
      RETURN
      END
C*MODULE TDDFT  *DECK ESUM
      SUBROUTINE ESUM(E,PMO,Z,LX,L7,NOCC,IVEC)
C
C     -- ADD (EA-EI)*ZAI
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DIMENSION E(*)
      DIMENSION PMO(L7,*)
      DIMENSION Z(L7,*)
C
      DO I=1,NOCC
         DO J=NOCC+1,LX
            JJ=J-NOCC
            IJ=(JJ-1)*NOCC+I
            PMO(IJ,IVEC)=PMO(IJ,IVEC)+(E(J)-E(I))*Z(IJ,IVEC)
         ENDDO
      ENDDO
      RETURN
      END
C*MODULE TDDFT  *DECK MOCVIR
      SUBROUTINE MOCVIR(VRO,ISTAT,NDIM,MOCC,MVIR,INXOV)
C
C     PICK UP THE MAIN TRANSITION : MOCC -> MVIR
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      LOGICAL GOPARR,DSKWRK,MASWRK
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      DIMENSION VRO(NDIM,*)
      DIMENSION INXOV(2,*)
C
C     SEARCH THE MOST DOMINANT COMPONENT : ITRA
C
      ITRA=0
      DUM=0.0D+00
      DO I=1,NDIM
      DUMI=ABS(VRO(I,ISTAT))
       IF(DUMI.GT.DUM) THEN
       DUM=DUMI
       ITRA=I
       ENDIF
      ENDDO
C
C     SEARCH THE MOST DOMINANT TRANSITION : MOCC->MVIR
C
      MOCC=INXOV(1,ITRA)
      MVIR=INXOV(2,ITRA)
      IF(MASWRK) WRITE(IW,*)'MOCC/MVIR IN TDDFT=',MOCC,MVIR
C
      RETURN
      END
C*MODULE TDDFT  *DECK TDPRNT
      SUBROUTINE TDPRNT(OS,LABSYM,TXYZ,INXOV,NDSR,EE,VRO,L7,IMAX)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      LOGICAL GOPARR,DSKWRK,MASWRK,TRIPLET,CHECK
      LOGICAL ABEL,ABELPT,SG1T
C
      PARAMETER (MXAO=8192, MXATM=2000, MXRT=100)
C
      COMMON /ENRGYS/ ENUCR,EELCT,ETOT,SZ,SZZ,ECORE,ESCF,EERD,E1,E2,
     *                VEN,VEE,EPOT,EKIN,ESTATE(MXRT),STATN,EDFT(2)
      COMMON /FUNCT / ENERGY,EG(3*MXATM)
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /INFOTD/ TDM(3,MXRT),CNVTOL,TRIPLET,MOCC(MXRT),MVIR(MXRT),
     *                NRADT,NTHET,NPHIT,NSTAT,NTRIAL,MAXVEC,NTHST,
     *                IRECTD,NLEBT,SG1T
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /SOSYM/  EULANG(4,48),GAM(48,48),IRMON(MXAO)
      COMMON /SYMMUL/ NIJREP(14,14),IJREP(2,14,14,14)
      COMMON /SYMREP/ IRPNAM(14),IPA(14),LAMBDA(14),LAMBD0(14),
     &                IADDR1(14),IADDR2(14),IADDR3(14)
C
      DIMENSION VRO(L7,*),INXOV(2,*)
      DIMENSION OS(NDSR),LABSYM(NDSR),TXYZ(3,NDSR),EE(NDSR)
C
      PARAMETER (TOEV=27.21138386D+00, THRESHOLD=3.0D-02)
      DATA IQMARK/4H????/
C-MWS DATA IBLANK/4H    /
C
      IF(.NOT.TRIPLET) THEN
       IF(MASWRK) WRITE(IW,8886)
      ELSE
       IF(MASWRK) WRITE(IW,8885)
      ENDIF
C
      ABEL = ABELPT()
      IF(.NOT.ABEL  .AND.  MASWRK) WRITE(IW,9050)
 9050 FORMAT(/1X,'THIS IS A NON-ABELIAN POINT GROUP, AS A RESULT,'/
     *       1X,'SOME STATE SYMMETRY LABELS MAY NOT BE CORRECTLY',
     *          ' PRINTED BELOW'/)
C
      DO 10 ISTAT=1,NDSR
      IF(ISTAT.LE.IMAX) GOTO 10
        OSI=OS(ISTAT)
        YDUM = TOEV*EE(ISTAT)
C
C      -- SYMMETRY OF EXCITED STATE
C     MIKE DOES NOT BELIEVE THAT THIS CODE ALWAYS PRINTS THE CORRECT
C     STATE SYMMETRY, AND HAS AN N2 EXAMPLE TO BACK THIS UP.  HE
C     ORIGINALLY USED THE LINES THAT YOU SEE COMMENTED OUT, PREFERRING
C     TO PRINT BLANKS FOR ALL STATES, RATHER THAN MISTAKES FOR SOME.
C
C-MWS   IF(ABEL) THEN
           CHECK = .FALSE.
           NAMSAV = -179
           DO IJ=1,L7
              I=INXOV(1,IJ)
              J=INXOV(2,IJ)
              XDUM=VRO(IJ,ISTAT)
              IF (ABS(XDUM).GT.THRESHOLD) THEN
                 IRP = IRMON(I)
                 JRP = IRMON(J)
                 IINAM = IJREP(2,1,IRP,JRP)
                 IF(CHECK) THEN
                    IF(IINAM.NE.NAMSAV) THEN
                       WRITE(IW,8892)
                       LABSYM(ISTAT) = IQMARK
                       GO TO 50
                    END IF
                 ELSE
                    CHECK = .TRUE.
                    LABSYM(ISTAT) = IRPNAM(IINAM)
                    NAMSAV = IINAM
                 END IF
              END IF
           ENDDO
C-MWS   ELSE
C-MWS      LABSYM(ISTAT)=IBLANK
C-MWS   END IF
C
   50   CONTINUE
        IF(MASWRK) THEN
           WRITE(IW,8888) ISTAT, YDUM
           WRITE(IW,8889) OSI
           WRITE(IW,8890) LABSYM(ISTAT)
           WRITE(IW,8887)
           WRITE(IW,8884)
        ENDIF
C
        DO IJ=1,L7
           I=INXOV(1,IJ)
           J=INXOV(2,IJ)
           XDUM=VRO(IJ,ISTAT)
           IF (ABS(XDUM).GT.THRESHOLD) THEN
               IF(MASWRK) WRITE(IW,8891) IJ,XDUM,I,J
           ENDIF
        ENDDO
   10 CONTINUE
C
C         PRINT SUMMARY TABLE
C
      IF(MASWRK) THEN
         IF(.NOT.ABEL) WRITE(IW,9050)
         WRITE(IW,9010)
         WRITE(IW,9020) 0,IRPNAM(1),ESCF,0.0D+00
         DO IST=1,NDSR
            WRITE(IW,9020) IST,LABSYM(IST),ESCF+EE(IST),TOEV*EE(IST),
     *                     (TXYZ(II,IST),II=1,3),OS(IST)
         ENDDO
         WRITE(IW,*) ' '
      END IF
C
C                ----- SAVE ENERGY RESULTS -----
C     CAN'T UPDATE K.E., ET CETERA WITHOUT A DENSITY MATRIX
C
      ENERGY    = ESCF + EE(NTHST)
      ETOT      = ENERGY
C     ESTATE(1) = ENERGY
      STATN     = MIN(1+NDSR,MXRT)
      DO I=1,MIN(NDSR,MXRT-1)
        ESTATE(I+1) = ESCF + EE(NTHST)
      ENDDO
      NSTMAX=NDSR
      DO I=1,NSTMAX
      ESTATE(I)       = EE(I)
      TDM(1,I)    = TXYZ(1,I)
      TDM(2,I)    = TXYZ(2,I)
      TDM(3,I)    = TXYZ(3,I)
      ENDDO
      CALL DAWRIT(IDAF,IODA,ENUCR,MXRT+15,2,0)
      IF(MASWRK) WRITE(IW,9030) NTHST,ENERGY
C
      RETURN
C
 8892 FORMAT(/1X,'WARNING: THIS STATE HAS BROKEN SYMMETRY, CHECK MOS')
 8890 FORMAT(1X,'SYMMETRY OF STATE   =',4X,A4)
 8889 FORMAT(1X,'OSCILLATOR STRENGTH =',F12.6)
 8888 FORMAT(/1X,'STATE #',I4,2X,'ENERGY =',F12.6,1X,'EV')
 8887 FORMAT(8X,'DRF',4X,'COEF',7X,'OCC',6X,'VIR')
 8886 FORMAT(/10X,19(1H-)/
     *        10X,'SINGLET EXCITATIONS'/
     *        10X,19(1H-))
 8885 FORMAT(/10X,19(1H-)/
     *        10X,'TRIPLET EXCITATIONS'/
     *        10X,19(1H-))
 8884 FORMAT(8X,3(1H-),4X,4(1H-),7X,3(1H-),6X,3(1H-))
 8891 FORMAT(7X,I4,1X,F9.6,3X,I4,2X,'->',1X,I4,2X)
 9010 FORMAT(/26X,'SUMMARY OF TDDFT RESULTS'//
     *         3X,'STATE',13X,'ENERGY',5X,'EXCITATION',
     *         6X,'TRANSITION DIPOLE',4X,'OSCILLATOR'/
     *        20X,'HARTREE',10X,'EV',9X,'X',7X,'Y',7X,'Z',5X,'STRENGTH')
 9020 FORMAT(1X,I3,2X,A4,1X,F20.10,1X,F8.3,3X,3F8.4,1X,F8.3)
 9030 FORMAT(1X,'SELECTING EXCITED STATE IROOT=',I4,' AT E=',F20.10/
     *       1X,'AS THE STATE OF INTEREST.'/)
      END
C*MODULE TDDFT  *DECK TDPCMDEN
      SUBROUTINE TDPCMDEN(DEN,V,L1)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      COMMON /IJPAIR/ IA(8192)
      DIMENSION DEN(*)
      DIMENSION V(L1,*)
C
      DO M=1,L1
        DO N=1,M
          MN=IA(M)+N
          DEN(MN)=V(M,N)+V(N,M)
        ENDDO
      ENDDO
C
      RETURN
      END
C*MODULE TDDFT  *DECK TDOSCALC
      SUBROUTINE TDOSCALC(OS,TXYZ,TDEN,SCR,VRO,VLO,EE,V,AX,AY,AZ,
     *                    L1,LX,L2,L7,NDSR,INXOV)
C
C     --- CALCULATE OSCILLATOR STRENGTH
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER (ZERO=0.0D+00,ONE=1.0D+00)
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /XYZPRP/ XP,YP,ZP,
     *                DMX,DMY,DMZ,
     *                QXX,QYY,QZZ,QXY,QXZ,QYZ,
     *                QMXX,QMYY,QMZZ,QMXY,QMXZ,QMYZ,
     *                OXXX,OXXY,OXXZ,OXYY,OYYY,OYYZ,
     *                OXZZ,OYZZ,OZZZ,OXYZ,
     *                OMXXX,OMXXY,OMXXZ,OMXYY,OMYYY,
     *                OMYYZ,OMXZZ,OMYZZ,OMZZZ,OMXYZ
      DIMENSION VRO(L7,NDSR)
      DIMENSION VLO(L7,NDSR)
      DIMENSION EE(*)
      DIMENSION V(L1,LX)
      DIMENSION TDEN(L1,L1),SCR(L1,L1)
      DIMENSION AX(L2),AY(L2),AZ(L2)
      DIMENSION TXYZ(3,NDSR)
      DIMENSION OS(NDSR)
      DIMENSION INXOV(2,*)
C
C     -- GET DIPOLE INTEGRALS AT CENTER OF MASS
C
      CALL CALCOM(XP,YP,ZP)
      CALL DIPINT(XP,YP,ZP,.FALSE.)
      CALL DAREAD(IDAF,IODA,AX,L2,95,0)
      CALL DAREAD(IDAF,IODA,AY,L2,96,0)
      CALL DAREAD(IDAF,IODA,AZ,L2,97,0)
C
      DO IST=1,NDSR
C
C     -- TRANSITION DENSITY
C
       CALL TRAD(TDEN,VRO,VLO,LX,L7,INXOV,NDSR,IST)
C
C     -- TRANSFORM TO AO-BASIS
C
       CALL DGEMM('N','N',L1,LX,LX,ONE,V,L1,TDEN,LX,ZERO,SCR,L1)
       CALL DGEMM('N','T',L1,L1,LX,ONE,SCR,L1,V,L1,ZERO,TDEN,L1)
C
C     -- TRANSITION DIPOLE
C
       CALL TRADIP(TXYZ,TDEN,AX,AY,AZ,L1,L2,NDSR,IST)
C
      ENDDO
C
C     -- OSCILLATOR STRENGTH, AND PRINTING
C
      CALL OSCALC(OS,TXYZ,EE,NDSR)
C
      RETURN
      END
C*MODULE TDDFT  *DECK TRAD
      SUBROUTINE TRAD(TDEN,VRO,VLO,LX,L7,INXOV,NDSR,IST)
C
C     -- TRANSITION DENSITY
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION TDEN(LX,LX)
      DIMENSION VRO(L7,NDSR)
      DIMENSION VLO(L7,NDSR)
      DIMENSION INXOV(2,*)
C     --CLEAR
      CALL VCLR(TDEN,1,LX*LX)
C
      DO IPR=1,L7
      I=INXOV(1,IPR)
      J=INXOV(2,IPR)
C     -- IA COMPONENTS
      TDEN(I,J)=VRO(IPR,IST)+VLO(IPR,IST)
C     -- AI COMPONENTS
      TDEN(J,I)=VRO(IPR,IST)-VLO(IPR,IST)
      ENDDO
C
      RETURN
      END
C*MODULE TDDFT  *DECK TRADIP
      SUBROUTINE TRADIP(TXYZ,TDEN,AX,AY,AZ,L1,L2,NDSR,IST)
C
C     --  TRANSITION DIPOLE MOMENT
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      COMMON /IJPAIR/ IA(8192)
      DIMENSION TDEN(L1,L1)
      DIMENSION TXYZ(3,NDSR)
      DIMENSION AX(L2),AY(L2),AZ(L2)
      DATA ZERO/0.0D+00/
C     --X
      DUMN=ZERO
      DO M=1,L1
      DO N=1,L1
      MN=IA(MAX(M,N))+MIN(M,N)
      DUM=TDEN(M,N)*AX(MN)
      DUMN=DUMN+DUM
      ENDDO
      ENDDO
      TXYZ(1,IST)=-DUMN/SQRT(2.0D+00)
C     --Y
      DUMN=ZERO
      DO M=1,L1
      DO N=1,L1
      MN=IA(MAX(M,N))+MIN(M,N)
      DUM=TDEN(M,N)*AY(MN)
      DUMN=DUMN+DUM
      ENDDO
      ENDDO
      TXYZ(2,IST)=-DUMN/SQRT(2.0D+00)
C     --Z
      DUMN=ZERO
      DO M=1,L1
      DO N=1,L1
      MN=IA(MAX(M,N))+MIN(M,N)
      DUM=TDEN(M,N)*AZ(MN)
      DUMN=DUMN+DUM
      ENDDO
      ENDDO
      TXYZ(3,IST)=-DUMN/SQRT(2.0D+00)
C
      RETURN
      END
C*MODULE TDDFT  *DECK OSCALC
      SUBROUTINE OSCALC(OS,TXYZ,EE,NDSR)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      DIMENSION OS(NDSR),TXYZ(3,NDSR),EE(NDSR)
C
      DATA ZERO /0.0D+00/
      DATA TWO  /2.0D+00/
      DATA THREE/3.0D+00/
C
      DO IST=1,NDSR
         OS(IST)=ZERO
         DO K=1,3
            DUM=TXYZ(K,IST)**2*EE(IST)*TWO/THREE
            OS(IST)=OS(IST)+DUM
         ENDDO
      ENDDO
C
      RETURN
      END
