C
C $Id: rdkeys.F,v 1.19 1999/03/02 20:02:25 arjan Exp arjan $
C
C------------------------------------------------------------------------
      SUBROUTINE RDKEYS(KEYWORDS,IERROR)
C
C     READS AND ECHOES KEYWORDS AND TITLE, AND EXTRACTS ANY NUMERICAL DATA
C     SUPPLIED VIA KEYWORDS.
C

      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
#include "divcon.dim"
#include "divcon.h"

      character*400 keywords
      logical push
CC---------------------------------------------------CC

      push = .false.
!     charge = 0.0
      xlog2 = log(2.0)

!      mais = .false.
!      pif = .false.
!      pm3 = .false.
!      pm3pif = .false.
c#ifdef MPI 
c      CALL OPNFIL(1,IERROR)
c      inpt=37
c#else
c      if (myid.eq.0) then
c         CALL OPNFIL(1,IERROR)
c      else
c         call opnpfil(11,ierror)
c      endif
c#endif
c      IF(IERROR.NE.0) RETURN

c      rewind(inpt)
c      lkeywrd = len(keywrd)
      lkeywrd = len(keywords)
c      DO 10 I=1,lkeywrd
c         KEYWRD(I:I) = ' '
c 10   CONTINUE

      CALL HEADER(ierror)
      if (ierror.ne.0) return

C
C     READ AND PROCESS KEYWORDS AND TITLE.
C
C     default is one line, but if an ampersand (&) is encountered,
C     a new line will be read.
C     if a line starts with a # on column 1, the line is treated as
C     comment

c      lstrip = 0
c      lkeywrd = lkeywrd-80
c 20   iold = lstrip
c      if (iold.gt.lkeywrd) then
         write(iout,'(/" ERROR: TOO MANY KEYWORDS SPECIFIED")')
c         ierror = 1
c         return
c      endif
c      read(inpt,'(A80)',END=500) keywrd(iold+1:iold+80)
c      if (keywrd(iold+1:iold+1).eq.'#') goto 20
      write(iout,'(A80)') keywrd(iold+1:iold+80)
c      call stripkeywd(keywrd,iold+80,lstrip)
      call stripkeywd(keywords,lkeywrd,lstrip)
      keywrd = keywords
c      i=index(keywrd(1:lstrip),'&')
c      i=index(keywords(1:lstrip),'&')
c      if (i.ne.0) then
c         keywords(i:i) = ' '
c         goto 20
c      endif
C
C     READ TITLE CARD.
C
c      READ(INPT,'(A80)',END=600) TITLE
C
C     ECHO TITLE BACK TO USER.
C
      WRITE(IOUT,'(/A80//)') TITLE
C
C     CONVERT KEYWORDS TO UPPER CASE THEN JOIN ANY WORDS SEPARATED BY
C     EQUALS SIGNS.
C
      CALL UPCASE1(KEYWRD,400)
      CALL WDJOIN(KEYWRD,400,'=',IERROR)
      IF(IERROR.NE.0)THEN
         WRITE(IOUT,'(/" IMPROPER USE OF EQUALS SIGN IN KEYWORDS")')
         RETURN
      ENDIF
C
C     LOCATE EACH WORD IN THE STRING KEYWRD AND TRY TO MATCH IT UP
C     WITH A RECOGNIZED KEYWORD.
C
C     note that if no precaution is taken to prevent printing keywords
C     longer than 15 characters (including the '=' and the value), bus
C     errors and core dumps occur. To prevent this, istopp is used.
C
      ISTART = 1
 100  CALL RDWORD(KEYWRD,ISTART,ISTOP)
      
      IF(ISTART.EQ.0) GO TO 400

      istopp = min(istart+14,istop)

      IF(KEYWRD(ISTART:ISTOP).EQ.'ADDMM')THEN
         WRITE(IOUT,'(" ADDMM",10X,
     .        "- ADD MM CORRECTION TO PEPTIDE TORSIONAL BARRIER")')

      elseif(keywrd(istart:istop).eq.'1SCF') then
         write(iout,'(" 1SCF",11X,"- PERFORM ONLY ONE SCF ITERATION")')
         onescf = .true.

#ifdef SCRF_IS_ON
      ELSEIF(INDEX(KEYWRD(ISTART:ISTOP),'SCRF-INSIDE').NE.0) THEN
         SCRF=.TRUE.
         scrfin = .true.
         WRITE(IOUT,'(" SCRF-INSIDE    - REACTION FIELD CALCULATION ",
     $        "FOR CAVITY ONLY")')

      elseif(index(keywrd(istart:istop),'DRAD=').ne.0) then
         call rdnum(keywrd,istart,DRAD,ierror)
         if (ierror.ne.0) goto 700
         write(iout,'(1x,15a1,"- PROBE RADIUS FOR INTERIOR ",
     $        "SURFACE CHARGES IS ",F12.5)')
     $        (KEYWRD(I:I),I=ISTART,ISTOPP),
     $        (' ',I=ISTOPP+1,ISTART+14),DRAD
         drad2 = drad*drad

      ELSEIF(INDEX(KEYWRD(ISTART:ISTOP),'SCRF').NE.0) THEN
         SCRF=.TRUE.
         WRITE(IOUT,'(" SCRF           - REACTION FIELD CALCULATION ")')
#endif

      elseif(index(keywrd(istart:istop),'GEOM=PDB').ne.0) then
         PDB=.True.
         write(iout,'(" GEOM=PDB       - COORDINATES IN PDB FORMAT ")')

      elseif(keywrd(istart:istop).eq.'NO-OVERLAP') then
         write(iout,'(" NO-OVERLAP",5X,"- NO OVERLAPPING BUFFERS",
     &        " FOR CLUSTSUB")')
         noovrlp = .true.

!      ELSEIF(keywrd(istart:istop).eq.'DIPOLE')THEN
!!        write(iout,'(" DIPOLE",9X,
!!    .        "- CALCULATE MOLECULAR DIPOLE")')
!         pole=.true.

      ELSEIF(keywrd(istart:istop).eq.'IP')THEN
         write(iout,'(" IP",13X,
     .        "- CALCULATE IONIZATION POTENTIAL")')
         ionpot=.true.

      elseif(keywrd(istart:istop).eq.'PUSH=') then
         CALL RDNUM(KEYWRD,ISTART,push0,IERROR)
         IF(IERROR.NE.0) GO TO 700
         WRITE(IOUT,'(1X,15A1,"- PUSH CLUSTER GROUPS APART BY ",
     .        E12.7, " A")')
     .        (KEYWRD(I:I),I=ISTART,ISTOPP),
     .        (' ',I=ISTOPP+1,ISTART+14),push0
         push = .true.

      elseif(keywrd(istart:istop).eq.'PUSH') THEN
         write(iout,'(" PUSH",11X,"- PUSH CLUSTER GROUPS APART")')
         push = .true.

      ELSEIF(keywrd(istart:istop).eq.'ZMAKE')THEN
         write(iout,'(" ZMAKE",10X,
     .        "- OUTPUT A Z-MATRIX")')
         zgen=.true.

      ELSEIF(index(keywrd(istart:istop),'ROTATE=').ne.0) THEN
         call rdnum(keywrd,istart,rotang,ierror)
         IF(ierror.ne.0) goto 700
         write(iout,'(1X,15A1,"- ROTATION ANGLE FOR BARRIER IS ",
     &        F8.3, " DEG")')
     &        (KEYWRD(I:I),I=ISTART,ISTOPP),
     &        (' ',I=ISTOPP+1,ISTART+14),rotang
         rotang=rotang*radian
         rotat=.true.

