C 20 AUG 07 - MWS - MEXINP: KILL ATTEMPT TO USE EFP PARTICLES
C 22 DEC 06 - DGF - SYNCHRONISE MP2PAR
C  6 NOV 06 - MWS - ADJUST WAVEFUNCTION COMMON BLOCK
C 13 MAR 06 - NM  - NEW MODULE TO SEARCH MINIMUM ENERGY CROSSING POINT
C
C*MODULE MEXING  *DECK MEXNGX
C======================
      SUBROUTINE MEXNGX
C======================
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      LOGICAL GOPARR,DSKWRK,MASWRK
C
      PARAMETER (MXATM=2000)
C
      COMMON /FMCOM / X(1)
      COMMON /INFOA / NAT,ICH,MUL,NUM,NQMT,NE,NA,NB,
     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
      COMMON /IOFILE/ IR,IW,IP,IJK,IJKT,IDAF,NAV,IODA(950)
      COMMON /MEXPAR/ SCF1,SCF2,TGMAX,TDE,TDXMAX,TDXRMS,TGRMS,STPSZ,
     *                MULT1,MULT2,NSTEP,NRDMOS,NMOS1,NMOS2,NPRT
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
C
C          ----- MINIMUM ENERGY CROSSING POINT SEARCH -----
C       LOCATE "DIABATIC TRANSITION STATE" BETWEEN TWO SURFACES,
C        WHICH MUST DIFFER IN SPIN OR SPACE SYMMETRY (OR BOTH).
C
C     AUTHORS:    JEREMEY HARVEY         NIKITA MATSUNAGA,
C                 SCHOOL OF CHEMISTRY    DEPARTMENT OF CHEMISTRY,
C                 UNIVERSITY OF BRISTOL  LONG ISLAND UNIVERSITY.
C                 CANTOCK'S CLOSE        1 UNIVERSITY PLAZA
C                 BRISTOL BS8 1TS, UK    BROOKLYN NY 11201
C
C     STANDALONE LATE '90S, 1ST EMBEDDING=10/2001, MCSCF ADDED 2005,
C     CODE INCORPORATED INTO PRODUCTION VERSION OF GAMESS IN 3/2006.
C
C     METHOD IS NEWTON-RAPHSON SEARCH WITH LAGRANGE MULTIPLIER
C     SIDE CONDITION TO CONSTRAIN TO EQUAL STATE ENERGIES.  SEE,
C     FOR EXAMPLE, A.FARAZDEL, M.DUPUIS J.COMPUT.CHEM. 12,276(1991)
C     BUT MANY OTHERS HAVE PUBLISHED SIMILAR ALGORITHMS, STARTING
C     WITH N.KOGA, K.MOROKUMA  CHEM.PHYS.LETT 119, 371(1985), SEE
C     J.N.HARVEY, M.ASCHI  PHYS.CHEM.CHEM.PHYS. 1, 5555(1999), AND
C     REFERENCES THERE IN.
C
      IF(MASWRK) WRITE(IW,9000)
C
C     ----- IF GUESS MOS ARE READ IN, DO IT HERE -----
C
      IF(NRDMOS.GT.0) THEN
         NBF = NUM
         CALL ONEEI
         CALL GESFMX(NBF,NE,MULT1,MULT2,SCF1,SCF2,NRDMOS,NMOS1,NMOS2)
      END IF
C
C     ----- GET MEMORY ALLOCATED FOR MEX RUN -----
C
      NBF   = NUM
      NATOM = NAT
      NX    = 3*NAT
C
      CALL VALFM(LOADFM)
      LHI1  = LOADFM + 1
      LHI2  = LHI1   + NX*NX
      LPRG  = LHI2   + NX*NX
      LPPG  = LPRG   + NX
      LGA2  = LPPG   + NX
      LGB2  = LGA2   + NX
      LX1   = LGB2   + NX
      LX2   = LX1    + NX
      LX3   = LX2    + NX
      LG1   = LX3    + NX
      LG2   = LG1    + NX
      LDG   = LG2    + NX
      LHDG  = LDG    + NX
      LCH   = LHDG   + NX
      LDX   = LCH    + NX
      LW    = LDX    + NX
      LORB1 = LW     + NX
      LORB2 = LORB1  + NBF*NBF
      LAST  = LORB2  + NBF*NBF
      NEED  = LAST - LOADFM - 1
C
      IF(MASWRK) WRITE(IW,9100) NEED
      CALL GETFM(NEED)
C
C     ----- GO SEARCH! ------
C
      CALL MEXDRV(NATOM,NX,X(LHI1),X(LHI2),X(LPRG),X(LPPG),
     *            X(LGA2),X(LGB2),X(LX1),X(LX2),X(LX3),
     *            X(LG1),X(LG2),MULT1,MULT2,SCF1,SCF2,
     *            X(LDG),X(LHDG),X(LCH),X(LDX),X(LW),NBF,X(LORB1),
     *            X(LORB2),NSTEP,TGMAX,TDE,TDXMAX,TDXRMS,
     *            TGRMS,NRDMOS,STPSZ,NPRT)
C
      CALL RETFM(NEED)
      IF(MASWRK) WRITE(IW,9200)
      CALL TIMIT(1)
      RETURN
C
 9000 FORMAT(/18X,36(1H-)/18X,'MINIMUM ENERGY CROSSING POINT SEARCH',
     *       /18X,36(1H-))
 9100 FORMAT(/5X,'MEMORY NEEDED FOR MEX POINT SEARCH =',I10,' WORDS')
 9200 FORMAT(1X,'.... DONE WITH MINIMUM ENERGY CROSSING',
     *          ' POINT SEARCH ....')
      END
C
C*MODULE MEXING  *DECK MEXDRV
C==================================================================
      SUBROUTINE MEXDRV(NATM,NX,HI1,HI2,PPG,PERPG,GA2,GB2,
     *                  X1,X2,X3,G1,G2,MULT1,MULT2,SCF1,SCF2,
     *                  DELG,HDELG,CHGEX,DELX,W,NBF,VEC1,VEC2,
     *                  NSTEP,TGMAX,TDE,TDXMAX,TDXRMS,TGRMS,NRDMOS,
     *                  STPSZ,NPRT)
C==================================================================
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      DOUBLE PRECISION MTHDMC
C
      LOGICAL MP2,GOPARR,DSKWRK,MASWRK,SKIPNO1,SKIPNO2,OK,LINEAR,
     *        CANONC,FCORE,FORSMC,NOCI,EKT,LINSER,GOTEG,GPSAVE
C
      DIMENSION HI1(NX,NX),HI2(NX,NX),PPG(NX),PERPG(NX),
     *          GA2(NX),X1(NX),X2(NX),X3(NX),G1(NX),G2(NX),GB2(NX),
     *          DELG(NX),HDELG(NX),CHGEX(NX),DELX(NX),W(NX),
     *          VEC1(NBF*NBF),VEC2(NBF*NBF)
C
      PARAMETER (MXATM=2000, MXAO=8192,ZERO=0.0D+00,ONE=1.0D+00)
      PARAMETER (TOANGS=0.52917724924D+00)
      PARAMETER (MXNORO=250)
C
      COMMON /ECP2  / CLP(400),ZLP(400),NLP(400),KFIRST(MXATM,6),
     *                KLAST(MXATM,6),LMAX(MXATM),LPSKIP(MXATM),
     *                IZCORE(MXATM)
      COMMON /FUNCT / E,EG(3*MXATM)
      COMMON /INFOA / NAT,ICH,MUL,NUM,NQMT,NE,NA,NB,
     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
      COMMON /IOFILE/ IR,IW,IP,IJK,IJKT,IDAF,NAV,IODA(950)
      COMMON /MCINP / MTHDMC,CISTEP,ACURCY,ENGTOL,DAMP,MICIT,NWRDMC,
     *                NORBMC,NOROT(2,MXNORO),MOFRZ(15),NPFLG(10),
     *                NOFO,MCFMO,CANONC,FCORE,FORSMC,NOCI,EKT,LINSER
      COMMON /MEXOPT/ MEXSKPGES,MEXSTATE
      COMMON /MP2PAR/ OSPT,CODEMP,SCSPT,TOL,METHOD,NWDMP2,MEMPRI,MPPROP,
     *                NACORE,NBCORE,NOA,NOB,NORB,NBFMP2,NOMIT,MOCPHF,
     *                MAXITC
      COMMON /OUTPUT/ NPRINT,ITOL,ICUT,NORMF,NORMP,NOPK
      COMMON /RUNLAB/ TITLE(10),ANAM(MXATM),BNAM(MXATM),BFLAB(MXAO)
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /RUNOPT/ RUNTYP,EXETYP,NEVALS,NGLEVL,NHLEVL
      COMMON /WFNOPT/ SCFTYP,VBTYP,DFTYPE,TDDFTYP,CITYP,CCTYP,
     *                MPLEVL,MPCTYP
      COMMON /ZMAT  / NZMAT,NZVAR,NVAR,NSYMC,LINEAR
C
      DATA RHF,UHF,ROHF,AMCSCF
     *    /8HRHF     ,8HUHF     ,8HROHF    ,8HMCSCF   /
      DATA ANONE,UMP,ZAPT
     *    /8HNONE    ,8HUMP     ,8HZAPT    /
      DATA CHECK /8HCHECK   /
      DATA GUGA,ALDET/8HGUGA    ,8HALDET   /
      DATA GRWRD1,GRWRD2/8H $GRAD1 ,8H $GRAD2 /
