      SUBROUTINE MOLDAT(MODE)
      implicit double precision (A-H,O-Z)
      include 'SIZES'
      include 'MOP6FILE'
      common /GEOKST/ NATOMS,LABELS(NUMATM),
     1                NA(NUMATM),NB(NUMATM),NC(NUMATM)
      common /MOLMEC/ HTYPE(4),NHCO(4,20),NNHCO,ITYPE
      common /MOLKST/ NUMAT,NAT(NUMATM),NFIRST(NUMATM),NMIDLE(NUMATM),
     1                NLAST(NUMATM), NORBS, NELECS,NALPHA,NBETA,
     2                NCLOSE,NOPEN,NDUMY,FRACT
     3       /KEYWRD/ KEYWRD
     4       /NATORB/ NATORB(107)
      common /CORE  / CORE(107)
     1       /BETAS / BETAS(107),BETAP(107),BETAD(107)
     2       /MOLORB/ USPD(MAXORB),PSPD(MAXORB)
     3       /VSIPS / VS(107),VP(107),VD(107)
     4       /ONELEC/ USS(107),UPP(107),UDD(107)
      common /ATHEAT/ ATHEAT
     1       /POLVOL/ POLVOL(107)
     2       /MULTIP/ DD(107),QQ(107),AM(107),AD(107),AQ(107)
     3       /TWOELE/ GSS(107),GSP(107),GPP(107),GP2(107),HSP(107)
     4                ,GSD(107),GPD(107),GDD(107)
     5       /IDEAA / GUESA1(107,10),GUESA2(107,10),GUESA3(107,10)
     6       /IDEAS / GUESS1(107,10),GUESS2(107,10),GUESS3(107,10)
     7       /IDEAP / GUESP1(107,10),GUESP2(107,10),GUESP3(107,10)
      common /ALPHA / ALP(107)
     1       /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)
      common /GEOM  / GEO(3,NUMATM)
      PARAMETER (MDUMY=MAXPAR**2-MPACK)
      common /SCRACH/ RXYZ(MPACK), XDUMY(MDUMY)
*
*  common BLOCKS FOR MINDO/3
*
      common /ONELE3 /  USS3(18),UPP3(18)
     1       /ATOMI3 /  EISOL3(18),EHEAT3(18)
     2       /EXPON3 /  ZS3(18),ZP3(18)
*
*  end OF MINDO/3 common BLOCKS
*
      common /EXPONT/ ZS(107),ZP(107),ZD(107)
      common /ATOMIC/ EISOL(107),EHEAT(107)
      dimension COORD(3,NUMATM), ISWAP(2,20), ESTORE(107)
      character*241 KEYWRD, OLDE(20)*6, ALLREF*80
      logical DEBUG, UHF,EXCI, TRIP, MINDO3, BIRAD, AM1, OPEN, LPM3,
     1LMNDO, HALFE, SLOW
      data ESTORE(1)/0.D0/
      if (ESTORE(1).EQ.0.D0) then
      do 9 I=1,107
   9  ESTORE(I)=EHEAT(I)
      end if
      do 8 I=1,107
  8   EHEAT(I)=ESTORE(I)
      DEBUG = (INDEX(KEYWRD,'MOLDAT').NE.0)
      LPM3  = (INDEX(KEYWRD,'PM3').NE.0)
      MINDO3= (INDEX(KEYWRD,'MINDO').NE.0)
      UHF=(INDEX(KEYWRD,'UHF') .NE. 0)
      AM1= (INDEX(KEYWRD,'AM1').NE.0)
      LMNDO=(.NOT.AM1.AND..NOT.LPM3)
      KHARGE=0
      I=INDEX(KEYWRD,'CHARGE')
      if (I.NE.0) KHARGE = READA(KEYWRD,I)
      ELECS=-KHARGE
      NDORBS=0
      ATHEAT=0.D0
      EAT=0.D0
      NUMAT=0
      if (MODE.EQ.1) go to 80
      if ( LMNDO ) then
*
*    SWITCH IN MNDO PARAMETERS
*
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
         GUESS1(6,1)=0.D0
         do 10 I=1,107
            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
      else if ( LPM3 ) then
*
*    SWITCH IN MNDO-PM3 PARAMETERS
*
         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)
            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
      else