!      ELSEIF(KEYWRD(ISTART:ISTOP).EQ.'AM1')THEN
!!        WRITE(IOUT,'(" AM1",12X,
!!    .        "- AM1 HAMILTONIAN TO BE USED")')
!         am1 = .true.
!
!      ELSEIF(KEYWRD(ISTART:ISTOP).EQ.'RM1')THEN
!!        WRITE(IOUT,'(" RM1",12X,
!!    .        "- RM1 HAMILTONIAN TO BE USED")')
!         rm1 = .true.
!
!      ELSEIF(KEYWRD(ISTART:ISTOP).EQ.'PM3PDDG')THEN
!!        WRITE(IOUT,'(" PM3PDDG",12X,
!!    .        "- PM3PDDG HAMILTONIAN TO BE USED")')
!         pm3pddg = .true.
!
!      ELSEIF(KEYWRD(ISTART:ISTOP).EQ.'AM1D')THEN
!!        WRITE(IOUT,'(" AM1D",12X,
!!    .        "- AM1/d-PhoT HAMILTONIAN TO BE USED")')
!         am1d = .true.

!      ELSEIF(KEYWRD(ISTART:ISTOP).EQ.'CLUSTER')THEN
!!        WRITE(IOUT,'(" CLUSTER",8X,
!!    .        "- USE CLUSTER BASED SUBSETTING")')
!         clust = .true.

!      ELSEIF(KEYWRD(ISTART:ISTOP).EQ.'RESCLUST')THEN
!!       WRITE(IOUT,'(" RESCLUST",8X,
!!    .        "- USE RESIDUE-CLUSTER BASED SUBSETTING")')
!         resclust = .true.

      ELSEIF(KEYWRD(ISTART:ISTOP).EQ.'ATMSUB')THEN
          atmsub = .true.

      ELSEIF(KEYWRD(ISTART:ISTOP).EQ.'RESSUB')THEN
          ressub = .true.

c      ELSEIF(KEYWRD(ISTART:ISTOP).EQ.'CARTESIAN')THEN
         WRITE(IOUT,'(" CARTESIAN",6X,
     .        "- CARTESIAN COORDINATE FORMAT SELECTED")')
c         cart = .true.

c      ELSEIF(KEYWRD(ISTART:ISTOP).EQ.'CENTRAL')THEN
         WRITE(IOUT,'(" CENTRAL",8X,
     .        "- USE CENTRAL DIFFERENCE IN GRADIENT CALCULATION")')
c         centrl = .true.

!      ELSEIF(INDEX(KEYWRD(ISTART:ISTOP),'CHARGE=').NE.0)THEN
!         CALL RDNUM(KEYWRD,ISTART,CHARGE,IERROR)
!         IF(IERROR.NE.0) GO TO 700
!!        WRITE(IOUT,'(1X,15A1,"- A NET CHARGE OF ",F8.2,
!!    .        " TO BE PLACED ON SYSTEM")')
!!    .        (KEYWRD(I:I),I=ISTART,ISTOPP),
!!    .        (' ',I=ISTOPP+1,ISTART+14),CHARGE

      ELSEIF(KEYWRD(ISTART:ISTOP).EQ.'CHKRES')THEN
        WRITE(IOUT,'(" CHKRES",9X,
     .        "- CHECK INTER-ATOMIC DISTANCES FOR EACH RESIDUE")')

      ELSEIF(KEYWRD(ISTART:ISTOP).EQ.'PRTCOORDS')THEN
         WRITE(IOUT,'(" PRTCOORDS",6X,
     .        "- PRINT ATOMIC COORDINATES")')

      ELSEIF(INDEX(KEYWRD(ISTART:ISTOP),'CUTBOND=').NE.0)THEN
         CALL RDNUM(KEYWRD,ISTART,BCUT,IERROR)
         IF(IERROR.NE.0) GO TO 700
         WRITE(IOUT,'(1X,15A1,"- CUTOFF BONDING BEYOND ",F7.2,
     .        " ANGSTROMS")')
     .        (KEYWRD(I:I),I=ISTART,ISTOPP),
     .        (' ',I=ISTOPP+1,ISTART+14),BCUT
         cutbnd = .true.
         bcut2 = bcut*bcut

      ELSEIF(INDEX(KEYWRD(ISTART:ISTOP),'DESCF=').NE.0)THEN
         CALL RDNUM(KEYWRD,ISTART,DESCF,IERROR)
         IF(IERROR.NE.0) GO TO 700
         WRITE(IOUT,'(1X,15A1,"- SCF ENERGY CRITERION = ",G8.2," eV")')
     .        (KEYWRD(I:I),I=ISTART,ISTOPP),
     .        (' ',I=ISTOPP+1,ISTART+14),DESCF

      ELSEIF(INDEX(KEYWRD(ISTART:ISTOP),'DPSCF=').NE.0)THEN
         CALL RDNUM(KEYWRD,ISTART,DPSCF,IERROR)
         IF(IERROR.NE.0) GO TO 700
         WRITE(IOUT,'(1X,15A1,"- SCF DENSITY MATRIX CRITERION = ",
     .        G8.2)')
     .        (KEYWRD(I:I),I=ISTART,ISTOPP),
     .        (' ',I=ISTOPP+1,ISTART+14),DPSCF

c      ELSEIF(INDEX(KEYWRD(ISTART:ISTOP),'DIRECT').NE.0)THEN
         WRITE(IOUT,'(" DIRECT",9X,
     .        "- DO NOT STORE 2-ELECTRON INTEGRALS IN SCF",
     .        " CALCULATIONS")')
c         direct = .true.

      ELSEIF(INDEX(KEYWRD(ISTART:ISTOP),'PDUMP=').NE.0)THEN
         CALL RDNUM(KEYWRD,ISTART,DUMP,IERROR)
         IF(IERROR.NE.0) GO TO 700
         NnDUMP = DUMP
         WRITE(IOUT,'(1X,15A1,"- WRITE RESTART DENSITY MATRIX ",
     .        "FILE EVERY ",I4,
     .        " STEPS")')
     .        (KEYWRD(I:I),I=ISTART,ISTOPP),
     .        (' ',I=ISTOPP+1,ISTART+14),NnDUMP

      ELSEIF(INDEX(KEYWRD(ISTART:ISTOP),'DUMP=').NE.0)THEN
         CALL RDNUM(KEYWRD,ISTART,DUMP,IERROR)
         IF(IERROR.NE.0) GO TO 700
         NDUMP = DUMP
         WRITE(IOUT,'(1X,15A1,"- WRITE RESTART FILE EVERY ",I4,
     .        " CYCLES")')
     .        (KEYWRD(I:I),I=ISTART,ISTOPP),
     .        (' ',I=ISTOPP+1,ISTART+14),NDUMP

      ELSEIF(INDEX(KEYWRD(ISTART:ISTOP),'ETEST=').NE.0)THEN
         CALL RDNUM(KEYWRD,ISTART,ETEST,IERROR)
         IF(IERROR.NE.0) GO TO 700
         WRITE(IOUT,'(1X,15A1,"- GEOMETRY OPTIMIZATION ENERGY",
     .        " CHANGE CRITERION =",F8.5)')
     .        (KEYWRD(I:I),I=ISTART,ISTOPP),
     .        (' ',I=ISTOPP+1,ISTART+14),ETEST

      ELSEIF(KEYWRD(ISTART:ISTOP).EQ.'FORCE-IT')THEN
         WRITE(IOUT,'(" FORCE-IT",7X,
     .        "- FORCE THE GEOMETRY OPTIMIZATION CRITERIA TO BE",
     .        " SATISFIED")')

!      ELSEIF(KEYWRD(ISTART:ISTOP).EQ.'FULLSCF')THEN
!!        WRITE(IOUT,'(" FULLSCF",8X,
!!    .        "- DO FULL DIAGONALIZATIONS IN SCF CALCULATION")')
!         fullscf = .true.