C
C     ----- STEPPER FOR MINIMUM ENERGY CROSSING SEARCH -----
C
C     E1,   E2  = ENERGIES OF FIRST AND SECOND STATES
C     GA2,  GB2 = GRADIENTS OF FIRST AND SECOND STATES
C     G1,   G2  = EFFECTIVE GRADIENT OF PREVIOUS STEP AND CURRENT STEP
C     PERPG,PPG = PERPENDICULAR AND PARALLEL GRADIENTS
C
C     NPRT = 0 DEFAULT PRINTING -- NO ORBITALS AFTER THE FIRST STEP
C     NPRT = 1 PRINT OUT ORBITALS IN EACH ITERATION
C
      GPSAVE = GOPARR
C
      IF(MPLEVL.EQ.2) MP2=.TRUE.
C
C     SO FAR, ONLY CISTEP=GUGA OR ALDET IS PERMITTED IN CASE OF MCSCF
C
      IF(SCF1.EQ.AMCSCF) THEN
         OK=.FALSE.
         IF(CISTEP.EQ.GUGA)  OK=.TRUE.
         IF(CISTEP.EQ.ALDET) OK=.TRUE.
         IF(.NOT.OK) THEN
            IF(MASWRK) WRITE(IW,9010) CISTEP
            CALL ABRT
         END IF
      END IF
C
      MEXSKPGES = 0
      SKIPNO1 = .FALSE.
      SKIPNO2 = .FALSE.
C
C     ----- COPY CURRENT GEOMETRY INTO ARRAY X2 -----
C
      K = 0
      DO I=1,NAT
         DO J=1,3
            K = K+1
            X2(K) = C(J,I)
         END DO
      END DO
C
C     ----- SET UP FOR TWO STATE CALCULATIONS -----
C
      NBFMP2 = NUM
      NORB = NUM
      MSPN = MULT1 - 1
      NB1 = (NE - MSPN)/2
      NA1 = NB1 + MSPN
      MSPN = MULT2 -1
      NB2 = (NE - MSPN)/2
      NA2 = NB2 + MSPN
C
C     ----- LOOP OVER NUMBER OF STEPS -----
C
      NBF2 = NBF*NBF
C
      DO 100 ISTEP = 0, NSTEP
C
      IF(NPRT.EQ.0) THEN
         IF(ISTEP.EQ.0) NPRINT= 7
         IF(ISTEP.GE.1) NPRINT=-5
      END IF
      IF(NPRT.EQ.1) NPRINT= 7
C
C     ----- 1ST STATE PREP FOR ENERGY AND GRADIENT -----
C
      MEXSTATE = 1
      SCFTYP = SCF1
      MUL = MULT1
      NA  = NA1
      NB  = NB1
      IF(MP2) THEN
        NOA = NA1
        NOB = NB1
        IF(SCFTYP.EQ.RHF)  OSPT=ANONE
        IF(SCFTYP.EQ.UHF)  OSPT=UMP
        IF(SCFTYP.EQ.ROHF) OSPT=ZAPT
      END IF
C
      IF(MASWRK) WRITE(IW,9020) 1,1,MUL
      CALL PRATM(IW,2)
      IF(NZMAT.GT.0) THEN
         CALL BANDBI
         CALL PZANDG(DUMMY,0)
      END IF
      IF(MASWRK) WRITE(IP,8000) ISTEP,' FIRST'
      CALL PRATM(IP,2)
C
      IF(MUL.NE.1) CALL SETCAN
C
C     ----- 1ST STATE ORBITALS -----
C           TONS OF INITIAL ORBITALS/GRADIENTS OPTIONS
C
      IF(ISTEP.EQ.0) THEN
         IF(NRDMOS.EQ.1 .OR. NRDMOS.EQ.3) THEN
            NEVALS = 0
            MEXSKPGES = 1
            CALL ORBMAN(MEXSTATE,VEC1,VEC2,NBF2,SCF1)
         ELSE
            MEXSKPGES = 0
            NEVALS = 0
         END IF
         CALL EGIN(EG,3*NAT,GOTEG,GRWRD1)
         IF (GOTEG) THEN
            SKIPNO1=.TRUE.
            IF(MASWRK) WRITE(IW,9030) MEXSTATE
         END IF
      END IF
      IF(ISTEP.EQ.1) THEN
         IF(SKIPNO1) THEN
            IF(NRDMOS.EQ.1.OR .NRDMOS.EQ.3) THEN
               NEVALS = 0
               MEXSKPGES = 1
               CALL ORBMAN(MEXSTATE,VEC1,VEC2,NBF2,SCF1)
               SKIPNO1 = .FALSE.
            ELSE
               MEXSKPGES = 0
               NEVALS = 0
            END IF
         ELSE
            NEVALS = 1
            MEXSKPGES = 1
            CALL ORBMAN(MEXSTATE,VEC1,VEC2,NBF2,SCF1)
         END IF
      END IF
C
      IF(ISTEP.GE.2) THEN
         NEVALS = 1
         MEXSKPGES = 1
         CALL ORBMAN(MEXSTATE,VEC1,VEC2,NBF2,SCF1)
      END IF
C
C     ----- MCSCF -----
C
      IF(ISTEP.EQ.0) THEN
         DO I=1,10
            NPFLG(I) = 0
         END DO
      END IF
C
C     ------ SKIP IF ALREADY HAVE GRADIENT -----
C
      IF(ISTEP.EQ.0) THEN
         IF(GOTEG) GO TO 2000
      END IF
C
C     ------ CRUNCH ENERGY ------
C
      CALL ENERGX
C
      IF(E.EQ.ZERO .AND. EXETYP.NE.CHECK) THEN
         IF(MASWRK) WRITE(IW,9040) 1
         CALL ABRT
      END IF
C
C     ----- 1ST STATE GRADIENT -----
C
      CALL HFGRAD
      GOPARR = GPSAVE
C
C     ----- BEFORE NEXT STEP STORE OLD ORBITALS -----
C
      NRECA = 406
      NRECB = 407
      CALL DAREAD(IDAF,IODA,VEC1,NBF2,15,0)
      CALL DAWRIT(IDAF,IODA,VEC1,NBF2,NRECA,0)
      IF(SCFTYP.EQ.UHF) THEN
         CALL DAREAD(IDAF,IODA,VEC2,NBF2,19,0)
         CALL DAWRIT(IDAF,IODA,VEC2,NBF2,NRECB,0)
      END IF
      IF(SCFTYP.EQ.AMCSCF) THEN
         CALL SEQCLO(11,'DELETE')
         CALL SEQCLO(12,'DELETE')
         CALL SEQCLO(14,'DELETE')
      END IF
C
 2000 CONTINUE
      IF(MASWRK) WRITE (IW,9050) ' FIRST'
      CALL EGPUN(EG,NAT,GRWRD1)
      DO IAT = 1,NAT
         N = 3*(IAT-1)
         ZNUC = ZAN(IAT) + IZCORE(IAT)
         IF(MASWRK) THEN
           WRITE(IW,9060) IAT,ANAM(IAT),BNAM(IAT),ZNUC,(EG(N+I),I=1,3)
         END IF
      END DO
C
      DO I=1,NX
         GA2(I) = EG(I)
      END DO
      E1 = E
C
C    ----- 2ND STATE PREP FOR ENERGY AND GRADIENT -----
C
      MEXSTATE = 2
      SCFTYP = SCF2
      MUL    = MULT2
      NA     = NA2
      NB     = NB2
      IF(MP2) THEN
        NOA = NA2
        NOB = NB2
        IF(SCFTYP.EQ.RHF)  OSPT=ANONE
        IF(SCFTYP.EQ.UHF)  OSPT=UMP
        IF(SCFTYP.EQ.ROHF) OSPT=ZAPT
      END IF
C
      IF(MASWRK) WRITE(IW,9020) 2,2,MUL
      IF(MASWRK) WRITE(IP,8000) ISTEP,'SECOND'
C
      IF(MUL.NE.1) CALL SETCAN
C
C     ----- 2ND STATE ENERGY -----
C
      IF(ISTEP.EQ.0) THEN
         IF(NRDMOS.EQ.2 .OR. NRDMOS.EQ.3) THEN
            NEVALS = 0
            MEXSKPGES = 1
            CALL ORBMAN(MEXSTATE,VEC1,VEC2,NBF2,SCF2)
         ELSE
            MEXSKPGES = 0
            NEVALS = 0
         END IF
         CALL EGIN(EG,3*NAT,GOTEG,GRWRD2)
         IF (GOTEG) THEN
            SKIPNO2=.TRUE.
            IF(MASWRK) WRITE(IW,9030) MEXSTATE
         END IF
      END IF
      IF(ISTEP.EQ.1) THEN
         IF(SKIPNO2) THEN
            IF(NRDMOS.EQ.2.OR .NRDMOS.EQ.3) THEN
               NEVALS = 0
               MEXSKPGES = 1
               CALL ORBMAN(MEXSTATE,VEC1,VEC2,NBF2,SCF2)
               SKIPNO2 = .FALSE.
            ELSE
               MEXSKPGES = 0
               NEVALS = 0
            END IF
         ELSE
            NEVALS = 1
            MEXSKPGES = 1
            CALL ORBMAN(MEXSTATE,VEC1,VEC2,NBF2,SCF2)
         END IF
      END IF
      IF(ISTEP.GE.2) THEN
         NEVALS = 1
         MEXSKPGES = 1
         CALL ORBMAN(MEXSTATE,VEC1,VEC2,NBF2,SCF2)
      END IF
