C 19 SEP 05 - MWS - ADD TRUE NUCLEAR CHARGE ARRAY TO INFOA COMMON
C  5 JUL 05 - MWS - SELECT NEW ATOM,BASIS,EFP,PCM,DAF DIMENSIONS
C 23 JUL 04 - MWS - MOLDAT,RDMOL: REMOVE RESIDUAL FTNCHEK COMMON ERRORS
C 16 JAN 04 - HL,MWS - MOLDAT: MAKE PEPTIDE LINKAGE MM TERM OPTIONAL
C  3 SEP 03 - MWS - MOLDAT: ENSURE PEPTIDE DIMENSION NOT EXCEEDED
C 16 JUN 03 - HL  - MOLDAT: INCLUDE MM CORRECTION FOR PEPTIDE
C 25 JUN 01 - MWS - ALTER COMMON BLOCK WFNOPT
C 13 JUN 96 - VAG - ADD VARIABLE FOR CI TYPE TO SCFOPT COMMON
C 30 OCT 95 - MWS - REMOVE MORE FTNCHEK ERRORS
C  6 APR 95 - MWS - MOLDAT: USE FORMAT STATEMENTS
C 11 NOV 94 - MWS - REMOVE MOST FTNCHEK ERRORS
C 10 AUG 94 - MWS - INCREASE NUMBER OF DAF RECORDS
C 21 JUL 94 - MWS - ADD DUMMY READA ROUTINE
C  9 DEC 93 - MWS - MPCHCK: CHANGE DAF RECORD OF OVERLAP INTEGRALS
C 23 JUN 93 - MWS - SETUPG: WRITE ERROR MESSAGE BEFORE ABRT CALL
C 13 APR 92 - MWS - MOLDAT: DELETE UHF ELECTRON/MULTIPLICITY TESTING.
C 15 MAR 92 - JHJ - NEW MODULE FOR MOPAC MOLECULE SETUP
C
C         THE FOLLOWING STATEMENT APPEARED IN MOPAC 6.0 AND
C         IS REPRODUCED HERE TO COMPLY WITH SAID STATEMENT.
C
C         NOTICE OF PUBLIC DOMAIN NATURE OF THIS PROGRAM
C
C      'THIS COMPUTER PROGRAM IS A WORK OF THE UNITED STATES
C       GOVERNMENT AND AS SUCH IS NOT SUBJECT TO PROTECTION BY
C       COPYRIGHT (17 U.S.C. # 105.)  ANY PERSON WHO FRAUDULENTLY
C       PLACES A COPYRIGHT NOTICE OR DOES ANY OTHER ACT CONTRARY
C       TO THE PROVISIONS OF 17 U.S. CODE 506(C) SHALL BE SUBJECT
C       TO THE PENALTIES PROVIDED THEREIN.  THIS NOTICE SHALL NOT
C       BE ALTERED OR REMOVED FROM THIS SOFTWARE AND IS TO BE ON
C       ALL REPRODUCTIONS.'
C
C*MODULE MPCMOL  *DECK MOLDAT
      SUBROUTINE MOLDAT(MODE)
C
C        ARGUMENT IS AN ARTIFACT, MODE=0 IN THIS ROUTINE'S ONLY CALL
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      CHARACTER*241 KEYWRD, OLDE(20)*6, ALLREF*80
      LOGICAL UHFFLG,MINDO3,LPM3,PEPTID,
     *        LMNDO,HALFE,SLOW,GOPARR,DSKWRK,MASWRK
C
      PARAMETER (MXATM=2000, MAXORB=4*MXATM)
C
      DIMENSION ISWAP(2,20), ESTORE(107)
C
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /INFOA / NNAT,ICH,MUL,NUM,NQMT,NELEC,NAE,NBE,
     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /WFNOPT/ SCFTYP,CITYP,DFTYPE,CCTYP,MPLEVL,MPCTYP
C
      COMMON /NATYPE/ NZTYPE(107),MTYPE(30),LTYPE
      COMMON /GEOKST/ NATOMS,LABELS(MXATM),
     1                NA(MXATM),NB(MXATM),NC(MXATM)
      COMMON /MOLMEC/ HTYPE(4),NHCO(4,200),NNHCO,ITYPE
      COMMON /MPCGEO/ GEO(3,MXATM)
      COMMON /N2ELCT/ N2EL
      COMMON /MOLKST/ NUMAT,NAT(MXATM),NFIRST(MXATM),NMIDLE(MXATM),
     1                NLAST(MXATM), NORBS, NELECS,NALPHA,NBETA,
     2                NCLOSE,NOPEN,NDUMY,FRACT
      COMMON /KEYWRD/ KEYWRD
      COMMON /NATORB/ NATORB(107)
      COMMON /CORE  / CORE(107)
      COMMON /BETAS / BETAS(107),BETAP(107),BETAD(107)
      COMMON /MOLORB/ USPD(MAXORB),PSPD(MAXORB)
      COMMON /ONELEC/ USS(107),UPP(107),UDD(107)
      COMMON /ATHEAT/ ATHEAT
C--   COMMON /POLVOL/ POLVOL(107)
      COMMON /MULTIP/ DD(107),QQ(107),AM(107),AD(107),AQ(107)
      COMMON /TWOELE/ GSS(107),GSP(107),GPP(107),GP2(107),HSP(107),
     *                GSD(107),GPD(107),GDD(107)
      COMMON /IDEAA / GUESA1(107,10),GUESA2(107,10),GUESA3(107,10)
      COMMON /IDEAS / GUESS1(107,10),GUESS2(107,10),GUESS3(107,10)
      COMMON /IDEAP / GUESP1(107,10),GUESP2(107,10),GUESP3(107,10)
      COMMON /ALPHA / ALP(107)
      COMMON /REFS/ ALLREF(107,4)
      COMMON /MNDO/  USSM(107), UPPM(107), UDDM(107), ZSM(107),
     1ZPM(107), ZDM(107), BETASM(107), BETAPM(107), BETADM(107),
     2ALPM(107), EISOLM(107), DDM(107), QQM(107), AMM(107),
     3ADM(107), AQM(107), GSSM(107), GSPM(107), GPPM(107),
     4GP2M(107), HSPM(107), POLVOM(107)
      COMMON /PM3 /  USSPM3(107), UPPPM3(107), UDDPM3(107), ZSPM3(107),
     1ZPPM3(107), ZDPM3(107), BETASP(107), BETAPP(107), BETADP(107),
     2ALPPM3(107), EISOLP(107), DDPM3(107), QQPM3(107), AMPM3(107),
     3ADPM3(107), AQPM3(107) ,GSSPM3(107), GSPPM3(107), GPPPM3(107),
     4GP2PM3(107), HSPPM3(107),POLVOP(107)
      COMMON /AM1BLO/USSAM1(107), UPPAM1(107), UDDAM1(107), ZSAM1(107),
     1ZPAM1(107), ZDAM1(107), BETASA(107), BETAPA(107), BETADA(107),
     2ALPAM1(107), EISOLA(107), DDAM1(107), QQAM1(107), AMAM1(107),
     3ADAM1(107), AQAM1(107) ,GSSAM1(107), GSPAM1(107), GPPAM1(107),
     4GP2AM1(107), HSPAM1(107),POLVOA(107)
C
C  COMMON BLOCKS FOR MINDO/3
C
      COMMON /ONELE3 /  USS3(18),UPP3(18)
      COMMON /ATOMI3 /  EISOL3(18),EHEAT3(18)
      COMMON /EXPON3 /  ZS3(18),ZP3(18)
C
C  END OF MINDO/3 COMMON BLOCKS
C
      COMMON /EXPONT/ ZS(107),ZP(107),ZD(107),ZF(107),ZG(107)
      COMMON /ATOMIC/ EISOL(107),EHEAT(107)
C
      PARAMETER (NNAM=1)
      DIMENSION QNAM(NNAM),KQNAM(NNAM)
      DATA GNAM/8HMOPAC   /
      DATA QNAM/8HPEPTID  /
      DATA KQNAM/0/
C
      DATA ESTORE(1)/0.0D+00/
C
      DATA KMINDO,KMNDO,KAM1,KPM3/4HMIND,4HMNDO,4HAM1 ,4HPM3 /
      DATA UHF,ROHF/8HUHF     ,8HROHF    /
C
C         THE KEYWORD ARRAY OF MOPAC IS COMPLETELY DEACTIVATED,
C         BUT A FEW PLACES STILL LOOK FOR OBSCURE OPTIONS HERE.
C         BLANK IT TO MAKE SURE NOTHING EXOTIC IS EVER SELECTED.
C
      DO I=1,241
         KEYWRD(I:I) = ' '
      ENDDO
C
C         ONLY INPUT OPTION ALLOWS TURNING OFF THE PEPTIDE FUDGE
C
      PEPTID=.TRUE.
      CALL NAMEIO(IR,JRET,GNAM,NNAM,QNAM,KQNAM,
     *            PEPTID,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,0,0,0,   0,0,0,0,0,   0,0,0,0,0)
      IF(JRET.EQ.2) THEN
         IF (MASWRK) WRITE (IW,*) 'ERROR READING $MOPAC GROUP'
         CALL ABRT
      END IF
C
      IF(ESTORE(1).EQ.0.0D+00) THEN
      DO 9 I=1,107
   9  ESTORE(I)=EHEAT(I)
      END IF
      DO 8 I=1,107
  8   EHEAT(I)=ESTORE(I)
C
      LPM3   = MPCTYP.EQ.KPM3
      MINDO3 = MPCTYP.EQ.KMINDO
      LMNDO  = MPCTYP.EQ.KMNDO  .OR.  MPCTYP.EQ.KMINDO
      UHFFLG = SCFTYP.EQ.UHF  .OR.  SCFTYP.EQ.ROHF
C
      KHARGE=0
      KHARGE = ICH
      ELECS=-KHARGE
      NDORBS=0
      ATHEAT=0.0D+00
      EAT=0.0D+00
      NUMAT=0
      IF(MODE.EQ.1) GO TO 80
C
C    SWITCH IN MNDO PARAMETERS
C
C
C       ZERO OUT GAUSSIAN 1 FOR CARBON.  THIS WILL BE USED IN
C       ROTATE TO DECIDE WHETHER OR NOT TO USE AM1-TYPE GAUSSIANS
C
      IF ( LMNDO ) THEN
         GUESS1(6,1)=0.0D+00
         DO 10 I=1,107
C--         IF(.NOT.MINDO3) POLVOL(I)=POLVOM(I)
            ZS(I)=ZSM(I)
            ZP(I)=ZPM(I)
            ZD(I)=ZDM(I)
            USS(I)=USSM(I)
            UPP(I)=UPPM(I)
            UDD(I)=UDDM(I)
            BETAS(I)=BETASM(I)
            BETAP(I)=BETAPM(I)
            BETAD(I)=BETADM(I)
            ALP(I)=ALPM(I)
            EISOL(I)=EISOLM(I)
            DD(I)=DDM(I)
            QQ(I)=QQM(I)
            AM(I)=AMM(I)
            AD(I)=ADM(I)
            AQ(I)=AQM(I)
            GSS(I)=GSSM(I)
            GPP(I)=GPPM(I)
            GSP(I)=GSPM(I)
            GP2(I)=GP2M(I)
            HSP(I)=HSPM(I)
   10    CONTINUE
C
C    SWITCH IN MNDO-PM3 PARAMETERS
C
      ELSE IF( LPM3 ) THEN
         DO 30 I=1,107
            DO 20 J=1,10
               GUESS1(I,J)=GUESP1(I,J)
               GUESS2(I,J)=GUESP2(I,J)
   20       GUESS3(I,J)=GUESP3(I,J)
C--         POLVOL(I)=POLVOP(I)
            ZS(I)=ZSPM3(I)
            ZP(I)=ZPPM3(I)
            ZD(I)=ZDPM3(I)
            USS(I)=USSPM3(I)
            UPP(I)=UPPPM3(I)
            UDD(I)=UDDPM3(I)
            BETAS(I)=BETASP(I)
            BETAP(I)=BETAPP(I)
            BETAD(I)=BETADP(I)
            ALP(I)=ALPPM3(I)
            EISOL(I)=EISOLP(I)
            DD(I)=DDPM3(I)
            QQ(I)=QQPM3(I)
            AM(I)=AMPM3(I)
            AD(I)=ADPM3(I)
            AQ(I)=AQPM3(I)
            GSS(I)=GSSPM3(I)
            GPP(I)=GPPPM3(I)
            GSP(I)=GSPPM3(I)
            GP2(I)=GP2PM3(I)
            HSP(I)=HSPPM3(I)
   30    CONTINUE
C
C    SWITCH IN AM1 PARAMETERS
C
      ELSE
         DO 50 I=1,107
            DO J=1,10
               GUESS1(I,J)=GUESA1(I,J)
               GUESS2(I,J)=GUESA2(I,J)
               GUESS3(I,J)=GUESA3(I,J)
            ENDDO
C--         POLVOL(I)=POLVOA(I)
            ZS(I)=ZSAM1(I)
            ZP(I)=ZPAM1(I)
            ZD(I)=ZDAM1(I)
            USS(I)=USSAM1(I)
            UPP(I)=UPPAM1(I)
            UDD(I)=UDDAM1(I)
            BETAS(I)=BETASA(I)
            BETAP(I)=BETAPA(I)
            BETAD(I)=BETADA(I)
            ALP(I)=ALPAM1(I)
            EISOL(I)=EISOLA(I)
            DD(I)=DDAM1(I)
            QQ(I)=QQAM1(I)
            AM(I)=AMAM1(I)
            AD(I)=ADAM1(I)
            AQ(I)=AQAM1(I)
            GSS(I)=GSSAM1(I)
            GPP(I)=GPPAM1(I)
            GSP(I)=GSPAM1(I)
            GP2(I)=GP2AM1(I)
            HSP(I)=HSPAM1(I)
   50    CONTINUE
      END IF
C
C        SWAP IN OLD PARAMETERS FOR ELEMENTS.  OLDE CONTAINS THE
C        CHARACTER NAME OF THE ELEMENT, AND ISWAP(1,1:NEWELE) CONTAINS
C        THE ATOMIC NUMBER OF THE ELEMENT. ISWAP(2,1:NEWELE) CONTAINS
C        THE STORAGE ADDRESS OF THE OLD SET OF PARAMETERS.
C
      NEWELE=2
      OLDE(1)=' S1978'
      ISWAP(1,1)=16
      ISWAP(2,1)=91
      OLDE(2)='SI1978'
      ISWAP(1,2)=14
      ISWAP(2,2)=90
C$DOIT ASIS
      DO 60 K=1,NEWELE
         IF(INDEX(KEYWRD,OLDE(K)).NE.0) THEN
            I=ISWAP(1,K)
            J=ISWAP(2,K)
            ALLREF(I,3)=ALLREF(J,1)
            ALLREF(I,1)=ALLREF(J,1)
            ZS(I)=ZS(J)
            ZP(I)=ZP(J)
            ZD(I)=ZD(J)
            USS(I)=USS(J)
            UPP(I)=UPP(J)
            UDD(I)=UDD(J)
            BETAS(I)=BETAS(J)
            BETAP(I)=BETAP(J)
            BETAD(I)=BETAD(J)
            ALP(I)=ALP(J)
            EISOL(I)=EISOL(J)
            DD(I)=DD(J)
            QQ(I)=QQ(J)
            AM(I)=AM(J)
            AD(I)=AD(J)
            AQ(I)=AQ(J)
            IF(GSS(J).NE.0)GSS(I)=GSS(J)
            IF(GPP(J).NE.0)GPP(I)=GPP(J)
            IF(GSP(J).NE.0)GSP(I)=GSP(J)
            IF(GP2(J).NE.0)GP2(I)=GP2(J)
            IF(HSP(J).NE.0)HSP(I)=HSP(J)
         END IF
   60 CONTINUE
      IF( MINDO3 ) THEN
         DO 70 I=1,17
            IF(I.NE.2.AND.I.NE.10) THEN
               USS(I)=USS3(I)
               UPP(I)=UPP3(I)
               EISOL(I)=EISOL3(I)
               EHEAT(I)=EHEAT3(I)
               ZS(I)=ZS3(I)
               ZP(I)=ZP3(I)
               GSS(I)=GSSM(I)
               GPP(I)=GPPM(I)
               GSP(I)=GSPM(I)
               GP2(I)=GP2M(I)
               HSP(I)=HSPM(I)
            END IF
   70    CONTINUE
      END IF
   80 CONTINUE
      IF(USS(1) .GT. -1.0D+00) THEN
         IF(MASWRK) WRITE(IW,9030)
         CALL ABRT
         STOP
      END IF
      IA=1
      IB=0
      NHEAVY=0
      DO 130 II=1,NATOMS
         IF(LABELS(II).EQ.99.OR.LABELS(II).EQ.107) GO TO 130
         NUMAT=NUMAT+1
         NAT(NUMAT)=LABELS(II)
         NFIRST(NUMAT)=IA
         NI=NAT(NUMAT)
         ATHEAT=ATHEAT+EHEAT(NI)
         EAT   =EAT   +EISOL(NI)
         ELECS=ELECS+CORE(NI)
         IB=IA+NATORB(NI)-1
         NMIDLE(NUMAT)=IB
         IF(NATORB(NI).EQ.9)NDORBS=NDORBS+5
         IF(NATORB(NI).EQ.9)NMIDLE(NUMAT)=IA+3
         NLAST(NUMAT)=IB
         IF(IA.GT.MAXORB) GO TO 810
         USPD(IA)=USS(NI)
         IF(IA.EQ.IB) GO TO 120
         K=IA+1
         K1=IA+3
C$DOIT ASIS
         DO 90 J=K,K1
            IF(J.GT.MAXORB) GO TO 810
            USPD(J)=UPP(NI)
   90    CONTINUE
         NHEAVY=NHEAVY+1
         IF(K1.EQ.IB) GO TO 120
         K=K1+1
C$DOIT ASIS
         DO 110 J=K,IB
  110    USPD(J)=UDD(NI)
  120    CONTINUE
  130 IA=IB+1
C
C     ---- FILL MTYPE AND NZTYPE ARRAYS -----
C
      LTYPE=0
      DO 133 I=1,NUMAT
         DO 131 J=1,LTYPE
  131    IF(NAT(I).EQ.MTYPE(J)) GO TO 132
         LTYPE=LTYPE+1
         MTYPE(LTYPE)=NAT(I)
         NZTYPE(NAT(I))=LTYPE
C
C       LTYPE = NUMBER OF TYPES OF REAL ATOM PRESENT
C       MTYPE = TYPES OF REAL ATOMS PRESENT
         J=LTYPE
  132    CONTINUE
  133 CONTINUE
C
C   PRINT THE APPROPRIATE REFERENCES FOR THE PARAMETERS
C
      IF(MODE.NE.1) CALL REFER
C
      ATHEAT=ATHEAT-EAT*23.061D+00
      NORBS=NLAST(NUMAT)
      IF(NORBS.GT.MAXORB .AND. MASWRK) THEN
         WRITE(IW,9020) MAXORB,NORBS
         CALL ABRT
         STOP
      END IF
C
      NLIGHT=NUMAT-NHEAVY
      N2ELEC = 50*MXATM*(MXATM-1)
      N2EL = 50*NHEAVY*(NHEAVY-1) + 10*NHEAVY*NLIGHT +
     *         (NLIGHT*(NLIGHT-1))/2
      IF(MASWRK) THEN
         WRITE(IW,9000) NHEAVY,NLIGHT,N2EL
         IF(N2EL.GT.N2ELEC) THEN
            WRITE(IW,9010) N2ELEC,N2EL
            CALL ABRT
            STOP
         END IF
      END IF
C
C NOW TO WORK OUT HOW MANY ELECTRONS ARE IN EACH TYPE OF SHELL
C
      NALPHA=0
      NBETA=0
C
C      PROTECT DUMB USERS FROM DUMB ERRORS!
C
      NELECS=INT(MAX(ELECS+0.5D+00,0.0D+00))
      NELECS=MIN(2*NORBS,NELECS)
      NCLOSE=0
      NOPEN=0
C
C UHF IS DEFAULTS FOR ODD NUMBER OF ELECTRONS
C
      IF( UHFFLG ) THEN
         NALPHA= NAE
         NBETA = NBE
         FRACT=1.0D+00
      ELSE
C
C   NOW TO DETERMINE OPEN AND CLOSED SHELLS
C
         IELEC = MUL - 1
         ILEVEL = MUL - 1
C
         I=INDEX(KEYWRD,'OPEN(')
         IF(I.NE.0) THEN
            IELEC = INT(READA(KEYWRD,I))
            ILEVEL= INT(READA(KEYWRD,I+7))
         END IF
         NCLOSE=NELECS/2
         NOPEN = NELECS-NCLOSE*2
         IF( IELEC.NE.0 ) THEN
            IF((NELECS/2)*2.EQ.NELECS .NEQV.
     1                  (IELEC/2)*2.EQ.IELEC) THEN
               CALL ABRT
               STOP
            END IF
            NCLOSE=NCLOSE-IELEC/2
            NOPEN=ILEVEL
            IF(NCLOSE+NOPEN.GT.NORBS) THEN
               CALL ABRT
               STOP
            END IF
            FRACT=IELEC*1.0D+00/ILEVEL
         END IF
         IF(INDEX(KEYWRD,'C.I.=(').NE.0) THEN
            I=INT(READA(KEYWRD,INDEX(KEYWRD,'C.I.=(')+5))
     1       -INT(READA(KEYWRD,INDEX(KEYWRD,'C.I.=(')+7))
            IF(NOPEN.GT.I) THEN
               CALL ABRT
               STOP
            END IF
         END IF
         IF(INDEX(KEYWRD,'C.I.').NE.0.AND.NOPEN.EQ.0) THEN
            NOPEN=1
            NCLOSE=NCLOSE-1
            FRACT=2.0D+00
         END IF
         NOPEN=NOPEN+NCLOSE
      END IF
C
C  WORK OUT IF DEFINED SPIN-STATE ALLOWED
C
      MSDEL=INDEX(KEYWRD,' MS')
      IF(MSDEL.NE.0) THEN
         MSDEL=INT(1.0001D+00*READA(KEYWRD,INDEX(KEYWRD,' MS')))
      ELSE
         IF(INDEX(KEYWRD,'TRIP')+INDEX(KEYWRD,'QUAR').GT.0)MSDEL=1
         IF(INDEX(KEYWRD,'QUIN')+INDEX(KEYWRD,'SEXT').GT.0)MSDEL=2
      END IF
      IF(MSDEL.NE.0.AND..NOT.UHFFLG) THEN
C
C   MSDEL = NUMBER OF ALPHA ELECTRONS - NUMBER OF BETA ELECTRONS
C
         NDOUBL=99
         IF(INDEX(KEYWRD,'C.I.=(').NE.0) THEN
            NDOUBL=INT(READA(KEYWRD,INDEX(KEYWRD,'C.I.=(')+7))
            NMOS=INT(READA(KEYWRD,INDEX(KEYWRD,'C.I.=(')+5))
         ELSE IF (INDEX(KEYWRD,'C.I.=').NE.0) THEN
            NMOS=INT(READA(KEYWRD,INDEX(KEYWRD,'C.I.=')+5))
         ELSE
            NMOS=NOPEN-NCLOSE
         END IF
         IF(NDOUBL.EQ.99) THEN
            J=MAX(MIN((NCLOSE+NOPEN+1)/2-(NMOS-1)/2,NORBS-NMOS+1),1)
         ELSE
            J=NCLOSE-NDOUBL+1
         END IF
         ELS=MAX(0.0D+00,(NCLOSE-J+1.0D+00))*2.0D+00+
     1       MAX(0.0D+00,(NOPEN-NCLOSE)*FRACT) + 0.5D+00
         NE = INT(ELS)
         NUPP=(NE+1)/2+MSDEL
         NDOWN=NE-NUPP
C
C  NUPP  = NUMBER OF ALPHA ELECTRONS IN ACTIVE SPACE
C  NDOWN = NUMBER OF BETA  ELECTRONS IN ACTIVE SPACE
C
         IF(NUPP*NDOWN.LT.0.OR.NUPP.GT.NMOS.OR.NDOWN.GT.NMOS) THEN
C           WRITE(IW,'(A)')
C    1' SPECIFIED SPIN COMPONENT NOT SPANNED BY ACTIVE SPACE'
            CALL ABRT
            STOP
         END IF
      END IF
C#      WRITE(IW,'(''  NOPEN,NCLOSE,NALPHA,NBETA,FRACT'',4I4,F12.5)')
C#     1 NOPEN, NCLOSE, NLAPHA, NBETA, FRACT
C
C   MAKE SURE ANALYT IS NOT USED WITH ANALYTICAL C.I. DERIVATIVES
C
      HALFE = (NOPEN.GT.NCLOSE.AND.FRACT.NE.2.0D+00.AND.FRACT.NE.0.0D+00
     1 .OR. INDEX(KEYWRD,'C.I.').NE.0)
      SLOW=(INDEX(KEYWRD,'EXCI').NE.0.OR.
     1INDEX(KEYWRD,'ROOT').NE.0.AND.INDEX(KEYWRD,'ROOT=1').EQ.0)
      IF(HALFE)HALFE=(.NOT.SLOW)
      IF(INDEX(KEYWRD,'NOANCI').EQ.0.AND.
     * INDEX(KEYWRD,'ANALYT').NE.0.AND.HALFE) THEN
C        WRITE(IW,*)
C        WRITE(IW,'(A)')' KEYWORD ''ANALYT'' CANNOT BE USED HERE: ',
C    1' ANALYICAL C.I. DERIVATIVES MUST USE FINITE DIFFERENCES',
C    2' TO CORRECT, REMOVE KEYWORD ''ANALYT'' OR ADD ''NOANCI'''
         CALL ABRT
         STOP
      END IF
      YY=KHARGE
      YY=YY/NORBS
      L=0
      DO 160 I=1,NUMAT
         NI=NAT(I)
         XX=1.0D+00/(NLAST(I)-NFIRST(I)+1+1.0D-10)
         W=CORE(NI)*XX-YY
         IA=NFIRST(I)
         IC=NMIDLE(I)
         IB=NLAST(I)
C$DOIT ASIS
         DO 140 J=IA,IC
            L=L+1
  140    PSPD(L)=W
C$DOIT ASIS
         DO 150 J=IC+1,IB
            L=L+1
  150    PSPD(L)=0.0D+00
  160 CONTINUE
C
C     SET UP MOLECULAR-MECHANICS CORRECTION TO -(C=O)-(NH)- LINKAGE
C
      NNHCO=0
      IF(.NOT.PEPTID) GO TO 300
C
      HTYPE(1)=6.1737D+00
      HTYPE(2)=3.3191D+00
      HTYPE(3)=7.1853D+00
      HTYPE(4)=1.7712D+00
C
                           ITYPE=1
      IF(MPCTYP.EQ.KAM1)   ITYPE=2
      IF(MPCTYP.EQ.KPM3)   ITYPE=3
      IF(MPCTYP.EQ.KMINDO) ITYPE=4
C
C     IDENTIFY O=C-N-H SYSTEMS VIA THE INTERATOMIC DISTANCES MATRIX
C     NOTE: GEO(*,*) IS CARTESIAN, NOT INTERNAL
C
      DO 220 I=1,NUMAT
         IF(NAT(I).NE.8) GO TO 220
         DO 210 J=1,NUMAT
            IF(NAT(J).NE.6) GO TO 210
            IF(ABS(GEO(1,I)-GEO(1,J)).GT.1.3D+00 .OR.
     *         ABS(GEO(2,I)-GEO(2,J)).GT.1.3D+00 .OR.
     *         ABS(GEO(3,I)-GEO(3,J)).GT.1.3D+00) GO TO 210
            RIJ=SQRT((GEO(1,I)-GEO(1,J))**2 +
     *               (GEO(2,I)-GEO(2,J))**2 +
     *               (GEO(3,I)-GEO(3,J))**2  )
            IF(RIJ.GT.1.3D+00) GO TO 210
            DO 200 K=1,NUMAT
               IF(NAT(K).NE.7) GO TO 200
               IF(ABS(GEO(1,J)-GEO(1,K)).GT.1.6D+00 .OR.
     *            ABS(GEO(2,J)-GEO(2,K)).GT.1.6D+00 .OR.
     *            ABS(GEO(3,J)-GEO(3,K)).GT.1.6D+00) GO TO 200
               RJK=SQRT((GEO(1,J)-GEO(1,K))**2 +
     *                  (GEO(2,J)-GEO(2,K))**2 +
     *                  (GEO(3,J)-GEO(3,K))**2  )
               IF(RJK.GT.1.6D+00) GO TO 200
               DO 190 L=1,NUMAT
                  IF(NAT(L).NE.1) GO TO 190
                  IF(ABS(GEO(1,K)-GEO(1,L)).GT.1.3D+00 .OR.
     *               ABS(GEO(2,K)-GEO(2,L)).GT.1.3D+00 .OR.
     *               ABS(GEO(3,K)-GEO(3,L)).GT.1.3D+00) GO TO 190
                  RKL=SQRT((GEO(1,K)-GEO(1,L))**2 +
     *                     (GEO(2,K)-GEO(2,L))**2 +
     *                     (GEO(3,K)-GEO(3,L))**2  )
                  IF(RKL.GT.1.3D+00) GO TO 190
C
C     WE HAVE A H-N-C=O SYSTEM.  THE ATOM NUMBERS ARE L-K-J-I
C     NOW SEARCH OUT ATOM ATTACHED TO NITROGEN, THIS SPECIFIES
C     THE SYSTEM X-N-C=O
C
                  DO 180 M=1,NUMAT
                     IF(M.EQ.K.OR.M.EQ.L.OR.M.EQ.J) GO TO 180
                     IF(ABS(GEO(1,K)-GEO(1,M)).GT.1.7D+00 .OR.
     *                  ABS(GEO(2,K)-GEO(2,M)).GT.1.7D+00 .OR.
     *                  ABS(GEO(3,K)-GEO(3,M)).GT.1.7D+00) GO TO 180
                     RKM=SQRT((GEO(1,K)-GEO(1,M))**2 +
     *                        (GEO(2,K)-GEO(2,M))**2 +
     *                        (GEO(3,K)-GEO(3,M))**2  )
                     IF(RKM.GT.1.7D+00) GO TO 180
C
C     PEPTIDE BOND FOUND.  TEST ON DIMENSION IN /MOLMEC/ EXCEEDED.
C
                     NNHCO=NNHCO+1
                     IF(NNHCO.GT.200) THEN
                        WRITE(6,*) 'TOO MANY PEPTIDE LINKAGES FOUND'
                        CALL ABRT
                     END IF
C
                     NHCO(1,NNHCO)=I
                     NHCO(2,NNHCO)=J
                     NHCO(3,NNHCO)=K
                     NHCO(4,NNHCO)=M
                     NNHCO=NNHCO+1
                     NHCO(1,NNHCO)=I
                     NHCO(2,NNHCO)=J
                     NHCO(3,NNHCO)=K
                     NHCO(4,NNHCO)=L
                     GO TO 190
  180             CONTINUE
  190          CONTINUE
  200       CONTINUE
  210    CONTINUE
  220 CONTINUE
C
      IF(MODE.NE.1.AND.NNHCO.NE.0) THEN
         WRITE(IW,9230) NNHCO
         DO I=1,NNHCO
            IF(NAT(NHCO(4,I)).EQ.1) THEN
               WRITE(IW,9235) (NHCO(K,I),K=1,4)
            ELSE
               WRITE(IW,9236) (NHCO(K,I),K=1,4)
            END IF
         ENDDO
      END IF
C
  300 CONTINUE
      RETURN
C
  810 CONTINUE
      IF(MASWRK) THEN
         WRITE(IW,*) 'MAXIMUM NO. OF ATOMIC ORBITALS EXCEEDED'
      END IF
      CALL ABRT
      STOP
C
 9000 FORMAT(/1X,'THERE ARE',I5,' HEAVY AND',I5,' LIGHT ATOMS,'/
     *        1X,'YIELDING A TOTAL OF',I10,' MOPAC 2E- INTEGRALS.')
 9010 FORMAT(//10X,'**** MAX. NUMBER OF TWO-ELECTRON INTEGRALS:',I8/
     *         10X,'NUMBER OF TWO ELECTRON INTEGRALS IN SYSTEM:',I8)
 9020 FORMAT(//10X,'**** MAX. NUMBER OF ORBITALS:',I4/
     *            10X,'NUMBER OF ORBITALS IN SYSTEM:',I4)
 9030 FORMAT(1X,'THE HAMILTONIAN REQUESTED IS NOT AVAILABLE IN',
     * ' THIS PROGRAM')
 9230 FORMAT(1X,'MOLECULAR MECHANICS CORRECTION APPLIED TO',I5,
     *          ' PEPTIDE LINKAGES:')
 9235 FORMAT(5X,'O=',I4,' C=',I4,' N=',I4,' H=',I4)
 9236 FORMAT(5X,'O=',I4,' C=',I4,' N=',I4,' X=',I4)
      END
C*MODULE MPCMOL  *DECK MPCHCK
      SUBROUTINE MPCHCK(H,EIG,SCR,VEC,SCFTYP,DBG,IWRK,L1,L2)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      LOGICAL DBG,GOPARR,DSKWRK,MASWRK
C
      PARAMETER (MXATM=2000, TWO=2.0D+00)
C
      DIMENSION H(L2),EIG(L1),SCR(L1,8),VEC(L1,L1),IWRK(L1)
      DIMENSION ENEG(2,86)
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 /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
C
      DATA ((ENEG(I,NUC),I=1,2),NUC=1,36) /
     *           0.500D+00,0.0D+00,   0.918D+00,0.0D+00,
     *           0.196D+00,0.05D+00,  0.309D+00,0.05D+00,
     *           0.495D+00,0.310D+00, 0.706D+00,0.433D+00,
     *           0.945D+00,0.568D+00, 1.244D+00,0.632D+00,
     *           1.573D+00,0.730D+00, 1.93D+00,0.850D+00,
     *           0.182D+00,0.05D+00,  0.253D+00,0.05D+00,
     *           0.393D+00,0.210D+00, 0.540D+00,0.297D+00,
     *           0.696D+00,0.392D+00, 0.880D+00,0.437D+00,
     *           1.073D+00,0.506D+00, 1.278D+00,0.591D+00,
     *           0.147D+00,0.050D+00, 0.196D+00,0.050D+00,
     *           18*0.0D+00,          0.293D+00,0.10D+00,
     *           0.424D+00,0.208D+00, 0.553D+00,0.287D+00,
     *           0.686D+00,0.369D+00, 0.838D+00,0.403D+00,
     *           0.993D+00,0.457D+00, 1.153D+00,0.524D+00/
      DATA ((ENEG(I,NUC),I=1,2),NUC=37,86) /
     *           0.138D+00,0.050D+00, 0.178D+00,0.050D+00,
     *           18*0.0D+00,          0.265D+00,0.10D+00,
     *           0.372D+00,0.197D+00, 0.476D+00,0.265D+00,
     *           0.582D+00,0.335D+00, 0.701D+00,0.360D+00,
     *           0.821D+00,0.403D+00, 0.944D+00,0.457D+00,
     *           0.124D+00,0.05D+00, 0.1575D+00,0.05D+00,
     *           46*0.0D+00,          0.2610D+00,0.10D+00,
     *           0.261D+00,0.192D+00, 0.459D+00,0.2398D+00,
     *           0.558D+00,0.286D+00, 0.660D+00,0.332D+00,
     *           0.765D+00,0.380D+00, 0.874D+00,0.428D+00/
C
      DATA UHF/8HUHF     /
C
C     ----- SET UP HUCKEL GUESS FOR MOPAC WAVEFUNCTIONS,
C           WRITTEN BY MIKE SCHMIDT, FEBRUARY 11, 1992 -----
C
C           READ OVERLAP MATRIX INTO HUCKEL OPERATOR
C
      CALL DAREAD(IDAF,IODA,H,L2,52,0)
C
C        CONSTRUCT HUCKEL MATRIX -H-, ON TOP OF THE OVERLAP MATRIX
C                 SET VALENCE ORBITAL ENERGIES ONLY.
C
      I0=0
      DO 390 N=1,NAT
         NUCZ = INT(ZAN(N))
         NVAL=1
         IF(NUCZ.GT.2) NVAL=2
         DO 360 J=1,NVAL
            IVAL = 1
            IF(J.EQ.2) IVAL=3
            DO 350 I=1,IVAL
               I0= I0+1
               II= (I0*I0+I0)/2
               H(II)=-ENEG(J,NUCZ)
  350       CONTINUE
  360    CONTINUE
  390 CONTINUE
C
C           FINISH UP THE OFF-DIAGONAL OF THE HUCKEL OPERATOR
C
      FUDGE = 1.75D+00/TWO
      II = 1
      IJ = 1
      DO 720 I=2,L1
         II = II + I
         HII=H(II)
         JJ = 0
         DO 710 J=1,I-1
            JJ = JJ + J
            IJ = IJ + 1
            SIJ = H(IJ)
            H(IJ)=FUDGE*SIJ*(HII+H(JJ))
  710    CONTINUE
         IJ = IJ + 1
  720 CONTINUE
C
      IF(DBG) THEN
         WRITE(IW,*) 'HUCKEL OPERATOR BEFORE DIAGONALIZATION'
         CALL PRTRI(H,L1)
      END IF
C
C           DIAGONALIZE THE HUCKEL OPERATOR
C
      IERR = 0
      CALL GLDIAG(L1,L1,L1,H,SCR,EIG,VEC,IERR,IWRK)
      IF(IERR.NE.0) THEN
         IF (MASWRK) WRITE(IW,*)
     *               'UNABLE TO DIAGONALIZE HUCKEL MATRIX'
         CALL ABRT
      END IF
C
      IF(DBG) THEN
         WRITE(IW,*) 'EXTENDED HUCKEL ORBITALS (IN MINI BASIS) ARE'
         CALL PREVNL(VEC,EIG,L1,L1,L1)
      END IF
C
      CALL DAWRIT(IDAF,IODA,VEC,L1*L1,15,0)
      IF(SCFTYP.EQ.UHF) CALL DAWRIT(IDAF,IODA,VEC,L1*L1,19,0)
C
      RETURN
      END
C*MODULE MPCMOL  *DECK OVRLAP
      SUBROUTINE OVRLAP(H,S,NORBS)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER (MXATM=2000, MAXORB=4*MXATM)
      DIMENSION H(*),S(*)
C
      DIMENSION IFACT(MAXORB)
C
      COMMON /MOLKST/ NUMAT,NAT(MXATM),NFIRST(MXATM),NMIDLE(MXATM),
     *                NLAST(MXATM), NORBX, NELECS,NALPHA,NBETA,
     *                NCLOSE,NOPEN,NDUMY,FRACT
      COMMON /BETAS / BETAS(107),BETAP(107),BETAD(107)
C
C THIS IS WHAT REMAINS OF THE MULLIK SUBROUTINE.  IT NOW CALCULATES
C OVERLAP MATRIX FROM THE ONE-ELECTRON MATRIX.
C
      DO 10 I=1,NORBS
   10 IFACT(I)=(I*(I-1))/2
      IFACT(NORBS+1)=(NORBS*(NORBS+1))/2
      DO 50 I=1,NUMAT
         IF=NFIRST(I)
         IL=NLAST(I)
         IM1=I-1
         BI=BETAS(NAT(I))
         DO 50 K=IF,IL
            II=(K*(K-1))/2
            DO 30 J=1,IM1
               JF=NFIRST(J)
               JL=NLAST(J)
               BJ=BETAS(NAT(J))
               DO 20 JJ=JF,JL
                  IJ=II+JJ
                  S(IJ)=2.0D+00*H(IJ)/(BI+BJ)     +1.0D-14
C  THE  +1.0D-14 IS TO PREVENT POSSIBLE ERRORS IN THE DIAGONALIZATION.
C                  STORE(IJ)=H(IJ)
   20          BJ=BETAP(NAT(J))
   30       CONTINUE
            DO 40 JJ=IF,K
               IJ=II+JJ
C               STORE(IJ)=0.0D+00
   40       S(IJ)=0.0D+00
   50 BI=BETAP(NAT(I))
      DO 60 I=1,NORBS
C         STORE(IFACT(I+1))=1.0D+00
   60 S(IFACT(I+1))=1.0D+00
C
      RETURN
      END
C*MODULE MPCMOL  *DECK READA
      DOUBLE PRECISION FUNCTION READA(STRING,ISTART)
      CHARACTER STRING*(*)
C         NOTE THAT THIS IS NEVER ACTUALLY INVOKED
      READA=1
      IF(STRING(1:1).EQ.'A'  .OR.  ISTART.EQ.0) READA=2
      RETURN
      END
C*MODULE MPCMOL  *DECK RDMOL
      SUBROUTINE RDMOL
C
      IMPLICIT DOUBLE PRECISION (A-H, O-Z)
C
      PARAMETER (MXATM=2000)
C
C MODULE TO READ IN GEOMETRY FILE, OUTPUT IT TO THE USER,
C AND CHECK THE DATA TO SEE IF IT IS REASONABLE.
C EXIT IF NECESSARY.
C
C  ON EXIT NATOMS    = NUMBER OF ATOMS PLUS DUMMY ATOMS (IF ANY).
C          KEYWRD    = KEYWORDS TO CONTROL CALCULATION
C          KOMENT    = COMMENT CARD
C          TITLE     = TITLE CARD
C          LABELS    = ARRAY OF ATOMIC LABELS INCLUDING DUMMY ATOMS.
C          GEO       = ARRAY OF INTERNAL COORDINATES.
C          LOPT      = FLAGS FOR OPTIMIZATION OF MOLECULE
C          NA        = ARRAY OF LABELS OF ATOMS, BOND LENGTHS.
C          NB        = ARRAY OF LABELS OF ATOMS, BOND ANGLES.
C          NC        = ARRAY OF LABELS OF ATOMS, DIHEDRAL ANGLES.
C          LATOM     = LABEL OF ATOM OF REACTION COORDINATE.
C          LPARAM    = RC: 1 FOR LENGTH, 2 FOR ANGLE, AND 3 FOR DIHEDRAL
C          REACT(200)= REACTION COORDINATE PARAMETERS
C          LOC(1,I)  = LABEL OF ATOM TO BE OPTIMIZED.
C          LOC(2,I)  = 1 FOR LENGTH, 2 FOR ANGLE, AND 3 FOR DIHEDRAL.
C          NVAR      = NUMBER OF PARAMETERS TO BE OPTIMIZED.
C          XPARAM    = STARTING VALUE OF PARAMETERS TO BE OPTIMIZED.
C
C *** INPUT THE TRIAL GEOMETRY  \IE.  KGEOM=0\
C   LABEL(I) = THE ATOMIC NUMBER OF ATOM\I\.
C            = 99, THEN THE I-TH ATOM IS A DUMMY ATOM USED ONLY TO
C              SIMPLIFY THE DEFINITION OF THE MOLECULAR GEOMETRY.
C   GEO(1,I) = THE INTERNUCLEAR SEPARATION \IN ANGSTROMS\ BETWEEN ATOMS
C              NA(I) AND (I).
C   GEO(2,I) = THE ANGLE NB(I):NA(I):(I) INPUT IN DEGREES; STORED IN
C              RADIANS.
C   GEO(3,I) = THE ANGLE BETWEEN THE VECTORS NC(I):NB(I) AND NA(I):(I)
C              INPUT IN DEGREES - STORED IN RADIANS.
C  LOPT(J,I) = -1 IF GEO(J,I) IS THE REACTION COORDINATE.
C            = +1 IF GEO(J,I) IS A PARAMETER TO BE OPTIMIZED
C            =  0 OTHERWISE.
C *** NOTE:    MUCH OF THIS DATA IS NOT INCLUDED FOR THE FIRST 3 ATOMS.
C     ATOM1  INPUT LABELS(1) ONLY.
C     ATOM2  INPUT LABELS(2) AND GEO(1,2) SEPARATION BETWEEN ATOMS 1+2
C     ATOM3  INPUT LABELS(3), GEO(1,3)    SEPARATION BETWEEN ATOMS 2+3
C              AND GEO(2,3)              ANGLE ATOM1 : ATOM2 : ATOM3
C
      PARAMETER (ONE=1.00D+00, UNITS=ONE/0.52917724924D+00)
C
      COMMON /INFOA / NAT,ICH,MUL,NUM,NQMT,NE,NA,NB,
     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
      COMMON /GEOKST/ NATOMS,LABELS(MXATM),INA(MXATM),INB(MXATM),
     *                INC(MXATM)
      COMMON /MPCGEO/ GEO(3,MXATM)
C--   COMMON /PATH  / LATOM,LPARAM,REACT(200)
C
C--   LATOM=0
C--   LPARAM=0
      NATOMS=NAT
      INA(1)=99
C
      DO 200 J = 1,NATOMS
        LABELS(J)=INT(ZAN(J))
        DO 100 I = 1,3
         GEO(I,J) = C(I,J)/UNITS
  100   CONTINUE
  200 CONTINUE
      RETURN
      END
C*MODULE MPCMOL  *DECK REFER
      SUBROUTINE REFER
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER (MXATM=2000)
C
      LOGICAL GOPARR,DSKWRK,MASWRK
C
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /WFNOPT/ SCFTYP,CITYP,DFTYPE,CCTYP,MPLEVL,MPCTYP
      COMMON
     1       /REFS/ ALLREF(107,4)
     2       /MOLKST/ NUMAT,NAT(MXATM),NFIRST(MXATM),NMIDLE(MXATM),
     3                NLAST(MXATM), NORBS, NELECS,NALPHA,NBETA,
     4                NCLOSE,NOPEN,NDUMY,FRACT
     5       /KEYWRD/ KEYWRD
      LOGICAL ALLOK, ELEMNS(107), MIXOK, MIX
      SAVE MIX
      CHARACTER KEYWRD*241, ALLREF*80
      DATA MIX/.FALSE./
C
      DATA KMINDO,KMNDO,KAM1,KPM3/4HMIND,4HMNDO,4HAM1 ,4HPM3 /
C
      MIXOK=(INDEX(KEYWRD,'PARASOK').NE.0)
      DO 10 I=1,102
   10 ELEMNS(I)=.FALSE.
C
      MODE=0
      IF(MPCTYP.EQ.KMNDO)  MODE=1
      IF(MPCTYP.EQ.KMINDO) MODE=2
      IF(MPCTYP.EQ.KAM1)   MODE=3
      IF(MPCTYP.EQ.KPM3)   MODE=4
      IF(MODE.EQ.0) THEN
         IF(MASWRK) WRITE(6,*) 'HAMILTONIAN CONFUSION'
         CALL ABRT
      END IF
C
      IF(MASWRK) WRITE(6,9000)
 9000 FORMAT(/1X,'THE PARAMETERS USED IN THIS CALCULATION ARE',
     *           ' DESCRIBED IN:'/)
      ALLREF(100,MODE)=' '
      DO 20 I=1,NUMAT
         J=NAT(I)
   20 ELEMNS(J)=.TRUE.
      ALLOK=.TRUE.
      DO 30 I=1,102
         IF(ELEMNS(I)) THEN
            IF(I.LT.99.AND..NOT.MIX.AND.MODE.EQ.3)
     1MIX=(INDEX(ALLREF(I,3),'MNDO').NE.0)
            IF(ALLREF(I,MODE)(1:1).NE.' ' .AND.MASWRK) THEN
               WRITE(6,'(A,I3)')' DATA ARE NOT AVAILABLE FOR ELEMENT NO.
     1',I
               ALLOK=.FALSE.
            ELSE IF(MASWRK) THEN
               WRITE(6,'(A)')ALLREF(I,MODE)
            END IF
         END IF
   30 CONTINUE
      IF(MIX.AND..NOT.MIXOK) THEN
         IF(MASWRK) THEN
         WRITE(6,40)
     1  'SOME ELEMENTS HAVE BEEN SPECIFIED FOR WHICH ONLY MNDO',
     2  'PARAMETERS ARE AVAILABLE.  SUCH MIXTURES OF METHODS ARE',
     3  'VERY RISKY AND HAVE NOT BEEN FULLY TESTED.  IF YOU FEEL',
     4  'THE RISK IS WORTH WHILE - CHECK THE MANUAL FIRST - THEN',
     5  'SPECIFY "PARASOK" IN THE KEYWORDS'
         END IF
         CALL ABRT
         STOP
      END IF
      IF(ALLOK)RETURN
      IF(MASWRK)
     *WRITE(6,40)
     1 'SOME ELEMENTS HAVE BEEN SPECIFIED FOR WHICH',
     2 'NO PARAMETERS ARE AVAILABLE.  CALCULATION STOPPED.'
      CALL ABRT
      STOP
   40 FORMAT(/////10X,A,4(/10X,A))
      END
C*MODULE MPCMOL  *DECK SETUPG
      SUBROUTINE SETUPG
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      COMMON /EXPONT/ ZS(107),ZP(107),ZD(107),ZF(107),ZG(107)
      COMMON /NATYPE/ NZTYPE(107),MTYPE(30),LTYPE
      COMMON /TEMP/  CC(60,6),ZZ(60,6)
      COMMON /STO6G/ ALLC(6,5,2),ALLZ(6,5,2)
C
C     SET-UP THE STEWART'S STO-6G EXPANSIONS
C                                            1S
      ALLZ(1,1,1) =2.310303149D+001
      ALLZ(2,1,1) =4.235915534D+00
      ALLZ(3,1,1) =1.185056519D+00
      ALLZ(4,1,1) =4.070988982D-01
      ALLZ(5,1,1) =1.580884151D-01
      ALLZ(6,1,1) =6.510953954D-02
C
      ALLC(1,1,1) =9.163596280D-03
      ALLC(2,1,1) =4.936149294D-02
      ALLC(3,1,1) =1.685383049D-01
      ALLC(4,1,1) =3.705627997D-01
      ALLC(5,1,1) =4.164915298D-01
      ALLC(6,1,1) =1.303340841D-01
C                                      2S
      ALLZ(1,2,1) =2.768496241D+01
      ALLZ(2,2,1) =5.077140627D+00
      ALLZ(3,2,1) =1.426786050D+00
      ALLZ(4,2,1) =2.040335729D-01
      ALLZ(5,2,1) =9.260298399D-02
      ALLZ(6,2,1) =4.416183978D-02
C
      ALLC(1,2,1) =-4.151277819D-03
      ALLC(2,2,1) =-2.067024148D-02
      ALLC(3,2,1) =-5.150303337D-02
      ALLC(4,2,1) =3.346271174D-01
      ALLC(5,2,1) =5.621061301D-01
      ALLC(6,2,1) =1.712994697D-01
C                                     2P
      ALLZ(1,2,2) =5.868285913D+00
      ALLZ(2,2,2) =1.530329631D+00
      ALLZ(3,2,2) =5.475665231D-01
      ALLZ(4,2,2) =2.288932733D-01
      ALLZ(5,2,2) =1.046655969D-01
      ALLZ(6,2,2) =4.948220127D-02
C
      ALLC(1,2,2) =7.924233646D-03
      ALLC(2,2,2) =5.144104825D-02
      ALLC(3,2,2) =1.898400060D-01
      ALLC(4,2,2) =4.049863191D-01
      ALLC(5,2,2) =4.012362861D-01
      ALLC(6,2,2) =1.051855189D-01
C                                      3S
      ALLZ(1,3,1) =3.273031938D+00
      ALLZ(2,3,1) =9.200611311D-01
      ALLZ(3,3,1) =3.593349765D-01
      ALLZ(4,3,1) =8.636686991D-02
      ALLZ(5,3,1) =4.797373812D-02
      ALLZ(6,3,1) =2.724741144D-02
      ALLC(1,3,1) =-6.775596947D-03
      ALLC(2,3,1) =-5.639325779D-02
      ALLC(3,3,1) =-1.587856086D-01
      ALLC(4,3,1) =5.534527651D-01
      ALLC(5,3,1) =5.015351020D-01
      ALLC(6,3,1) =7.223633674D-02
C                                     3P
      ALLZ(1,3,2) =5.077973607D+00
      ALLZ(2,3,2) =1.340786940D+00
      ALLZ(3,3,2) =2.248434849D-01
      ALLZ(4,3,2) =1.131741848D-01
      ALLZ(5,3,2) =6.076408893D-02
      ALLZ(6,3,2) =3.315424265D-02
      ALLC(1,3,2) =-3.329929840D-03
      ALLC(2,3,2) =-1.419488340D-02
      ALLC(3,3,2) =1.639395770D-01
      ALLC(4,3,2) =4.485358256D-01
      ALLC(5,3,2) =3.908813050D-01
      ALLC(6,3,2) =7.411456232D-02
C                                     4S
      ALLZ(1,4,1) = 1.365346D+00
      ALLZ(2,4,1) = 4.393213D-01
      ALLZ(3,4,1) = 1.877069D-01
      ALLZ(4,4,1) = 9.360270D-02
      ALLZ(5,4,1) = 5.052263D-02
      ALLZ(6,4,1) = 2.809354D-02
      ALLC(1,4,1) = 3.775056D-03
      ALLC(2,4,1) =-5.585965D-02
      ALLC(3,4,1) =-3.192946D-01
      ALLC(4,4,1) =-2.764780D-02
      ALLC(5,4,1) = 9.049199D-01
      ALLC(6,4,1) = 3.406258D-01
C                                   4P
      ALLC(1,4,2) =-7.052075D-03
      ALLC(2,4,2) =-5.259505D-02
      ALLC(3,4,2) =-3.773450D-02
      ALLC(4,4,2) = 3.874773D-01
      ALLC(5,4,2) = 5.791672D-01
      ALLC(6,4,2) = 1.221817D-01
      ALLZ(1,4,2) = 1.365346D+00
      ALLZ(2,4,2) = 4.393213D-01
      ALLZ(3,4,2) = 1.877069D-01
      ALLZ(4,4,2) = 9.360270D-02
      ALLZ(5,4,2) = 5.052263D-02
      ALLZ(6,4,2) = 2.809354D-02
C                                     5S
      ALLZ(1,5,1) = 7.701420258D-01
      ALLZ(2,5,1) = 2.756268915D-01
      ALLZ(3,5,1) = 1.301847480D-01
      ALLZ(4,5,1) = 6.953441940D-02
      ALLZ(5,5,1) = 4.002545502D-02
      ALLZ(6,5,1) = 2.348388309D-02
      ALLC(1,5,1) = 1.267447151D-02
      ALLC(2,5,1) = 3.266734789D-03
      ALLC(3,5,1) =-4.307553999D-01
      ALLC(4,5,1) =-3.231998963D-01
      ALLC(5,5,1) = 1.104322879D+00
      ALLC(6,5,1) = 4.368498703D-01
C                                      5P
      ALLZ(1,5,2) = 7.701420258D-01
      ALLZ(2,5,2) = 2.756268915D-01
      ALLZ(3,5,2) = 1.301847480D-01
      ALLZ(4,5,2) = 6.953441940D-02
      ALLZ(5,5,2) = 4.002545502D-02
      ALLZ(6,5,2) = 2.348388309D-02
      ALLC(1,5,2) =-1.105673292D-03
      ALLC(2,5,2) =-6.243132446D-02
      ALLC(3,5,2) =-1.628476766D-01
      ALLC(4,5,2) = 3.210328714D-01
      ALLC(5,5,2) = 6.964579592D-01
      ALLC(6,5,2) = 1.493146125D-01
      DO 30 I=1,10
         IF(MTYPE(I).EQ.0) GO TO 30
         NI=MTYPE(I)
         XI=ZS(NI)
         IA=I*4-3
         IB=IA+3
         IF(NI.LT.2) THEN
            NQN=1
         ELSE IF(NI.LT.10) THEN
            NQN=2
         ELSE IF(NI.LT.18) THEN
            NQN=3
         ELSE IF(NI.LT.36) THEN
            NQN=4
         ELSE IF(NI.LT.54) THEN
            NQN=5
         ELSE
            WRITE(6,*) 'SETUPG: NO PQN=6 STO-6G BASIS IS AVAILABLE'
            CALL ABRT
            STOP
         END IF
         DO 20 K=IA,IB
            L=1
            IF(K.GT.IA) L=2
            IF(K.GT.IA) XI=ZP(NI)
            DO 10 J=1,6
               CC(K,J)=ALLC(J,NQN,L)
   10       ZZ(K,J)=ALLZ(J,NQN,L)*XI**2
   20    CONTINUE
   30 CONTINUE
      RETURN
      END