c      ELSEIF(INDEX(KEYWRD(ISTART:ISTOP),'GRADIENT').NE.0)THEN
!        WRITE(IOUT,'(" GRADIENT",7X,
!    .        "- OUTPUT FINAL GRADIENT")')
c 	  gradient = .true.

      ELSEIF(INDEX(KEYWRD(ISTART:ISTOP),'GTEST=').NE.0)THEN
         CALL RDNUM(KEYWRD,ISTART,GTEST,IERROR)
         IF(IERROR.NE.0) GO TO 700
         WRITE(IOUT,'(1X,15A1,"- GEOMETRY OPTIMIZATION GRADIENT",
     .        " COMPONENT CRITERION =",F6.3)')
     .        (KEYWRD(I:I),I=ISTART,ISTOPP),
     .        (' ',I=ISTOPP+1,ISTART+14),GTEST

c      ELSEIF(KEYWRD(ISTART:ISTOP).EQ.'INTERNAL')THEN
         WRITE(IOUT,'(" INTERNAL",7X,
     .        "- INTERNAL COORDINATE FORMAT SELECTED")')

      elseif(index(keywrd(istart:istop),'MAXIT=').ne.0) then
         call rdinum(keywrd,istart,maxitscf,ierror)
         if(ierror.ne.0) goto 700
         write(iout,'(1X,15A1,"- MAXIMUM NUMBER OF SCF ITERATIONS IS ",
     $        I5)')
     &        (KEYWRD(I:I),I=ISTART,ISTOPP),
     &        (' ',I=ISTOPP+1,ISTART+14),maxitscf

      ELSEIF(INDEX(KEYWRD(ISTART:ISTOP),'MAXOPT=').NE.0)THEN
         CALL RDNUM(KEYWRD,ISTART,OPTMAX,IERROR)
         IF(IERROR.NE.0) GO TO 700
         MAXOPT = OPTMAX
         WRITE(IOUT,'(1X,15A1,"- DO A MAXIMUM OF ",I6," CYCLES",
     .        " OF GEOMETRY OPTIMIZATION")')
     .        (KEYWRD(I:I),I=ISTART,ISTOPP),
     .           (' ',I=ISTOPP+1,ISTART+14),MAXOPT

!      ELSEIF(KEYWRD(ISTART:ISTOP).EQ.'MNDO')THEN
!         WRITE(IOUT,'(" MNDO",11X,
!     .        "- MNDO HAMILTONIAN TO BE USED")')
!         mndo = .true.

      ELSEIF(KEYWRD(ISTART:ISTOP).EQ.'NODIIS')THEN
         WRITE(IOUT,'(" NODIIS",9X,
     .        "- DO NOT USE DIIS IN GEOMETRY OPTIMIZATION")')

      ELSEIF(KEYWRD(ISTART:ISTOP).EQ.'NOMM')THEN
         WRITE(IOUT,'(" NOMM",11X,
     .        "- DO NOT USE MOLECULAR MECHANICS CORRECTION FOR",
     .        " PEPTIDE TORSIONS")')

      ELSEIF(KEYWRD(ISTART:ISTOP).EQ.'PBC')THEN
         WRITE(IOUT,'(" PBC",12X,
     .        "- PERIODIC BOUNDARY CONDITIONS IN EFFECT")')
         pbc = .true.

!      ELSEIF(KEYWRD(ISTART:ISTOP).EQ.'PIF2')THEN
!         WRITE(IOUT,'(" PIF",12X,
!     .        "- PM3-PIF HAMILTONIAN TO BE USED")')
!         pif = .true.
!         pm3 = .true.
!         pm3pif = .true.

!      ELSEIF(KEYWRD(ISTART:ISTOP).EQ.'PM3')THEN
!!        WRITE(IOUT,'(" PM3",12X,
!!    .        "- PM3 HAMILTONIAN TO BE USED")')
!         pm3 = .true.

!      ELSEIF(KEYWRD(ISTART:ISTOP).EQ.'MAIS')THEN
!         WRITE(IOUT,'(" MAIS",11X,
!     .        "- PM3-MAIS HAMILTONIAN TO BE USED")')
!         mais = .true.
!         pm3 = .true.

!      ELSEIF(KEYWRD(ISTART:ISTOP).EQ.'PIF')THEN
!         WRITE(IOUT,'(" PIF",12X,
!     .        "- PM3-PIF HAMILTONIAN TO BE USED")')
!         pif = .true.
!         pm3 = .true.

!      ELSEIF(KEYWRD(ISTART:ISTOP).EQ.'PM3PIF')THEN
!         WRITE(IOUT,'(" PIF",12X,
!     .        "- PM3-PIF HAMILTONIAN TO BE USED")')
!         pif = .true.
!         pm3 = .true.
!         pm3pif = .true.

!      ELSEIF(KEYWRD(ISTART:ISTOP).EQ.'SOLUTE')THEN
!!        WRITE(IOUT,'(" PIF",12X,
!!    .        "- DIFFERENT INTERACTONS SOLUTE/WATER")')
!         solute = .true.

!      ELSEIF(INDEX(KEYWRD(ISTART:ISTOP),'NSOL=').NE.0)THEN
!         CALL RDNUM(KEYWRD,ISTART,ANSOL,IERROR)
!         NSOL = NINT(ANSOL)
!         IF(IERROR.NE.0) GO TO 700
      ELSEIF(KEYWRD(ISTART:ISTOP).EQ.'DEBUG')THEN
         WRITE(IOUT,'(" DEBUG",12X,
     .        "- DEBUGGING INFO")')
         debug = .true.

      ELSEIF(KEYWRD(ISTART:ISTOP).EQ.'OPT=BFGS')THEN
         WRITE(IOUT,'(" OPT=BFGS",7X,
     .        "- CARRY OUT A BFGS GEOMETRY OPTIMIZATION")')

      ELSEIF(KEYWRD(ISTART:ISTOP).EQ.'OPT=CONJGRAD')THEN
         WRITE(IOUT,'(" OPT=CONJGRAD",3X,
     .        "- CARRY OUT A CONJUGATE-GRADIENT GEOMETRY",
     .        " OPTIMIZATION")')

      ELSEIF(KEYWRD(ISTART:ISTOP).EQ.'OPT=STEEP')THEN
         WRITE(IOUT,'(" OPT=STEEP",6X,
     .        "- CARRY OUT A STEEPEST-DESCENT GEOMETRY",
     .        " OPTIMIZATION")')

!      ELSEIF(KEYWRD(ISTART:ISTOP).EQ.'PRTSUB')THEN
!         WRITE(IOUT,'(" PRTSUB",9X,
!     .        "- PRINT SUBSYSTEM ATOM LISTS")')
!        prtsub = .true.

      ELSEIF(KEYWRD(ISTART:ISTOP).EQ.'PRTVEC')THEN
        WRITE(IOUT,'(" PRTVEC",9X,
     .        "- PRINT FINAL EIGENVECTORS")')
        prtvec = .true.

      ELSEIF(KEYWRD(ISTART:ISTOP).EQ.'DOS')THEN
        WRITE(IOUT,'(" DOS",12x,
     .        "- PRINT DENSITY OF STATES")')
        dos = .true.

      ELSEIF(KEYWRD(ISTART:ISTOP).EQ.'RESTART')THEN
         WRITE(IOUT,'(" RESTART",8X,
     .        "- READ IN RESTART FILE FROM PREVIOUS RUN")')

      elseif (keywrd(istart:istop).eq.'GUESS') then
         write (iout,'(" GUESS",10X,
     .   "- CONSTRUCT INITIAL DENSITY MATRIX FROM OLD ",
     &        "DENSITY MATRICES")')
         guess = .true.

!!      elseif (keywrd(istart:istop).eq.'IPOLYN') then
!!        write (iout,'(" IPOLYN",10X,
!!    .   "- CONSTRUCT INITIAL DENSITY OF THE CURRENT MD ",
!!    &        "STEP FROM THE THREE PREVIOUS STEPS")')
!         ipolyn = .true.