*
*    SWITCH IN AM1 PARAMETERS
*
         do 50 I=1,107
            do 40 J=1,10
               GUESS1(I,J)=GUESA1(I,J)
               GUESS2(I,J)=GUESA2(I,J)
   40       GUESS3(I,J)=GUESA3(I,J)
            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.D0) then
         write(MIOB,
     %     '(''  THE HAMILTONIAN REQUESTED IS NOT AVAILABLE IN'',
     %       '' THIS PROGRAM'')')
         call filclo ! J/M
         stop ' MOLDAT - nonexistent Hamiltonian used'
      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 270
         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 270
            USPD(J)=UPP(NI)
   90    continue
         NHEAVY=NHEAVY+1
  100    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
      if (NUMAT.EQ.1) then
         if (INDEX(KEYWRD,'FORCE').NE.0) then
            write(MIOB,
     %       '(///,A)')'      A SINGLE ATOM HAS NO VIBRATIONAL MODES'
c
            call filclo ! J/M
            stop ' MOLDAT - no vibrations for single atom!'
c
         end if
      end if
      if (MODE.NE.1) call REFER
      ATHEAT=ATHEAT-EAT*23.061D0
      NORBS=NLAST(NUMAT)
      if (NORBS.GT.MAXORB) then
         write(MIOB,'(//10X,''**** MAX. NUMBER OF ORBITALS:'',I4,/
     %                   10X,''NUMBER OF ORBITALS IN SYSTEM:'',I4)')
     %       MAXORB,NORBS
c
         call filclo ! J/M
         stop ' MOLDAT - too many orbitals'
c
      end if
      NLIGHT=NUMAT-NHEAVY
      N2EL=50*NHEAVY*(NHEAVY-1)+10*NHEAVY*NLIGHT+(NLIGHT*(NLIGHT-1))/2
      if (N2EL.GT.N2ELEC) then
         write(MIOB,
     %     '(//10X,''**** MAX. NUMBER OF TWO-ELECTRON INTEGRALS:'',I8/
     %         10X,''NUMBER OF TWO ELECTRON INTEGRALS IN SYSTEM:'',
     %                                                            I8)')
     %     N2ELEC,N2EL
c
         call filclo ! J/M
         stop ' MOLDAT - too many integrals'
c
      end if
C
C   NOW TO CALCULATE THE NUMBER OF LEVELS OCCUPIED
      TRIP=(INDEX(KEYWRD,'TRIP').NE.0)
      EXCI=(INDEX(KEYWRD,'EXCI').NE.0)
      BIRAD=(EXCI.OR.INDEX(KEYWRD,'BIRAD').NE.0)
      if (INDEX(KEYWRD,'C.I.') .NE. 0 .AND. UHF ) then
         write(MIOB,'(//10X,''C.I. NOT ALLOWED WITH UHF '')')
c
         call filclo ! J/M
         stop ' MOLDAT - C.I. and UHF not allowed in same run'
c
      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=MAX(ELECS+0.5D0,0.D0)
      NELECS=MIN(2*NORBS,NELECS)
      NCLOSE=0
      NOPEN=0
      if ( UHF ) then
         FRACT=1.D0
         NBETA=NELECS/2
         if ( TRIP ) then
            if (NBETA*2 .NE. NELECS) then
               write(MIOB,'(//10X,''TRIPLET SPECIFIED WITH ODD NUMBER'',
     1                            '' OF ELECTRONS, CORRECT FAULT '')')
c
               call filclo ! J/M
               stop ' MOLDAT - triplet with odd # of electrons?'
c
            else
               if (MODE.NE.1)
     1           write(MIOB,'(//'' TRIPLET STATE CALCULATION'')')
               NBETA=NBETA-1
            end if
         end if
         if (INDEX(KEYWRD,'QUAR').NE.0) then
            if (NBETA*2 .EQ. NELECS) then
              write(MIOB,'(//10X,''QUARTET SPECIFIED WITH EVEN NUMBER'',
     1                           '' OF ELECTRONS, CORRECT FAULT '')')
c
               call filclo ! J/M
               stop ' MOLDAT - quartet wit even # of electrons?'
c
            else
               if (MODE.NE.1)
     1           write(MIOB,'(//'' QUARTET STATE CALCULATION'')')
               NBETA=NBETA-1
            end if
         end if
         if (INDEX(KEYWRD,'QUIN').NE.0) then
            if (NBETA*2 .NE. NELECS) then
               write(MIOB,'(//10X,''QUINTET SPECIFIED WITH ODD NUMBER'',
     %                            '' OF ELECTRONS, CORRECT FAULT '')')