C
C     ----- MCSCF -----
C
      IF(ISTEP.EQ.0) THEN
         DO I=1,10
            NPFLG(I) = 0
         END DO
      END IF
C
C     ----- IF GRADIENT IS ALREADY KNOWN SKIP-----
C
      IF(ISTEP.EQ.0) THEN
         IF(GOTEG) GO TO 3000
      END IF
C
      CALL ENERGX
C
      IF(E.EQ.ZERO .AND. EXETYP.NE.CHECK) THEN
         IF(MASWRK) WRITE(IW,9040) 2
         CALL ABRT
      END IF
C
C     ----- 2ND STATE GRADIENT -----
C
      CALL HFGRAD
      GOPARR = GPSAVE
C
C     ----- BEFORE GOING TO NEXT STEP, STORE ORBITALS -----
C
      NRECA = 408
      NRECB = 409
      CALL DAREAD(IDAF,IODA,VEC1,NBF2,15,0)
      CALL DAWRIT(IDAF,IODA,VEC1,NBF2,NRECA,0)
      IF(SCFTYP.EQ.UHF) THEN
         CALL DAREAD(IDAF,IODA,VEC2,NBF2,19,0)
         CALL DAWRIT(IDAF,IODA,VEC2,NBF2,NRECB,0)
      END IF
      IF(SCFTYP.EQ.AMCSCF) THEN
         CALL SEQCLO(11,'DELETE')
         CALL SEQCLO(12,'DELETE')
         CALL SEQCLO(14,'DELETE')
      END IF
C
 3000 CONTINUE
      IF(MASWRK) WRITE(IW,9050) 'SECOND'
      CALL EGPUN(EG,NAT,GRWRD2)
      DO IAT = 1,NAT
         N = 3*(IAT-1)
         ZNUC = ZAN(IAT) + IZCORE(IAT)
         IF(MASWRK) THEN
            WRITE(IW,9060) IAT,ANAM(IAT),BNAM(IAT),ZNUC,(EG(N+I),I=1,3)
         END IF
      END DO
C
      DO I=1,NX
         GB2(I) = EG(I)
      END DO
      E2 = E
C
C     ----- CONSTRUCTION OF EFFECTIVE GRADIENTS -----
C
      IF(ISTEP.EQ.0) THEN
         CALL VCLR(G1,   1,NX)
         CALL VCLR(G2,   1,NX)
         CALL VCLR(PERPG,1,NX)
         CALL VCLR(PPG,  1,NX)
         CALL VCLR(X1,   1,NX)
         CALL VCLR(X3,   1,NX)
         CALL VCLR(HI1,  1,NX*NX)
         CALL VCLR(HI2,  1,NX*NX)
         DO I=1,NX
            HI1(I,I) = ONE
         END DO
      END IF
      IF(MASWRK) WRITE(IW,9070) ISTEP
      CALL EFFGRD(NX,E1,E2,GA2,GB2,PERPG,PPG,G2)
C
C     ----- UPDATE GEOMETRY (BFGS STEP) -----
C
      CALL UPDTX(NX,ISTEP,STPSZ,X1,X2,X3,G1,G2,HI1,HI2,
     *           DELG,HDELG,CHGEX,DELX,W)
C
C     ----- CONVERGENCE TEST -----
C
      CALL TSTCNV(TGMAX,TDE,TDXMAX,TDXRMS,TGRMS,
     *            NX,NATM,E1,E2,X2,X3,PERPG,PPG,G2,CHGEX,ICONV,
     *            ISTEP)
C
C     ----- PRINT NEW GEOMETRY AND UPDATE OTHER ARRAYS  -----
C
      DO I = 1, NAT
         C(1,I) = X3(3*(I-1)+1)
         C(2,I) = X3(3*(I-1)+2)
         C(3,I) = X3(3*(I-1)+3)
      END DO
      DO I=1,NX
         X1(I) = X2(I)
         X2(I) = X3(I)
         G1(I) = G2(I)
         DO J=1,NX
            HI1(I,J) = HI2(I,J)
         END DO
      END DO
C
      IF(ICONV .EQ. 1) THEN
C                CARTESIAN COORDINATES
         IF(MASWRK) THEN
            WRITE(IW,9100) ISTEP
            WRITE(IW,9110)
            DO I=1,NAT
               ZNUC = ZAN(I) + IZCORE(I)
               X = TOANGS*C(1,I)
               Y = TOANGS*C(2,I)
               Z = TOANGS*C(3,I)
               WRITE(IW,9120) ANAM(I),BNAM(I),ZNUC,X,Y,Z
            END DO
         END IF
C                 INTERNAL COORDINATES AND DISTANCE MATRIX
         IF(NZMAT.GT.0) CALL PZANDG(DUMMY,0)
         CALL INTR
C
         IF(MASWRK) WRITE(IW,9200)
C
         GO TO 1000
      END IF
C
      IF(EXETYP.EQ.CHECK) GO TO 1000
  100 CONTINUE
C         THIS ENDS THE LOOP OVER POINTS IN THE SEARCH,
C         FALLING OUT MEANS THE SEARCH DID NOT CONVERGE.
C
      IF(MASWRK) THEN
         WRITE(IW,9105) NSTEP
         WRITE(IW,9110)
      END IF
      DO I=1,NAT
         ZNUC = ZAN(I) + IZCORE(I)
         X = TOANGS*C(1,I)
         Y = TOANGS*C(2,I)
         Z = TOANGS*C(3,I)
         IF(MASWRK) WRITE(IW,9120) ANAM(I),BNAM(I),ZNUC,X,Y,Z
      END DO
      IF(MASWRK) WRITE(IW,*)
C
 1000 CONTINUE
      RETURN
C
 9010 FORMAT(/1X,'*** YOUR RUNTYP=MEX WITH SCFTYP OF MCSCF REQUESTS',
     *           ' CISTEP=',A8/
     *        1X,'ONLY CISTEP=GUGA OR CISTEP=ALDET IS PERMITTED NOW')
 9020 FORMAT(/5X,9(1H-),/5X,'STATE # ',I1,
     *       20X,'MULTIPLICITY OF STATE ',I1,' IS:',I2,
     *       /5X,9(1H-)/)
 9030 FORMAT(/10X,'TAKING NUCLEAR GRADIENT FOR STATE #',I1,
     *            ' FROM THE INPUT FILE...')
 9040 FORMAT(/1X,'SCF FOR STATE #',I1,' DID NOT CONVERGE.',
     *       '  ABORTING RUN.'/)
 9050 FORMAT(/18X,39(1H-)/18X,'GRADIENT OF ',A6,' STATE ',
     *         '(HARTREE/BOHR)'/18X,39(1H-)/
     *         8X,'ATOM     ZNUC',7X,'DE/DX',9X,'DE/DY',9X,'DE/DZ'/
     *         1X,62(1H-))
 9060 FORMAT(I3,2X,A8,A2,F6.1,3F14.7)
 9070 FORMAT(/18X,19(1H-)/18X,'EFFECTIVE GRADIENTS'/18X,19(1H-),
     *       /1X,'NSERCH = ',I3)
C
 9100 FORMAT(/2X,34(1H=)/2X,'CONVERGENCE ACHIEVED IN',I5,' STEPS'/
     *        2X,34(1H=)/
     *       2X,'CONVERGED GEOMETRY IN ANGS.'/)
 9105 FORMAT(2X,'MINIMUM ENERGY CROSSING POINT SEARCH DID NOT ',
     *       'CONVERGE IN ',I3,' STEPS.'/'  RESTART FROM THE',
     *       ' FOLLOWING GEOMETRY (IN ANGSTROM).')
 9110 FORMAT(3X,'ATOM',3X,'CHARGE',7X,'X',14X,'Y',14X,'Z'/
     *       1X,60(1H-))
 9120 FORMAT(1X,A8,A2,F5.1,3F15.10)
 9200 FORMAT(1X,'THE ORBITALS $VEC1 AND $VEC2 AT THE LAST GEOMETRY',
     *          ' WERE PUNCHED.'//
     *       1X,'IT IS A GOOD IDEA TO USE THEM IN RUNTYP=ENERGY JOBS,',
     *       1X,'AT THE MEX POINT,'/
     *       1X,'TO BE SURE THE ELECTRONIC STATES REMAIN AS YOU WANTED',
     *       1X,'THEM TO BE,'/
     *       1X,'TO DOUBLE CHECK ENERGIES, AND TO SEE THE PROPERTIES',
     *          ' OF EACH SURFACE.'/)
C
 8000 FORMAT('--------------- MEX SEARCH NSERCH=',I4,' ---',A6,
     *      ' STATE ---------------')
      END