c      ELSEIF(KEYWRD(ISTART:ISTOP).EQ.'RESIDUE')THEN
!        WRITE(IOUT,'(" RESIDUE",8X,
!    .        "- STORE RESIDUE POINTERS")')
c         resdue = .true.

!      ELSEIF(KEYWRD(ISTART:ISTOP).EQ.'SCREEN')THEN
!!        WRITE(IOUT,'(" SCREEN",9X,
!!    .        "- OUTPUT VITAL INFORMATION TO SCREEN")')
!         screen = .true.

      ELSEIF(INDEX(KEYWRD(ISTART:ISTOP),'SHIFT=').NE.0)THEN
         CALL RDNUM(KEYWRD,ISTART,SHIFT,IERROR)
         IF(IERROR.NE.0) GO TO 700
         WRITE(IOUT,'(1X,15A1,"- INITIAL DYNAMIC LEVEL SHIFT",
     .        " PARAMETER =",F6.1," eV")')
     .        (KEYWRD(I:I),I=ISTART,ISTOPP),
     .        (' ',I=ISTOPP+1,ISTART+14),SHIFT

!      ELSEIF(KEYWRD(ISTART:ISTOP).EQ.'STANDARD')THEN
!!        WRITE(IOUT,'(" STANDARD",7X,
!!    .        "- STANDARD CLOSED-SHELL CALCULATION (NO D&C)")')
!         stand = .true.

      ELSEIF(INDEX(KEYWRD(ISTART:ISTOP),'TEMPK=').NE.0)THEN
         CALL RDNUM(KEYWRD,ISTART,TEMPK,IERROR)
         IF(IERROR.NE.0) GO TO 700
         IF(TEMPK.LT.99999.9)THEN
            WRITE(IOUT,'(1X,15A1,"- A D&C TEMPERATURE OF ",
     .           F7.1," KELVIN REQUESTED")')
     .           (KEYWRD(I:I),I=ISTART,ISTOPP),
     .           (' ',I=ISTOPP+1,ISTART+14),TEMPK
         ELSE
            WRITE(IOUT,'(1X,15A1,"- A D&C TEMPERATURE OF ",
     .           G10.3," KELVIN REQUESTED")')
     .           (KEYWRD(I:I),I=ISTART,ISTOPP),
     .           (' ',I=ISTOPP+1,ISTART+14),TEMPK
         ENDIF

      ELSEIF(KEYWRD(ISTART:ISTOP).EQ.'TESTRUN')THEN
         WRITE(IOUT,'(" TESTRUN",8X,
     .        "- DO SETUP WORK AND STOP BEFORE FIRST ENERGY",
     .        " EVALUATION")')

      ELSEIF(INDEX(KEYWRD(ISTART:ISTOP),'TMAX=').NE.0)THEN
         CALL RDNUM(KEYWRD,ISTART,TMAX,IERROR)
         IF(IERROR.NE.0) GO TO 700
         IF(TMAX.LT.99999.9)THEN
            WRITE(IOUT,'(1X,15A1,"- A MAXIMUM CPU TIME OF ",
     .           F7.1," SECONDS REQUESTED")')
     .           (KEYWRD(I:I),I=ISTART,ISTOPP),
     .           (' ',I=ISTOPP+1,ISTART+14),TMAX
         ELSE
            WRITE(IOUT,'(1X,15A1,"- A MAXIMUM CPU TIME OF ",
     .           G10.3," SECONDS REQUESTED")')
     .           (KEYWRD(I:I),I=ISTART,ISTOPP),
     .           (' ',I=ISTOPP+1,ISTART+14),TMAX
         ENDIF
      ELSEIF(KEYWRD(ISTART:ISTOP).EQ.'TRAJECTORY')THEN
         WRITE(IOUT,'(" TRAJECTORY",5X,
     .        "- DUMP COORDINATES TO TRAJECTORY FILE AT",
     .        " RESTART POINTS")')
      ELSEIF(KEYWRD(ISTART:ISTOP).EQ.'WRTPDB')THEN
         WRITE(IOUT,'(" WRTPDB",9X,"- WRITE FINAL COORDINATES IN PDB",
     .        " FORMAT")')
      ELSEIF(INDEX(KEYWRD(ISTART:ISTOP),'XTEST=').NE.0)THEN
         CALL RDNUM(KEYWRD,ISTART,XTEST,IERROR)
         IF(IERROR.NE.0) GO TO 700
         WRITE(IOUT,'(1X,15A1,"- GEOMETRY OPTIMIZATION COORDINATE",
     .        " CHANGE CRITERION =",F8.5)')
     .        (KEYWRD(I:I),I=ISTART,ISTOPP),
     .        (' ',I=ISTOPP+1,ISTART+14),XTEST
      ELSEIF(KEYWRD(ISTART:ISTOP).EQ.'XYZSPACE')THEN
         WRITE(IOUT,'(" XYZSPACE",7X,
     .        "- DO ALL OPERATIONS IN XYZ SPACE")')
         xyzspace = .true.

      elseif(index(keywrd(istart:istop),'INTER').ne.0) then
         write(iout,'(1X,15A1,"- INCLUDE ONLY INTERMOLECULAR ",
     &        "INTERACTIONS FOR",/18X,"THE GRADIENT (DEFAULT ",
     &        "FOR MC SIMULATIONS)")')
     &        (KEYWRD(I:I),I=ISTART,ISTOPP),
     &        (' ',I=ISTOPP+1,ISTART+14)
         inter = .true.

      elseif(index(keywrd(istart:istop),'RECIPINTRA').ne.0) then
         write(iout,'(1X,15A1,"- INCLUDE INTRAMOLECULAR ",
     &        "CONTRIBUTIONS OF",/18X,"THE RECIPROCAL ENERGY ",
     &        "TO THE GRADIENT")')
     &        (KEYWRD(I:I),I=ISTART,ISTOPP),
     &        (' ',I=ISTOPP+1,ISTART+14)
         recipintr = .true.

      elseif(index(keywrd(istart:istop),'RMIN=').ne.0) then
         call rdnum(keywrd,istart,rijmin,ierror)
         if(ierror.ne.0) goto 700
         dangle = abs(dangle)
         write(iout,'(1X,15A1,"- MIMIMUM DISTANCE BETWEEN ",
     &        "ATOMS IS ",F8.3, " A")')
     &        (KEYWRD(I:I),I=ISTART,ISTOPP),
     &        (' ',I=ISTOPP+1,ISTART+14),rijmin

C     charge data
      elseif(index(keywrd(istart:istop),'MULLIKEN').ne.0) then
         write(iout,'(1X,15A1,"- USE MULLIKEN CHARGES FOR ",
     $        "PME/SCRF")')
     &        (KEYWRD(I:I),I=ISTART,ISTOPP),
     &        (' ',I=ISTOPP+1,ISTART+14)
         cm1 = .false.
         cm2 = .false.

      elseif(index(keywrd(istart:istop),'CM1').ne.0) then
         write(iout,'(1X,15A1,"- USE CM1 CHARGES FOR PME/SCRF")')
     &        (KEYWRD(I:I),I=ISTART,ISTOPP),
     &        (' ',I=ISTOPP+1,ISTART+14)
         cm1 = .true.
         cm2 = .false.

      elseif(index(keywrd(istart:istop),'CM2').ne.0) then
         write(iout,'(1X,15A1,"- USE CM2 CHARGES FOR PME/SCRF")')
     &        (KEYWRD(I:I),I=ISTART,ISTOPP),
     &        (' ',I=ISTOPP+1,ISTART+14)
         cm2 = .true.
         cm1 = .false.