c
               call filclo ! J/M
               stop ' MOLDAT - quintet with odd # of electrons?'
c
            else
               if (MODE.NE.1)
     1           write(MIOB,'(//'' QUINTET STATE CALCULATION'')')
               NBETA=NBETA-2
            end if
         end if
         if (INDEX(KEYWRD,'SEXT').NE.0) then
            if (NBETA*2 .EQ. NELECS) then
               write(MIOB,'(//10X,''SEXTET SPECIFIED WITH EVEN NUMBER'',
     %                            '' OF ELECTRONS, CORRECT FAULT '')')
c
               call filclo ! J/M
               stop ' MOLDAT - sextet with even # of electrons?'
c
            else
               if (MODE.NE.1)
     1           write(MIOB,'(//'' SEXTET STATE CALCULATION'')')
               NBETA=NBETA-2
            end if
         end if
         NALPHA=NELECS-NBETA
         if (MODE.NE.1)
     1     write(MIOB,'(//10X,
     2            ''UHF CALCULATION, NO. OF ALPHA ELECTRONS ='',I3/
     3        27X,''NO. OF BETA  ELECTRONS ='',I3)')
     4                                          NALPHA,NBETA
      else
C
C   NOW TO DETERMINE OPEN AND CLOSED SHELLS
C
         OPEN=.FALSE.
         IELEC=0
         ILEVEL=0
         if ( TRIP .OR. EXCI .OR. BIRAD ) then
            if ( (NELECS/2)*2 .NE. NELECS) then
               write(MIOB,'(//10X,''SYSTEM SPECIFIED WITH ODD NUMBER'',
     1                            '' OF ELECTRONS, CORRECT FAULT '')')
c
               call filclo ! J/M
               stop ' MOLDAT - odd # of electrons and trip,exci,birad?'
c
            end if
            if (MODE.NE.1) then
               if (BIRAD) write(MIOB,'(//'' SYSTEM IS A BIRADICAL'')')
               if (TRIP )
     1               write(MIOB,'(//'' TRIPLET STATE CALCULATION'')')
               if (EXCI )
     1               write(MIOB,'(//'' EXCITED STATE CALCULATION'')')
            end if
            IELEC=2
            ILEVEL=2
         else if ((NELECS/2)*2.NE.NELECS) then
            IELEC=1
            ILEVEL=1
         end if
         if (INDEX(KEYWRD,'QUAR').NE.0) then
            if (MODE.NE.1)
     1         write(MIOB,'(//'' QUARTET STATE CALCULATION'')')
            IELEC=3
            ILEVEL=3
         end if
         if (INDEX(KEYWRD,'QUIN').NE.0) then
            if (MODE.NE.1)
     1         write(MIOB,'(//'' QUINTET STATE CALCULATION'')')
            IELEC=4
            ILEVEL=4
         end if
         if (INDEX(KEYWRD,'SEXT').NE.0) then
            if (MODE.NE.1)
     1         write(MIOB,'(//'' SEXTET STATE CALCULATION'')')
            IELEC=5
            ILEVEL=5
         end if
         I=INDEX(KEYWRD,'OPEN(')
         if (I.NE.0) then
            IELEC = READA(KEYWRD,I)
            ILEVEL = 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
               write(MIOB,
     %           '('' IMPOSSIBLE NUMBER OF OPEN SHELL ELECTRONS'')')
c
               call filclo ! J/M
               stop ' MOLDAT - strange # of open-shell electrons'
c
            end if
            NCLOSE=NCLOSE-IELEC/2
            NOPEN=ILEVEL
            if (NCLOSE+NOPEN.GT.NORBS) then
               write(MIOB,'(A)')
     %           ' NUMBER OF DOUBLY FILLED PLUS PARTLY FILLED LEVELS'
               write(MIOB,'(A)')
     %           ' GREATER THAN TOTAL NUMBER OF ORBITALS.'
c
               call filclo ! J/M
               stop ' MOLDAT - check sum of orbitals'