C
C*MODULE MEXING  *DECK EFFGRD
C==================================================
      SUBROUTINE EFFGRD(N,E1,E2,GA,GB,PERPG,PARG,G)
C==================================================
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      DIMENSION GA(N),GB(N),G(N),PARG(N),PERPG(N)
C
      COMMON /RUNOPT/ RUNTYP,EXETYP,NEVALS,NGLEVL,NHLEVL
C
      PARAMETER (FACPP=30.0D+00,FACP=1.0D+00,ZERO=0.0D+00)
C           THESE FACTORS ARE ONLY REALLY IMPORTANT FOR THE FIRST STEP
C
      DATA CHECK/8HCHECK   /
C
C      COMPUTES THE PARALLEL AND PERPENDICULAR COMPENENTS OF
C            THE EFFECTIVE GRADIENT
C
      IF(EXETYP.EQ.CHECK) THEN
         CALL VCLR(PERPG,1,N)
         CALL VCLR(PARG ,1,N)
         CALL VCLR(G    ,1,N)
         RETURN
      END IF
      PG = ZERO
      PP = ZERO
      DO I = 1, N
          PERPG(I) = GA(I) - GB(I)
          PG = PG + PERPG(I)*PERPG(I)
          PP = PP + GA(I) * PERPG(I)
      END DO
      PG = SQRT(PG)
      PP = PP / PG
      DO I = 1, N
          PARG(I) = GA(I) - PERPG(I) / PG * PP
          G(I) = (E1 - E2) * FACPP * PERPG(I) + FACP * PARG(I)
      END DO
C
      RETURN
      END
C
C*MODULE MEXING  *DECK TSTCNV
C=====================================================
      SUBROUTINE TSTCNV(TGMAX,TDE,TDXMAX,TDXRMS,TGRMS,
     *           N,NATOM,EA,EB,X2,X3,
     *           PARG,PERPG,G,DELTAX,ICONV,ISTEP)
C=====================================================
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
C     CHECKS CONVERGENCE, AND UPDATES REPORT FILE
C     THERE ARE 5 CRITERIA FOR TESTING CONVERGENCE
C     AV.DELTAX, MAX.DELTAX, AV.GRAD., MAX.GRAD., DELTAE
C
      LOGICAL GOPARR,DSKWRK,MASWRK
C
      PARAMETER (MXATM=2000, MXAO=8192)
      PARAMETER (ZERO=0.0D+00)
C
      COMMON /ECP2  / CLP(400),ZLP(400),NLP(400),KFIRST(MXATM,6),
     *                KLAST(MXATM,6),LMAX(MXATM),LPSKIP(MXATM),
     *                IZCORE(MXATM)
      COMMON /INFOA / NAT,ICH,MUL,NUM,NQMT,NE,NA,NB,
     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
      COMMON /IOFILE/ IR,IW,IP,IJK,IJKT,IDAF,NAV,IODA(950)
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /RUNLAB/ TITLE(10),ANAM(MXATM),BNAM(MXATM),BFLAB(MXAO)
C
      DIMENSION X2(N),PARG(N),PERPG(N),X3(N),G(N)
      DIMENSION DELTAX(N)
C
      DE     = ABS(EA - EB)
      DXMAX  = ZERO
      DXRMS  = ZERO
      GMAX   = ZERO
      GRMS   = ZERO
      PPGRMS = ZERO
      PGRMS  = ZERO
C
      DO I = 1, N
          DELTAX(I) = X3(I) - X2(I)
          IF (ABS(DELTAX(I)) .GT. DXMAX) DXMAX = ABS(DELTAX(I))
          IF (ABS(G(I))      .GT. GMAX)   GMAX = ABS(G(I))
          DXRMS  = DXRMS  + DELTAX(I)*DELTAX(I)
          GRMS   = GRMS   + G(I)*G(I)
          PPGRMS = PPGRMS + PERPG(I)*PERPG(I)
          PGRMS  = PGRMS  + PARG(I)*PARG(I)
      END DO
C
      DXRMS  = SQRT(DXRMS  / N)
      GRMS   = SQRT(GRMS   / N)
      PPGRMS = SQRT(PPGRMS / N)
      PGRMS  = SQRT(PGRMS  / N)
C
      ICONV = 0
      IF( (GMAX  .LT. TGMAX)  .AND. (GRMS  .LT. TGRMS)   .AND.
     *    (DXMAX .LT. TDXMAX) .AND. (DXRMS .LT. TDXRMS)  .AND.
     *    (DE    .LT. TDE)   ) ICONV = 1
C
      IF(MASWRK) WRITE (IW,9155)
      IF(MASWRK) WRITE (IW,9090)
      DO I = 1, NATOM
          K = 3 * (I - 1) + 1
          ZNUC = ZAN(I) + IZCORE(I)
          IF(MASWRK) THEN
             WRITE (IW,9100) ANAM(I),BNAM(I),ZNUC,(G(J),J=K,K+2)
          END IF
      END DO
      IF(MASWRK) WRITE (IW,9020)
      IF(MASWRK) WRITE (IW,9110) PPGRMS
      DO I = 1, NATOM
          K = 3 * (I - 1) + 1
          ZNUC = ZAN(I) + IZCORE(I)
          IF(MASWRK) THEN
             WRITE (IW,9100) ANAM(I),BNAM(I),ZNUC,(PERPG(J),J=K,K+2)
          END IF
      END DO
      IF(MASWRK) WRITE (IW,9020)
      IF(MASWRK) WRITE(IW,9120) PGRMS
      DO I = 1, NATOM
          K = 3 * (I - 1) + 1
          ZNUC = ZAN(I) + IZCORE(I)
          IF(MASWRK) THEN
             WRITE (IW,9100) ANAM(I),BNAM(I),ZNUC,(PARG(J),J=K,K+2)
          END IF
      END DO
      IF(MASWRK) WRITE (IW,9021) ISTEP
      IF(MASWRK) WRITE (IW,9000) EA
      IF(MASWRK) WRITE (IW,9010) EB
      IF(MASWRK) WRITE (IW,9015) DE,TDE
      IF(MASWRK) WRITE (IW,9040) GMAX,TGMAX
      IF(MASWRK) WRITE (IW,9050) GRMS,TGRMS
      IF(MASWRK) WRITE (IW,9060) DXMAX,TDXMAX
      IF(MASWRK) WRITE (IW,9070) DXRMS,TDXRMS
      IF(MASWRK) WRITE (IW,9020)
C
      IF (ICONV .EQ. 1) THEN
          IF(MASWRK) THEN
             WRITE (IW,9130)
             WRITE (IW,9140)
          END IF
      ELSE
          IF(MASWRK) WRITE (IW,9150)
      END IF
      RETURN
C
 9020 FORMAT(' ')
 9021 FORMAT(/2X,'NSERCH=',I3,8X,
     *      'ACTUAL',17X,'THRESHOLD  '/5X,36(1H-),2X,10(1H-))
 9000 FORMAT(5X,'ENERGY OF FIRST STATE  =',F12.6)
 9010 FORMAT(5X,'ENERGY OF SECOND STATE =',F12.6)
 9015 FORMAT(5X,'ENERGY DIFFERENCE =     ',F12.6,2X,F10.6)
 9040 FORMAT(5X,'MAX EFFECTIVE GRADIENT =',F12.6,2X,F10.6)
 9050 FORMAT(5X,'RMS EFFECTIVE GRADIENT =',F12.6,2X,F10.6)
 9060 FORMAT(5X,'MAX CHANGE OF X =       ',F12.6,2X,F10.6)
 9070 FORMAT(5X,'RMS CHANGE OF X =       ',F12.6,2X,F10.6)
 9090 FORMAT(5X,'OVERALL EFFECTIVE GRADIENT:')
 9100 FORMAT(1X,A8,A2,F5.1,3F15.10)
 9110 FORMAT(5X,'DIFFERENCE GRADIENT: (RMS * DE:',F11.6,')')
 9120 FORMAT(5X,'PARALLEL GRADIENT: (RMS:',F11.6,')')
 9130 FORMAT(5X,'MEX OPTIMIZATION HAS CONVERGED.')
 9140 FORMAT(5X,'BYE, BYE, BUH-BYE.  THANK YOU, AND FLY WITH US ',
     *        'AGAIN.')
 9150 FORMAT(5X,'CONVERGENCE NOT ACHIEVED YET',/5X,28(1H=))
 9155 FORMAT(/3X,'ATOM',3X,'CHARGE',7X,'X',14X,'Y',14X,'Z'/
     *       1X,60(1H-))
      END
C*MODULE MEXING  *DECK UPDTX
C===========================================================
      SUBROUTINE UPDTX(N,NSTEP,STPMX,X1,X2,X3,G1,G2,HI1,HI2,
     *                   DELG,HDELG,CHGEX,DELX,W)
C===========================================================
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
C     ADAPTED BFGS ROUTINE FROM NUMERICAL RECIPES
C
      DIMENSION X1(N),X2(N),G1(N), G2(N),
     *           HI1(N,N),X3(N),HI2(N,N)
      DIMENSION DELG(N),HDELG(N),CHGEX(N),DELX(N),W(N)
      PARAMETER (ZERO  = 0.0D+00,ONE=1.0D+00)