C     external parameter data
      elseif(index(keywrd(istart:istop),'EXTERNAL').ne.0) then
         i = index(keywrd(istart:istop),'=')
         ii = min(20,istop-i-istart)
         fname(26) = keywrd(i+istart:i+ii+istart)
         extrn = .true.
         write(iout,'(1X,15A1,"- READ PARAMETERS FROM ",a20)')
     &        (KEYWRD(I:I),I=ISTART,ISTOPP),
     &        (' ',I=ISTOPP+1,ISTART+14), fname(26)

      elseif(index(keywrd(istart:istop),
     &        'NOTNAMED=AM1').ne.0) then
         notnamed = 2
         write(iout,'(1X,15A1,"- INITIALIZE PARAMETERS TO AM1 ",
     &        "PARAMETERS")')
     &        (KEYWRD(I:I),I=ISTART,ISTOPP),
     &        (' ',I=ISTOPP+1,ISTART+14)

      elseif (index(keywrd(istart:istop),'PRTPAR').ne.0) then
         prtpar = .true.
         write(iout,'(1X,15A1,"- PRINT SEMI-EMPIRICAL PARAMETERS")')
     &        (KEYWRD(I:I),I=ISTART,ISTOPP),
     &        (' ',I=ISTOPP+1,ISTART+14)

      elseif(index(keywrd(istart:istop),
     &        'NOTNAMED=PM3').ne.0) then
         notnamed = 3
         write(iout,'(1X,15A1,"- INITIALIZE PARAMETERS TO PM3 ",
     &        "PARAMETERS")')
     &        (KEYWRD(I:I),I=ISTART,ISTOPP),
     &        (' ',I=ISTOPP+1,ISTART+14)

      elseif(index(keywrd(istart:istop),
     &        'NOTNAMED=MNDO').ne.0) then
         notnamed = 1
         write(iout,'(1X,15A1,"- INITIALIZE PARAMETERS TO MNDO ",
     &        "PARAMETERS")')
     &        (KEYWRD(I:I),I=ISTART,ISTOPP),
     &        (' ',I=ISTOPP+1,ISTART+14)

C     MC data

      elseif(index(keywrd(istart:istop),'FROZENMC').ne.0) then
         write(iout,'(1X,15A1,"- FROZEN DENSITY MATRIX ",
     &        "APPROACH FOR MC")')
     &        (KEYWRD(I:I),I=ISTART,ISTOPP),
     &        (' ',I=ISTOPP+1,ISTART+14)
         frozenmc = .true.

      elseif(index(keywrd(istart:istop),'SMARTMC=').ne.0) then
         call rdinum(keywrd,istart,nmc,ierror)
         if(ierror.ne.0) goto 700
         smartmc = .true.
         write(iout,'(1X,15A1,"- NUMBER OF SMART MC CYCLES IS ",i5)')
     &        (KEYWRD(I:I),I=ISTART,ISTOPP),
     &        (' ',I=ISTOPP+1,ISTART+14),nmc
         im = mod(nmc, nproc)
         nmc = nmc / nproc
C     .  mc calculation must always be a direct calculation
C     .  because of the virial calculation in gvir
         direct = .true.
         if ((nmc .le. 0) .or. (im .ne. 0)) then
            write(iout,'(/" THE NUMBER OF MONTE CARLO CYCLES ",
     &           "SHOULD BE DIVIDABLE BY",
     &           /" THE NUMBER OF PROCESSORS (",I4,").")') nproc
            ierror = 1
            return
         endif

      elseif(index(keywrd(istart:istop),'MC=').ne.0) then
         call rdinum(keywrd,istart,nmc,ierror)
         if(ierror.ne.0) goto 700
         write(iout,'(1X,15A1,"- NUMBER OF MC CYCLES IS ",i5)')
     &        (KEYWRD(I:I),I=ISTART,ISTOPP),
     &        (' ',I=ISTOPP+1,ISTART+14),nmc
         im = mod(nmc, nproc)
         nmc = nmc / nproc