c
            end if
            FRACT=IELEC*1.D0/ILEVEL
            if (MODE.NE.1)
     1        write(MIOB,
     2         '('' THERE ARE'',I3,'' DOUBLY FILLED LEVELS'')') NCLOSE
         end if
         if (MODE.NE.1)
     1        write(MIOB,'(//10X,''RHF CALCULATION, NO. OF '',
     2                   ''DOUBLY OCCUPIED LEVELS ='',I3)')     NCLOSE
         if (MODE.NE.1 .AND. NOPEN.NE.0.AND.ABS(FRACT-1.D0).LT.1.D-4)
     1        write(MIOB,
     2         '(/27X,''NO. OF SINGLY OCCUPIED LEVELS ='',I3)')  NOPEN
         if (MODE.NE.1 .AND. NOPEN.NE.0.AND.ABS(FRACT-1.D0).GT.1.D-4)
     1      write(MIOB,
     2       '(/27X,''NO. OF LEVELS WITH OCCUPANCY'',F6.3,''  ='',I3)')
     3                                                     FRACT,NOPEN
         if (INDEX(KEYWRD,'C.I.=(').NE.0) then
            I = READA(KEYWRD,INDEX(KEYWRD,'C.I.=(')+5)-
     1      READA(KEYWRD,INDEX(KEYWRD,'C.I.=(')+7)
            if (NOPEN.GT.I) then
               write(MIOB,
     1          '(//,'' NUMBER OF OPEN-SHELLS ALLOWED IN C.I. IS LESS''
     2              /''    THAN THAT SPECIFIED BY OTHER KEYWORDS'')')
c
               call filclo ! J/M
               stop ' MOLDAT - inconsistent # of open shells'
c
            end if
         end if
         if (INDEX(KEYWRD,'C.I.').NE.0.AND.NOPEN.EQ.0) then
            NOPEN=1
            NCLOSE=NCLOSE-1
            FRACT=2.D0
         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=1.0001D0*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.UHF) then
C
C   MSDEL = NUMBER OF ALPHA ELECTRONS - NUMBER OF BETA ELECTRONS
C
         NDOUBL=99
         if (INDEX(KEYWRD,'C.I.=(').NE.0) then
            NDOUBL = READA(KEYWRD,INDEX(KEYWRD,'C.I.=(')+7)
            NMOS = READA(KEYWRD,INDEX(KEYWRD,'C.I.=(')+5)
         else if (INDEX(KEYWRD,'C.I.=').NE.0) then
            NMOS = 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
         NE=MAX(0.D0,(NCLOSE-J+1.D0))*2.D0+
     1     MAX(0.D0,(NOPEN-NCLOSE)*FRACT) + 0.5D0
         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
            write(MIOB,'(A)')
     1' SPECIFIED SPIN COMPONENT NOT SPANNED BY ACTIVE SPACE'
c
            call filclo ! J/M
            stop ' MOLDAT - check requsted spin component'
c
         end if
      end if
C#      write(MIOB,'(''  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.D0.AND.FRACT.NE.0.D0
     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.
     1   INDEX(KEYWRD,'ANALYT').NE.0.AND.HALFE) then
         write(MIOB,*)
         write(MIOB,'(A/)')' KEYWORD ''ANALYT'' CANNOT BE USED HERE: ',
     %     ' ANALYTICAL C.I. DERIVATIVES MUST USE FINITE DIFFERENCES',
     %     ' TO CORRECT, REMOVE KEYWORD ''ANALYT'' OR ADD ''NOANCI'''
c
         call filclo ! J/M
         stop ' MOLDAT - dont use ANALYT or add NOANCI'