CNM   PARAMETER (STPMX = 0.1D0)
CNM   PARAMETER (STPMX = 0.03D0)
CNM   PARAMETER (STPMX = 0.2D0)
C
      STPMAX = STPMX * N
      IF (NSTEP .EQ. 0) THEN
          DO I = 1, N
             CHGEX(I) = -G2(I)
             DO J = 1, N
                 HI2(I,J) = HI1(I,J)
             END DO
         END DO
      ELSE
         DO I = 1, N
            DELG(I) = G2(I) - G1(I)
            DELX(I) = X2(I) - X1(I)
         END DO
         DO I = 1, N
            HDELG(I) = ZERO
            DO J = 1, N
               HDELG(I) = HDELG(I) + HI1(I,J) * DELG(J)
            END DO
         END DO
         FAC = ZERO
         FAE = ZERO
         SUMDG = ZERO
         SUMDX = ZERO
         DO I = 1, N
            FAC = FAC + DELG(I) * DELX(I)
            FAE = FAE + DELG(I) * HDELG(I)
            SUMDG = SUMDG + DELG(I)*DELG(I)
            SUMDX = SUMDX + DELX(I)*DELX(I)
         END DO
         FAC = ONE / FAC
         FAD = ONE / FAE
         DO I = 1, N
            W(I) = FAC * DELX(I) - FAD * HDELG(I)
         END DO
         DO I = 1, N
            DO J = 1, N
               HI2(I,J) = HI1(I,J) + FAC * DELX(I) * DELX(J) -
     *         FAD * HDELG(I) * HDELG(J) + FAE * W(I) * W(J)
            END DO
         END DO
         DO I = 1, N
            CHGEX(I) = ZERO
            DO J = 1, N
               CHGEX(I) = CHGEX(I) - HI2(I,J) * G2(J)
            END DO
         END DO
      END IF
C
      FAC = ZERO
      DO I = 1, N
         FAC = FAC + CHGEX(I)*CHGEX(I)
      END DO
      STPL = SQRT(FAC)
      IF(STPL .GT. STPMAX) THEN
         DO I = 1, N
            CHGEX(I) = CHGEX(I) / STPL * STPMAX
         END DO
      END IF
      GSTST = ZERO
      DO I = 1, N
         IF(ABS(CHGEX(I)) .GT. GSTST) GSTST = ABS(CHGEX(I))
      END DO
      IF(GSTST .GT. STPMX) THEN
         DO I = 1, N
            CHGEX(I) = CHGEX(I) / GSTST * STPMX
         END DO
      END IF
C
C           ENFORCE POINT GROUP SYMMETRY.
C        THIS PROJECTS THE TOTALLY SYMMETRIC COMPONENT OF THE STEP.
C        IF ALL IS GOING WELL, THE GRADIENTS SHOULD BE TOTALLY
C        SYMMETRIC TO MACHINE PRECISION, AND HENCE THE STEP TOO,
C        SO MOST OF THE TIME THIS IS JUST BEING REALLY CAREFUL.
C
      CALL SYMDR(CHGEX)
C
      DO I = 1, N
         X3(I) = X2(I) + CHGEX(I)
      END DO
C
      RETURN
      END
C
C*MODULE MEXING  *DECK MEXINP
C======================
      SUBROUTINE MEXINP
C======================
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      LOGICAL GOPARR,DSKWRK,MASWRK,MP2
C
      PARAMETER (MXATM=2000, MXFRG=1050, MXDFG=5,
     *           MXDPPT=MXFRG*MXDFG*12)
      PARAMETER (ZERO=0.0D+00)
C
      COMMON /DFTPAR/ DFTTYP(20),EXENA,EXENB,EXENC,IDFT34,NAUXFUN,
     *                                                    NAUXSHL
      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 /IOFILE/ IR,IW,IP,IJK,IJKT,IDAF,NAV,IODA(950)
      COMMON /MEXPAR/ SCF1,SCF2,TGMAX,TDE,TDXMAX,TDXRMS,TGRMS,STPSZ,
     *                MULT1,MULT2,NSTEP,NRDMOS,NMOS1,NMOS2,NPRT
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /RUNOPT/ RUNTYP,EXETYP,NEVALS,NGLEVL,NHLEVL
      COMMON /WFNOPT/ SCFTYP,VBTYP,DFTYPE,TDDFTYP,CITYP,CCTYP,
     *                MPLEVL,MPCTYP
C
      PARAMETER (NNAM=15)
      DIMENSION QNAM(NNAM),KQNAM(NNAM)
      DATA WRDMEX /8HMEX     /
      DATA QNAM   /8HMULT1   ,8HMULT2   ,8HSCF1    ,8HSCF2    ,
     *             8HNSTEP   ,8HTDE     ,8HTDXMAX  ,
     *             8HTDXRMS  ,8HTGMAX   ,8HTGRMS   ,8HNRDMOS  ,
     *             8HNMOS1   ,8HNMOS2   ,8HSTPSZ   ,8HNPRT    /
      DATA KQNAM  /1,1,5,5,   1,3,3,   3,3,3,1,   1,1,3,1/
C
      DATA RHF,UHF,ROHF,GVB,RMC
     *     /8HRHF     ,8HUHF     ,8HROHF    ,8HGVB     ,8HMCSCF   /
      DATA RNONE/8HNONE    /
C
      SCF1    = RHF
      SCF2    = ROHF
      MULT1   = 1
      MULT2   = 3
      NSTEP   = 50
      TDE     = 1.0D-05
      TDXMAX  = 2.0D-03
      TDXRMS  = 1.5D-03
      TGMAX   = 5.0D-04
      TGRMS   = 3.0D-04
      STPSZ   = 0.1D+00
      NRDMOS  = 0
      NMOS1   = 0
      NMOS2   = 0
      NPRT    = 0
C
C     ----- GET NAMELIST $MEX -----
C
      CALL NAMEIO(IR,JRET,WRDMEX,NNAM,QNAM,KQNAM,
     *            MULT1,MULT2,SCF1,SCF2,
     *            NSTEP,TDE,TDXMAX,TDXRMS,TGMAX,
     *            TGRMS,NRDMOS,NMOS1,NMOS2,STPSZ,NPRT,
     *            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)
C
      IF(JRET.EQ.2) THEN
         IF (MASWRK) WRITE(IW,*) 'ERROR READING $MEX GROUP'
         CALL ABRT
      END IF
C
C     ----- KILL IF ...
C
      NERR=0
C
C        THE USE OF NUMERICAL GRADIENTS HAS TWO PITFALLS:
C        ONE: POSSIBLE COLLAPSE OF THE ELECTRONIC STATE TO A LOWER
C             STATE, ALTHOUGH DISPLACEMENTS ARE POINT GROUP CONSERVING.
C        TWO: NEW PAIRINGS OF METHODS WILL NOT HAVE BEEN SOUNDED OUT
C             TO SEE IF THEY SHARE DAF RECORD NUMBERS OR OTHER ISSUES.
C        THE 2ND IS THE REAL REASON FOR TRAPPING NUMERICAL GRADIENTS!
C
      IF(NGLEVL.NE.0) THEN
         IF(MASWRK) WRITE(IW,9100)
         NERR = NERR+1
      END IF
 9100 FORMAT(/1X,'*** RUNTYP=MEX REQUIRES ANALYTIC GRADIENTS ***')
C
C        CERTAINLY -EFP- COMPUTATIONS DON'T WORK, FOR THERE IS
C        ABSOLUTELY NO CODE HERE TO MOVE EFP PARTICLES.  THE
C        USE OF -PCM- IS PERMITTED, BUT PERHAPS DODGY, SINCE THE
C        PROGRAMMERS HAVEN'T PUT ANY THOUGHT INTO WHETHER THE
C        PCM MODULE IS ABLE TO DEAL WITH TWO DIFFERENT STATES
C        IN A SINGLE CAVITY.
C
      IF(NFRG.GT.0) THEN
         IF(MASWRK) WRITE(IW,9105)
         NERR = NERR+1
      END IF
 9105 FORMAT(/1X,'*** RUNTYP=MEX LACKS ANY CODING FOR -EFP- RUNS ***')
C
C        THE WAVEFUNCTIONS SHOULD BE SIMILAR IN NATURE, SO
C        MIXING ROHF WITH MCSCF, OR DFT WITH MP2, IS TRAPPED.
C
C        BEYOND COMMON SENSE, IT IS ALSO NECESSARY TO TRY BOTH
C        WAVEFUNCTIONS IN THE SAME JOB, TO MAKE SURE THAT THE
C        DIRECT ACCESS FILE RECORDS DON'T CLASH, AND THAT THE
C        ENERGIES IN RUNTYP=MEX MATCH SEPARATE RUNTYP=ENERGY.
C        THIS IS HARDER THAN IT SOUNDS, AND NOT ALL THE CASES
C        WILL HAVE BEEN TESTED.
C
C        THERE ARE SOME FUNDAMENTAL ISSUES, E.G. GVB WITH GVB
C        MIGHT REALLY BE OF USE, BUT PRESENTLY, ONLY ONE $SCF
C        GROUP IS READ.  IN THE CASE OF MCSCF, WE DEAL WITH THIS
C        BY SPECIFYING THE ACTIVE SPACE IN DISTINCT $DET1/$DET2
C        OR $DRT1/$DRT2 GROUPS.
C
C        WHAT IS MISSING?  CITYP=CIS WOULD CLEARLY BE USEFUL,
C        BUT THE APPROACH HERE (1 SET OF INTEGRALS, THEN 2
C        DISTINCT SCF, EACH FOLLOWED BY A GRADIENT) WILL NEED
C        QUITE A BIT OF WORK: CIS CODE SHOULD RUN 1 INTEGRALS,
C        THEN 1 SCF, AND TRY TO SET UP 2 GRADIENTS IN THE
C        COURSE OF FINDING 2 OR MORE ROOTS DURING 1 CIS STEP.
C
      MP2 = MPLEVL.EQ.2