C     .  mc calculation must always be a direct calculation
C     .  because of the virial calculation in gvir
         direct = .true.
         if ((nmc .le. 0) .or. (im .ne. 0)) then
            write(iout,'(/" THE NUMBER OF MONTE CARLO CYCLES ",
     &           "SHOULD BE DIVIDABLE BY",
     &           /" THE NUMBER OF PROCESSORS (",I4,").")') nproc
            ierror = 1
            return
         endif

      elseif(index(keywrd(istart:istop),'SEED=').ne.0) then
         call rdinum(keywrd,istart,iseed,ierror)
         if(ierror.ne.0) goto 700
         write(iout,'(1X,15A1,"- SEED FOR RANDOM NUMBER GENERATOR IS ",
     &        i10)') (KEYWRD(I:I),I=ISTART,ISTOPP),
     &        (' ',I=ISTOPP+1,ISTART+14),iseed
         iseed = -abs(iseed)

      elseif(index(keywrd(istart:istop),'SMARTA=').ne.0) then
         call rdnum(keywrd,istart,smarta,ierror)
         if(ierror.ne.0) goto 700
         write(iout,'(1X,15A1,"- HALVE THE VARIANCE OF THE ",
     &        "GAUSSIAN DISTRIBUTION",
     &        /,18X,"FOR SMART MC IS ",
     &        F8.3, " A^2")')
     &        (KEYWRD(I:I),I=ISTART,ISTOPP),
     &        (' ',I=ISTOPP+1,ISTART+14),smarta

      elseif(index(keywrd(istart:istop),'DANGLE=').ne.0) then
         call rdnum(keywrd,istart,dangle,ierror)
         if(ierror.ne.0) goto 700
         dangle = abs(dangle)
         write(iout,'(1X,15A1,"- MAXIMUM ROTATION ANGLE FOR MC IS ",
     &        F8.3, " DEG")')
     &        (KEYWRD(I:I),I=ISTART,ISTOPP),
     &        (' ',I=ISTOPP+1,ISTART+14),dangle
         dangle = dangle*radian

      elseif(index(keywrd(istart:istop),'DBOX=').ne.0) then
         call rdnum(keywrd,istart,drboxmc,ierror)
         if(ierror.ne.0) goto 700
         drboxmc = abs(drboxmc)
         write(iout,'(1X,15A1,"- MAXIMUM DISPLACEMENT BOX (1 DIM.)",
     &        " FOR MC IS ",F8.3, " A")')
     &        (KEYWRD(I:I),I=ISTART,ISTOPP),
     &        (' ',I=ISTOPP+1,ISTART+14),drboxmc

       elseif(index(keywrd(istart:istop),'DR=').ne.0) then
         call rdnum(keywrd,istart,drmc,ierror)
         if(ierror.ne.0) goto 700
         drmc = abs(drmc)
         write(iout,'(1X,15A1,"- MAXIMUM DISPLACEMENT RESIDUES",
     &        " FOR MC IS ",F8.3, " A")')
     &        (KEYWRD(I:I),I=ISTART,ISTOPP),
     &        (' ',I=ISTOPP+1,ISTART+14),drmc

      elseif(index(keywrd(istart:istop),'NMOVE=').ne.0) then
         call rdinum(keywrd,istart,nmove,ierror)
         if(ierror.ne.0) goto 700
         nmove = abs(nmove)
         write(iout,'(1X,15A1,"- NUMBER OF MC MOVES A TIME IS ",
     &        i3)')
     &        (KEYWRD(I:I),I=ISTART,ISTOPP),
     &        (' ',I=ISTOPP+1,ISTART+14),nmove

      elseif(index(keywrd(istart:istop),'NSUBMOVE=').ne.0) then
         call rdinum(keywrd,istart,nsubmove,ierror)
         if(ierror.ne.0) goto 700
         nsubmove = abs(nsubmove)
         write(iout,'(1X,15A1,"- MOVE THE CORE RESIDUES OF ",I4,
     &        "SUBSYSTEMS PER MC STEP")')
     &        (KEYWRD(I:I),I=ISTART,ISTOPP),
     &        (' ',I=ISTOPP+1,ISTART+14),nsubmove

      elseif(index(keywrd(istart:istop),'NSUBSTEP=').ne.0) then
         call rdinum(keywrd,istart,nsubstep,ierror)
         if(ierror.ne.0) goto 700
         nsubstep = abs(nsubstep)
         write(iout,'(1X,15A1,"- MOVE THE RESIDUES SUBSYSTEM-WISE",
     &        " IN ",I4," STEPS")')
     &        (KEYWRD(I:I),I=ISTART,ISTOPP),
     &        (' ',I=ISTOPP+1,ISTART+14),nsubstep

      elseif(index(keywrd(istart:istop),'UPD=').ne.0) then
         call rdinum(keywrd,istart,nupdate,ierror)
         if(ierror.ne.0) goto 700
         write(iout,'(1X,15A1,"- UPDATE NEIGHBOURLIST EVERY ",
     &        i3, "-TH MC CYCLE")')
     &        (KEYWRD(I:I),I=ISTART,ISTOPP),
     &        (' ',I=ISTOPP+1,ISTART+14),nupdate

      elseif(index(keywrd(istart:istop),'ADJ=').ne.0) then
         call rdinum(keywrd,istart,nadjst,ierror)
         if(ierror.ne.0) goto 700
         write(iout,'(1X,15A1,"- ADJUST DX TO RATIO EVERY ",
     &        i8, "-TH MC CYCLE")')
     &        (KEYWRD(I:I),I=ISTART,ISTOPP),
     &        (' ',I=ISTOPP+1,ISTART+14),nadjst

      elseif(index(keywrd(istart:istop),'RATIO=').ne.0) then
         call rdnum(keywrd,istart,ratio,ierror)
         if(ierror.ne.0) goto 700
         write(iout,'(1X,15A1,"- DESIRED RATIO ACCEPTED/REJECTED IS ",
     &        F8.3)')
     &        (KEYWRD(I:I),I=ISTART,ISTOPP),
     &        (' ',I=ISTOPP+1,ISTART+14),ratio

      elseif(index(keywrd(istart:istop),'NWRITE=').ne.0) then
         call rdinum(keywrd,istart,nwrt,ierror)
         if(ierror.ne.0) goto 700
         write(iout,'(1X,15A1,"- WRITE COORDINATES EVERY ",i5,
     &        "-TH MC CYCLE")')
     &        (KEYWRD(I:I),I=ISTART,ISTOPP),
     &        (' ',I=ISTOPP+1,ISTART+14),nwrt

      elseif(index(keywrd(istart:istop),'MCECRIT=').ne.0) then
         call rdnum(keywrd,istart,eecrtmc,ierror)
         if(ierror.ne.0) goto 700
         write(iout,'(1X,15A1,"- CONVERGENCE CRITERION FOR ENERGY ",
     &        /18X,"(INTERMEDIATE MC STEPS) IS ",E10.3)')
     &        (KEYWRD(I:I),I=ISTART,ISTOPP),
     &        (' ',I=ISTOPP+1,ISTART+14),eecrtmc

      elseif(index(keywrd(istart:istop),'MCDCRIT=').ne.0) then
         call rdnum(keywrd,istart,dencrtmc,ierror)
         if(ierror.ne.0) goto 700
         write(iout,'(1X,15A1,"- CONVERGENCE CRITERION FOR DENSITY ",
     &        "MATRIX ",/18X,"(INTERMEDIATE MC STEPS) IS ",E10.3)')
     &        (KEYWRD(I:I),I=ISTART,ISTOPP),
     &        (' ',I=ISTOPP+1,ISTART+14),dencrtmc

       elseif(keywrd(istart:istop).eq.'NVT') then
          iensemb = 0
          write(iout,'(" NVT            - SAMPLE NVT ENSEMBLE")')

       elseif(keywrd(istart:istop).eq.'NPT') then
          iensemb = 1
          write(iout,'(" NPT            - SAMPLE NPT ENSEMBLE")')


      elseif(index(keywrd(istart:istop),'MCTEMP=').ne.0) then
         call rdnum(keywrd,istart,tempmc,ierror)
         if(ierror.ne.0) goto 700
         write(iout,'(1X,15A1,"- TEMPERATURE FOR MC SIMULATION IS ",
     &        F8.3," K")')
     &        (KEYWRD(I:I),I=ISTART,ISTOPP),
     &        (' ',I=ISTOPP+1,ISTART+14),tempmc
C     .  boltzt = 1/kT
         boltzt = 8.61733D-5*tempmc
         boltzt = 1.0/boltzt

      elseif(index(keywrd(istart:istop),'MCPRES=').ne.0) then
         call rdnum(keywrd,istart,press,ierror)
         if(ierror.ne.0) goto 700
         write(iout,'(1X,15A1,"- PRESSURE FOR MC SIMULATION IS ",
     &        F8.3," BAR")')
     &        (KEYWRD(I:I),I=ISTART,ISTOPP),
     &        (' ',I=ISTOPP+1,ISTART+14),press
C     .  convert press so press(bar)*V(A^3) is in eV
         press = 6.241457D-7*press

C     convergence criteria:

      elseif(index(keywrd(istart:istop),'ECRIT=').ne.0) then
         call rdnum(keywrd,istart,eecrt,ierror)
         if(ierror.ne.0) goto 700
         write(iout,'(1X,15A1,"- CONVERGENCE CRITERION FOR ENERGY IS ",
     &        E10.3," EV")')
     &        (KEYWRD(I:I),I=ISTART,ISTOPP),
     &        (' ',I=ISTOPP+1,ISTART+14),eecrt

      elseif(index(keywrd(istart:istop),'DCRIT=').ne.0) then
         call rdnum(keywrd,istart,dencrt,ierror)
         if(ierror.ne.0) goto 700
         write(iout,'(1X,15A1,"- CONVERGENCE CRITERION FOR DENSITY ",
     &        "MATRIX IS ",E10.3)')
     &        (KEYWRD(I:I),I=ISTART,ISTOPP),
     &        (' ',I=ISTOPP+1,ISTART+14),dencrt

#ifdef CUTREPUL_IS_ON

C     cutoff repulsion integrals
      elseif (index(keywrd(istart:istop),'CUTREPUL=').ne.0) then
         call rdnum(keywrd,istart,cutrepul,ierror)
         if(ierror.ne.0) goto 700
         cutrepul = abs(cutrepul)
         write(iout,'(1X,15A1,"- CUTOFF FOR REPULSION INTEGRALS IS ",
     &        E10.3)')
     &        (KEYWRD(I:I),I=ISTART,ISTOPP),
     &        (' ',I=ISTOPP+1,ISTART+14),cutrepul

#endif

C     grid subsetting
      elseif(index(keywrd(istart:istop),'ATGRID').ne.0) then
         write(iout,'(" ATGRID",9X,
     &        "- ATOM WISE, GRID BASED SUBSETTING")')
         atgr = .true.

      elseif(index(keywrd(istart:istop),'MIXGRID').ne.0) then
         write(iout,'(" MIXGRID",8X,
     &        "- RESIDUE WISE, GRID BASED SUBSETTING FOR CORE,"/
     &        18X,"ATOM WISE, GRID BASED SUBSETTING FOR BUFFERS")')
         mixgr = .true.

      elseif(index(keywrd(istart:istop),'RESGRID').ne.0) then
         write(iout,'(" RESGRID",8X,
     &        "- RESIDUE WISE, GRID BASED SUBSETTING")')
         resgr = .true.

C     combination subsetting
      elseif(index(keywrd(istart:istop),'COMBSUB').ne.0) then
         write(iout,'(1x,15a1,
     &        "- DO A COMBINATION SUBSETTING")')
     &        (KEYWRD(I:I),I=ISTART,ISTOPP),
     &        (' ',I=ISTOPP+1,ISTART+14)
         combsub = .true.

C     PME-stuff