c
      end if
      YY=FLOAT(KHARGE)/(NORBS+1.D-10)
      L=0
      do 160 I=1,NUMAT
         NI=NAT(I)
         XX=1.D0/(NLAST(I)-NFIRST(I)+1+1.D-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.D0
  160 continue
C
C   WRITE OUT THE INTERATOMIC DISTANCES
C
      call GMETRY(GEO,COORD)
      RMIN=100.D0
      L=0
      do 170 I=1,NUMAT
         do 170 J=1,I
            L=L+1
            RXYZ(L)=SQRT((COORD(1,I)-COORD(1,J))**2+
     1                     (COORD(2,I)-COORD(2,J))**2+
     2                     (COORD(3,I)-COORD(3,J))**2)
            if (RMIN.GT.RXYZ(L) .AND. I .NE. J .AND.
     1         (NAT(I).LT.103 .OR. NAT(J).LT.103)) then
               IMINR=I
               JMINR=J
               RMIN=RXYZ(L)
            end if
  170 continue
      NNHCO=0
C
C   SET UP MOLECULAR-MECHANICS CORRECTION TO -(C=O)-(NH)- LINKAGE
C   THIS WILL BE USED IF MMOK HAS BEEN SPECIFIED.
C
      ITYPE=1
      if (INDEX(KEYWRD,'AM1').NE.0)   ITYPE=2
      if (INDEX(KEYWRD,'PM3').NE.0)   ITYPE=3
      if (INDEX(KEYWRD,'MINDO').NE.0) ITYPE=4
C
C   IDENTIFY O=C-N-H SYSTEMS VIA THE INTERATOMIC DISTANCES MATRIX
      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
            IJ=MAX(I,J)
            JI=I+J-IJ
            if (RXYZ((IJ*(IJ-1))/2+JI).GT.1.3) go to 210
            do 200 K=1,NUMAT
               if (NAT(K).NE.7) go to 200
               JK=MAX(J,K)
               KJ=J+K-JK
               if (RXYZ((JK*(JK-1))/2+KJ).GT.1.6) go to 200
               do 190 L=1,NUMAT
                  if (NAT(L).NE.1) go to 190
                  KL=MAX(K,L)
                  LK=K+L-KL
                  if (RXYZ((KL*(KL-1))/2+LK).GT.1.3) 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
                     MK=MAX(M,K)
                     KM=M+K-MK
                     if (RXYZ((MK*(MK-1))/2+KM).GT.1.7) go to 180
                     NNHCO=NNHCO+1
                     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
      if (MODE.NE.1.AND.NNHCO.NE.0) then
         if (INDEX(KEYWRD,'MMOK').NE.0) then
           write(MIOB,'(A)')
     %      ' MOLECULAR MECHANICS CORRECTION APPLIED TO PEPTIDE LINKAGE'
         else if (INDEX(KEYWRD,'NOMM').NE.0) then
           write(MIOB,'(A,I2,2A)')
     %      ' THERE ARE ',NNHCO/2,' PEPTIDE LINKAGES',
     %      ' IDENTIFIED IN THIS SYSTEM'
           write(MIOB,'(A)')
     %      ' IF YOU WANT MM CORRECTION TO THE CONH BARRIER,'//
     %      ' ADD THE KEY-WORD "MMOK"'
            NNHCO=0
         else
            write(MIOB,'(A)')' THIS SYSTEM CONTAINS -HNCO- GROUPS.'
            write(MIOB,'(A)')
     %       ' YOU MUST SPECIFY "NOMM" OR "MMOK" REGARDING'//
     %       ' THE MOLECULAR MECHANICS CORRECTION'
c
            call filclo ! J/M
            stop ' MOLDAT - HNCO groups use NOMM or MMOK is mandatory'
c
         end if
      end if
      if (MODE.NE.1.AND.INDEX(KEYWRD,'NOINTER') .EQ. 0) then
         write(MIOB,'(//10X,''  INTERATOMIC DISTANCES'')')
         call VECPRT(RXYZ,NUMAT)
      end if
      if (RMIN.LT.0.8D0.AND.INDEX(KEYWRD,'GEO-OK') .EQ.0) then
         write(MIOB,230) IMINR,JMINR,RMIN
  230    format(//,'   ATOMS',I3,' AND',I3,' ARE SEPARATED BY',F8.4,
     1   ' ANGSTROMS.',/'   TO CONTINUE CALCULATION SPECIFY "GEO-OK"')
c
         call filclo ! J/M
         stop ' MOLDAT - atoms too close, use GEO-OK'
c
      end if
      if (.NOT. DEBUG) return
      write(MIOB,240) NUMAT,NORBS,NDORBS,NATOMS
  240 format('   NUMBER OF REAL ATOMS:',I4,/
     1      ,'   NUMBER OF ORBITALS:  ',I4,/
     2      ,'   NUMBER OF D ORBITALS:',I4,/
     3      ,'   TOTAL NO. OF ATOMS:  ',I4)
      write(MIOB,250) (USPD(I),I=1,NORBS)
  250 format('   ONE-ELECTRON DIAGONAL TERMS',/,10(/,10F8.3))
      write(MIOB,260) (PSPD(I),I=1,NORBS)
  260 format('   INITIAL P FOR ALL ATOMIC ORBITALS',/,10(/,10F8.3))
      return
  270 write(MIOB,
     %  '(//10X,'' MAXIMUM NUMBER OF ATOMIC ORBITALS EXCEEDED'')')
      write(MIOB,'(  10X,'' MAXIMUM ALLOWED ='',I4)') MAXORB
c
      call filclo ! J/M
      stop ' MOLDAT - too many orbitals'
c
      end
