SUBROUTINE ESP
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      INCLUDE 'SIZES'
C***********************************************************************
C
C     THIS IS A DRIVER ROUTINE FOR ELECTROSTATIC POTENTIAL GENERATION
C     WRITTEN BY K.M.MERZ FEB. 1989 AT UCSF
C
C***********************************************************************
      COMMON /KEYWRD/ KEYWRD
      CHARACTER*241 KEYWRD
C
C     SET STANDARD PARAMETERS FOR THE SURFACE GENERATION
C
      IF(INDEX(KEYWRD,'SCALE=') .NE. 0)THEN
         SCALE = READA(KEYWRD,INDEX(KEYWRD,'SCALE='))
      ELSE
         SCALE = 1.4D0
      ENDIF
C
      IF(INDEX(KEYWRD,'DEN=') .NE. 0)THEN
         DEN = READA(KEYWRD,INDEX(KEYWRD,'DEN='))
      ELSE
         DEN = 1.0D0
      ENDIF
C
      IF(INDEX(KEYWRD,'SCINCR=') .NE. 0)THEN
         SCINCR = READA(KEYWRD,INDEX(KEYWRD,'SCINCR='))
      ELSE
         SCINCR = 0.20D0
      ENDIF
C
      IF(INDEX(KEYWRD,'NSURF=') .NE. 0)THEN
         N = READA(KEYWRD,INDEX(KEYWRD,'NSURF='))
      ELSE
         N = 4
      ENDIF
C
      TIME1=SECOND()
C
C     NOW CALCULATE THE SURFACE POINTS
C
      IF(INDEX(KEYWRD,'WILLIAMS') .NE. 0) THEN
         CALL PDGRID
      ELSE
         DO 10 I = 1,N
            CALL SURFAC(SCALE,DEN,I)
            SCALE = SCALE + SCINCR
   10    CONTINUE
      ENDIF
C
C     NEXT CALCULATE THE ESP AT THE POINTS CALCULATED BY SURFAC
C
      CALL POTCAL
C
C     END OF CALCULATION
C
      TIME1=SECOND()-TIME1
      WRITE(6,20) 'TIME TO CALCULATE ESP:',TIME1,' SECONDS'
   20 FORMAT(/9X,A,F8.2,A)
      RETURN
      END
      SUBROUTINE PDGRID
C
C     ROUTINE TO CALCULATE WILLIAMS SURFACE
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      INCLUDE 'SIZES'
      DIMENSION IZ(100),XYZ(3,100),VDERW(53),DIST(100)
      DIMENSION XMIN(3),XMAX(3),COORD(3,NUMATM)
      COMMON /GEOM/   GEO(3,NUMATM)
      COMMON /GEOKST/ NATOMS,LABELS(NUMATM), NABC(3*NUMATM)
C
      COMMON /ABC/    CO(3,NUMATM),IAN(NUMATM),NATOM
      COMMON /WORK1/  POTPT(3,MESP),WORK1D(4*MESP)
      COMMON /POTESP/ XC,YC,ZC,ESPNUC,ESPELE,NESP
C
      DATA VDERW/53*0.0D0/
      VDERW(1)=2.4D0
      VDERW(5)=3.0D0
      VDERW(6)=2.9D0
      VDERW(7)=2.7D0
      VDERW(8)=2.6D0
      VDERW(9)=2.55D0
      VDERW(15)=3.1D0
      VDERW(16)=3.05D0
      VDERW(17)=3.0D0
      VDERW(35)=3.15D0
      VDERW(53)=3.35D0
      SHELL=1.2D0
      NESP=0
      GRID=0.8D0
      CLOSER=0.D0
C     CHECK IF VDERW IS DEFINED FOR ALL ATOMS
C
C     CONVERT INTERNAL TO CARTESIAN COORDINATES
C
      CALL GMETRY(GEO,COORD)
C
C     STRIP COORDINATES AND ATOM LABEL FOR DUMMIES (I.E. 99)
C
      ICNTR = 0
      DO 20 I=1,NATOMS
         DO 10 J=1,3
   10    CO(J,I) = COORD(J,I)
         IF(LABELS(I) .EQ. 99) GOTO 20
         ICNTR = ICNTR + 1
         IAN(ICNTR) = LABELS(I)
   20 CONTINUE
      NATOM=ICNTR
C
      DO 30 I=1,NATOM
         J=IAN(I)
         IF (VDERW(J).EQ.0.0D0) GO TO 40
   30 CONTINUE
      GO TO 50
   40 CONTINUE
      WRITE(6,*) 'VAN DER WAALS'' RADIUS NOT DEFINED FOR ATOM',I
      WRITE(6,*) 'IN WILLIAMS SURFACE ROUTINE PDGRID!'
      STOP
C     NOW CREATE LIMITS FOR A BOX
   50 DO 100 IX = 1,3
         XMIN(IX)= 100000.0D0
         XMAX(IX)=-100000.0D0
         DO 90 IA = 1,NATOM
            IF (CO(IX,IA)-XMIN(IX))60,70,70
   60       XMIN(IX)=CO(IX,IA)
   70       IF (CO(IX,IA)-XMAX(IX))90,90,80
   80       XMAX(IX)=CO(IX,IA)
   90    CONTINUE
  100 CONTINUE
C     ADD (OR SUBTRACT) THE MAXIMUM VDERW PLUS SHELL
      VDMAX=0.0D0
      DO 110 I=1,53
         IF (VDERW(I).GT.VDMAX) VDMAX=VDERW(I)
  110 CONTINUE
      DO 120 I=1,3
         XMIN(I)=XMIN(I)-VDMAX-SHELL
  120 XMAX(I)=XMAX(I)+VDMAX+SHELL
C STEP GRID BACK FROM ZERO TO FIND STARTING POINTS
      XSTART=0.0D0
  130 XSTART=XSTART-GRID
      IF (XSTART.GT.XMIN(1)) GO TO 130
      YSTART=0.0D0
  140 YSTART=YSTART-GRID
      IF (YSTART.GT.XMIN(2)) GO TO 140
      ZSTART=0.0D0
  150 ZSTART=ZSTART-GRID
      IF (ZSTART.GT.XMIN(3)) GO TO 150
      NPNT=0
      ZGRID=ZSTART
  160 YGRID=YSTART
  170 XGRID=XSTART
  180 DO 190 L=1,NATOM
         JZ=IAN(L)
         DIST(L)=SQRT((CO(1,L)-XGRID)**2+(CO(2,L)-YGRID)**2+
     1 (CO(3,L)-ZGRID)**2)
C     REJECT GRID POINT IF ANY ATOM IS TOO CLOSE
         IF(DIST(L).LT.(VDERW(JZ)-CLOSER)) GO TO 220
  190 CONTINUE
C BUT AT LEAST ONE ATOM MUST BE CLOSE ENOUGH
      DO 200 L=1,NATOM
         JZ=IAN(L)
         IF(DIST(L).GT.(VDERW(JZ)+SHELL)) GO TO 200
         GO TO 210
  200 CONTINUE
      GO TO 220
  210 NPNT=NPNT+1
      NESP=NESP+1
      POTPT(1,NESP)=XGRID
      POTPT(2,NESP)=YGRID
      POTPT(3,NESP)=ZGRID
  220 XGRID=XGRID+GRID
      IF (XGRID.LE.XMAX(1)) GO TO 180
      YGRID=YGRID+GRID
      IF (YGRID.LE.XMAX(2)) GO TO 170
      ZGRID=ZGRID+GRID
      IF (ZGRID.LE.XMAX(3)) GO TO 160
      RETURN
      END
C***********************************************************************
      SUBROUTINE SURFAC(SCALE,DENS,IPT)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      INCLUDE 'SIZES'
C***********************************************************************
C
C      THIS SUBROUTINE CALCULATES THE MOLECULAR SURFACE OF A MOLECULE
C      GIVEN THE COORDINATES OF ITS ATOMS.  VAN DER WAALS' RADII FOR
C      THE ATOMS AND THE PROBE RADIUS MUST ALSO BE SPECIFIED.
C
C      ON INPUT    SCALE = INITIAL VAN DER WAALS' SCALE FACTOR
C                  DENS  = DENSITY OF POINTS PER UNIT AREA
C
C      THIS SUBROUTINE WAS LIFTED FROM MICHAEL CONNOLLY'S SURFACE
C      PROGRAM FOR UCSF GRAPHICS SYSTEM BY U.CHANDRA SINGH AND
C      P.A.KOLLMAN AND MODIFIED FOR USE IN QUEST. K.M.MERZ
C      ADAPTED AND CLEANED UP THIS PROGRAM FOR USE IN AMPAC/MOPAC
C      IN FEB. 1989 AT UCSF.
C
C***********************************************************************
      COMMON /GEOM/   GEO(3,NUMATM)
      COMMON /GEOKST/ NATOMS,LABELS(NUMATM),
     1                NA(NUMATM),NB(NUMATM),NC(NUMATM)
      COMMON /KEYWRD/ KEYWRD
C
      COMMON /ABC/    CO(3,NUMATM),IAN(NUMATM),NATOM
      COMMON /WORK1/  POTPT(3,MESP),PAD1(2*MESP),RAD(MESP),IAS(MESP)
      COMMON /POTESP/ XC,YC,ZC,ESPNUC,ESPELE,NESP
C
      CHARACTER*241 KEYWRD
C
C     CARTESIAN COORDINATE AND ATOM LABELS
C
      DIMENSION COORD(3,NUMATM),VANDER(100)
      DIMENSION CON(3,1000),ROT(3,3)
C
C     NEIGHBOR ARRAYS
C
C     THIS SAME DIMENSION FOR THE MAXIMUM NUMBER OF NEIGHBORS
C     IS USED TO DIMENSION ARRAYS IN THE LOGICAL FUNCTION COLLID
C
      DIMENSION INBR(200),CNBR(3,200),RNBR(200)
      LOGICAL SNBR(200),MNBR(200)
C
C     ARRAYS FOR ALL ATOMS
C
C     IATOM, JATOM AND KATOM COORDINATES
C
      DIMENSION CI(3), IELDAT(56), TEMP0(3)
C
C     GEOMETRIC CONSTRUCTION VECTORS
C
      DIMENSION CW(3,2)
C
C     LOGICAL VARIABLES
C
      LOGICAL SI
C
C     LOGICAL FUNCTIONS
C
      LOGICAL COLLID
C
C     DATA FOR VANDER VALL RADII
C
      CHARACTER MARKER*3, MARKSS*3, MYNAM*3, IELDAT*4, NAMATM*4
      DATA VANDER/1.20D0,1.20D0,1.37D0,1.45D0,1.45D0,1.50D0,1.50D0,
     1            1.40D0,1.35D0,1.30D0,1.57D0,1.36D0,1.24D0,1.17D0,
     2            1.80D0,1.75D0,1.70D0,17*0.0D0,2.3D0,65*0.0D0/
      DATA MARKER/'A  '/,MARKSS/'SS0'/,MYNAM/'UC '/
C
      DATA IELDAT/'  BQ','  H ','  HE','  LI','  BE','  B ',
     1            '  C ','  N ','  O ','  F ','  NE','  NA',
     2            '  MG','  AL','  SI','  P ','  S ','  CL',
     3            '  AR','  K ','  CA','  SC','  TI','  V ',
     4            '  CR','  MN','  FE','  CO','  NI','  CU',
     5            '  ZN','  GA','  GE','  AS','  SE','  BR',
     6            '  KR','  RB','  SR','   Y','  ZR','  NB',
     7            '  MO','  TC','  RU','  RH','  PD','  AG',
     8            '  CD','  IN','  SN','  SB','  TE','   I',
     9            '   X','  CS'/
      PI=4.D0*ATAN(1.D0)
C     INSERT VAN DER WAAL RADII FOR ZINC
      VANDER(30)=1.00D0
C
C     CONVERT INTERNAL TO CARTESIAN COORDINATES
C
      CALL GMETRY(GEO,COORD)
C
C     STRIP COORDINATES AND ATOM LABEL FOR DUMMIES (I.E. 99)
C
      ICNTR = 0
      DO 20 I=1,NATOMS
         DO 10 J=1,3
   10    CO(J,I) = COORD(J,I)
         IF(LABELS(I) .EQ. 99) GOTO 20
         ICNTR = ICNTR + 1
         IAN(ICNTR) = LABELS(I)
   20 CONTINUE