!      elseif(index(keywrd(istart:istop),'PME').ne.0) then
!         write(iout,'(1x,15a1,
!     &        "- PME-QM IS USED ")')
!     &        (KEYWRD(I:I),I=ISTART,ISTOPP),
!     &        (' ',I=ISTOPP+1,ISTART+14)
!         pme = .true.

!      elseif(index(keywrd(istart:istop),'EWALDQM').ne.0) then
!         write(iout,'(1x,15a1,
!     &        "- EWALD IS USED ")')
!     &        (KEYWRD(I:I),I=ISTART,ISTOPP),
!     &        (' ',I=ISTOPP+1,ISTART+14)
!         pmeqm = .true.

!      elseif(index(keywrd(istart:istop),'MULLEWALD').ne.0) then
!         write(iout,'(1x,15a1,
!     &        "- WITH MULLIKEN CHARGES ")')
!     &        (KEYWRD(I:I),I=ISTART,ISTOPP),
!     &        (' ',I=ISTOPP+1,ISTART+14)
!         mullewald = .true.

!      elseif(index(keywrd(istart:istop),'CMEWALD').ne.0) then
!         write(iout,'(1x,15a1,
!     &        "- WITH CM1 CHARGES      ")')
!     &        (KEYWRD(I:I),I=ISTART,ISTOPP),
!     &        (' ',I=ISTOPP+1,ISTART+14)
!         cmewald = .true.
!
!      elseif(index(keywrd(istart:istop),'CHEWALD').ne.0) then
!c       use constant charges for ewald calculation
!         chewald = .true.

      elseif(index(keywrd(istart:istop),'BETAPME=HBOX').ne.0) then
         write(iout,'(1x,15a1,
     &        "- SET BETA FOR PME TO 3.5/(0.5*BOXLENGTH)")')
     &        (KEYWRD(I:I),I=ISTART,ISTOPP),
     &        (' ',I=ISTOPP+1,ISTART+14)
         setbpme = .true.
         pme = .true.

      elseif(index(keywrd(istart:istop),'BETAPME=').ne.0) then
         call rdnum(keywrd,istart,betapme,ierror)
         if (ierror.ne.0) goto 700
         write(iout,'(1x,15a1,
     &        "- BETA FOR PME IS ",F8.5)')
     &        (KEYWRD(I:I),I=ISTART,ISTOPP),
     &        (' ',I=ISTOPP+1,ISTART+14),betapme
         betapme2 = betapme*betapme
         pme = .true.

      elseif(index(keywrd(istart:istop),'NSPLINE=').ne.0) then
         call rdinum(keywrd,istart,nspline,ierror)
         if (ierror.ne.0) goto 700
         write(iout,'(1x,15a1,
     &        "- ORDER OF CARDINAL B-SPLINE FOR PME IS ",I4)')
     &        (KEYWRD(I:I),I=ISTART,ISTOPP),
     &        (' ',I=ISTOPP+1,ISTART+14),nspline

         if ((nspline.gt.maxspline).or.(nspline.le.2)) then
            write(iout,'(/"ERROR: NSPLINE SHOULD BE BETWEEN 2 AND ",I4,
     &           ".",/"        DECREASE NSPLINE OR INCREASE MAXSPLINE",
     &           " IN divcon.h")') maxspline
            ierror = 1
            return
         elseif (mod(nspline,2).ne.0) then
            write(iout,'(/"ERROR: NSPLINE SHOULD BE AN EVEN NUMBER")')
            ierror = 1
            return
         endif
         pme = .true.

      elseif(index(keywrd(istart:istop),'K1PME=').ne.0) then
         call rdinum(keywrd,istart,k1pme,ierror)
         if (ierror.ne.0) goto 700
         write(iout,'(1x,15a1,
     &        "- NUMBER OF GRIDS IN X DIRECTION FOR PME IS ",I4)')
     &        (KEYWRD(I:I),I=ISTART,ISTOPP),
     &        (' ',I=ISTOPP+1,ISTART+14),
     &        k1pme
         xk1pme = k1pme+0.1
         n1 = log(xk1pme)/xlog2
         k1 = 2**n1
         if (k1.ne.k1pme) then
            write(iout,'(/" ERROR: K1PME IS NOT A POWER OF 2")')
            ierror = 1
            return
         elseif (k1pme.gt.maxkpme) then
            write(iout,'(/
     &           " ERROR: K1PME SHOULD BE EQUAL TO OR SMALLER ",
     &           "THAN ",I4,".",
     &           /"        DECREASE K1PME OR INCREASE MAXKPME IN",
     &           " divcon.h")') maxkpme
            ierror = 1
            return
         endif
         pme = .true.

      elseif(index(keywrd(istart:istop),'K2PME=').ne.0) then
         call rdinum(keywrd,istart,k2pme,ierror)
         if (ierror.ne.0) goto 700
         write(iout,'(1x,15a1,
     &        "- NUMBER OF GRIDS IN Y DIRECTION FOR PME IS ",I4)')
     &        (KEYWRD(I:I),I=ISTART,ISTOPP),
     &        (' ',I=ISTOPP+1,ISTART+14),
     &        k2pme
         xk2pme = k2pme+0.1
         n2 = log(xk2pme)/xlog2
         k2 = 2**n2
         if (k2.ne.k2pme) then
            write(iout,'(/" ERROR: K2PME IS NOT A POWER OF 2")')
            ierror = 1
            return
         elseif (k2pme.gt.maxkpme) then
            write(iout,'(
     &          /" ERROR: K2PME SHOULD BE EQUAL TO OR SMALLER ",
     &           "THAN ",I4,".",
     &           /"        DECREASE K2PME OR INCREASE MAXKPME IN",
     &           " divcon.h")') maxkpme
            ierror = 1
            return
         endif
         pme = .true.

      elseif(index(keywrd(istart:istop),'K3PME=').ne.0) then
         call rdinum(keywrd,istart,k3pme,ierror)
         if (ierror.ne.0) goto 700
         write(iout,'(1x,15a1,
     &        "- NUMBER OF GRIDS IN Z DIRECTION FOR PME IS ",I4)')
     &        (KEYWRD(I:I),I=ISTART,ISTOPP),
     &        (' ',I=ISTOPP+1,ISTART+14),
     &        k3pme
         xk3pme = k3pme+0.1
         n3 = log(xk3pme)/xlog2
         k3 = 2**n3
         if (k3.ne.k3pme) then
            write(iout,'(/" ERROR: K3PME IS NOT A POWER OF 2")')
            ierror = 1
            return
         elseif (k3pme.gt.maxkpme) then
            write(iout,'(
     &           /" ERROR: K3PME SHOULD BE EQUAL TO OR SMALLER ",
     &           "THAN ",I4,".",
     &           /"        DECREASE K3PME OR INCREASE MAXKPME IN",
     &           " divcon.h")') maxkpme
            ierror = 1
            return
         endif
         pme = .true.

      elseif(index(keywrd(istart:istop),'PME').ne.0) then
         pme = .true.

      ELSE
         IERROR = 1
         WRITE(IOUT,'(/" UNRECOGNIZED KEYWORD: ",15A1)')
     .        (KEYWRD(I:I),I=ISTART,ISTOP)
         RETURN
      ENDIF
      ISTART = ISTOP+1
 400  IF(ISTART.EQ.0)THEN
C
C     THERE ARE NO MORE WORDS IN THE STRING KEYWRD.  DO SOME CONSISTENCY
C     CHECKING, THEN RETURN.
C
C     .  convert smarta to A*beta and include the energy conversion
C     .  term for the forces
         smarta2 = smarta*2.0
         smarta = smarta*boltzt/23.061