C
      IF(MP2  .AND.  DFTTYP(1).NE.ZERO) THEN
         IF(MASWRK) WRITE(IW,9110)
         NERR = NERR+1
      END IF
 9110 FORMAT(/1X,'CAN NOT MIX DFT WITH MP2!  QUITTING... '/)
C
      IF(MP2) THEN
         IF(SCF1.EQ.GVB .OR. SCF2.EQ.GVB) THEN
            IF(MASWRK) WRITE(IW,9120)
            NERR=NERR+1
         END IF
         IF(SCF1.EQ.RMC .OR. SCF2.EQ.RMC) THEN
            IF(MASWRK) WRITE(IW,9130)
            NERR=NERR+1
         END IF
      END IF
 9120 FORMAT(/1X,'*** THERE IS NO MP2 THEORY FOR GVB ***')
 9130 FORMAT(/1X,'*** MULTI-REFERENCE PERTURBATION THEORY LACKS',
     *          ' ANALYTIC GRADIENTS ***')
C
      IF(SCF1.EQ.RMC .AND. SCF2.NE.RMC .OR.
     *   SCF1.NE.RMC .AND. SCF2.EQ.RMC) THEN
         IF(MASWRK) WRITE(IW,9140)
         NERR=NERR+1
      END IF
 9140 FORMAT(/1X,'*** FOR CONSISTENCY, BOTH STATES SHOULD BE',
     *          ' DESCRIBED BY MCSCF ***')
C
      IF(SCF1.EQ.RMC .OR. SCF2.EQ.RMC) THEN
         IF(NRDMOS.EQ.0) THEN
            IF(MASWRK) WRITE(IW,9150)
            NERR=NERR+1
         END IF
      END IF
 9150 FORMAT(/1X,'*** YOU SHOULD READ IN ORBITALS FOR MCSCF ***'/
     *  ' PLEASE CHANGE NRDMOS=0 TO =2, AND GIVE NMOS1 AND NMOS2')
C
      KERR = 0
      IF(SCF1.EQ.SCF2) THEN
         KERR=KERR+1
         IF(MULT1.NE.MULT2  .OR.  SCF1.EQ.RMC) KERR=KERR-1
      END IF
      IF(KERR.GT.0) THEN
         IF(MASWRK) WRITE(IW,9160) SCF1,MULT1,SCF2,MULT2
         NERR=NERR+1
      END IF
 9160 FORMAT(/1X,'*** EQUIVALENT WAVEFUNCTIONS WERE CHOSEN,',
     *           ' SCF1=',A8,' SCF2=',A8/
     *        1X,'THESE CANNOT BE DISTINCT ELECTRONIC STATES.')
C
      IF((SCF1.EQ.ROHF  .AND.  SCF2.EQ.UHF) .OR.
     *   (SCF2.EQ.ROHF  .AND.  SCF1.EQ.UHF)) THEN
         IF(MASWRK) WRITE(IW,9170)
         NERR=NERR+1
      END IF
 9170 FORMAT(/1X,'*** SCF1=',A8,' SCF2=',A8,' IS NOT PERMISSIBLE.'/
     *        1X,'USE BOTH AS ROHF, OR BOTH AS UHF.')
C
      IF(CITYP.NE.RNONE  .OR.  CCTYP.NE.RNONE) THEN
         IF(MASWRK) WRITE(IW,9180) CITYP,CCTYP
         NERR=NERR+1
      END IF
 9180 FORMAT(/1X,'*** CITYP=',A8,' AND/OR CCTYP=',A8,
     *           ' CANNOT BE USED HERE.')
C
C        TRY TO SET THE WAVEFUNCTION SPECIFICATIONS SO AS TO
C        READ THE APPROPRIATE INPUT GROUPS TO SET THEM UP.
C        FOR EXAMPLE, IF -UHF- IS TO BE USED, IT IS GOOD FOR THE
C        $SCF TO DECIDE THAT DIIS IS THE APPROPRIATE CONVERGER,
C        SO KEEP UHF AFTER RHF/ROHF.
C        GVB HAS TO COME AFTER THOSE THREE, SO THAT WE GET THE
C        WAVEFUNCTION DEFINITIONS IN $SCF PROCESSED PROPERLY.
C        SO, GENERALLY SPEAKING, CHOOSE THE TRICKER WAVEFUNCTION.
C
      SCFTYP=RHF
      IF(SCF1.EQ.ROHF .OR.  SCF2.EQ.ROHF) SCFTYP=ROHF
      IF(SCF1.EQ.UHF  .OR.  SCF2.EQ.UHF)  SCFTYP=UHF
      IF(SCF1.EQ.GVB  .OR.  SCF2.EQ.GVB)  SCFTYP=GVB
      IF(SCF1.EQ.RMC  .OR.  SCF2.EQ.RMC)  SCFTYP=RMC
                      MUL=MAX(MULT1,MULT2)
      IF(SCF1.EQ.GVB) MUL=MULT1
      IF(SCF2.EQ.GVB) MUL=MULT2
C
      IF(MASWRK) THEN
         WRITE(IW,9000) SCF1,SCF2,
     *                  MULT1,MULT2,
     *                  MPLEVL,MPLEVL,
     *                  DFTYPE,DFTYPE,
     *                  NMOS1,NMOS2,NRDMOS,
     *                  NPRT,NSTEP,STPSZ,
     *                  TDE,TDXMAX,TDXRMS,TGMAX,TGRMS
      END IF
C
      IF(NERR.GT.0) THEN
         IF(MASWRK) WRITE(IW,9010)
         CALL ABRT
      ENDIF
      RETURN
C
 9000 FORMAT(/5X,12(1H-)/5X,'$MEX OPTIONS'/5X,12(1H-),
     *       //5X,'SCF1    =',A8,5X,'SCF2   =',A8/
     *         5X,'MULT1   =',I8,5X,'MULT2  =',I8/
     *         5X,'MPLEVL  =',I8,5X,'MPLEVL =',I8/
     *         5X,'DFTTYP  =',A8,5X,'DFTTYP =',A8/
     *         5X,'NMOS1   =',I8,5X,'NMOS2  =',I8,5X,'NRDMOS=',I8/
     *         5X,'NPRT    =',I8/
     *         5X,'NSTEP   =',I8,5X,'STPSZ   =',1P,E9.2/
     *         5X,'ENERGY DIFFERENCE THRESHOLD,    TDE =',1P,E9.2/
     *         5X,'MAX. CHANGE IN DISPLACEMENT, TDXMAX =',1P,E9.2/
     *         5X,'RMS CHANGE IN DISPLACEMENT,  TDXRMS =',1P,E9.2/
     *         5X,'MAX. CHANGE IN GRADIENT,      TGMAX =',1P,E9.2/
     *         5X,'RMS CHANGE IN GRADIENT,       TGRMS =',1P,E9.2)
 9010 FORMAT(/1X,'*** PLEASE FIX THE ERROR QUOTED ABOVE THE INPUT,',
     *      ' AND RESUBMIT YOUR MEX RUN'/)
      END
C
C*MODULE MEXING  *DECK SETCAN
C======================
      SUBROUTINE SETCAN
C======================
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      LOGICAL COUPLE,GOPARR,DSKWRK,MASWRK
C
      PARAMETER (MXATM=2000, MXAO=8192)
      PARAMETER (ZERO=0.0D+00, HALF=0.5D+00, ONE=1.0D+00, TWO=2.0D+00)
C
      COMMON /GVBWFN/ CICOEF(2,12),F(25),ALPHA(325),BETA(325),NO(10),
     *                NCO,NSETO,NOPEN,NPAIR,NORB,NCONF(MXAO),NHAM
      COMMON /INFOA / NAT,ICH,MUL,NUM,NQMT,NE,NA,NB,
     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
      COMMON /IOFILE/ IR,IW,IP,IJK,IJKT,IDAF,NAV,IODA(950)
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /SCFWFN/ AROHF(3),BROHF(3),PACAVO(6),IACAVO,IUHFNO,MVOQ
      COMMON /WFNOPT/ SCFTYP,VBTYP,DFTYPE,TDDFTYP,CITYP,CCTYP,
     *                MPLEVL,MPCTYP
C
      DATA ROHF,GVB/8HROHF    ,8HGVB     /
      DATA RHF,UHF,RMC/8HRHF     ,8HUHF     ,8HMCSCF   /
C
C     ----- SET UP DEFAULT COUPLING CONSTANTS -----
C
      IF(SCFTYP.EQ.RHF) RETURN
      IF(SCFTYP.EQ.UHF) RETURN
      IF(SCFTYP.EQ.RMC) RETURN