C
C     ONLY VAN DER WAALS' TYPE SURFACE IS GENERATED
C
      IOP = 1
      RW =0.0D0
      NATOM = ICNTR
      DEN = DENS
      DO 30 I=1,NATOM
         IPOINT = IAN(I)
         RAD(I) = VANDER(IPOINT)*SCALE
         IF (RAD(I) .LT. 0.01D0) THEN
            WRITE(6,'(T2,''VAN DER WAALS'''' RADIUS FOR ATOM '',I3,
     1         '' IS ZERO, SUPPLY A VALUE IN SUBROUTINE SURFAC)''
     2         )')
         ENDIF
         IAS(I) = 2
   30 CONTINUE
C
C     BIG LOOP FOR EACH ATOM
C
      DO 110 IATOM = 1, NATOM
         IF (IAS(IATOM) .EQ. 0) GO TO 110
C
C     TRANSFER VALUES FROM LARGE ARRAYS TO IATOM VARIABLES
C
         NAMATM =IELDAT(IAN(IATOM)+1)
         RI = RAD(IATOM)
         SI = IAS(IATOM) .EQ. 2
         DO 40 K = 1,3
            CI(K) = CO(K,IATOM)
   40    CONTINUE
C
C     GATHER THE NEIGHBORING ATOMS OF IATOM
C
         NNBR = 0
         DO 60 JATOM = 1, NATOM
            IF (IATOM .EQ. JATOM .OR. IAS(JATOM) .EQ. 0) GO TO 60
            D2 = DIST2(CI,CO(1,JATOM))
            IF (D2 .GE. (2*RW+RI+RAD(JATOM)) ** 2) GO TO 60
C
C     WE HAVE A NEW NEIGHBOR
C     TRANSFER ATOM COORDINATES, RADIUS AND SURFACE REQUEST NUMBER
C
            NNBR = NNBR + 1
            IF (NNBR .GT. 200)THEN
               WRITE (6,'(''ERROR'',2X,''TOO MANY NEIGHBORS:'',I5)')NNBR
               STOP
            ENDIF
            INBR(NNBR) = JATOM
            DO 50 K = 1,3
               CNBR(K,NNBR) = CO(K,JATOM)
   50       CONTINUE
            RNBR(NNBR) = RAD(JATOM)
            SNBR(NNBR) = IAS(JATOM) .EQ. 2
   60    CONTINUE
C
C     CONTACT SURFACE
C
         IF (.NOT. SI) GO TO 110
         NCON = (4 * PI * RI ** 2) * DEN
         IF (NCON .GT. 1000) NCON = 1000
C
C     THIS CALL MAY DECREASE NCON SOMEWHAT
C
         IF ( NCON .EQ. 0) THEN
            WRITE(6,'(T2,''VECTOR LENGTH OF ZERO IN SURFAC'')')
            STOP
         ENDIF
         CALL GENUN(CON,NCON)
         AREA = (4 * PI * RI ** 2) / NCON
C
C     CONTACT PROBE PLACEMENT LOOP
C
         DO 100 I = 1,NCON
            DO 70 K = 1,3
               CW(K,1) = CI(K) + (RI + RW) * CON(K,I)
   70       CONTINUE
C
C     CHECK FOR COLLISION WITH NEIGHBORING ATOMS
C
            IF (COLLID(CW(1,1),RW,CNBR,RNBR,MNBR,NNBR,1,
     1      JNBR,KNBR)) GO TO 100
            DO 80 KK=1,3
               TEMP0(KK) =CI(KK)+RI*CON(KK,I)
   80       CONTINUE
C
C     STORE POINT IN POTPT AND INCREMENT NESP
C
            NESP = NESP + 1
            IF (NESP .GT. MESP) THEN
               WRITE(6,90)
   90          FORMAT(/'ERROR - TO MANY POINTS GENERATED IN SURFAC')
               WRITE(6,'(''    REDUCE NSURF, SCALE, DEN, OR SCINCR'')')
               STOP
            ENDIF
            POTPT(1,NESP) = TEMP0(1)
            POTPT(2,NESP) = TEMP0(2)
            POTPT(3,NESP) = TEMP0(3)
  100    CONTINUE
  110 CONTINUE
      RETURN
      END
C****************************************************************
      FUNCTION DIST2(A,B)
C
C     DETERMINE DISTANCES BETWEEN NEIGHBORING ATOMS
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION A(3)
      DIMENSION B(3)
      DIST2 = (A(1)-B(1))**2 + (A(2)-B(2))**2 + (A(3)-B(3))**2
      RETURN
      END
C****************************************************************
      LOGICAL FUNCTION COLLID(CW,RW,CNBR,RNBR,MNBR,NNBR,ISHAPE,
     1JNBR,KNBR)
C****************************************************************
C
C     COLLISION CHECK OF PROBE WITH NEIGHBORING ATOMS
C     USED BY SURFAC ONLY.
C
C****************************************************************
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION CW(3)
      DIMENSION CNBR(3,200)
      DIMENSION RNBR(200)
      LOGICAL MNBR(200)
      IF (NNBR .LE. 0) GO TO 20
C
C     CHECK WHETHER PROBE IS TOO CLOSE TO ANY NEIGHBOR
C
      DO 10 I = 1, NNBR
         IF (ISHAPE .GT. 1 .AND. I .EQ. JNBR) GO TO 10
         IF (ISHAPE .EQ. 3 .AND. (I .EQ. KNBR .OR. .NOT. MNBR(I)))
     1   GO TO 10
         SUMRAD = RW + RNBR(I)
         VECT1 = DABS(CW(1) - CNBR(1,I))
         IF (VECT1 .GE. SUMRAD) GO TO 10
         VECT2 = DABS(CW(2) - CNBR(2,I))
         IF (VECT2 .GE. SUMRAD) GO TO 10
         VECT3 = DABS(CW(3) - CNBR(3,I))
         IF (VECT3 .GE. SUMRAD) GO TO 10
         SR2 = SUMRAD ** 2
         DD2 = VECT1 ** 2 + VECT2 ** 2 + VECT3 ** 2
         IF (DD2 .LT. SR2) GO TO 30
   10 CONTINUE
   20 CONTINUE
      COLLID = .FALSE.
      GO TO 40
   30 CONTINUE
      COLLID = .TRUE.
   40 CONTINUE
      RETURN
      END
C****************************************************************
      SUBROUTINE GENUN(U,N)
C****************************************************************
C
C     GENERATE UNIT VECTORS OVER SPHERE. USED BY SURFAC ONLY.
C
C****************************************************************
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION U(3,N)
      PI=4.D0*ATAN(1.D0)
      NEQUAT = SQRT(N * PI)
      NVERT = NEQUAT/2
      NU = 0
      DO 20 I = 1,NVERT+1
         FI = (PI * (I-1)) / NVERT
         Z = COS(FI)
         XY = SIN(FI)
         NHOR = NEQUAT * XY
         IF (NHOR .LT. 1) NHOR = 1
         DO 10 J = 1,NHOR
            FJ = (2.D0 * PI * (J-1)) / NHOR
            X = DCOS(FJ) * XY
            Y = DSIN(FJ) * XY
            IF (NU .GE. N) GO TO 30
            NU = NU + 1
            U(1,NU) = X
            U(2,NU) = Y
            U(3,NU) = Z
   10    CONTINUE
   20 CONTINUE
   30 CONTINUE
      N = NU
      RETURN
      END
C***********************************************************************
      SUBROUTINE POTCAL
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      INCLUDE 'SIZES'
C***********************************************************************
C
C     THIS SUBROUTINE CALCULATES THE TOTAL ELECTROSTATIC POTENTIAL
C     THE NUCLEAR CONTRIBUTION IS EVALUATED BY NUCPOT
C     THE ELECTRONIC CONTRIBUTION IS EVALUATED BY ELESP
C     ESPFIT FITS THE QUANTUM POTENTIAL TO A CLASSICAL POINT CHARGE
C     MODEL.
C     THIS SUBROUTINE WAS WRITTEN BY B.H.BESLER AND K.M.MERZ IN FEB.
C     1989 AT UCSF
C
C***********************************************************************
      COMMON /KEYWRD/ KEYWRD
      COMMON /CORE/ TORE(107)
      COMMON /ELEMTS/ ELEMNT(107)
      COMMON /DENSTY/ P(MPACK),PA(MPACK),PB(MPACK)
      COMMON /POTESP/ XC,YC,ZC,ESPNUC,ESPELE,NESP
      COMMON /WORK1/  POTPT(3,MESP), ES(MESP), ESP(MESP), WORK1D(2*MESP)
      COMMON /ABC/    CO(3,NUMATM),IAN(NUMATM),NATOM
      COMMON /DIPSTO/ UX,UY,UZ,CH(NUMATM)
      COMMON /ESPF/  AL((NUMATM+4)**2),A(NUMATM,NUMATM),B(NUMATM),
     1Q(NUMATM+4),QSC(NUMATM+4),CF, ESPFD(MAXORB**2-NUMATM-5)
      CHARACTER*241 KEYWRD
      CHARACTER *2  ELEMNT
      LOGICAL DEBUG,WRTESP,CEQUIV(NUMATM,NUMATM)
C
C     DEBUG PRINTING - RESULTS IN COPIOUS OUTPUT
C
      DEBUG = (INDEX(KEYWRD,'DEBUG') .NE. 0)
C
C
      CALL ELESP
      BOHR = 0.529167D00
C
C     NOW FIT THE ELECTROSTATIC POTENTIAL
C
      WRITE(6,'(//12X,''ELECTROSTATIC POTENTIAL CHARGES'',/)')
      IZ=0
      IF(INDEX(KEYWRD,'CHARGE=') .NE. 0) IZ=READA(KEYWRD,INDEX(KEYWRD,
     1'CHARGE='))
C
C     DIPOLAR CONSTRAINTS IF DESIRED
C
      IF(INDEX(KEYWRD,'DIPOLE') .NE. 0) THEN
         IDIP = 1
         IF(IZ .NE. 0)THEN
            IDIP = 0
            WRITE(6,'(/12X,''  DIPOLE CONSTRAINTS NOT USED'')')
            WRITE(6,'(12X,''        CHARGED MOLECULE'',/)')
         ENDIF
      ELSE
         IDIP = 0
      ENDIF
      IF (IDIP .EQ. 1) THEN
         WRITE(6,'(/12X,''DIPOLE CONSTRAINTS WILL BE USED'',/)')
      ENDIF
C
C     GET X,Y,Z DIPOLE COMPONENTS IF DESIRED
C
      IF(INDEX(KEYWRD,'DIPX=') .NE. 0) THEN
         DX = READA(KEYWRD,INDEX(KEYWRD,'DIPX='))
      ELSE
         DX = UX
      ENDIF
      IF(INDEX(KEYWRD,'DIPY=') .NE. 0) THEN
         DY = READA(KEYWRD,INDEX(KEYWRD,'DIPY='))
      ELSE
         DY = UY
      ENDIF
      IF(INDEX(KEYWRD,'DIPZ=') .NE. 0) THEN
         DZ = READA(KEYWRD,INDEX(KEYWRD,'DIPZ='))
      ELSE
         DZ = UZ
      ENDIF
      CALL ESPFIT(IDIP,NATOM,NESP,IZ,ESP,POTPT,CO,DX,DY,DZ,RMS,RRMS)
C
C     WRITE OUT OUR RESULTS TO CHANNEL 6
C     THE CHARGES ARE SCALED TO REPRODUCE 6-31G* CHARGES FOR MNDO ONLY
C     AM1 AND MINDO/3 CHARGES ARE NOT SCALED DUE TO THE LOW COORELATION
C     COEFFICIENT. SEE BESLER,MERZ,KOLLMAN IN J. COMPUT. CHEM.
C     (IN PRESS)
C
      IF((INDEX(KEYWRD,'AM1') .NE. 0) .OR.
     1(INDEX(KEYWRD,'MINDO') .NE. 0) .OR.
     2(INDEX(KEYWRD,'PM3') .NE. 0))THEN
         WRITE(6,'(15X,''ATOM NO.    TYPE    CHARGE'')')
         DO 10 I=1,NATOM
            WRITE(6,'(17X,I2,9X,A2,1X,F10.4)')I,ELEMNT(IAN(I)),Q(I)
   10    CONTINUE
      ELSE
C
C     MNDO CALCULATION-SCALE THE CHARGES. TEST FOR SLOPE KEYWORD
C
         IF(INDEX(KEYWRD,'SLOPE=') .NE. 0) THEN
            SLOPE = READA(KEYWRD,INDEX(KEYWRD,'SLOPE='))
         ELSE
            SLOPE = 1.422D0
         ENDIF
         DO 20 I=1,NATOM
            QSC(I) = SLOPE*Q(I)
   20    CONTINUE
         WRITE(6,'(7X,''ATOM NO.    TYPE    CHARGE   SCALED CHARGE'')')
         DO 30 I=1,NATOM
            WRITE(6,'(9X,I2,9X,A2,1X,F10.4,2X,F10.4)')I,ELEMNT(IAN(I
     1)),   Q(I),QSC(I)
   30    CONTINUE
      ENDIF
      WRITE(6,'(/12X,A,4X,I6)') 'THE NUMBER OF POINTS IS:',NESP
      WRITE(6,'(12X,A,4X,F9.4)') 'THE RMS DEVIATION IS:',RMS
      WRITE(6,'(12X,A,3X,F9.4)') 'THE RRMS DEVIATION IS:',RRMS
C
C     CALCULATE DIPOLE MOMENT IF NEUTRAL MOLECULE
C
      IF (IZ .NE. 0) THEN
         GO TO 60
      ELSE
         WRITE(6,40)
   40    FORMAT (//5X,'DIPOLE MOMENT EVALUATED FROM '
     1,'THE POINT CHARGES',/)
         DO 50 I=1,NATOM
            DIPX=DIPX+CO(1,I)*Q(I)/BOHR
            DIPY=DIPY+CO(2,I)*Q(I)/BOHR
            DIPZ=DIPZ+CO(3,I)*Q(I)/BOHR
   50    CONTINUE
         DIP=SQRT(DIPX**2+DIPY**2+DIPZ**2)
         WRITE(6,'(12X,'' X        Y        Z       TOTAL'')')
         WRITE(6,'(8X,4F9.4)')DIPX*CF,DIPY*CF,DIPZ*CF,DIP*CF
      ENDIF
   60 CONTINUE
C     DETERMINE WHICH CHARGES SHOULD BE EQUIVALENT BY SYMMETRY AND
C     AVERAGE THEM IF DESIRED
      IF(INDEX(KEYWRD,'SYMAVG') .NE. 0) THEN
         DO 70 I=1,NATOM
            DO 70 J=1,NATOM
               CEQUIV(I,J)=.FALSE.
               IF(ABS(ABS(CH(I))-ABS(CH(J))) .LT. 1.D-5)  CEQUIV(I,J)=.T
     1RUE.
   70    CONTINUE
         DO 90 I=1,NATOM
            IEQ=0
            QSC(I)=0.D0
            DO 80 J=1,NATOM
               IF(CEQUIV(I,J)) THEN
                  QSC(I)=QSC(I)+ABS(Q(J))
                  IEQ=IEQ+1
               ENDIF
   80       CONTINUE
            CH(I)=Q(I)/ABS(Q(I))*QSC(I)/IEQ
   90    CONTINUE
         WRITE(6,*) ' '
         WRITE(6,*)'   ELECTROSTATIC POTENTIAL CHARGES AVERAGED FOR'
         WRITE(6,*)'   SYMMETRY EQUIVALENT ATOMS'
         WRITE(6,*) ' '
         IF((INDEX(KEYWRD,'AM1') .NE. 0) .OR.
     1(INDEX(KEYWRD,'MINDO') .NE. 0) .OR.
     2(INDEX(KEYWRD,'PM3') .NE. 0))THEN
            WRITE(6,'(7X,''ATOM NO.    TYPE    CHARGE'')')
            DO 100 I=1,NATOM
               WRITE(6,'(9X,I2,9X,A2,1X,F10.4)')I,ELEMNT(IAN(I)),
     1   CH(I)
  100       CONTINUE
         ELSE
            WRITE(6,'(7X,''ATOM NO.    TYPE    CHARGE   SCALED CHARGE'')
     1')
            DO 110 I=1,NATOM
               WRITE(6,'(9X,I2,9X,A2,1X,F10.4,2X,F10.4)')I,ELEMNT(IA
     1N(I)),   CH(I),CH(I)*SLOPE
  110       CONTINUE
         ENDIF
      ENDIF
      RETURN
      END
      SUBROUTINE ELESP
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C***********************************************************************
C     ELESP LOADS THE STO-6G BASIS SET ONTO THE ATOMS, PERFOMS THE
C     DEORTHOGONALIZATION OF THE COEFFICIENTS AND EVALUATES THE
C     ELECTRONIC CONTRIBUTION TO THE ESP. IT WAS WRITTEN BY B.H.BESLER
C     AND K.M.MERZ IN FEB. 1989 AT UCSF.
C
C***********************************************************************
      CHARACTER*241 KEYWRD
      DOUBLE PRECISION NORM,OVL
      LOGICAL CALLED,POTWRT,RST,STO3G
      INCLUDE 'SIZES'
      COMMON/ESPF/ AL((NUMATM+4)**2),A(NUMATM,NUMATM),B(NUMATM),
     1Q(NUMATM+4),CESPM(MAXORB,MAXORB)
      COMMON /DENSTY/ P(MPACK),PA(MPACK),PB(MPACK)
      COMMON /POTESP/ XC,YC,ZC,ESPNUC,ESPELE,NESP
      COMMON /ABC/    CO(3,NUMATM),IAN(NUMATM),NATOM
      COMMON /WORK1/  POTPT(3,MESP), ES(MESP), ESP(MESP), WORK1D(2*MESP)
      COMMON /STO6G/  ALLC(6,6,2),ALLZ(6,6,2)
      COMMON /VECTOR/ C(MORB2*2+MAXORB*2)
      COMMON /MOLKST/ NUMAT,NAT(NUMATM),NFIRST(NUMATM),NMIDLE(NUMATM),
     1                NLAST(NUMATM), NORBS, NELECS,NALPHA,NBETA,
     2                NCLOSE,NOPEN,NDUMY,FRACT
      COMMON /KEYWRD/ KEYWRD
      COMMON /ESPC/  CC(MAXPR),CEN(MAXPR,3),IAM(MAXPR,2),IND(MAXPR),
     1                EX(MAXPR),ESPI(MAXORB,MAXORB),
     2                FV(0:8,821),FAC(0:7),
     3                DEX(-1:96),TF(0:2),TEMP(MAXPR),ITEMP(MAXPR),
     4                OVL(MAXORB,MAXORB),FC(MAXPR*6)
     6       /CORE  / TORE(107)
     7       /EXPONT/ ZS(107),ZP(107),ZD(107)
*
*  END OF MINDO/3 COMMON BLOCKS
*
      COMMON /INDX/   INDC(MAXORB)
      DIMENSION CESPM2(MAXORB,MAXORB),SLA(10)
      DIMENSION CESPML(MAXORB*MAXORB),CESP(MAXORB*MAXORB)
      DATA BOHR/0.529167D0/
      PI=4.D0*ATAN(1.D0)
C
C     PUT STO-6G BASIS SET ON ATOM CENTERS
C
      DO 10 I=-1,10
         DEX(I)=DEX2(I)
   10 CONTINUE
      DO 20   I=0,7
         FAC(I)=1.D0/FAC(I)
   20 CONTINUE
      DO 30 M=0,8
         K=1
         FV(M,1)=1.D0/(2.D0*M+1.D0)
         DO 30 T=0.05D0,41.D0,0.05D0
            K=K+1
            CALL FSUB(M,T,FVAL)
            FV(M,K)=FVAL
   30 CONTINUE
C
C     LOAD BASIS FUNCTIONS INTO ARRAYS
C
      STO3G=(INDEX(KEYWRD,'STO3G') .NE. 0)
      IF(STO3G) THEN
         ICD=3
         CALL SETUP3
      ELSE
         ICD=6
         CALL SETUPG
      ENDIF
      NC=0
      NPR=0
      DO 80 I=1,NATOM
         IF (IAN(I) .LE. 2) THEN
            DO 40 J=1,ICD
               CC(NPR+J)=ALLC(J,1,1)
               EX(NPR+J)=ALLZ(J,1,1)*ZS(1)**2
               CEN(NPR+J,1)=CO(1,I)/BOHR
               CEN(NPR+J,2)=CO(2,I)/BOHR
               CEN(NPR+J,3)=CO(3,I)/BOHR
               IAM(NPR+J,1)=0
               IAM(NPR+J,2)=0
               FC(NPR+J)=I
   40       CONTINUE
            NC=NC+1
            NPR=NPR+ICD
         ELSE
C        DETERMINE PRINCIPAL QUANTUM NUMBER(NQN)
C        OF ORBITALS TO BE USED
C
            NQN=2
            IF(IAN(I) .GT. 10 .AND. IAN(I) .LE. 18) NQN=3
            IF(IAN(I) .GT. 18 .AND. IAN(I) .LE. 36) NQN=4
            IF(IAN(I) .GT. 36 .AND. IAN(I) .LE. 54) NQN=5
C
            DO 50 J=1,ICD
               CC(NPR+J)=ALLC(J,NQN,1)
               EX(NPR+J)=ALLZ(J,NQN,1)*ZS(IAN(I))**2
               CEN(NPR+J,1)=CO(1,I)/BOHR
               CEN(NPR+J,2)=CO(2,I)/BOHR
               CEN(NPR+J,3)=CO(3,I)/BOHR
               IAM(NPR+J,1)=0
               IAM(NPR+J,2)=0
   50       CONTINUE
            NC=NC+1
            NPR=NPR+ICD
            DO 70 K=1,3
               DO 60  J=1,ICD
                  CC(NPR+J)=ALLC(J,NQN,2)
                  EX(NPR+J)=ALLZ(J,NQN,2)*ZP(IAN(I))**2
                  CEN(NPR+J,1)=CO(1,I)/BOHR
                  CEN(NPR+J,2)=CO(2,I)/BOHR
                  CEN(NPR+J,3)=CO(3,I)/BOHR
                  IAM(NPR+J,1)=1
                  IAM(NPR+J,2)=K
   60          CONTINUE
               NC=NC+1
               NPR=NPR+ICD
   70       CONTINUE
         ENDIF
   80 CONTINUE
C
C     CALCULATE NORMALIZATION CONSTANTS AND INCLUDE
C     THEM IN THE CONTRACTION COEFFICIENTS
C
      DO 90 I=1,NPR
         NORM=(2.D0*EX(I)/PI)**0.75D0*(4.D0*EX(I))**(IAM(I,1)/2.D0)/
     1   SQRT(DEX(2*IAM(I,1)-1))
         CC(I)=CC(I)*NORM
   90 CONTINUE
      IPR=0
C
C     PERFORM SORT OF PRIMITIVES BY ANGULAR MOMENTUM
C
      IS=0
      IP=0
      IPC=0
      ISC=0
      J=0
      DO 100 I=1,NPR
         IF (IAM(I,1) .EQ. 0) THEN
            IS=IS+1
            IND(IS)=I
         ENDIF
  100 CONTINUE
      IP=IS
      DO 110 I=1,NPR
         IF (IAM(I,1) .EQ. 1 .AND. IAM(I,2) .EQ. 1) THEN
            IP=IP+1
            IND(IP)=I
         ENDIF
  110 CONTINUE
      DO 120 I=1,NPR
         IF (IAM(I,1) .EQ. 1 .AND. IAM(I,2) .EQ. 2) THEN
            IP=IP+1
            IND(IP)=I
         ENDIF
  120 CONTINUE
      DO 130 I=1,NPR
         IF (IAM(I,1) .EQ. 1 .AND. IAM(I,2) .EQ. 3) THEN
            IP=IP+1
            IND(IP)=I
         ENDIF
  130 CONTINUE
      DO 140 I=1,NC
         IN=I*ICD-ICD+1
         IF (IAM(IN,1) .EQ. 0) THEN
            ISC=ISC+1
            INDC(ISC)=I
         ENDIF
  140 CONTINUE
      IPC=ISC
      DO 150 I=1,NC
         IN=I*ICD-ICD+1
         IF (IAM(IN,1) .EQ. 1 .AND. IAM(IN,2) .EQ. 1) THEN
            IPC=IPC+1
            INDC(IPC)=I
         ENDIF
  150 CONTINUE
      DO 160 I=1,NC
         IN=I*ICD-ICD+1
         IF (IAM(IN,1) .EQ. 1 .AND. IAM(IN,2) .EQ. 2) THEN
            IPC=IPC+1
            INDC(IPC)=I
         ENDIF
  160 CONTINUE
      DO 170 I=1,NC
         IN=I*ICD-ICD+1
         IF (IAM(IN,1) .EQ. 1 .AND. IAM(IN,2) .EQ. 3) THEN
            IPC=IPC+1
            INDC(IPC)=I
         ENDIF
  170 CONTINUE
      DO 180 I=1,NPR
         TEMP(I)=CC(IND(I))
  180 CONTINUE
      DO 190 I=1,NPR
         CC(I)=TEMP(I)
  190 CONTINUE
      DO 200 I=1,NPR
         TEMP(I)=EX(IND(I))
  200 CONTINUE
      DO 210 I=1,NPR
         EX(I)=TEMP(I)
  210 CONTINUE
      DO 220 I=1,NPR
         TEMP(I)=CEN(IND(I),1)
  220 CONTINUE
      DO 230 I=1,NPR
         CEN(I,1)=TEMP(I)
  230 CONTINUE
      DO 240 I=1,NPR
         TEMP(I)=CEN(IND(I),2)
  240 CONTINUE
      DO 250 I=1,NPR
         CEN(I,2)=TEMP(I)
  250 CONTINUE
      DO 260 I=1,NPR
         TEMP(I)=CEN(IND(I),3)
  260 CONTINUE
      DO 270 I=1,NPR
         CEN(I,3)=TEMP(I)
  270 CONTINUE
      DO 280 I=1,NPR
         ITEMP(I)=IAM(IND(I),1)
  280 CONTINUE
      DO 290 I=1,NPR
         IAM(I,1)=ITEMP(I)
  290 CONTINUE
      DO 300 I=1,NPR
         ITEMP(I)=IAM(IND(I),2)
  300 CONTINUE
      DO 310 I=1,NPR
         IAM(I,2)=ITEMP(I)
  310 CONTINUE
C     CALCULATE OVERLAP MATRIX OF STO-6G FUNCTIONS
C
      DO 320 J=1,NC
         CALL OVLP(J,1,IS,IP,NPR,NC,ICD)
  320 CONTINUE
C
      DO 330 J=1,NC
         DO 330 K=1,NC
            CESPM2(INDC(J),INDC(K))=OVL(J,K)
  330 CONTINUE
      DO 340 J=1,NC
         DO 340 K=1,NC
            OVL(J,K)=CESPM2(J,K)
  340 CONTINUE
      L=0
      DO 350 I=1,NC
         DO 350 J=1,I
            L=L+1
            CESP(L)=OVL(I,J)
  350 CONTINUE
C
C     DEORTHOGONALIZE THE COEFFICIENTS AND REFORM THE DENSITY MATRIX
C
      CALL RSP(CESP,NC,1,TEMP,CESPML)
      DO 360 I=1,NC
         DO 360 J=1,I
            SUM=0.D0
            DO 360 K=1,NC
               SUM=SUM+CESPML(I+(K-1)*NC)/SQRT(TEMP(K))*CESPML(J+(K-1)*N
     1C)
               CESP(I+(J-1)*NC)=SUM
               CESP(J+(I-1)*NC)=SUM
  360 CONTINUE
      CALL MULT(C,CESP,CESPML,NC)
      CALL DENSIT(CESPML,NC,NC,NCLOSE,NOPEN,FRACT,CESP,2)
C
C     NOW CALCULATE THE ELECTRONIC CONTRIBUTION TO THE ELECTROSTATIC POT
C
      L=0
      DO 370 I=1,NC
         DO 370 J=1,I
            L=L+1
            CESPM(I,J)=CESP(L)
            CESPM(J,I)=CESP(L)
  370 CONTINUE
      IPX=(NPR-IS)/3
      IPE=IS+IPX
      DO 380 I=1,NESP
         ES(I)=0.D0
  380 CONTINUE
      CALL NAICAS(ISC,IS,IP,NPR,NC,IPE,IPX,ICD)
      CALL NAICAP(ISC,IS,IP,NPR,NC,IPE,IPX,ICD)
C     CALCULATE TOTAL ESP AND FORM ARRAYS FOR ESPFIT
      DO 400 I=1,NESP
         ESP(I)=0.D0
         DO 390 J=1,NATOM
            RA=SQRT((CO(1,J)-POTPT(1,I))**2+(CO(2,J)-POTPT(2,I))**2+(CO(
     13,J)-POTPT(3,I))**2)
            ESP(I)=ESP(I)+TORE(IAN(J))/(RA/BOHR)
  390    CONTINUE
         ESP(I)=ESP(I)-ES(I)
         DO 400  J=1,NATOM
            RIJ=SQRT((CO(1,J)-POTPT(1,I))**2+(CO(2,J)-POTPT(2,I))**2
     1+(CO(3,J)-POTPT(3,I))**2)/BOHR
            B(J)=B(J)+ESP(I)*1.D0/RIJ
  400 CONTINUE
C
C     IF REQUESTED WRITE OUT ELECTRIC POTENTIAL DATA TO
C     UNIT 21
C
      POTWRT=(INDEX(KEYWRD,'POTWRT') .NE. 0)
      IF(POTWRT) THEN
         OPEN(UNIT=21)
         WRITE(21,'(I5)') NESP
         DO 410 I=1,NESP
  410    WRITE(21,420) ESP(I),POTPT(1,I)/BOHR,POTPT(2,I)/BOHR,
     1POTPT(3,I)
      ENDIF
  420 FORMAT(1X,4E16.7)
      RETURN
      END
      DOUBLE PRECISION FUNCTION DEX2(M)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      IF(M .LT. 2) THEN
         DEX2=1
      ELSE
         DEX2=1
         DO 10 I=1,M,2
   10    DEX2=DEX2*I
      ENDIF
      RETURN
      END
      BLOCK DATA ESPBLO
      IMPLICIT DOUBLE PRECISION (A-H, O-Z)
      INCLUDE 'SIZES'
      COMMON /ESPC/  CC(MAXPR),CEN(MAXPR,3),IAM(MAXPR,2),IND(MAXPR),
     1                EX(MAXPR),ESPI(MAXORB,MAXORB),
     2                FV(0:8,821),FAC(0:7),
     3                DEX(-1:96),TF(0:2),TEMP(MAXPR),ITEMP(MAXPR),
     4                OVL(MAXORB,MAXORB),FC(MAXPR*6)
      DATA TF/33.D0,37.D0,41.D0/
      DATA FAC/1.D0,1.D0,2.D0,6.D0,24.D0,120.D0,720.D0,5040.D0/
      END
C***********************************************************************
      SUBROUTINE ESPFIT(IDIP,NATOM,NESP,IZ,ESP,POTPT,CO,
     1DX,DY,DZ,RMS,RRMS)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      INCLUDE 'SIZES'
C***********************************************************************
C
C     THIS ROUTINE FITS THE ELECTROSTATIC POTENTIAL TO A MONOPOLE
C     EXPANSION. FITTING TO THE DIPOLE MONENT CAN ALSO BE DONE.
C     THIS ROUTINE WAS WRITTEN BY B.H.BESLER AND K.M.MERZ
C     IN FEB. 1989 AT UCSF.
C
C     ON INPUT:  IDIP = FLAG TO INDICATE IF THE DIPOLE IS FIT
C                NATOM = NUMBER OF ATOMS
C                NESP = NUMBER OF ESP POINTS
C                IZ = MOLECULAR CHARGE
C                ESP = TOTAL ESP AT THE POINTS
C                POTPT = ESP POINTS
C                CO = COORDINATES
C                DX = X COMPONENT OF THE DIPOLE
C                DY = Y COMPONENT OF THE DIPOLE
C                DZ = Z COMPONENT OF THE DIPOLE
C
C     ON OUTPUT: Q = ESP CHARGES
C                RMS = ROOT MEAN SQUARE FIT
C                RRMS = RELATIVE ROOT MEAN SQUARE FIT
C
C     FOR MORE DETAILS SEE: BESLER,MERZ,KOLLMAN J. COMPUT. CHEM.
C     (IN PRESS)
C***********************************************************************
      COMMON/ESPF/ AL((NUMATM+4)**2),A(NUMATM,NUMATM),B(NUMATM),
     1Q(NUMATM+4),QSC(NUMATM+4),CF, ESPFD(MAXORB**2-NUMATM-5)
      DIMENSION CO(3,*),ESP(*),POTPT(3,*)
      BOHR = 0.529167D00
C     CONVERSION FACTOR FOR DEBYE TO ATOMIC UNITS
      CF=5.2917715D-11*1.601917D-19/3.33564D-30
C
C     THE FOLLOWING SETS UP THE LINEAR EQUATION A*Q=B
C     SET UP THE A(J,K) ARRAY
C
      DO 20  K=1,NATOM
         DO 10  J=1,NATOM
            DO 10  I=1,NESP
               RIK=SQRT((CO(1,K)-POTPT(1,I))**2+(CO(2,K)-POTPT(2,I))**2
     1      +(CO(3,K)-POTPT(3,I))**2)/BOHR
               RIJ=SQRT((CO(1,J)-POTPT(1,I))**2+(CO(2,J)-POTPT(2,I))**2
     1      +(CO(3,J)-POTPT(3,I))**2)/BOHR
               A(J,K)=A(J,K)+1.D0/RIK*1.D0/RIJ
   10    CONTINUE
         A(NATOM+1,K)=1.D0
         A(K,NATOM+1)=1.D0
         A(NATOM+1,NATOM+1)=0.D0
         IF(IDIP .EQ. 1) THEN
            A(NATOM+2,K)=CO(1,K)/BOHR
            A(K,NATOM+2)=CO(1,K)/BOHR
            A(NATOM+2,NATOM+2)=0.D0
            A(NATOM+3,K)=CO(2,K)/BOHR
            A(K,NATOM+3)=CO(2,K)/BOHR
            A(NATOM+3,NATOM+3)=0.D0
            A(NATOM+4,K)=CO(3,K)/BOHR
            A(K,NATOM+4)=CO(3,K)/BOHR
            A(NATOM+4,NATOM+4)=0.D0
         ENDIF
   20 CONTINUE
      B(NATOM+1)=FLOAT(IZ)
      B(NATOM+2)=DX/CF
      B(NATOM+3)=DY/CF
      B(NATOM+4)=DZ/CF
C
C     INSERT CHARGE AND DIPOLAR (IF DESIRED) CONSTRAINTS
C
      IF(IDIP .EQ. 1) THEN
         L=0
         DO 30 I=1,NATOM+4
            DO 30 J=1,NATOM+4
               L=L+1
   30    AL(L)=A(I,J)
      ELSE
         L=0
         DO 40 I=1,NATOM+1
            DO 40 J=1,NATOM+1
               L=L+1
   40    AL(L)=A(I,J)
      ENDIF
      IF (IDIP .EQ. 1) THEN
         CALL OSINV(AL,NATOM+4,DET)
      ELSE
         CALL OSINV(AL,NATOM+1,DET)
      ENDIF
      IF(IDIP .EQ. 1) THEN
         L=0
         DO 50 I=1,NATOM+4
            DO 50 J=1,NATOM+4
               L=L+1
   50    A(I,J)=AL(L)
      ELSE
         L=0
         DO 60 I=1,NATOM+1
            DO 60 J=1,NATOM+1
               L=L+1
   60    A(I,J)=AL(L)
      ENDIF
C
C     SOLVE FOR THE CHARGES
C
      IF(IDIP .EQ. 1) THEN
         DO 70 I=1,NATOM+4
            DO 70 J=1,NATOM+4
               Q(I)=Q(I)+A(I,J)*B(J)
   70    CONTINUE
      ELSE
         DO 80 I=1,NATOM+1
            DO 80 J=1,NATOM+1
               Q(I)=Q(I)+A(I,J)*B(J)
   80    CONTINUE
      ENDIF
C
C     CALCULATE ROOT MEAN SQUARE FITS AND RELATIVE ROOT MEAN SQUARE FITS
C
      CTOT=0.0
      DO 100 I=1,NESP
         ESPC=0.D0
         DO 90 J=1,NATOM
            RIJ=SQRT((CO(1,J)-POTPT(1,I))**2+(CO(2,J)-POTPT(2,I))**2
     1      +(CO(3,J)-POTPT(3,I))**2)/BOHR
   90    ESPC=ESPC+Q(J)/RIJ
         RMS=RMS+(ESPC-ESP(I))**2
  100 RRMS=RRMS+ESP(I)**2
      RMS=SQRT(RMS/NESP)
      RRMS=RMS/SQRT(RRMS/NESP)
      RMS=RMS*627.51D0
      RETURN
      END
C***********************************************************************
      SUBROUTINE FSUB(N,X,FVAL)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C***********************************************************************
C
C     CALCULATE THE FM(T). KINDLY SUPPLIED BY RUS PITZER AND CLEANED UP
C     BY K.M.MERZ IN FEB. 1989 AT UCSF.
C
C     ON INPUT:  N = INDEX
C                X = EXPONENT
C     ON OUTPUT: FVAL = VALUE OF THE FUNCTION
C
C     FOR MORE DETAILS SEE: OBARA AND SAIKA J. CHEM. PHYS. 1986,84,3963
C***********************************************************************
      DIMENSION FF(21),TERM(200),A(10),RT(10)
      DATA A0, A1S2, PIE4, A1
     1   /0.0D0,0.5D0,0.7853981633974483096156608D0,1.0D0/
      DATA XSW /24.0D0/
      E=A1S2*EXP(-X)
      FAC0=N
      FAC0=FAC0+A1S2
      IF(X.GT.XSW) GO TO 50
C
C     USE POWER SERIES
C
   10 FAC=FAC0
      TERM0=E/FAC
      SUM=TERM0
      KU=(X-FAC0)
      IF(KU.LT.1) GO TO 30
C
C     SUM INCREASING TERMS FORWARDS
C
      DO 20 K=1,KU
         FAC=FAC+A1
         TERM0=TERM0*X/FAC
         SUM=SUM+TERM0
   20 CONTINUE
   30 I=1
      FAC=FAC+A1
      TERM(1)=TERM0*X/FAC
      SUMA=SUM+TERM(1)
      IF(SUM.EQ.SUMA) GO TO 90
   40 I=I+1
      FAC=FAC+A1
      TERM(I)=TERM(I-1)*X/FAC
      SUM1=SUMA
      SUMA=SUMA+TERM(I)
      IF(SUM1-SUMA) 40,90,40
C
C     USE ASYMPTOTIC SERIES
C
   50 SUM=SQRT(PIE4/X)
      IF(N.EQ.0) GO TO 70
      FAC=-A1S2
      DO 60 K=1,N
         FAC=FAC+A1
         SUM=SUM*FAC/X
   60 CONTINUE
   70 I=1
      TERM(1)=-E/X
      SUMA=SUM+TERM(1)
      IF(SUM.EQ.SUMA) GO TO 90
      FAC=FAC0
      KU=(X+FAC0-A1)
      DO 80 I=2,KU
         FAC=FAC-A1
         TERM(I)=TERM(I-1)*FAC/X
         SUM1=SUMA
         SUMA=SUMA+TERM(I)
         IF(SUM1.EQ.SUMA) GO TO 90
   80 CONTINUE
C
C     XSW SET TOO LOW. USE POWER SERIES.
C
      GO TO 10
C
C     SUM DECREASING TERMS BACKWARDS
C
   90 SUM1=A0
      DO 100 K=1,I
         SUM1=SUM1+TERM(I+1-K)
  100 CONTINUE
      FF(N+1)=SUM+SUM1
C
C     USE RECURRENCE RELATION
C
      IF(N.EQ.0) GOTO 120
      DO 110 K=1,N
         FAC0=FAC0-A1
         FF(N+1-K)=(E+X*FF(N+2-K))/FAC0
  110 CONTINUE
  120 FVAL=FF(N+1)
      RETURN
      END
      SUBROUTINE SETUP3
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      INCLUDE 'SIZES'
      COMMON /NATYPE/ NZTYPE(107),MTYPE(30),LTYPE
      COMMON /STO6G/ ALLC(6,6,2),ALLZ(6,6,2)
C     SET-UP THE STEWART'S STO-3G EXPANSIONS
C     FROM J. CHEM. PHYS. 52 431.
C                                            1S
      ALLZ(1,1,1) =2.227660584D00
      ALLZ(2,1,1) =4.057711562D-01
      ALLZ(3,1,1) =1.098175104D-01
C
      ALLC(1,1,1) =1.543289673D-01
      ALLC(2,1,1) =5.353281423D-01
      ALLC(3,1,1) =4.446345422D-01
C                                      2S
      ALLZ(1,2,1) =2.581578398D00
      ALLZ(2,2,1) =1.567622104D-01
      ALLZ(3,2,1) =6.018332272D-02
C
      ALLC(1,2,1) =-5.994474934D-02
      ALLC(2,2,1) =5.960385398D-01
      ALLC(3,2,1) =4.581786291D-01
C                                     2P
      ALLZ(1,2,2) =9.192379002D-01
      ALLZ(2,2,2) =2.359194503D-01
      ALLZ(3,2,2) =8.009805746D-02
C
      ALLC(1,2,2) =1.623948553D-01
      ALLC(2,2,2) =5.661708862D-01
      ALLC(3,2,2) =4.223071752D-01
C                                      3S
      ALLZ(1,3,1) =5.641487709D-01
      ALLZ(2,3,1) =6.924421391D-02
      ALLZ(3,3,1) =3.269529097D-02
C
      ALLC(1,3,1) =-1.782577972D-01
      ALLC(2,3,1) =8.612761663D-01
      ALLC(3,3,1) =2.261841969D-01
C                                     3P
      ALLZ(1,3,2) =2.692880368D00
      ALLZ(2,3,2) =1.489359592D-01
      ALLZ(3,3,2) =5.739585040D-02
C
      ALLC(1,3,2) =-1.061945788D-02
      ALLC(2,3,2) =5.218564264D-01
      ALLC(3,3,2) =5.450015143D-01
C                                      4S
      ALLZ(1,4,1) =2.267938753D-01
      ALLZ(2,4,1) =4.448178019D-02
      ALLZ(3,4,1) =2.195294664D-02
C
      ALLC(1,4,1) =-3.349048323D-01
      ALLC(2,4,1) =1.056744667D00
      ALLC(3,4,1) =1.256661680D-01
C                                     4P
      ALLZ(1,4,2) =4.859692220D-01
      ALLZ(2,4,2) =7.430216918D-02
      ALLZ(3,4,2) =3.653340923D-02
C
      ALLC(1,4,2) =-6.147823411D-02
      ALLC(2,4,2) =6.604172234D-01
      ALLC(3,4,2) =3.932639495D-01
C                                      5S
      ALLZ(1,5,1) =1.080198458D-01
      ALLZ(2,5,1) =4.408119382D-02
      ALLZ(3,5,1) =2.610811810D-02
C
      ALLC(1,5,1) =-6.617401158D-01
      ALLC(2,5,1) =7.467595004D-01
      ALLC(3,5,1) =7.146490945D-01
C                                     5P
      ALLZ(1,5,2) =2.127482317D-01
      ALLZ(2,5,2) =4.729648620D-02
      ALLZ(3,5,2) =2.604865324D-02
C
      ALLC(1,5,2) =-1.389529695D-01
      ALLC(2,5,2) =8.076691064D-01
      ALLC(3,5,2) =2.726029342D-01
C
      RETURN
      END
      SUBROUTINE OVLP(IC,IESP,IS,IP,NPR,NC,ICD)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C***********************************************************************
C
C     OVLP CALCULATES THE OVERLAP INTEGRALS FOR A STO-6G BASIS SET.
C     THE RESULTING INTEGRALS ARE USED IN THE DEORTHOGONALIZATION
C     PROCESS.
C     THE CODE WAS WRITTEN BY B.H.BESLER AND K.M.MERZ IN FEB. 1989
C     AT UCSF.
C
C     ON INPUT:  IC = LOOP INDEX
C                IESP = LOOP INDEX
C                IS = NUMBER OF S ORBITALS
C                IP = NUMBER OF P ORBITALS
C                NPR = NUMBER OF PRIMITIVES
C                NC = NUMBER OF CONTRACTED FUNCTIONS
C
C     ON OUTPUT: OVL IS FILLED WITH THE OVERLAP INTEGRAL VALUE
C
C     FOR FURTHER INFO SEE: OBARA & SAIKA J.CHEM.PHYS. 1986,84,3963
C***********************************************************************
      LOGICAL CALLED
      DOUBLE PRECISION NAI,NAI1,NAI2
      INCLUDE 'SIZES'
      COMMON /DENSTY/ P(MPACK),PA(MPACK),PB(MPACK)
      COMMON /POTESP/ XC,YC,ZC,ESPNUC,ESPELE,NESP
      COMMON /ABC/    CO(3,NUMATM),IAN(NUMATM),NATOM
      COMMON /WORK1/  POTPT(3,MESP), ES(MESP), ESP(MESP), WORK1D(2*MESP)
      COMMON /EXPONT/ ZS(107),ZP(107),ZD(107)
      COMMON /STO6G/  ALLC(6,6,2),ALLZ(6,6,2)
      COMMON /ESPC/ CC(MAXPR),CEN(MAXPR,3),IAM(MAXPR,2),IND(MAXPR),
     1EX(MAXPR),ESPI(MAXORB,MAXORB),FV(0:8,821),
     2FAC(0:7),DEX(-1:96),TF(0:2),
     3TEMP(MAXPR),ITEMP(MAXPR),OVL(MAXORB,MAXORB),XDMY(MAXPR*6)
      COMMON/X/ DX(MAXPR),DY(MAXPR),DZ(MAXPR),F1(MAXPR,6),F2(MAXPR,6),
     1TD(MAXPR),CE(MAXPR,6),U(MAXPR,6),EXS(MAXPR,6),EXPN(MAXPR,6),
     2NAI(MAXPR,6),EWCX(MAXPR,6),EWCY(MAXPR,6),EWCZ(MAXPR,6),F0(MAXPR,6)
     3,NAI1(MAXPR,6),NAI2(MAXPR,6)
      DATA BOHR/0.529167D0/
C
C     CALCULATE DISTANCE ARRAYS
C
      PI=4.D0*ATAN(1.D0)
      IPR=IC*ICD-ICD+1
      ISTART=IPR
      DO 10 I=ISTART,NPR
         DX(I)=CEN(IPR,1)-CEN(I,1)
         DY(I)=CEN(IPR,2)-CEN(I,2)
         DZ(I)=CEN(IPR,3)-CEN(I,3)
         TD(I)=DX(I)**2+DY(I)**2+DZ(I)**2
   10 CONTINUE
C
C     CALCULATE EXPONENT SUM
C
      DO 20 I=ISTART,NPR
         DO 20 J=1,ICD
            EXS(I,J)=1.D0/(EX(IPR+J-1)+EX(I))
            CE(I,J)=EX(IPR+J-1)*EX(I)*EXS(I,J)
   20 CONTINUE
C
C     CALCULATE EXPONENT WEIGHTED CENTERS
C
      DO 30 I=ISTART,NPR
         DO 30 J=1,ICD
            EWCX(I,J)=(EX(I)*CEN(I,1)+EX(IPR+J-1)
     1*CEN(IPR+J-1,1))*EXS(I,J)
            EWCY(I,J)=(EX(I)*CEN(I,2)+EX(IPR+J-1)
     1*CEN(IPR+J-1,2))*EXS(I,J)
            EWCZ(I,J)=(EX(I)*CEN(I,3)+EX(IPR+J-1)
     1*CEN(IPR+J-1,3))*EXS(I,J)
   30 CONTINUE
      DO 40 I=1,NPR
         DO 40 J=1,ICD
            EXPN(I,J)=EXP(-CE(I,J)*TD(I))
            NAI(I,J)=(PI*EXS(I,J))**1.5D0*EXPN(I,J)
            EXPN(I,J)=NAI(I,J)
   40 CONTINUE
C
C     CALCULATE (S||P) ESP INTEGRALS
C
      IF((IAM(IPR,1) .EQ. 0) .AND. (IS .NE. IP)) THEN
         NP=IS+1
         DO 80 I=NP,NPR
            DO 80 J=1,ICD
               GO TO (50,60,70),IAM(I,2)
   50          NAI(I,J)=(EWCX(I,J)-CEN(I,1))*EXPN(I,J)
               go TO 80
   60          NAI(I,J)=(EWCY(I,J)-CEN(I,2))*EXPN(I,J)
               GO TO 80
   70          NAI(I,J)=(EWCZ(I,J)-CEN(I,3))*EXPN(I,J)
   80    CONTINUE
      ENDIF
C
C     CALCULATE (P||S) ESP INTEGRALS
C
      IF((IAM(IPR,1) .EQ. 1) .AND. (IS .NE. IP)) THEN
         NP=IS+1
         DO 120 I=	ISTART,NPR
            DO 120 J=1,ICD
               GO TO (90,100,110),IAM(IPR+J-1,2)
   90          NAI(I,J)=(EWCX(I,J)-CEN(IPR+J-1,1))*EXPN(I,J)
               GO TO 120
  100          NAI(I,J)=(EWCY(I,J)-CEN(IPR+J-1,2))*EXPN(I,J)
               GO TO 120
  110          NAI(I,J)=(EWCZ(I,J)-CEN(IPR+J-1,3))*EXPN(I,J)
  120    CONTINUE
      ENDIF
C
C     CALCULATE (P||P) ESP INTEGRALS
C
      IF((IAM(IPR,1) .EQ. 1) .AND. (IS .NE. IP)) THEN
         DO 160 I=ISTART,NPR
            DO 160 J=1,ICD
               GO TO (130,140,150),IAM(I,2)
  130          NAI(I,J)=(EWCX(I,J)-CEN(I,1))*NAI(I,J)
               IF(IAM(IPR+J-1,2) .EQ. IAM(I,2))
     1NAI(I,J)=NAI(I,J)+EXS(I,J)*0.5D0
     2      *EXPN(I,J)
               GO TO 160
  140          NAI(I,J)=(EWCY(I,J)-CEN(I,2))*NAI(I,J)
               IF(IAM(IPR+J-1,2) .EQ. IAM(I,2))
     1NAI(I,J)=NAI(I,J)+EXS(I,J)*0.5D0
     2      *EXPN(I,J)
               GO TO 160
  150          NAI(I,J)=(EWCZ(I,J)-CEN(I,3))*NAI(I,J)
               IF(IAM(IPR+J-1,2) .EQ. IAM(I,2))
     1NAI(I,J)=NAI(I,J)+EXS(I,J)*0.5D0
     2      *EXPN(I,J)
  160    CONTINUE
      ENDIF
      IPS=IC*ICD-ICD+1
      DO 180 I=IC,NC
         JPS=I*ICD-ICD+1
         OVL(IC,I)=0.D0
         DO 170 J=JPS,JPS+ICD-1
            DO 170 K=IPS,IPS+ICD-1
               OVL(IC,I)=OVL(IC,I)+CC(J)*CC(K)*NAI(J,K-IPS+1)
  170    CONTINUE
         OVL(I,IC)=OVL(IC,I)
  180 CONTINUE
      RETURN
      END
      SUBROUTINE NAICAS(ISC,IS,IP,NPR,NC,IPE,IPX,ICD)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C***********************************************************************
C
C     THIS SUBROUTINE EVALUATES (S|S) , (S|P) TYPE NUCLEAR ATTRACTION
C     INTEGRALS FOR A STO-NG BASIS SET
C     WRITTEN BY B.H. BESLER AT FORD SCIENTIFIC RESEARCH LABS IN
C     DECEMBER 1989.
C
C     ON INPUT:  IC = LOOP INDEX OF THE GAUSSIAN
C                IESP = LOOP INDEX OF THE ESP POINT
C                IPE = INDEX OF LAST Px PRIMITIVE
C                IPX = NUMBER OF Px PRIMITIVES
C                IS = NUMBER OS S ORBITALS
C                ISC = NUMBER OF CONTRACTED S ORBITALS
C                IP = NUMBER OF P ORBITALS
C                NPR = NUMBER OF PRIMITIVES
C                NC = NUMBER OF CONTRACTED FUNCTIONS
C
C
C     FOR MORE INFO SEE: OBARA&SAIKA J.CHEM.PHYS. 1986,84,3963.
C***********************************************************************
      INCLUDE 'SIZES'
      DOUBLE PRECISION NAI,NAI1,NAI2
      CHARACTER*241 KEYWRD
      COMMON/KEYWRD/ KEYWRD
      COMMON/ESPF/ AL((NUMATM+4)**2),A(NUMATM,NUMATM),B(NUMATM),
     1Q(NUMATM+4),CESPM(MAXORB,MAXORB)
      COMMON /INDX/ INDC(MAXORB)
      COMMON /DENSTY/ P(MPACK),PA(MPACK),PB(MPACK)
      COMMON /POTESP/ XC,YC,ZC,ESPNUC,ESPELE,NESP
      COMMON /ABC/    CO(3,NUMATM),IAN(NUMATM),NATOM
      COMMON /WORK1/  POTPT(3,MESP), ES(MESP), ESP(MESP), WORK1D(2*MESP)
      COMMON /EXPONT/ ZS(107),ZP(107),ZD(107)
      COMMON /STO6G/  ALLC(6,6,2),ALLZ(6,6,2)
      COMMON /ESPC/ CC(MAXPR),CEN(MAXPR,3),IAM(MAXPR,2),IND(MAXPR),
     1EX(MAXPR),ESPI(MAXORB,MAXORB),FV(0:8,821),
     2FAC(0:7),DEX(-1:96),TF(0:2),
     3TEMP(MAXPR),ITEMP(MAXPR),OVL(MAXORB,MAXORB),EXSR(MAXPR,6)
      COMMON/X/ DX(MAXPR),DY(MAXPR),DZ(MAXPR),F1(MAXPR,6),F2(MAXPR,6),
     1TD(MAXPR),CE(MAXPR,6),U(MAXPR,6),EXS(MAXPR,6),EXPN(MAXPR,6),
     2NAI(MAXPR,6),EWCX(MAXPR,6),EWCY(MAXPR,6),EWCZ(MAXPR,6),F0(MAXPR,6)
     3,NAI1(MAXPR,6),NAI2(MAXPR,6)
      DATA BOHR/0.529167D0/
C
C     CALCULATE DISTANCE ARRAYS
C
C *** it seems that this is not necessary...
C      WRITE(6,*)
      PI=4.D0*ATAN(1.D0)
      IPX2=2*IPX
C     IF THIS IS A RESTART RUN, READ IN RESTART INFO
C *** skip all restart stuff, we don't need that...
C      IF(INDEX(KEYWRD,'ESPRST') .NE. 0) THEN
C         OPEN(UNIT=15,FILE='ESP.DUMP',STATUS='OLD',FORM='UNFORMATTED')
C         READ(15) JSTART,IESPS
C         IF(JSTART .EQ. ISC*2) THEN
C            CLOSE(15)
C            RETURN
C         ENDIF
C         DO 10 I=1,NESP
C            READ(15) ES(I)
C   10    CONTINUE
C         CLOSE(15)
CC
C         JSTART=JSTART+1
C      ELSE
         JSTART=1
C      ENDIF
      NP=IS+1
      DO 200 IC=JSTART,ISC
         IPR=IC*ICD-ICD+1
         ISTART=IPR
         DO 20 I=ISTART,IPE
            DX(I)=CEN(IPR,1)-CEN(I,1)
            DY(I)=CEN(IPR,2)-CEN(I,2)
            DZ(I)=CEN(IPR,3)-CEN(I,3)
            TD(I)=DX(I)**2+DY(I)**2+DZ(I)**2
   20    CONTINUE
C
C     CALCULATE EXPONENT SUM
C
         DO 30 I=ISTART,IPE
            DO 30 J=1,ICD
               EXSR(I,J)=EX(IPR+J-1)+EX(I)
               EXS(I,J)=1.D0/EXSR(I,J)
               CE(I,J)=EX(IPR+J-1)*EX(I)*EXS(I,J)
               EXPN(I,J)=EXP(-CE(I,J)*TD(I))
   30    CONTINUE
C
C     CALCULATE EXPONENT WEIGHTED CENTERS
C
         DO 40 I=ISTART,IPE
            DO 40 J=1,ICD
               EWCX(I,J)=(EX(I)*CEN(I,1)+EX(IPR+J-1)
     1*CEN(IPR+J-1,1))*EXS(I,J)
               EWCY(I,J)=(EX(I)*CEN(I,2)+EX(IPR+J-1)
     1*CEN(IPR+J-1,2))*EXS(I,J)
               EWCZ(I,J)=(EX(I)*CEN(I,3)+EX(IPR+J-1)
     1*CEN(IPR+J-1,3))*EXS(I,J)
   40    CONTINUE
C
C     BEGIN LOOP OVER ESP POINTS
C
         DO 180 IESP=1,NESP
            POTP1=POTPT(1,IESP)/BOHR
            POTP2=POTPT(2,IESP)/BOHR
            POTP3=POTPT(3,IESP)/BOHR
C
C     BEGIN LOOP OVER COMPONENTS OF CONTRACTED FUNCTION IC
C
            DO 150 J=1,ICD
C
C     CALCULATE DISTANCE BETWEEN EXPONENT WEIGHTED AND PROBE POINT
C
               DO 50 I=ISTART,IPE
                  U(I,J)=((EWCX(I,J)-POTP1)**2+(EWCY(I,J)-POTP2)**2+
     1      (EWCZ(I,J)-POTP3)**2)*EXSR(I,J)
                  NAI(I,J)=SQRT(PI/U(I,J))
   50          CONTINUE
C
C     CALCULATE ESP INTEGRALS
C
               DO 70 I=ISTART,IPE
                  IF(U(I,J) .LE. TF(0)) THEN
                     IREF=DNINT(U(I,J)*20.D0)
                     REF=0.05D0*IREF
                     RES=U(I,J)-REF
                     TERM=1.D0
                     F0(I,J)=0.D0
                     DO 60 K=0,6
                        F=FV(K,IREF+1)
                        TS=F*TERM*FAC(K)
                        TERM=-TERM*RES
                        F0(I,J)=F0(I,J)+TS
   60                CONTINUE
                  ELSE
                     F0(I,J)=NAI(I,J)*0.5D0
                  ENDIF
   70          CONTINUE
               DO 90 I=NP,IPE
                  IF(U(I,J) .LE. TF(1)) THEN
                     IREF=DNINT(U(I,J)*20.D0)
                     REF=0.05D0*IREF
                     RES=U(I,J)-REF
                     TERM1=1.D0
                     F1(I,J)=0.D0
                     DO 80 K=0,6
                        FI=FV(K+1,IREF+1)
                        TS1=FI*TERM1*FAC(K)
                        TERM1=-TERM1*RES
                        F1(I,J)=F1(I,J)+TS1
   80                CONTINUE
                  ELSE
                     F1(I,J)=NAI(I,J)*0.25D0/U(I,J)
                  ENDIF
   90          CONTINUE
               DO 100 I=ISTART,IS
  100          U(I,J)=2.D0*PI*EXS(I,J)*EXPN(I,J)*F0(I,J)
               NP=IS+1
               DO 110 I=NP,IPE
                  NAI(I,J)=2.D0*PI*EXS(I,J)*EXPN(I,J)*F0(I,J)
                  NAI1(I,J)=2.D0*PI*EXS(I,J)*EXPN(I,J)*F1(I,J)
  110          CONTINUE
C
C     CALCULATE (S||P) ESP INTEGRALS
C
               IF((IAM(IPR,1) .EQ. 0) .AND. (IS .NE. IP)) THEN
                  DO 120 I=NP,IPE
  120             U(I,J)=(EWCX(I,J)-CEN(I,1))*NAI(I,J)
     1-(EWCX(I,J)-POTP1)*NAI1(I,J)
                  DO 130 I=IPE+1,IPE+1+IPX
  130             U(I,J)=(EWCY(I-IPX,J)-CEN(I-IPX,2))*NAI(I-IPX,J)
     1-(EWCY(I-IPX,J)-POTP2)*NAI1(I-IPX,J)
                  DO 140 I=IPE+1+IPX,NPR
  140             U(I,J)=(EWCZ(I-IPX2,J)-CEN(I-IPX2,3))*NAI(I-IPX2,J)
     1-(EWCZ(I-IPX2,J)-POTP3)*NAI1(I-IPX2,J)
               ENDIF
  150       CONTINUE
            IPS=IC*ICD-ICD+1
            DO 170 I=IC,NC
               JPS=I*ICD-ICD+1
               ESPI(I,IC)=0.D0
               DO 160 J=JPS,JPS+ICD-1
                  DO 160 K=IPS,IPS+ICD-1
                     ESPI(I,IC)=ESPI(I,IC)+CC(J)*CC(K)*U(J,K-IPS+1)
  160          CONTINUE
               ES(IESP)=ES(IESP)+2.D0*CESPM(INDC(I),INDC(IC))*ESPI(I,IC)
  170       CONTINUE
            ES(IESP)=ES(IESP)-CESPM(INDC(IC),INDC(IC))*ESPI(IC,IC)
  180    CONTINUE
C     WRITE OUT RESTART INFORMATION
C *** no dumps please...
C *** no dumps please...
C *** no dumps please...
C         OPEN(UNIT=15,FILE='ESP.DUMP',STATUS='UNKNOWN',FORM='UNFORMATTED
C     1')
C         IESPS=0
C         WRITE(15) IC,IESPS
C         DO 190 I=1,NESP
C            WRITE(15) ES(I)
C  190    CONTINUE
C         CLOSE(15)
C
C         WRITE(6,'(A,F6.2,A)')
C     1'NAICAS DUMPED: ',100.D0/ISC*IC,' PERCENT COMPLETE'
  200 CONTINUE
      RETURN
      END
      SUBROUTINE NAICAP(ISC,IS,IP,NPR,NC,IPE,IPX,ICD)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C***********************************************************************
C     THIS ROUTINE EVALUATES (P|P) NUCLEAR ATTRACTION INTEGRALS OVER
C
C     A STO-NG BASIS SET.
C     WRITTEN BY B.H. BESLER AT FORD SCIENTIFIC RESEARCH LABS IN
C     SEPT. 1989
C
C     ON INPUT:  IC = LOOP INDEX OF THE GAUSSIAN
C                ICD = CONTRACTION DEPTH OF BASIS SET
C                IESP = LOOP INDEX OF THE ESP POINT
C                IS = NUMBER OS S PRIMITIVES
C                IPE = INDEX OF LAST PX PRIMITIVE
C                IPX = NUMBER OF PX PRIMITIVES
C                IS = NUMBER OS S PRIMITIVES
C                ISC = NUMBER OF CONTRACTED
C                NPR = NUMBER OF PRIMITIVES
C                NC = NUMBER OF CONTRACTED FUNCTIONS
C
C
C     FOR MORE INFO SEE: OBARA&SAIKA J.CHEM.PHYS. 1986,84,3963.
C***********************************************************************
      INCLUDE 'SIZES'
      DOUBLE PRECISION NAI,NAI1,NAI2
      CHARACTER*241 KEYWRD
      COMMON /KEYWRD/ KEYWRD
      COMMON/ESPF/ AL((NUMATM+4)**2),A(NUMATM,NUMATM),B(NUMATM),
     1Q(NUMATM+4),CESPM(MAXORB,MAXORB)
      COMMON /INDX/ INDC(MAXORB)
      COMMON /DENSTY/ P(MPACK),PA(MPACK),PB(MPACK)
      COMMON /POTESP/ XC,YC,ZC,ESPNUC,ESPELE,NESP
      COMMON /ABC/    CO(3,NUMATM),IAN(NUMATM),NATOM
      COMMON /WORK1/  POTPT(3,MESP), ES(MESP), ESP(MESP), WORK1D(2*MESP)
      COMMON /EXPONT/ ZS(107),ZP(107),ZD(107)
      COMMON /STO6G/  ALLC(6,6,2),ALLZ(6,6,2)
      COMMON /ESPC/ CC(MAXPR),CEN(MAXPR,3),IAM(MAXPR,2),IND(MAXPR),
     1EX(MAXPR),ESPI(MAXORB,MAXORB),FV(0:8,821),
     2FAC(0:7),DEX(-1:96),TF(0:2),
     3TEMP(MAXPR),ITEMP(MAXPR),OVL(MAXORB,MAXORB),EXSR(MAXPR,6)
      COMMON/X/ DX(MAXPR),DY(MAXPR),DZ(MAXPR),F1(MAXPR,6),F2(MAXPR,6),
     1TD(MAXPR),CE(MAXPR,6),U(MAXPR,6),EXS(MAXPR,6),EXPN(MAXPR,6),
     2NAI(MAXPR,6),EWCX(MAXPR,6),EWCY(MAXPR,6),EWCZ(MAXPR,6),F0(MAXPR,6)
     3,NAI1(MAXPR,6),NAI2(MAXPR,6)
      COMMON/FP/ PF0(MAXHES),PF1(MAXHES),PF2(MAXHES),ID(MAXPAR),
     1PEXS(MAXHES),PCE(MAXHES),PEXPN(MAXHES),PTD(MAXHES),
     2PEWCX(MAXHES),PEWCY(MAXHES),PEWCZ(MAXHES),IRD(MAXHES)
      DATA BOHR/0.529167D0/
C     SET NUMBER OF EQUALLY SPACED DUMPS
      IDN=10
C
      IDC=0
C *** it seems that this is not necessary...
C      WRITE(6,*)
      IPX2=2*IPX
      PI=4.D0*ATAN(1.D0)
      NP=IS+1
C     SETUP INDEX ARRAY
      DO 10 I=NP,IPE
         IRD(I)=I-IS
         IRD(I+IPX)=I-IS
         IRD(I+IPX2)=I-IS
   10 CONTINUE
C
C     CALCULATE QUANTITIES INVARIANT WITH ESP POINT FOR
C     (P|P) ESP INTEGRALS
C
      IL=L
      L=0
      DO 30 I=NP,IPE
         DO 20 J=I,IPE
            L=L+1
            PTD(L)=(CEN(I,1)-CEN(J,1))**2+(CEN(I,2)-CEN(J,2))**2+
     1(CEN(I,3)-CEN(J,3))**2
            PEXS(L)=1.d0/(EX(I)+EX(J))
            PCE(L)=EX(I)*EX(J)*PEXS(L)
            PEXPN(L)=EXP(-PCE(L)*PTD(L))
            PEWCX(L)=(EX(I)*CEN(I,1)+EX(J)*CEN(J,1))*PEXS(L)
            PEWCY(L)=(EX(I)*CEN(I,2)+EX(J)*CEN(J,2))*PEXS(L)
            PEWCZ(L)=(EX(I)*CEN(I,3)+EX(J)*CEN(J,3))*PEXS(L)
   20    CONTINUE
C
C     SET UP OTHER INDEX ARRAY FOR PACKED SYMMETRIC ARRAY
C     STORAGE
C
         ID(I-IS)=L-IPX
   30 CONTINUE
C
C     READ IN RESTART INFORMATION IF THIS IS A RESTART
C
C *** skip all restart stuff, we don't need that...
C      IF(INDEX(KEYWRD,'ESPRST') .NE. 0) THEN
C         OPEN(UNIT=15,FILE='ESP.DUMP',STATUS='UNKNOWN',FORM='UNFORMATTED
C     1')
C         READ(15) JSTART,IESPS
C         IF(JSTART .NE. ISC*2) THEN
C            IESPS=0
C            CLOSE(15)
C            GOTO 50
C         ENDIF
C         DO 40 I=1,NESP
C            READ(15) ES(I)
C   40    CONTINUE
C         CLOSE(15)
C         IDC=FLOAT(IESPS)/FLOAT(NESP)*10
C      ELSE
         IESPS=0
C      ENDIF
   50 CONTINUE
C
C     LOOP OVER ESP PROBE POINTS
C
      DO 250 IESP=IESPS+1,NESP
         POTP1=POTPT(1,IESP)/BOHR
         POTP2=POTPT(2,IESP)/BOHR
         POTP3=POTPT(3,IESP)/BOHR
C     CALCULATE QUANTITY U
C
         L=0
         DO 60 I=NP,IPE
            DO 60 J=I,IPE
               L=L+1
               PTD(L)=((PEWCX(L)-POTP1)**2+(PEWCY(L)-POTP2)**2+
     1      (PEWCZ(L)-POTP3)**2)/PEXS(L)
               PCE(L)=SQRT(PI/PTD(L))
   60    CONTINUE
C
C     CALCULATE F0, F1, AND F2(U) USING TAYLOR SERIES
C     OR ASYMPTOTIC EXPANSION
C
         IL=L
         L=0
         DO 100 I=1,IL
            IF(PTD(I) .LE. TF(0)) THEN
               IREF=DNINT(PTD(I)*20.D0)
               REF=0.05D0*IREF
               RES=PTD(I)-REF
               TERM=1.D0
               PF0(I)=0.D0
               DO 70 K=0,6
                  F=FV(K,IREF+1)
                  TS=F*TERM*FAC(K)
                  TERM=-TERM*RES
                  PF0(I)=PF0(I)+TS
   70          CONTINUE
            ELSE
               PF0(I)=PCE(I)*0.5D0
            ENDIF
            IF(PTD(I) .LE. TF(1)) THEN
               IREF=DNINT(PTD(I)*20.D0)
               REF=0.05D0*IREF
               RES=PTD(I)-REF
               TERM1=1.D0
               PF1(I)=0.D0
               DO 80 K=0,6
                  FI=FV(K+1,IREF+1)
                  TS1=FI*TERM1*FAC(K)
                  TERM1=-TERM1*RES
                  PF1(I)=PF1(I)+TS1
   80          CONTINUE
            ELSE
               PF1(I)=PCE(I)*0.25D0/PTD(I)
            ENDIF
            IF(PTD(I) .LE. TF(2)) THEN
               IREF=DNINT(PTD(I)*20.D0)
               REF=0.05D0*IREF
               RES=PTD(I)-REF
               TERM2=1.D0
               PF2(I)=0.D0
               DO 90 K=0,6
                  FII=FV(K+2,IREF+1)
                  TS2=FII*TERM2*FAC(K)
                  TERM2=-TERM2*RES
                  PF2(I)=PF2(I)+TS2
   90          CONTINUE
            ELSE
               PF2(I)=PCE(I)*0.375D0/(PTD(I)*PTD(I))
            ENDIF
  100    CONTINUE
C
C     CALCULATE (S||S) TYPE INTEGRALS
C
         DO 110 I=1,IL
            PF0(I)=2.D0*PI*PEXS(I)*PEXPN(I)*PF0(I)
            PTD(I)=PF0(I)
            PF1(I)=2.D0*PI*PEXS(I)*PEXPN(I)*PF1(I)
            PF2(I)=2.D0*PI*PEXS(I)*PEXPN(I)*PF2(I)
  110    CONTINUE
C
         DO 230 IC=ISC+1,NC
            IPR=IC*ICD-ICD+1
            ISTART=IPR
            DO 200 J=1,ICD
C
C     CALCULATE (P||S) ESP INTEGRALS
C
               IF((IAM(IPR,1) .EQ. 1) .AND. (IS .NE. IP)) THEN
                  DO 150 I=ISTART,NPR
                     IN=IPR+J-1
                     IR=IRD(I)+ID(IRD(IN))
                     IR2=ID(IRD(I))+IRD(IN)
                     IF(IR2 .LE. IR ) IR=IR2
                     GO TO (120,130,140),IAM(IN,2)
  120                NAI2(I,J)=(PEWCX(IR)-CEN(IN,1))*PF1(IR)-PF2(IR)*
     1      (PEWCX(IR)-POTP1)
                     NAI(I,J)=(PEWCX(IR)-CEN(IN,1))*PF0(IR)-PF1(IR)*
     1      (PEWCX(IR)-POTP1)
                     GO TO 150
  130                NAI2(I,J)=(PEWCY(IR)-CEN(IN,2))*PF1(IR)-PF2(IR)*
     1      (PEWCY(IR)-POTP2)
                     NAI(I,J)=(PEWCY(IR)-CEN(IN,2))*PF0(IR)-PF1(IR)*
     1      (PEWCY(IR)-POTP2)
                     GO TO 150
  140                NAI2(I,J)=(PEWCZ(IR)-CEN(IN,3))*PF1(IR)-PF2(IR)*
     1      (PEWCZ(IR)-POTP3)
                     NAI(I,J)=(PEWCZ(IR)-CEN(IN,3))*PF0(IR)-PF1(IR)*
     1      (PEWCZ(IR)-POTP3)
  150             CONTINUE
               ENDIF
C
C     CALCULATE (P||P) ESP INTEGRALS
C
               IF((IAM(IPR,1) .EQ. 1) .AND. (IS .NE. IP)) THEN
                  DO 190 I=ISTART,NPR
                     IN=IPR+J-1
                     IR=IRD(I)+ID(IRD(IN))
                     IR2=ID(IRD(I))+IRD(IN)
                     IF(IR2 .LE. IR ) IR=IR2
                     GO TO (160,170,180),IAM(I,2)
  160                NAI(I,J)=(PEWCX(IR)-CEN(I,1))*NAI(I,J)-(PEWCX(IR)-P
     1OTP1)*      NAI2(I,J)
                     IF(IAM(IN,2) .EQ. IAM(I,2)) NAI(I,J)=NAI(I,J)+PEXS(
     1IR)*      0.5D0*(PTD(IR)-PF1(IR))
                     GO TO 190
  170                NAI(I,J)=(PEWCY(IR)-CEN(I,2))*NAI(I,J)-(PEWCY(IR)-P
     1OTP2)*      NAI2(I,J)
                     IF(IAM(IN,2) .EQ. IAM(I,2)) NAI(I,J)=NAI(I,J)+PEXS(
     1IR)*      0.5D0*(PTD(IR)-PF1(IR))
                     GO TO 190
  180                NAI(I,J)=(PEWCZ(IR)-CEN(I,3))*NAI(I,J)-(PEWCZ(IR)-P
     1OTP3)*      NAI2(I,J)
                     IF(IAM(IN,2) .EQ. IAM(I,2)) NAI(I,J)=NAI(I,J)+PEXS(
     1IR)*      0.5D0*(PTD(IR)-PF1(IR))
  190             CONTINUE
               ENDIF
  200       CONTINUE
C
C     FORM INTEGRALS OVER CONTRACTED FUNCTIONS
C
            IPS=IC*ICD-ICD+1
            DO 220 I=IC,NC
               JPS=I*ICD-ICD+1
               ESPI(I,IC)=0.D0
               DO 210 J=JPS,JPS+ICD-1
                  DO 210 K=IPS,IPS+ICD-1
                     ESPI(I,IC)=ESPI(I,IC)+CC(J)*CC(K)*NAI(J,K-IPS+1)
  210          CONTINUE
               ES(IESP)=ES(IESP)+2.D0*CESPM(INDC(I),INDC(IC))*ESPI(I,IC)
  220       CONTINUE
            ES(IESP)=ES(IESP)-CESPM(INDC(IC),INDC(IC))*ESPI(IC,IC)
  230    CONTINUE
C
C     WRITE OUT RESTART INFORMATION EVERY NESP/10 POINTS
C
C *** no dumps please...
C *** no dumps please...
C *** no dumps please...
C         IF(MOD(IESP,NESP/IDN) .EQ. 0) THEN
C            OPEN(UNIT=15,FILE='ESP.DUMP',STATUS='UNKNOWN',FORM='UNFORMAT
C     1TED')
C            JSTART=ISC*2
C            WRITE(15) JSTART,IESP
C            DO 240 I=1,NESP
C               WRITE(15) ES(I)
C  240       CONTINUE
C            CLOSE(15)
C            IDC=IDC+1
C            WRITE(6,'(A,F6.2,A)')
C     1'NAICAP DUMPED: ',100.D0/IDN*IDC,' PERCENT COMPLETE'
C         ENDIF
  250 CONTINUE
      RETURN
      END
C *** extensions for "miniMOPAC" plotting start here...
C *** extensions for "miniMOPAC" plotting start here...
C *** extensions for "miniMOPAC" plotting start here...
      SUBROUTINE GETGEOM
C *** this is a start of PDGRID subroutine with small modifications.
C *** this will just copy the geometry data for orginal ELESP.
C
C     ROUTINE TO CALCULATE WILLIAMS SURFACE
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      INCLUDE 'SIZES'
      DIMENSION IZ(100),XYZ(3,100),VDERW(53),DIST(100)
      DIMENSION XMIN(3),XMAX(3),COORD(3,NUMATM)
      COMMON /GEOM/   GEO(3,NUMATM)
      COMMON /GEOKST/ NATOMS,LABELS(NUMATM), NABC(3*NUMATM)
C
      COMMON /ABC/    CO(3,NUMATM),IAN(NUMATM),NATOM
      COMMON /WORK1/    POTPT(3,MESP), WORK1D(4*MESP)
      COMMON /POTESP/ XC,YC,ZC,ESPNUC,ESPELE,NESP
C
      DATA VDERW/53*0.0D0/
      VDERW(1)=2.4D0
      VDERW(5)=3.0D0
      VDERW(6)=2.9D0
      VDERW(7)=2.7D0
      VDERW(8)=2.6D0
      VDERW(9)=2.55D0
      VDERW(15)=3.1D0
      VDERW(16)=3.05D0
      VDERW(17)=3.0D0
      VDERW(35)=3.15D0
      VDERW(53)=3.35D0
      SHELL=1.2D0
C     NESP=0
      GRID=0.8D0
      CLOSER=0.D0
C     CHECK IF VDERW IS DEFINED FOR ALL ATOMS
C
C     CONVERT INTERNAL TO CARTESIAN COORDINATES
C
      CALL GMETRY(GEO,COORD)
C
C     STRIP COORDINATES AND ATOM LABEL FOR DUMMIES (I.E. 99)
C
      ICNTR = 0
      DO 20 I=1,NATOMS
         DO 10 J=1,3
   10    CO(J,I) = COORD(J,I)
         IF(LABELS(I) .EQ. 99) GOTO 20
         ICNTR = ICNTR + 1
         IAN(ICNTR) = LABELS(I)
   20 CONTINUE
      NATOM=ICNTR
      RETURN
      END
      SUBROUTINE LM7INIPLT
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C *** this is a modification to ELESP. it initializes the ELESP
C *** calculation and also stores some extra data for other plots.
C***********************************************************************
C     ELESP LOADS THE STO-6G BASIS SET ONTO THE ATOMS, PERFOMS THE
C     DEORTHOGONALIZATION OF THE COEFFICIENTS AND EVALUATES THE
C     ELECTRONIC CONTRIBUTION TO THE ESP. IT WAS WRITTEN BY B.H.BESLER
C     AND K.M.MERZ IN FEB. 1989 AT UCSF.
C
C***********************************************************************
      CHARACTER*241 KEYWRD
      DOUBLE PRECISION NORM,OVL
      LOGICAL CALLED,POTWRT,RST,STO3G
      INCLUDE 'SIZES'
      COMMON/ESPF/ AL((NUMATM+4)**2),A(NUMATM,NUMATM),B(NUMATM),
     1Q(NUMATM+4),CESPM(MAXORB,MAXORB)
      COMMON /DENSTY/ P(MPACK),PA(MPACK),PB(MPACK)
      COMMON /POTESP/ XC,YC,ZC,ESPNUC,ESPELE,NESP
      COMMON /ABC/    CO(3,NUMATM),IAN(NUMATM),NATOM
      COMMON /WORK1/  POTPT(3,MESP), ES(MESP), ESP(MESP), WORK1D(2*MESP)
      COMMON /STO6G/  ALLC(6,6,2),ALLZ(6,6,2)
      COMMON /VECTOR/ C(MORB2*2+MAXORB*2)
      COMMON /MOLKST/ NUMAT,NAT(NUMATM),NFIRST(NUMATM),NMIDLE(NUMATM),
     1                NLAST(NUMATM), NORBS, NELECS,NALPHA,NBETA,
     2                NCLOSE,NOPEN,NDUMY,FRACT
      COMMON /KEYWRD/ KEYWRD
      COMMON /ESPC/  CC(MAXPR),CEN(MAXPR,3),IAM(MAXPR,2),IND(MAXPR),
     1                EX(MAXPR),ESPI(MAXORB,MAXORB),
     2                FV(0:8,821),FAC(0:7),
     3                DEX(-1:96),TF(0:2),TEMP(MAXPR),ITEMP(MAXPR),
     4                OVL(MAXORB,MAXORB),FC(MAXPR*6)
     6       /CORE  / TORE(107)
     7       /EXPONT/ ZS(107),ZP(107),ZD(107)
*
*  END OF MINDO/3 COMMON BLOCKS
*
      COMMON /INDX/   INDC(MAXORB)
C *** an additional common block that carries variables for plotting routines.
      COMMON /PLOTS/  CESPM2(MAXORB,MAXORB),SLA(10),
     1                CESPML(MAXORB*MAXORB),CESP(MAXORB*MAXORB),
     2                INC(MAXPR),NC,NPR,IS,IP,IPC,ISC,ICD,IORB
C *** old arrays that are no longer needed are here, commented out.
C     DIMENSION CESPM2(MAXORB,MAXORB),SLA(10)
C     DIMENSION CESPML(MAXORB*MAXORB),CESP(MAXORB*MAXORB)
      DATA BOHR/0.529167D0/
C *** now we call our GETGEOM subroutine here...
      CALL GETGEOM
      PI=4.D0*ATAN(1.D0)
C
C     PUT STO-6G BASIS SET ON ATOM CENTERS
C
      DO 10 I=-1,10
         DEX(I)=DEX2(I)
   10 CONTINUE
      DO 20   I=0,7
         FAC(I)=1.D0/FAC(I)
   20 CONTINUE
      DO 30 M=0,8
         K=1
         FV(M,1)=1.D0/(2.D0*M+1.D0)
         DO 30 T=0.05D0,41.D0,0.05D0
            K=K+1
            CALL FSUB(M,T,FVAL)
            FV(M,K)=FVAL
   30 CONTINUE
C
C     LOAD BASIS FUNCTIONS INTO ARRAYS
C
      STO3G=(INDEX(KEYWRD,'STO3G') .NE. 0)
      IF(STO3G) THEN
         ICD=3
         CALL SETUP3
      ELSE
         ICD=6
         CALL SETUPG
      ENDIF
C *** NC is number of contractions
C *** NPR is number of primitives
      NC=0
      NPR=0
C *** the new array INC() will store the contraction indices...
C *** the new array INC() will store the contraction indices...
C *** the new array INC() will store the contraction indices...
      DO 80 I=1,NATOM
         IF (IAN(I) .LE. 2) THEN
            NC=NC+1
            DO 40 J=1,ICD
               CC(NPR+J)=ALLC(J,1,1)
               EX(NPR+J)=ALLZ(J,1,1)*ZS(1)**2
               CEN(NPR+J,1)=CO(1,I)/BOHR
               CEN(NPR+J,2)=CO(2,I)/BOHR
               CEN(NPR+J,3)=CO(3,I)/BOHR
               IAM(NPR+J,1)=0
               IAM(NPR+J,2)=0
	       INC(NPR+J)=NC
               FC(NPR+J)=I
   40       CONTINUE
            NPR=NPR+ICD
         ELSE
C        DETERMINE PRINCIPAL QUANTUM NUMBER(NQN)
C        OF ORBITALS TO BE USED
C
            NQN=2
            IF(IAN(I) .GT. 10 .AND. IAN(I) .LE. 18) NQN=3
            IF(IAN(I) .GT. 18 .AND. IAN(I) .LE. 36) NQN=4
            IF(IAN(I) .GT. 36 .AND. IAN(I) .LE. 54) NQN=5
C
            NC=NC+1
            DO 50 J=1,ICD
               CC(NPR+J)=ALLC(J,NQN,1)
               EX(NPR+J)=ALLZ(J,NQN,1)*ZS(IAN(I))**2
               CEN(NPR+J,1)=CO(1,I)/BOHR
               CEN(NPR+J,2)=CO(2,I)/BOHR
               CEN(NPR+J,3)=CO(3,I)/BOHR
               IAM(NPR+J,1)=0
               IAM(NPR+J,2)=0
	       INC(NPR+J)=NC
   50       CONTINUE
            NPR=NPR+ICD
            DO 70 K=1,3
               NC=NC+1
               DO 60  J=1,ICD
                  CC(NPR+J)=ALLC(J,NQN,2)
                  EX(NPR+J)=ALLZ(J,NQN,2)*ZP(IAN(I))**2
                  CEN(NPR+J,1)=CO(1,I)/BOHR
                  CEN(NPR+J,2)=CO(2,I)/BOHR
                  CEN(NPR+J,3)=CO(3,I)/BOHR
                  IAM(NPR+J,1)=1
                  IAM(NPR+J,2)=K
	          INC(NPR+J)=NC
   60          CONTINUE
               NPR=NPR+ICD
   70       CONTINUE
         ENDIF
   80 CONTINUE
C
C     CALCULATE NORMALIZATION CONSTANTS AND INCLUDE
C     THEM IN THE CONTRACTION COEFFICIENTS
C
      DO 90 I=1,NPR
         NORM=(2.D0*EX(I)/PI)**0.75D0*(4.D0*EX(I))**(IAM(I,1)/2.D0)/
     1   SQRT(DEX(2*IAM(I,1)-1))
         CC(I)=CC(I)*NORM
   90 CONTINUE
      IPR=0
C
C     PERFORM SORT OF PRIMITIVES BY ANGULAR MOMENTUM
C
C *** IS is count of S primitives???
C *** IP is count of P primitives???
      IS=0
      IP=0
      IPC=0
      ISC=0
      J=0
      DO 100 I=1,NPR
         IF (IAM(I,1) .EQ. 0) THEN
            IS=IS+1
            IND(IS)=I
         ENDIF
  100 CONTINUE
      IP=IS
      DO 110 I=1,NPR
         IF (IAM(I,1) .EQ. 1 .AND. IAM(I,2) .EQ. 1) THEN
            IP=IP+1
            IND(IP)=I
         ENDIF
  110 CONTINUE
      DO 120 I=1,NPR
         IF (IAM(I,1) .EQ. 1 .AND. IAM(I,2) .EQ. 2) THEN
            IP=IP+1
            IND(IP)=I
         ENDIF
  120 CONTINUE
      DO 130 I=1,NPR
         IF (IAM(I,1) .EQ. 1 .AND. IAM(I,2) .EQ. 3) THEN
            IP=IP+1
            IND(IP)=I
         ENDIF
  130 CONTINUE
      DO 140 I=1,NC
         IN=I*ICD-ICD+1
         IF (IAM(IN,1) .EQ. 0) THEN
            ISC=ISC+1
            INDC(ISC)=I
         ENDIF
  140 CONTINUE
      IPC=ISC
      DO 150 I=1,NC
         IN=I*ICD-ICD+1
         IF (IAM(IN,1) .EQ. 1 .AND. IAM(IN,2) .EQ. 1) THEN
            IPC=IPC+1
            INDC(IPC)=I
         ENDIF
  150 CONTINUE
      DO 160 I=1,NC
         IN=I*ICD-ICD+1
         IF (IAM(IN,1) .EQ. 1 .AND. IAM(IN,2) .EQ. 2) THEN
            IPC=IPC+1
            INDC(IPC)=I
         ENDIF
  160 CONTINUE
      DO 170 I=1,NC
         IN=I*ICD-ICD+1
         IF (IAM(IN,1) .EQ. 1 .AND. IAM(IN,2) .EQ. 3) THEN
            IPC=IPC+1
            INDC(IPC)=I
         ENDIF
  170 CONTINUE
      DO 180 I=1,NPR
         TEMP(I)=CC(IND(I))
  180 CONTINUE
      DO 190 I=1,NPR
         CC(I)=TEMP(I)
  190 CONTINUE
      DO 200 I=1,NPR
         TEMP(I)=EX(IND(I))
  200 CONTINUE
      DO 210 I=1,NPR
         EX(I)=TEMP(I)
  210 CONTINUE
      DO 220 I=1,NPR
         TEMP(I)=CEN(IND(I),1)
  220 CONTINUE
      DO 230 I=1,NPR
         CEN(I,1)=TEMP(I)
  230 CONTINUE
      DO 240 I=1,NPR
         TEMP(I)=CEN(IND(I),2)
  240 CONTINUE
      DO 250 I=1,NPR
         CEN(I,2)=TEMP(I)
  250 CONTINUE
      DO 260 I=1,NPR
         TEMP(I)=CEN(IND(I),3)
  260 CONTINUE
      DO 270 I=1,NPR
         CEN(I,3)=TEMP(I)
  270 CONTINUE
      DO 280 I=1,NPR
         ITEMP(I)=IAM(IND(I),1)
  280 CONTINUE
      DO 290 I=1,NPR
         IAM(I,1)=ITEMP(I)
  290 CONTINUE
      DO 300 I=1,NPR
         ITEMP(I)=IAM(IND(I),2)
  300 CONTINUE
      DO 310 I=1,NPR
         IAM(I,2)=ITEMP(I)
  310 CONTINUE
C *** also arrange our new array INC() like the others...
C *** also arrange our new array INC() like the others...
C *** also arrange our new array INC() like the others...
      DO 315 I=1,NPR
         ITEMP(I)=INC(IND(I))
  315 CONTINUE
      DO 316 I=1,NPR
         INC(I)=ITEMP(I)
  316 CONTINUE
C     CALCULATE OVERLAP MATRIX OF STO-6G FUNCTIONS
C
      DO 320 J=1,NC
         CALL OVLP(J,1,IS,IP,NPR,NC,ICD)
  320 CONTINUE
C
      DO 330 J=1,NC
         DO 330 K=1,NC
            CESPM2(INDC(J),INDC(K))=OVL(J,K)
  330 CONTINUE
      DO 340 J=1,NC
         DO 340 K=1,NC
            OVL(J,K)=CESPM2(J,K)
  340 CONTINUE
      L=0
      DO 350 I=1,NC
         DO 350 J=1,I
            L=L+1
            CESP(L)=OVL(I,J)
  350 CONTINUE
C
C     DEORTHOGONALIZE THE COEFFICIENTS AND REFORM THE DENSITY MATRIX
C
      CALL RSP(CESP,NC,1,TEMP,CESPML)
      DO 360 I=1,NC
         DO 360 J=1,I
            SUM=0.D0
            DO 360 K=1,NC
               SUM=SUM+CESPML(I+(K-1)*NC)/SQRT(TEMP(K))*CESPML(J+(K-1)*N
     1C)
               CESP(I+(J-1)*NC)=SUM
               CESP(J+(I-1)*NC)=SUM
  360 CONTINUE
      CALL MULT(C,CESP,CESPML,NC)
      CALL DENSIT(CESPML,NC,NC,NCLOSE,NOPEN,FRACT,CESP,2)
C *** does CESPML now contain the eigenvectors??? and TEMP the eigenvalues???
C *** does CESPML now contain the eigenvectors??? and TEMP the eigenvalues???
C *** does CESPML now contain the eigenvectors??? and TEMP the eigenvalues???
      RETURN
      END
      SUBROUTINE GETESP
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C *** this is the end part of ELESP subroutine.
C *** this is the end part of ELESP subroutine.
C *** this is the end part of ELESP subroutine.
      CHARACTER*241 KEYWRD
      DOUBLE PRECISION NORM,OVL
      LOGICAL CALLED,POTWRT,RST,STO3G
      INCLUDE 'SIZES'
      COMMON/ESPF/ AL((NUMATM+4)**2),A(NUMATM,NUMATM),B(NUMATM),
     1Q(NUMATM+4),CESPM(MAXORB,MAXORB)
      COMMON /DENSTY/ P(MPACK),PA(MPACK),PB(MPACK)
      COMMON /POTESP/ XC,YC,ZC,ESPNUC,ESPELE,NESP
      COMMON /ABC/    CO(3,NUMATM),IAN(NUMATM),NATOM
      COMMON /WORK1/  POTPT(3,MESP), ES(MESP), ESP(MESP), WORK1D(2*MESP)
      COMMON /STO6G/  ALLC(6,6,2),ALLZ(6,6,2)
      COMMON /VECTOR/ C(MORB2*2+MAXORB*2)
      COMMON /MOLKST/ NUMAT,NAT(NUMATM),NFIRST(NUMATM),NMIDLE(NUMATM),
     1                NLAST(NUMATM), NORBS, NELECS,NALPHA,NBETA,
     2                NCLOSE,NOPEN,NDUMY,FRACT
      COMMON /KEYWRD/ KEYWRD
      COMMON /ESPC/  CC(MAXPR),CEN(MAXPR,3),IAM(MAXPR,2),IND(MAXPR),
     1                EX(MAXPR),ESPI(MAXORB,MAXORB),
     2                FV(0:8,821),FAC(0:7),
     3                DEX(-1:96),TF(0:2),TEMP(MAXPR),ITEMP(MAXPR),
     4                OVL(MAXORB,MAXORB),FC(MAXPR*6)
     6       /CORE  / TORE(107)
     7       /EXPONT/ ZS(107),ZP(107),ZD(107)
*
*  END OF MINDO/3 COMMON BLOCKS
*
      COMMON /INDX/   INDC(MAXORB)
C *** an additional common block that carries variables for plotting routines.
      COMMON /PLOTS/  CESPM2(MAXORB,MAXORB),SLA(10),
     1                CESPML(MAXORB*MAXORB),CESP(MAXORB*MAXORB),
     2                INC(MAXPR),NC,NPR,IS,IP,IPC,ISC,ICD,IORB
C *** old arrays that are no longer needed are here, commented out.
C     DIMENSION CESPM2(MAXORB,MAXORB),SLA(10)
C     DIMENSION CESPML(MAXORB*MAXORB),CESP(MAXORB*MAXORB)
      DATA BOHR/0.529167D0/
C *** end of ELESP starts here...
C
C     NOW CALCULATE THE ELECTRONIC CONTRIBUTION TO THE ELECTROSTATIC POT
C
      L=0
      DO 370 I=1,NC
         DO 370 J=1,I
            L=L+1
            CESPM(I,J)=CESP(L)
            CESPM(J,I)=CESP(L)
  370 CONTINUE
      IPX=(NPR-IS)/3
      IPE=IS+IPX
      DO 380 I=1,NESP
         ES(I)=0.D0
  380 CONTINUE
      CALL NAICAS(ISC,IS,IP,NPR,NC,IPE,IPX,ICD)
      CALL NAICAP(ISC,IS,IP,NPR,NC,IPE,IPX,ICD)
C     CALCULATE TOTAL ESP AND FORM ARRAYS FOR ESPFIT
      DO 400 I=1,NESP
         ESP(I)=0.D0
         DO 390 J=1,NATOM
            RA=SQRT((CO(1,J)-POTPT(1,I))**2+(CO(2,J)-POTPT(2,I))**2+(CO(
     13,J)-POTPT(3,I))**2)
            ESP(I)=ESP(I)+TORE(IAN(J))/(RA/BOHR)
  390    CONTINUE
         ESP(I)=ESP(I)-ES(I)
         DO 400  J=1,NATOM
            RIJ=SQRT((CO(1,J)-POTPT(1,I))**2+(CO(2,J)-POTPT(2,I))**2
     1+(CO(3,J)-POTPT(3,I))**2)/BOHR
            B(J)=B(J)+ESP(I)*1.D0/RIJ
  400 CONTINUE
C
C     IF REQUESTED WRITE OUT ELECTRIC POTENTIAL DATA TO
C     UNIT 21
C
      POTWRT=(INDEX(KEYWRD,'POTWRT') .NE. 0)
      IF(POTWRT) THEN
         OPEN(UNIT=21)
         WRITE(21,'(I5)') NESP
         DO 410 I=1,NESP
  410    WRITE(21,420) ESP(I),POTPT(1,I)/BOHR,POTPT(2,I)/BOHR,
     1POTPT(3,I)
      ENDIF
  420 FORMAT(1X,4E16.7)
      RETURN
      END
      SUBROUTINE GETORB
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C *** this will calculate values for orbital plots...
C *** this will calculate values for orbital plots...
C *** this will calculate values for orbital plots...
      CHARACTER*241 KEYWRD
      DOUBLE PRECISION NORM,OVL
      LOGICAL CALLED,POTWRT,RST,STO3G
      INCLUDE 'SIZES'
      COMMON/ESPF/ AL((NUMATM+4)**2),A(NUMATM,NUMATM),B(NUMATM),
     1Q(NUMATM+4),CESPM(MAXORB,MAXORB)
      COMMON /DENSTY/ P(MPACK),PA(MPACK),PB(MPACK)
      COMMON /POTESP/ XC,YC,ZC,ESPNUC,ESPELE,NESP
      COMMON /ABC/    CO(3,NUMATM),IAN(NUMATM),NATOM
      COMMON /WORK1/  POTPT(3,MESP), ES(MESP), ESP(MESP), WORK1D(2*MESP)
      COMMON /STO6G/  ALLC(6,6,2),ALLZ(6,6,2)
      COMMON /VECTOR/ C(MORB2*2+MAXORB*2)
      COMMON /MOLKST/ NUMAT,NAT(NUMATM),NFIRST(NUMATM),NMIDLE(NUMATM),
     1                NLAST(NUMATM), NORBS, NELECS,NALPHA,NBETA,
     2                NCLOSE,NOPEN,NDUMY,FRACT
      COMMON /KEYWRD/ KEYWRD
      COMMON /ESPC/  CC(MAXPR),CEN(MAXPR,3),IAM(MAXPR,2),IND(MAXPR),
     1                EX(MAXPR),ESPI(MAXORB,MAXORB),
     2                FV(0:8,821),FAC(0:7),
     3                DEX(-1:96),TF(0:2),TEMP(MAXPR),ITEMP(MAXPR),
     4                OVL(MAXORB,MAXORB),FC(MAXPR*6)
     6       /CORE  / TORE(107)
     7       /EXPONT/ ZS(107),ZP(107),ZD(107)
*
*  END OF MINDO/3 COMMON BLOCKS
*
      COMMON /INDX/   INDC(MAXORB)
C *** an additional common block that carries variables for plotting routines.
      COMMON /PLOTS/  CESPM2(MAXORB,MAXORB),SLA(10),
     1                CESPML(MAXORB*MAXORB),CESP(MAXORB*MAXORB),
     2                INC(MAXPR),NC,NPR,IS,IP,IPC,ISC,ICD,IORB
C *** old arrays that are no longer needed are here, commented out.
C     DIMENSION CESPM2(MAXORB,MAXORB),SLA(10)
C     DIMENSION CESPML(MAXORB*MAXORB),CESP(MAXORB*MAXORB)
      DATA BOHR/0.529167D0/
      ESP(1)=0.D0
C *** variable I loops over all gaussian primitives.
C *** we calculate value of the primitive to PRIM and weight it according to the eigenvector.
C *** eigenvector contains weights for contracted functions; the array INC() contains contraction indices.
      DO 500 I=1,NPR
         DX=POTPT(1,1)-CEN(I,1)
         DY=POTPT(2,1)-CEN(I,2)
         DZ=POTPT(3,1)-CEN(I,3)
         TD=DX**2+DY**2+DZ**2
         PRIM=CC(I)*EXP(-EX(I)*TD)
         IF(IAM(I,2) .EQ. 1) THEN
           PRIM=PRIM*DX
         ENDIF
         IF(IAM(I,2) .EQ. 2) THEN
           PRIM=PRIM*DY
         ENDIF
         IF(IAM(I,2) .EQ. 3) THEN
           PRIM=PRIM*DZ
         ENDIF
         ESP(1)=ESP(1)+CESPML(INC(I)+(IORB-1)*NC)*PRIM
500   CONTINUE
      RETURN
      END
      SUBROUTINE GETELDENS
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C *** this will calculate values for the electron density plot...
C *** this will calculate values for the electron density plot...
C *** this will calculate values for the electron density plot...
      CHARACTER*241 KEYWRD
      DOUBLE PRECISION NORM,OVL
      LOGICAL CALLED,POTWRT,RST,STO3G
      INCLUDE 'SIZES'
      COMMON/ESPF/ AL((NUMATM+4)**2),A(NUMATM,NUMATM),B(NUMATM),
     1Q(NUMATM+4),CESPM(MAXORB,MAXORB)
      COMMON /DENSTY/ P(MPACK),PA(MPACK),PB(MPACK)
      COMMON /POTESP/ XC,YC,ZC,ESPNUC,ESPELE,NESP
      COMMON /ABC/    CO(3,NUMATM),IAN(NUMATM),NATOM
      COMMON /WORK1/  POTPT(3,MESP), ES(MESP), ESP(MESP), WORK1D(2*MESP)
      COMMON /STO6G/  ALLC(6,6,2),ALLZ(6,6,2)
      COMMON /VECTOR/ C(MORB2*2+MAXORB*2)
      COMMON /MOLKST/ NUMAT,NAT(NUMATM),NFIRST(NUMATM),NMIDLE(NUMATM),
     1                NLAST(NUMATM), NORBS, NELECS,NALPHA,NBETA,
     2                NCLOSE,NOPEN,NDUMY,FRACT
      COMMON /KEYWRD/ KEYWRD
      COMMON /ESPC/  CC(MAXPR),CEN(MAXPR,3),IAM(MAXPR,2),IND(MAXPR),
     1                EX(MAXPR),ESPI(MAXORB,MAXORB),
     2                FV(0:8,821),FAC(0:7),
     3                DEX(-1:96),TF(0:2),TEMP(MAXPR),ITEMP(MAXPR),
     4                OVL(MAXORB,MAXORB),FC(MAXPR*6)
     6       /CORE  / TORE(107)
     7       /EXPONT/ ZS(107),ZP(107),ZD(107)
*
*  END OF MINDO/3 COMMON BLOCKS
*
      COMMON /INDX/   INDC(MAXORB)
C *** an additional common block that carries variables for plotting routines.
      COMMON /PLOTS/  CESPM2(MAXORB,MAXORB),SLA(10),
     1                CESPML(MAXORB*MAXORB),CESP(MAXORB*MAXORB),
     2                INC(MAXPR),NC,NPR,IS,IP,IPC,ISC,ICD,IORB
C *** old arrays that are no longer needed are here, commented out.
C     DIMENSION CESPM2(MAXORB,MAXORB),SLA(10)
C     DIMENSION CESPML(MAXORB*MAXORB),CESP(MAXORB*MAXORB)
      DATA BOHR/0.529167D0/
      ESP(1)=0.D0
C *** this is quite similar to GETORB, we just loop over all occupied orbitals here...
C *** here we assume that we have an open-shell RHF model...
      ILOOP=NELECS/2
C *** variable I loops over all gaussian primitives.
C *** we calculate value of the primitive to PRIM and weight it according to the eigenvector.
C *** eigenvector contains weights for contracted functions; the array INC() contains contraction indices.
      DO 500 I=1,NPR
         DX=POTPT(1,1)-CEN(I,1)
         DY=POTPT(2,1)-CEN(I,2)
         DZ=POTPT(3,1)-CEN(I,3)
         TD=DX**2+DY**2+DZ**2
         DO 600 J=1,ILOOP
            PRIM=CC(I)*EXP(-EX(I)*TD)
            IF(IAM(I,2) .EQ. 1) THEN
              PRIM=PRIM*DX
            ENDIF
            IF(IAM(I,2) .EQ. 2) THEN
              PRIM=PRIM*DY
            ENDIF
            IF(IAM(I,2) .EQ. 3) THEN
              PRIM=PRIM*DZ
            ENDIF
            ORB=CESPML(INC(I)+(J-1)*NC)*PRIM
C *** here we assume that we have an open-shell RHF model...
            ESP(1)=ESP(1)+ORB*ORB*2.0D0
600      CONTINUE
500   CONTINUE
      RETURN
      END