C     .  Check for redundancy of keywords in z-matrix creation cases.

         IF (zgen) THEN
            IF (.not.cart) THEN
               write(iout,'(" UNNECESSARY COMBINATION",
     .              " - INTERNAL MEANS ZMAKE IS NOT REQUIRED.")')
               zgen=.false.
            ELSEIF(rotat) then
               write(iout,'(" UNNECESSARY COMBINATION - ROTATE",
     .              " AUTOMATICALLY FORMS A Z-MATRIX.")')
               zgen=.false.
            ENDIF
         ENDIF

         if (prtvec.and.(.not.frozenmc)
     $        .and.(.not.stand)) then
            pvec = .true.
         else
            pvec = .false.
         endif

         if (index(keywrd,'MC=').ne.0) mcsim = .true.
         if (mcsim) inter = .true.
C     .  set recipintra true in order to save time, if inter is false
         if (.not.inter) recipintr = .true.
         if ((index(keywrd,'GRID').ne.0).or.combsub) gridsub = .true.
         submv = (nsubmove.ne.0).and.(nsubstep.ne.0)
         if (submv) then
            write(iout,'(/" ERROR: USE EITHER NSUBMOVE OR NSUBSTEP")')
            ierror = 1
            return
         endif
         submv = (nsubmove.ne.0).or.(nsubstep.ne.0)
         if (((nmove.ne.0).and.submv.and.mcsim)
     &        .or.(mcsim.and.((nmove.eq.0).and.(nsubmove.eq.0)
     &        .and.(nsubstep.eq.0)))) then
            write(iout,'(/" ERROR: USE EITHER NMOVE, NSUBMOVE",
     &           " OR NSUBSTEP TO DEFINE HOW THE MONTE CARLO",
     &           /"        CONFIGURATIONS ARE GENERATED.")')
            ierror = 1
            return
         endif
         if (submv.and.stand) then
            write(iout,'(/" ERROR: NSUBMOVE OR NSUBSTEP CAN ONLY ",
     &           "BE USED IN A DIVIDE AND CONQUER CALCULATION")')
            ierror = 1
            return
         endif
         if (mcsim.and.gridsub) then

C     .     in this case only a residue-wise grid-based subsetting
C     .     scheme is allowed.
C     .     Reason for this is rather subtle:
C
C     .     Imagine that during the MC-simulation a molecule would
C     .     penetrate the box, such that the geometric center is still
C     .     inside the box, but some atoms are outside the box.
C     .     If an atombased subsetting was performed, the atoms outside
C     .     the box wouldn't be included in any subsystem.
C     .     Making the "grid-subsetting"-box artificially larger than
C     .     the pbc-box wouldn't work either: in that case there's is
C     .     an artificially larger distance between the molecules and
C     .     the (virtual, pbc) images of other molecules. This would mean
C     .     that some atoms will be skipped in making the buffer regions:
C     .     atoms that, according to their pbc-image should be included.
C     .     This will lead to non-optimal subsettings and can have a
C     .     rather drastic effect on energies as was found experimentally.

            if (index(keywrd,'RESGRID').eq.0) then
               write(iout,'(/" ERROR: A MONTE-CARLO SIMULATION CAN ",
     &              "ONLY BE DONE BY EITHER A STANDARD ",
     &              /"        CALCULATION OR A DIVIDE AND CONQUER ",
     &              "CALCULATION USING CLUSTER-",
     &              /"        BASED SUBSETTING OR RESIDUE-WISE ",
     &              "GRID-BASED SUBSETTING. NO ",
     &              /"        OTHER GRID-GASED SUBSETTING SCHEMES ",
     &              "ARE ALLOWED.")')
               ierror = 1
               return
            endif
         endif

C     .  warning if gradient calculation is requested while doing a
C     .  MC simulation
         if (mcsim.and.gradient) then
            write(iout,'(/" WARNING: THE GRADIENT FOR A MC SIMULATION",
     &           " ONLY CONTAINS INTERMOLECULAR TERMS",
     &           /"          (NO INTRAMOLECULAR TERMS)"/)')
         endif

C     .  mc-angle should be less than 180 deg because of accuracy
         if (dangle .ge. pi) then
            write(iout,'(/" DANGLE SHOULD BE LESS THAN 180 DEG")')
            ierror = 1
            return
         endif

         IF(INDEX(KEYWRD,'OPT=BFGS').NE.0.AND.
     .        INDEX(KEYWRD,'XYZSPACE').NE.0)THEN
C
C     .     BFGS OPTIMIZATION IN XYZ COORDINATES WOULD RESULT IN SINGULAR
C     .     HESSIAN MATRIX.
C
            IERROR = 1
            WRITE(IOUT,'(/" CANNOT DO A BFGS OPTIMIZATION ON XYZ",
     .           " COORDINATES")')
         ENDIF
         IF(INDEX(KEYWRD,'CHGCUT=').NE.0.AND.
     .        INDEX(KEYWRD,'RESIDUE').EQ.0)THEN
            IERROR = 1
            WRITE(IOUT,'(/" MUST SPECIFY ''RESIDUE'' KEYWORD IN ORDER",
     .           " TO USE ''CHGCUT='' OPTION")')
         ENDIF


         if ((.not.stand).and.(.not.combsub).and.(.not.clust)
     &        .and.(.not.gridsub).and.(.not.resclust)
     &        .and.(.not.atmsub).and.(.not.ressub)) handsub = .true.

         if ((k1pme.lt.nspline).or.(k2pme.lt.nspline).or.
     &        (k3pme.lt.nspline)) then
            write(iout,'(/" ERROR: THE NUMBER OF GRIDS IN ANY ",
     &           "DIRECTION (K1PME ETC.) SHOULD BE GREATER THAN",
     &           /" OR EQUAL TO NSPLINE (",I4,").")') nspline
            ierror = 1
         endif

         if (push) then
            if (mcsim) then
               write(iout,'(/" ERROR: PUSH CANNOT BE USED IN ",
     &              "COMBINATION WITH MC")')
               ierror = 1
               return
            elseif (pme) then
               write(iout,'(/" ERROR: PUSH CANNOT BE USED IN ",
     &              "COMBINATION WITH PME")')
               ierror = 1
               return
            elseif(gridsub.or.combsub) then
               write(iout,'(/" ERROR: PUSH CAN ONLY BE USED IN ",
     &              "COMBINATION WITH CLUSTER")')
               ierror = 1
               return
            elseif ((index(keywrd,'OPT=').ne.0).and.push) then
               write(iout,'(/" ERROR: PUSH CANNOT BE USED FOR ",
     &              "GEOMETRY OPTIMIZATIONS")')
               ierror = 1
               return
            endif
         endif

!         xyzspc = xyzspace.or.cart
!         mndo = .not.am1.and.
!     &          .not.pm3.and.
!     &          .not.rm1.and.
!     &          .not.am1d.and.
!     &          .not.pm3pddg
!         autosub = clust.or.resclust.or.atmsub.or.ressub
!     &             .or.(index(keywrd,'GRID').ne.0)

!         if (screen) then
!            if (mcsim) then
!               wrtmc = .true.
!            else
!               wrtscr = .true.
!            endif
!         endif
         RETURN
      ELSE
C
C     MORE WORDS MAY REMAIN.  TRY TO READ ANOTHER ONE.
C
         GO TO 100
      ENDIF
 500  IERROR = 1
      WRITE(IOUT,'(/" ERROR READING KEYWORDS FROM FILE ",A20,
     .     " NODE ", I5,
     .     /8X,"*** SUBROUTINE RDKEYS ***")') FNAME(7), myid
      RETURN
 600  IERROR = 1
      WRITE(IOUT,'(/" ERROR READING TITLE CARD FROM FILE ",A20,
     .     " NODE ", I5,
     .     /10X,"*** SUBROUTINE RDKEYS ***")') FNAME(7), myid
      RETURN
 700  IERROR = 1
      I1 = ISTART
      I2 = MIN(ISTOP,ISTART+19)
      WRITE(IOUT,'(/" IMPROPER KEYWORD SPECIFICATION ",20A1)')
     .     (KEYWRD(I:I),I=I1,I2)
      return

      END