C
C     NOTE THAT ALL -GVB- CASES ALLOW PAIRS TO BE COUPLED BEHIND
C     THE CLOSED AND OPEN SHELL SPACES, IF DESIRED.
C
      COUPLE = .FALSE.
C
C     ----- HIGH SPIN CASE (ROHF MODULE) -----
C     THE DEFAULT IS THE ROOTHAAN SINGLE MATRIX CASE.
C     FOR OTHER INTERESTING CASES, SEE ROUTINE ROFOCK.
C
      IF(SCFTYP.EQ.ROHF) THEN
         ALPHA(1)= -HALF
         ALPHA(2)=  HALF
         ALPHA(3)=  ONE+HALF
         BETA(1) =  ONE+HALF
         BETA(2) =  HALF
         BETA(3) = -HALF
         COUPLE = .TRUE.
      END IF
C
C     ----- CLOSED SHELL SINGLET (GVB MODULE) -----
C
      IF(NSETO.EQ.0  .AND.  MUL.EQ.1  .AND.  SCFTYP.EQ.GVB) THEN
         F(1) = ONE
         ALPHA(1) = TWO
         BETA(1) = -ONE
         COUPLE=.TRUE.
      END IF
C
C     ----- OPEN SHELL DOUBLET (GVB MODULE) -----
C
      IF(NSETO.EQ.1  .AND.  NO(1).EQ.1  .AND.
     *     MUL.EQ.2  .AND.  SCFTYP.EQ.GVB) THEN
         IF(NCO.GT.0) THEN
            F(1) = ONE
            F(2) = HALF
            ALPHA(1) = TWO
            ALPHA(2) = ONE
            ALPHA(3) = HALF
            BETA(1) = -ONE
            BETA(2) = -HALF
            BETA(3) = -HALF
         ELSE
            F(1) = HALF
            ALPHA(1) = HALF
            BETA(1) = -HALF
         END IF
         COUPLE=.TRUE.
      END IF
C
C     ----- OPEN SHELL SINGLET (GVB MODULE) -----
C
      IF(NSETO.EQ.2  .AND.  NO(1).EQ.1  .AND.  NO(2).EQ.1   .AND.
     *     MUL.EQ.1  .AND.  SCFTYP.EQ.GVB) THEN
         IF(NCO.GT.0) THEN
            F(1) = ONE
            F(2) = HALF
            F(3) = HALF
            ALPHA(1) = TWO
            ALPHA(2) = ONE
            ALPHA(3) = ZERO
            ALPHA(4) = ONE
            ALPHA(5) = HALF
            ALPHA(6) = ZERO
            BETA(1) = -ONE
            BETA(2) = -HALF
            BETA(3) =  ZERO
            BETA(4) = -HALF
            BETA(5) =  HALF
            BETA(6) =  ZERO
         ELSE
            F(1) = HALF
            F(2) = HALF
            ALPHA(1) = ZERO
            ALPHA(2) = HALF
            ALPHA(3) = ZERO
            BETA(1) =  ZERO
            BETA(2) =  HALF
            BETA(3) =  ZERO
         END IF
         COUPLE=.TRUE.
      END IF
C
C     ----- OPEN SHELL TRIPLET (GVB MODULE) -----
C
      IF(NSETO.EQ.2  .AND.  NO(1).EQ.1  .AND.  NO(2).EQ.1   .AND.
     *     MUL.EQ.3  .AND.  SCFTYP.EQ.GVB) THEN
         IF(NCO.GT.0) THEN
            F(1) = ONE
            F(2) = HALF
            F(3) = HALF
            ALPHA(1) = TWO
            ALPHA(2) = ONE
            ALPHA(3) = HALF
            ALPHA(4) = ONE
            ALPHA(5) = HALF
            ALPHA(6) = HALF
            BETA(1) = -ONE
            BETA(2) = -HALF
            BETA(3) = -HALF
            BETA(4) = -HALF
            BETA(5) = -HALF
            BETA(6) = -HALF
         ELSE
            F(1) = HALF
            F(2) = HALF
            ALPHA(1) = HALF
            ALPHA(2) = HALF
            ALPHA(3) = HALF
            BETA(1) = -HALF
            BETA(2) = -HALF
            BETA(3) = -HALF
         END IF
         COUPLE=.TRUE.
      END IF
C
C     IF VALUES WERE NOT DEFAULTS AND WERE NOT INPUT, BOMB THE JOB.
C
      IF(.NOT.COUPLE) THEN
         IF (MASWRK) WRITE (IW,9228) SCFTYP,MUL,NSETO
         CALL ABRT
      END IF
C
C     ----- THE ROHF-IN-UHF CODE WORKS WITH AROHF AND BROHF -----
C     WE SET F,ALPHA,BETA FOR THE OPEN SHELL CPHF CODE'S USE.
C
      IF(SCFTYP.EQ.ROHF) THEN
         AROHF(1) = ALPHA(1)
         AROHF(2) = ALPHA(2)
         AROHF(3) = ALPHA(3)
         BROHF(1) =  BETA(1)
         BROHF(2) =  BETA(2)
         BROHF(3) =  BETA(3)
C
         IF(NCO.GT.0) THEN
            IF(NA.EQ.NB) THEN
               F(1) = ONE
               ALPHA(1) = TWO
               BETA(1)  = -ONE
            ELSE
               F(1) = ONE
               F(2) = HALF
               ALPHA(1) = TWO
               ALPHA(2) = ONE
               ALPHA(3) = HALF
               BETA(1) = -ONE
               BETA(2) = -HALF
               BETA(3) = -HALF
            END IF
         ELSE
            F(1) = HALF
            ALPHA(1) = HALF
            BETA(1) = -HALF
         END IF
      END IF
C
 9228 FORMAT(//1X,'**** ERROR ****'/1X,'SCFTYP=',A8,' MULT=',I4,
     *            ' NSETO=',I4,' HAS NO DEFAULT COUPLING PARAMETERS'/
     *      1X,'READ THE ROHF-GVB SECTION IN THE REFS.DOC CHAPTER!')
      RETURN
      END
C
C*MODULE MEXING  *DECK GUSMX
C==========================================================
      SUBROUTINE GUSMX(NORB,NBF,NA,NB,IWRK,EIG,SCR,S,Q,VEC,
     *                 L1,L2,L3,SCFTYP,STVEC,NRECA,NRECB)
C==========================================================
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      LOGICAL GOPARR,DSKWRK,MASWRK
C
      CHARACTER*8 STVEC
C
      COMMON /IOFILE/ IR,IW,IP,IJK,IJKT,IDAF,NAV,IODA(950)
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
C
      DATA UHF/8HUHF     /
C
      DIMENSION IWRK(L1),EIG(L1),SCR(L1,8),S(L2),Q(L3),VEC(L1,L1)
C
      IPASS = 1
      NSAV  = NRECA
      TOLZ=1.0D-08
      TOLE=1.0D-05
      NUM = NBF
C
C     ----- READ THE INITIAL ORBITALS FROM THE INPUT FILE -----
C
      IZERO = 0
      IONE  = 1
      ITWO  = 2
C
C            GET CANONICAL ORBITALS
C
      CALL DAREAD(IDAF,IODA,S,L2,12,0)
      CALL QMTSYM(S,VEC,Q,EIG,SCR,IWRK,L0,L1,L2,L3,DBG)
      CALL DAWRIT(IDAF,IODA,Q,L3,45,0)
C
C            DEFAULT FOR NMOS IS THE NUMBER OF OCCUPIED ORBITALS
C            OR WHAT THE USER INPUT (NORB)
C
      NUM1 = NUM
      NMOS = NA
      IF(NORB.GT.0)     NMOS = NORB
C
      IF(NMOS.EQ.0.AND.MASWRK) WRITE(IW,9080)
C
C               POSITION THE INPUT FILE TO $VEC
C
      CALL SEQREW(IR)
      JJ = 0
      ICC = 0
      CALL FNDGRP(IR,STVEC,IEOF)
      IF (IEOF.EQ.0) GO TO 200
C
C               $VEC CARD NOT FOUND, KILL JOB
C
      IF (MASWRK) WRITE(IW,9020) STVEC
      CALL ABRT
      STOP
C
C               $VEC CARD FOUND, READ THE ORBITALS
C
  200 CONTINUE
      IF (MASWRK) THEN
      DO 280 J = 1,NMOS
         IMAX = 0
         IC = 0
  240    CONTINUE
            IMIN = IMAX+1
            IMAX = IMAX+5
            IC = IC+1
            IF(IMAX .GT. NUM1) IMAX = NUM1
            READ(IR,9040,END=300,ERR=300) JJ,ICC,(VEC(I,J),I=IMIN,IMAX)
            MODJ=MOD(J,100)
            IF(JJ.EQ.MODJ . AND.  ICC.EQ.IC) GO TO 260
               WRITE(IW,9060) J,IC,JJ,ICC
               IF (GOPARR) CALL DDI_BCAST(350,'I',ITWO,1,MASTER)
               CALL ABRT
               STOP
  260       CONTINUE
         IF(IMAX .LT. NUM1) GO TO 240
  280 CONTINUE
      IF (GOPARR) CALL DDI_BCAST(350,'I',IZERO,1,MASTER)
C
C SLAVE WORK
C
      ELSE
         IF (GOPARR) CALL DDI_BCAST(350,'I',IEND,1,MASTER)
         IF (IEND.EQ.IONE) GO TO 300
         IF (IEND.EQ.ITWO) CALL ABRT
      END IF
C
C  GIVE VECTORS TO ALL PROCESSES
C
      IF (GOPARR) CALL DDI_BCAST(351,'F',VEC,NMOS*L1,MASTER)
C
C            ORTHONORMALIZE THE ORBITALS
C            BACK-TRANSFORM THE MO'S
C
      CALL CLENMO(VEC,L1,NMOS,TOLZ,TOLE,IW,DBG)
C
C     IF MFRZ, REARRANGE THE VECTORS BEFORE TRANSFORMING THEM
C
      CALL ORTHO(Q,S,VEC,SCR,NMOS,L0,L1,L2,L1)
      CALL TFSQB(VEC,Q,SCR,L0,L1,L1)
      CALL STFASE(VEC,L1,L1,L1)
      CALL CLENMO(VEC,L1,L1,TOLZ,TOLE,IW,.FALSE.)
      CALL DAWRIT(IDAF,IODA,VEC,L3,NSAV,0)
C
C            GO BACK AND GET BETA MO S FOR UHF
C
      IF(SCFTYP.NE.UHF) RETURN
      IF(NB.EQ.0) RETURN
      IF(IPASS.EQ.2) RETURN
      IPASS = 2
      NMOS = NB
CNM   IF(NORB.GT.0) NMOS = NORB
      NSAV = NRECB
      GO TO 200
C
C        PREMATURE END OF $VEC INPUT ENCOUNTERED
C
  300 CONTINUE
      IF (MASWRK) THEN
         IF (GOPARR) CALL DDI_BCAST(350,'I',IONE,1,MASTER)
         WRITE(IW,9100) J,IC
      END IF
      CALL ABRT
      STOP
C
 9020 FORMAT(' NO ',A8,' INPUT FOUND - STOP')
 9040 FORMAT(I2,I3,5E15.8)
 9060 FORMAT(' *** ERROR READING THE MO COEFFICIENTS AT',2I4/
     *       '     THE INPUT VALUES WERE',2I4)
 9080 FORMAT(1X,'YOU FORGOT TO GIVE NORB IN $GUESS')
 9100 FORMAT(' *** ERROR: PREMATURE END OF ORBITAL INPUT ENCOUNTERED'/
     *       '            LOOKING FOR ORBITAL',I4,' ELEMENT',I4)
      END
C
C*MODULE MEXING  *DECK GESFMX
C=======================================================================
      SUBROUTINE GESFMX(NBF,NE,MULT1,MULT2,SCF1,SCF2,NRDMOS,NMOS1,NMOS2)
C=======================================================================
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      CHARACTER*8 STVEC
      LOGICAL GOPARR,DSKWRK,MASWRK
C
      COMMON /FMCOM / X(1)
      COMMON /IOFILE/ IR,IW,IP,IJK,IJKT,IDAF,NAV,IODA(950)
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
C
      DATA UHF/8HUHF     /
C
C     ----- PREP FOR READING $VECS -----
C
      L1   = NBF
      L3   = NBF*NBF
      L2   = (L3+NBF)/2
      LL1  = L1
      LL2  = L2
      LL3  = L3
      MSPN = MULT1 - 1
      NB1  = (NE - MSPN)/2
      NA1  = NB1 + MSPN
      MSPN = MULT2 -1
      NB2  = (NE - MSPN)/2
      NA2  = NB2 + MSPN
C
      CALL VALFM(LOADFM)
      IWRK = LOADFM + 1
      IEIG = IWRK   + L1
      ISCR = IEIG   + LL1
      IH   = ISCR   + 8*L1
      IQ   = IH     + LL2
      IVEC = IQ     + LL3
      LAST = IVEC   + LL3
      NEED = LAST - LOADFM -1
      CALL GETFM(NEED)
C
      IF(NRDMOS.EQ.3) THEN
         DO I=1,2
            IF(I.EQ.1) THEN
               STVEC  = ' $VEC1  '
               SCFTYP = SCF1
               NRECA  = 406
               NRECB  = 407
               SCFTYP = SCF1
               NORB   = NMOS1
               NA     = NA1
               NB     = NB1
            END IF
            IF(I.EQ.2) THEN
               STVEC  = ' $VEC2  '
               SCFTYP = SCF2
               NRECA  = 408
               NRECB  = 409
               SCFTYP = SCF2
               NORB   = NMOS2
               NA     = NA2
               NB     = NB2
            END IF
            IF(SCFTYP.EQ.UHF) THEN
               IF(MASWRK) WRITE(IW,9001) NORB,I,I,NRECA,NRECB,SCFTYP
            ELSE
               IF(MASWRK) WRITE(IW,9000) NORB,I,I,NRECA,SCFTYP
            END IF
            CALL GUSMX(NORB,NBF,NA,NB,X(IWRK),X(IEIG),X(ISCR),X(IH),
     *                 X(IQ),X(IVEC),L1,L2,L3,SCFTYP,STVEC,NRECA,NRECB)
         END DO
      ELSE IF(NRDMOS.EQ.1) THEN
         STVEC  = ' $VEC1  '
         NRECA  = 406
         NRECB  = 407
         SCFTYP = SCF1
         NORB   = NMOS1
         NA     = NA1
         NB     = NB1
         IF(SCFTYP.EQ.UHF) THEN
            IF(MASWRK) WRITE(IW,9001) NORB,1,1,NRECA,NRECB,SCFTYP
         ELSE
            IF(MASWRK) WRITE(IW,9000) NORB,1,1,NRECA,SCFTYP
         END IF
         CALL GUSMX(NORB,NBF,NA,NB,X(IWRK),X(IEIG),X(ISCR),X(IH),X(IQ),
     *              X(IVEC),L1,L2,L3,SCFTYP,STVEC,NRECA,NRECB)
      ELSE IF(NRDMOS.EQ.2) THEN
         STVEC  = ' $VEC2  '
         NRECA  = 408
         NRECB  = 409
         SCFTYP = SCF2
         NORB   = NMOS2
         NA     = NA2
         NB     = NB2
         IF(SCFTYP.EQ.UHF) THEN
            IF(MASWRK) WRITE(IW,9001) NORB,2,2,NRECA,NRECB,SCFTYP
         ELSE
            IF(MASWRK) WRITE(IW,9000) NORB,2,2,NRECA,SCFTYP
         END IF
         CALL GUSMX(NORB,NBF,NA,NB,X(IWRK),X(IEIG),X(ISCR),X(IH),X(IQ),
     *              X(IVEC),L1,L2,L3,SCFTYP,STVEC,NRECA,NRECB)
      END IF
C
      CALL RETFM(NEED)
      RETURN
 9000 FORMAT(/1X,'READING AND STORING',I4,' ORBITALS FROM ',
     *  '$VEC',I1,' FOR STATE ',I1,/1X,57(1H-),/20X,'RECORD WRITTEN TO',
     *  I4,' FOR ',A8)
 9001 FORMAT(/1X,'READING AND STORING',I4,' ORBITALS FROM ',
     *  '$VEC',I1,' FOR STATE ',I1,/1X,57(1H-),/20X,'RECORD WRITTEN TO',
     *  I4,' AND ',I4,' FOR ',A8)
      END
C
C*MODULE MEXING  *DECK ORBMAN
C======================================================
      SUBROUTINE ORBMAN(MEXSTATE,VEC1,VEC2,NBF2,SCFTYP)
C======================================================
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
C     PUT ORBITALS IN RIGHT PLACES
C
      LOGICAL GOPARR,DSKWRK,MASWRK
      DIMENSION VEC1(NBF2),VEC2(NBF2)
      COMMON /IOFILE/ IR,IW,IP,IJK,IJKT,IDAF,NAV,IODA(950)
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
C
      DATA UHF /8HUHF     /
C
      IF(MEXSTATE.EQ.1) THEN
         NRECA = 406
         NRECB = 407
      END IF
      IF(MEXSTATE.EQ.2) THEN
         NRECA = 408
         NRECB = 409
      END IF
C
      IF(SCFTYP.EQ.UHF) THEN
         CALL DAREAD(IDAF,IODA,VEC1,NBF2,NRECA,0)
         CALL DAWRIT(IDAF,IODA,VEC1,NBF2,15,0)
         CALL DAREAD(IDAF,IODA,VEC2,NBF2,NRECB,0)
         CALL DAWRIT(IDAF,IODA,VEC2,NBF2,19,0)
         IF(MASWRK) WRITE(IW,9001) MEXSTATE,15,19
      ELSE
         CALL DAREAD(IDAF,IODA,VEC1,NBF2,NRECA,0)
         CALL DAWRIT(IDAF,IODA,VEC1,NBF2,15,0)
         IF(MASWRK) WRITE(IW,9000) MEXSTATE,15
      END IF
      RETURN
C
 9000 FORMAT(/5X,'COPIED ORBITALS FOR STATE #',I1,' TO RECORD ',I4)
 9001 FORMAT(/5X,'COPIED ORBITALS FOR STATE #',I1,' TO RECORDS ',I4,
     *       ' AND ',I4)
      END
C
