C
C       ***********************************************************
C       ********************  PSI/88 - PART 1  ********************
C       ***********************************************************
C
C       Version 1.0  Any questions to the author should specify
C                    the version being used.
C
      PROGRAM PSI1
      IMPLICIT REAL (A-H,O-Z)
C
C      DANIEL L. SEVERANCE
C      WILLIAM L. JORGENSEN
C      DEPARTMENT OF CHEMISTRY
C      YALE UNIVERSITY
C      NEW HAVEN, CT 06511
C
C      THIS CODE DERIVED FROM THE PSI1 PORTION OF THE ORIGINAL PSI77
C      PROGRAM WRITTEN BY WILLIAM L. JORGENSEN, PURDUE.
C      IT HAS BEEN REWRITTEN TO ADD SPEED AND BASIS FUNCTIONS. DLS
C
C      THE CONTOURING CODE HAS BEEN MOVED TO A SEPARATE PROGRAM TO ALLOW
C      MULTIPLE CONTOURS TO BE PLOTTED WITHOUT RECOMPUTING THE
C      ORBITAL VALUE MATRIX.
C
C Redistribution and use in source and binary forms are permitted
C provided that the above paragraphs and this one are duplicated in 
C all such forms and that any documentation, advertising materials,
C and other materials related to such distribution and use acknowledge 
C that the software was developed by Daniel Severance at Purdue University
C The name of the University or Daniel Severance may not be used to endorse 
C or promote products derived from this software without specific prior 
C written permission.  The authors are now at Yale University.
C THIS SOFTWARE IS PROVIDED ``AS IS'' AND WITHOUT ANY EXPRESS OR
C IMPLIED WARRANTIES, INCLUDING, WITHOUT LIMITATION, THE IMPLIED
C WARRANTIES OF MERCHANTIBILITY AND FITNESS FOR A PARTICULAR PURPOSE.
C
C      THIS PROGRAM READS A WAVEFUNCTION FROM A SPECIFIED BASIS SET AND
C      COMPUTES A 3 DIMENSIONAL MATRIX OF ORBITAL VALUE.
C      THIS MATRIX IS WRITTEN TO A BINARY FILE FOR INPUT TO THE
C      CONTOURING PROGRAM.  THE BINARY FORMAT IS CHOSEN TO REDUCE I/O
C      TIME AS WELL AS DISK STORAGE SPACE.
C      CARD INPUT IS DESCRIBED IN SUBROUTINE MOGRID.
C
      COMMON /IO/ IRD,ILST
      IRD = 5
      ILST = 6
C
C       MAKE 3-D ORBITAL VALUE OR DENSITY GRID AND SAVE TO DISK.
C       3-D CONTOURS ( WHICH USED TO BE GENERATED IN THREED )
C       WILL BE CREATED IN A SEPARATE PROGRAM ( MOCON )
C
      CALL MOGRID
      STOP
      END
C
C
      SUBROUTINE MOGRID
      IMPLICIT REAL (A-H,O-Z)
C
C       PROGRAM TO PLOT SEMI EMPIRICAL, STO-3G, 3-21[++]G[(*)] AND
C       6-31[++]G[(**)] WAVEFUNCTIONS IN 3 DIMENSIONS FOR ATOMS H - AR
C       W.L. JORGENSEN - 12/19/71, MAJOR REVISION 5/15/77
C       D.L. SEVERANCE - COMPLETE REWRITE OF PSI1, SPLIT OFF CONTOURING
C           AND MODIFICATION OF PSI2 H.L. ALGORITHM 12/87
C
C       THE CONTOURING ROUTINES HAVE BEEN REMOVED COMPLETELY AND PLACED
C       IN A STAND-ALONE PROGRAM USING THE OUTPUT FROM THIS PROGRAM.
C       CODES DEALING WITH COMPUTATION OF THE ORBITAL VALUE OR DENSITY
C       ARE COMPLETELY DIFFERENT AND EACH IS PLACED IN A SEPARATE
C       SUBROUTINE. STOMO, MOSEMI, MO321, AND MO631.
C       PRESENTLY DOES STO-3G,3-21[++]G[(*)],6-31[++]G[**] THE PARAMETER
C       BRACKETS ARE OPTIONAL, I.E. 6-31+G**, 6-31G*, ETC.
C       SEMI-EMPIRICAL WF'S NOW SUPPORTED BY MOSEMI.
C       THE MO631 ROUTINE IS EASILY MODIFIED TO SUPPORT OTHER BASIS
C       FUNCTIONS OF THE K-LMG VARIETY, THE 3-21G ROUTINE WAS DERIVED
C       FROM IT DIRECTLY.
C
C       THIS ROUTINE AND THE CONTOURING ROUTINES SHARE THE SAME INPUT
C       FILE, THIS ROUTINE IS A STRIPPED DOWN VERSION OF THE ORIGINAL
C       THREED ROUTINE.
C
      PARAMETER (MXPTS=51)
      PARAMETER (MAXATM=50)
      PARAMETER (MAXORB=200)
      COMMON /DENS/ DENSIT(MXPTS,MXPTS,MXPTS),V(MAXORB,MAXORB),C(MAXATM,
     *   3),OCMO(MAXORB),IAN(MAXATM)
      COMMON /SPLIT/ XMI,YMI,ZMI,XINC,YINC,ZINC,MONE,MOLAST,NAT,
     *   IONEMO
      COMMON /IO/ IRD,ILST
      COMMON /BASIS/ CALC
      CHARACTER CALC*20,AUTO*4,TITLE*80
      DIMENSION CO(3),CMIN(3),CMAX(3)
      DATA AU / 0.5291670E+0 /
C
C       V IS THE EIGENVECTOR MATRIX. DENSIT IS THE CHARGE DENSITY OR
C       ORBITAL VALUE MATRIX. IAN IS THE ARRAY OF ATOMIC NUMBERS
C       FOR THE ATOMS  WHOSE COORDINATES, READ IN IN ANGSTROMS, ARE IN
C       C. 
C
C       CARD INPUT IS AS FOLLOWS - THIS IS NOT ALL USED IN THIS PROGRAM
C               BUT THE SAME INPUT DECK IS USED BY THE CONTOURING
C               PACKAGE
C
C       CARD  1 - TYPE OF WAVEFUNCTION = STO3G , 3-21G, 3-21+G, 3-21G*,
C                   6-31G TO 6-31++G**. (A20)
C
C       CARD  2 - COLUMN 1-4 = 'AUTO' FOR AUTOMATIC DETERMINATION OF THE
C                    LIMITS WITHIN WHICH THE COMPUTATION IS DONE.
C                    IF USED, SKIP TO CARD 5. THE OTHER PARAMETER
C                    ON THIS LINE IS STILL ACTIVE. (A4)
C                 COLUMN 5 IONEMO = 1 IF THE MO IS TO BE CARD INPUT.
C                    SEE BELOW FOR AO ORDER.
C
C       THE NEXT 3 CARDS ONLY IF AUTO NOT SPECIFIED:
C
C       CARD  2 - XMIN, XMAX OF THE PLOTTING REGION, 2F10.6
C       CARD  3 - YMIN, YMAX OF THE PLOTTING REGION, 2F10.6
C       CARD  4 - ZMIN, ZMAX OF THE PLOTTING REGION, 2F10.6
C
C       *** FOR THE PRESENT, USE ONLY MONE=MOLAST ON CARD 5 ***
C
C       CARD  5 - FIRST MO, LAST MO, SCALE (FOR AUTO SCALE)- 2I2, F10.6
C                 FOR AN ORBITAL VALUE PLOT, MONE=MOLAST=THE NUMBER OF
C                 THE DESIRED MO. IF IONEMO=1, MONE=MOLAST=01.
C
C       ***IF IONEMO=1, THE GEOMETRY CARDS AND A 99 CARD FOLLOWED BY
C       THE ONE MO TO BE PLOTTED (8F10.6) ARE INSERTED HERE ***
C
C       FORMAT FOR GEOMETRY CARDS AND WAVEFUNCTION INPUT:
C
C       CARD  1 - CHARGE,TITLE  (I2,A80)
C       CARD  2 - ATOMIC NUMBER,X COORD,Y COORD, AND Z COORD OF ATOM 1
C                 (I2,8X,3F10.6) CC 1-2 FOR AN, CC 11-20,21-30,31-40
C                 FOR CARTESIAN COORDINATES
C
C       REPEAT CARD 2 FOR EACH ATOM IN THE MOLECULE - END GEOMETRY
C       INPUT WITH 99 FOR THE ATOMIC NUMBER
C
C       WAVEFUNCTION INPUT:
C               BEGINNING IMMEDIATELY AFTER THE 99 CARD OF THE
C               GEOMETRY INPUT, PLACE THE SINGLE MO TO BE PLOTTED
C               TO BE READ WITH 8F10.6 FORMAT (8 AO'S PER LINE,
C               CC 1-10, 11-20, 21-30, ETC.
C
C       ***IF WAVEFUNCTION IS **SEMI** THEN
C       CARD  3 - ZETA VALUES, ONE ATOM PER LINE S,P FOR EACH ATOM
C                 IN THE MOLECULE.  2F10.6 - D ORBITALS NOT IMPLEMENTED
C
      READ (IRD,10) CALC,AUTO,IONEMO
   10 FORMAT (A/A,I1)
C
C     CHECK EARLY FOR A VALID BASIS SET
C
      CALL UPCASE (CALC)
      IF (INDEX(CALC,'STO-3G').EQ.0) THEN
         IF (INDEX(CALC,'SEMI').EQ.0) THEN
            IF (INDEX(CALC,'6-31').EQ.0) THEN
               IF (INDEX(CALC,'3-21').EQ.0) THEN
                  WRITE (ILST,*)
     *               'PSI1 IS NOT SET UP TO DO THE CALCULATIONS ON'
      WRITE (ILST,*) 'THIS BASIS SET. YOU CAN USE STO-3G, 3-21G, 3-21+G,
     *'
                  WRITE (ILST,*)
     *               '6-31G THROUGH 6-31++G** AND SEMI_EMPIRICAL.'
                  STOP
               ENDIF
            ENDIF
         ENDIF
      ENDIF
      CALL UPCASE (AUTO)
      IF (AUTO.NE.'AUTO') THEN
C
C       READ ABSCISSA AND ORDINATE MIN AND MAX IF NOT AUTO
C
         READ (IRD,20) XMIN,XMAX
         READ (IRD,20) YMIN,YMAX
         READ (IRD,20) ZMIN,ZMAX
      ENDIF
   20 FORMAT (2F10.6)
C
C       READ IN MONE AND MOLAST
C
      READ (IRD,30) MONE,MOLAST,SCALE
   30 FORMAT (2I2,F10.6)
C
C       OBTAIN PLOTTING DATA FROM DISK
C
      READ (IRD,40) ICHG,TITLE
   40 FORMAT (I2,A)
      NOEL = -ICHG
      NAT = 1
C
C       READ COORDINATES AND DETERMINE NORMALIZING FACTORS FOR AOS.
C       THE FORM OF THE SLATER ORBITALS MAY BE FOUND IN I.G. CSIZMADIA,
C       THEORY AND PRACTICE OF MO CALCULATIONS,ELSEVIER,1976, P 313.
C
   50 READ (IRD,60) IN,(C(NAT,J),J=1,3)
   60 FORMAT (I2,8X,3F10.6)
C
C       COMPUTATION OF THE NUMBER OF AO'S AND READING IN THE
C       WAVEFUNCTION MOVED TO THE INDIVIDUAL SUBROUTINES
C       IT IS EASIER TO ADD NEW WAVEFUNCTION ROUTINES WITHOUT
C       HAVING TO MODIFY THE MAINLINE WITH SPAGHETTI.
C       THANKS TO JIM BRIGGS - PURDUE FOR NOTING THE PROBLEM.
C
C       DLS 1-17-87
C
      IF (IN.NE.99) THEN
         IAN(NAT) = IN
         NAT = NAT+1
         GO TO 50
      ENDIF
      NAT = NAT-1
      CALL DRAMNP (C,NAT,CO,ICM,CMIN,CMAX)
C
C       CONVERT COORDINATES TO ATOMIC UNITS
C
      AUINV = 1.0E+0/AU
      DO 80 J = 1, 3
         DO 70 I = 1, NAT
            C(I,J) = C(I,J)*AUINV
   70    CONTINUE
   80 CONTINUE
C
C       NAT=NO. OF ATOMS, NOEL=NO. OF ELECTRONS, N= NO.
C       OF BASIS FUNCTIONS.
C
      NMO = (NOEL+1)/2
C
C       DEFAULT VALUES
C
      IF (MONE.EQ.0) MONE = 1
      IF (MOLAST.EQ.0) MOLAST = NMO
      IF (IONEMO.EQ.1) THEN
         MONE = 1
         MOLAST = 1
      ENDIF
      IF (SCALE.EQ.0.0) SCALE = 1.0E+0
C
C       GREATER RESOLUTION IS OBTAINABLE BY INCREASING THE
C       SIZE OF MXPTS WHICH IS THE DIMENSIONS OF THE DENSITY
C       OR ORBITAL VALUE ARRAY.
C       MXPTS MUST BE ODD.
C
C       SEE THE PARAMETER STATEMENT AT TOP FOR THIS VALUE
C
      SPACES = FLOAT(MXPTS-1)*AU
C
C       DETERMINE DEFAULT RANGES, IF NOT AUTO.  INCREASING SCALE
C       INCREASES THE RANGE OF THE PLOT.
C
      IF (AUTO.EQ.'AUTO') THEN
         R = 1.30E+0*SCALE
         XMIN = CMIN(1)-R
         XMAX = CMAX(1)+R
         YMIN = CMIN(2)-R
         YMAX = CMAX(2)+R
         ZMIN = CMIN(3)-R
         ZMAX = CMAX(3)+R
      ENDIF
      XMI = XMIN/AU
      YMI = YMIN/AU
      ZMI = ZMIN/AU
      XINC = (XMAX-XMIN)/SPACES
      YINC = (YMAX-YMIN)/SPACES
      ZINC = (ZMAX-ZMIN)/SPACES
C
C       CALL THE APPROPRIATE ROUTINE TO GENERATE THE 3-D ORBITAL
C       VALUE OR DENSITY MATRIX TO WHICH THE CONTOURS WILL BE APPLIED.
C
      IF (INDEX(CALC,'STO-3G').NE.0) THEN
         CALL STOMO
C
C       SEND ALL 3-21... VARIANTS TO MO321G, 3-21G,3-21+G,3-21++G,
C
      ELSEIF (INDEX(CALC,'3-21').NE.0) THEN
         CALL MO321G
C
C       SEND ALL 6-31... VARIANTS TO MO631G, 6-31G,6-31+G,6-31++G,
C       AS WELL AS ALL COMBINATIONS WITH *, **
C
      ELSEIF (INDEX(CALC,'6-31').NE.0) THEN
         CALL MO631G
      ELSEIF (INDEX(CALC,'SEMI').NE.0) THEN
C
C       GENERIC ROUTINE REQUIRING THAT THE ZETA VALUES FOR THE SHELL
C       BE AT THE END OF THE WAVEFUNCTION READ IN. 8F10.6, ONE ATOM
C       PER LINE.  IF THE SEMI-EMPIRICAL METHOD IS OTHER THAN S=P,
C       OR INCLUDES D FUNCTIONS, MODIFICATION WILL BE NECESSARY.
C
         CALL MOSEMI
      ELSE
         WRITE (ILST,*) ' NO COMPUTATIONS DONE - UNKNOWN BASIS? ',CALC
         STOP
      ENDIF
C
C       WRITE OUT THE COMPUTED DENSITY MATRIX
C
      OPEN (17,FILE='FOR017',FORM='UNFORMATTED',STATUS='UNKNOWN')
      MAXPTS = MXPTS
C
C     COMPUTE MIN AND MAX DENSITY VALUES COMPUTED, SHOW TO THE USER
C
      DMIN = 1000.0E+0
      DMAX = -1000.0E+0
      DO 110 K = 1, MAXPTS
         DO 100 J = 1, MAXPTS
            DO 90 I = 1, MAXPTS
               DMIN = MIN(DMIN,DENSIT(I,J,K))
               DMAX = MAX(DMAX,DENSIT(I,J,K))
   90       CONTINUE
  100    CONTINUE
  110 CONTINUE
      WRITE (ILST,*) 'MIN, MAX DENSITY(VALUE) COMPUTED IS ',DMIN,DMAX
      WRITE (17) MAXPTS,NAT
      WRITE (17) (IAN(I),I=1,NAT)
      WRITE (17) ((C(I,J),J=1,3),I=1,NAT)
      WRITE (17) (((DENSIT(NX,NY,NZ),NX=1,MAXPTS),NY=1,MAXPTS),NZ=1,
     *   MAXPTS)
      WRITE (17) XMIN,XMAX,YMIN,YMAX,ZMIN,ZMAX
      RETURN
      END
C
C
      SUBROUTINE DRAMNP (C,NAT,CO,ICM,CMIN,CMAX)
      IMPLICIT REAL (A-H,O-Z)
      PARAMETER (MAXATM=50)
      DIMENSION C(MAXATM,3),CO(3),CMIN(3),CMAX(3)
C
C       THIS ROUTINE DETERMINES CO AND CM WHICH ARE USED
C       FOR AUTOMATIC SCALING OF PLOTTING AREAS
C
C       THE ROUTINE WAS ADAPTED FROM THE ROUTINE CALLED
C       DRAMOL WHICH IS USED IN THE HIDDEN LINE PART OF
C       THE PROGRAM
C
      CM = -100.0E+0
      DO 20 I = 1, 3
         CMI = 100.0E+0
         CMA = -100.0E+0
         DO 10 J = 1, NAT
            P = C(J,I)
            CMI = MIN(CMI,P)
            CMA = MAX(CMA,P)
   10    CONTINUE
         CO(I) = (CMA+CMI)*0.50E+0
         CMIN(I) = CMI
         CMAX(I) = CMA
         P = CMA-CMI
         IF (P.GT.CM) THEN
            ICM = I
            CM = P
         ENDIF
   20 CONTINUE
      RETURN
      END
C
C
      SUBROUTINE UPCASE (STRNG)
C
C     CONVERT LOWER CASE CHARACTERS TO UPPER CASE
C
      CHARACTER*(*) STRNG
      INTEGER LNGTH,I,CHRVAL
C
      LNGTH = LEN(STRNG)
      DO 10 , I = 1, LNGTH
         CHRVAL = ICHAR(STRNG(I:I))
         IF (CHRVAL.GE.97.AND.CHRVAL.LE.122) THEN
            CHRVAL = CHRVAL-32
            STRNG(I:I) = CHAR(CHRVAL)
         ENDIF
   10 CONTINUE
      RETURN
      END
C
C
      SUBROUTINE MO631G
      IMPLICIT REAL (A-H,O-Z)
C
C       THIS DETERMINES THE DENSITY OF PLANES COMPUTED
C
      PARAMETER (MXPTS=51)
      PARAMETER (MAXATM=50)
      PARAMETER (MAXORB=200)
C
C       WHEN GENERATING A NEW WAVEFUNCTION ROUTINE, DEFINE THE
C       PARAMETERS KD,LD,AND MD FOR THE K-LMG WAVEFUNCTION
C       BEING USED - 6-31G WILL BE 6,3, AND 1 RESPECTIVELY.  REPLACE
C       THE DATA STATEMENTS FOR A,S,P,APL, AND ASTAR AS NECESSARY,
C       AND CHANGE ALL OCCURANCES OF 6-31 TO WHATEVER BASIS SET
C       YOU ARE DEFINING.   YOU WILL ALSO WANT TO REMOVE ANY CODE
C       NOT APPLICABLE TO THE BASIS (** FOR 3-21, FOR EXAMPLE), OR
C       NOT IMPLIMENTED AT THE PRESENT TIME.
C
C       DAN SEVERANCE 12/24/87
C
      PARAMETER (KD=6,LD=3,MD=1)
C
C     THESE PARAMETERS ARE GENERATED FROM THOSE ABOVE
C
      PARAMETER (K3=KD*3,K2=KD*2)
      PARAMETER (KLM=(KD+LD+MD),LM=(MD+LD),LM3=LM*3,LM2=LM*2)
C
C       WRITTEN AS A STAND-ALONE SUBROUTINE FOR CALCULATING ORBITAL
C       VALUES. WORKED TO REDUCE REDUNDANT COMPUTATION OF POWERS,
C       SQUARE ROOTS, AND EXPONENTIALS.  ALSO WRITTEN TO ALLOW
C       EASY VECTORIZATION BY COMPILERS ON VECTOR MACHINES.
C
C       DAN SEVERANCE 12/22/87-->2/88 MANY MODIFICATIONS
C
      COMMON /SPLIT/ XMI,YMI,ZMI,XINC,YINC,ZINC,MONE,MOLAST,NAT,
     *   IONEMO
      COMMON /DENS/ DENSIT(MXPTS,MXPTS,MXPTS),V(MAXORB,MAXORB),C(MAXATM,
     *   3),OCMO(MAXORB),IAN(MAXATM)
      COMMON /IO/ IRD,ILST
      DIMENSION XDEL(MXPTS),YDEL(MXPTS),ZDEL(MXPTS)
      DIMENSION XDELSQ(MXPTS),YDELSQ(MXPTS),ZDELSQ(MXPTS)
      DIMENSION XYZ(3,50)
      DIMENSION X(MXPTS),Y(MXPTS),Z(MXPTS)
      DIMENSION CNSXX(MXPTS),CNSYY(MXPTS),CNSZZ(MXPTS),AO1S(KD)
      DIMENSION ANEG(50),CNSYZ(MXPTS),CNSXY(MXPTS),CNSXZ(MXPTS)
      DIMENSION AO2S(KD),AO2P(KD),AO2PX(KD),AO2PY(KD),AO2PZ(KD)
      DIMENSION AO3S(KD),AO3P(KD),AO3PX(KD),AO3PY(KD),AO3PZ(KD)
      DIMENSION CNS2PX(MXPTS,KD),CNS2PY(MXPTS,KD),CNS2PZ(MXPTS,KD)
      DIMENSION CNS3PX(MXPTS,KD),CNS3PY(MXPTS,KD),CNS3PZ(MXPTS,KD)
      DIMENSION E1SX(MXPTS,KD),E1SY(MXPTS,KD),E1SZ(MXPTS,KD)
      DIMENSION E2SPX(MXPTS,KD),E2SPY(MXPTS,KD),E2SPZ(MXPTS,KD)
      DIMENSION E3SPX(MXPTS,KD),E3SPY(MXPTS,KD),E3SPZ(MXPTS,KD)
      DIMENSION EXPDX(MXPTS),EXPDY(MXPTS),EXPDZ(MXPTS)
C
C       THE DATA STATEMENTS CONTAINING THE EXPONENTS, S,P AND
C       "STAR" COEFFS WERE READ FROM THE 631.GBS, 631S.GBS, AND 631SS.GB
C       THE LATTER TWO CONTAINING THE STAR COEFFS.  THE APL DATA
C       STATEMENTS WERE TAKEN FROM MO321G ROUTINE, AS THE DIFFUSE
C       EXPONENTS ARE THE SAME.
C
C       D. SEVERANCE 12/22/87
C
C     FORM OF S,P, AND D GAUSSIAN FUNCTIONS:
C       J.COMP.CHEM.,7,359,(1986)
C
C       (^ DENOTES "RAISED TO THE POWER" I.E. R^2 IS R SQUARED)
C       S    TYPE : (   8 * ALPHA^3 /PI^3 )^(1/4)       * E^(-ALPHA*R^2)
C       PX   TYPE : ( 128 * ALPHA^5 /PI^3 )^(1/4) * X   * E^(-ALPHA*R^2)
C       DYZ  TYPE : (2048 * ALPHA^7 /PI^3 )^(1/4) * Y*Z * E^(-ALPHA*R^2)
C
C       FOR PY, REPLACE X WITH Y
C       FOR DX^2 REPLACE Y AND Z WITH X AND X, I.E. Y*Z-->X*X=X^2
C
C       THE FOLLOWING COMMENT LINES (EDITED) COURTESY AOSV:
C       WRITTEN BY JIM BRIGGS  NOV. 1986
C       SUPERCEDED ALONG WITH THE SUPPORT CODE BY MO321G.
C
C----------------------------------------------------------------------
C     -------------------------------------------------------
C     | FOR THE FORM OF THE SPLIT VALENCE AO FUNCTIONS SEE: |
C     -------------------------------------------------------
C     JACS,102,939,(1980)  - INNER SHELL H-NE. | ALSO CONTAINS THE FORM
C     JACS,104,2798,(1982) - INNER SHELL NA-AR | OF THE FUNCTIONS
C     J. CHEM. PHYS.,51,2657(1962) - GAUSSIANS | INCLUDING GAUSSIANS.
C     J. CHEM. PHYS.,80,3265(1984) - | GENERAL DESCRIPTION OF DIFFUSE
C                                    | FUNCTIONS  AND THEIR COEFS.
C
C     NUMERICAL VALUES FOR THE DIFFUSE EXPONENTS WERE TAKEN FROM:
C     "AB INITIO MOLECULAR ORBITAL THEORY",HEHRE,RADOM,SCHLEYER,
C     POPLE, JOHN WILEY & SONS, 1986. PG. 87.
C
C     NORMALIZATION - FOR HYDROGEN AND HELIUM (ZETA = 1.0 FOR OTHERS)
C
C     AO1S = AO1S*(ZETA**1.5)*((2./PI)**0.750E+0)
C     ((2./PI)**0.750E+0) = 0.712705470E+0
C     ZETA = 1.100E+0
C     ZETA*SQRT(ZETA) = 1.153689730E+0
C     1.153689730E+0*0.712705470E+0 = 0.822240980E+0
C
C----------------------------------------------------------------------
C
      COMMON /BASIS/ CALC
      CHARACTER CALC*20
C
C     THE DIMENSIONS OF THE FOLLOWING ARRAYS ARE:
C     16 IS (NROW-1)*KD+LM, WHERE H,HE IS FIRST, ETC. NROW=3
C     18 IS THE MAXIMUM ATOMIC NUMBER IMPLEMENTED
C     10 IS (NROW-2)*KD+LM
C
      DIMENSION A(16,18),S(16,18),P(10,18),APL(18),ASTAR(18)
      DATA (A(I,1),I=1,4) / 0.1873110E+2,0.2825390E+1,0.6401220E+0,
     *   0.1612780E+0 /
      DATA (A(I,2),I=1,4) / 0.3842160E+2,0.5778030E+1,0.1241770E+1,
     *   0.2979640E+0 /
      DATA (A(I,3),I=1,10) / 0.6424190E+3,0.9679850E+2,0.2209110E+2,
     *   0.6201070E+1,0.1935120E+1,0.6367360E+0,0.2324920E+1,0.632430E+
     *   00,0.7905340E-1,0.3596200E-1 /
      DATA (A(I,4),I=1,10) / 0.1264590E+4,0.1899370E+3,0.4315910E+2,
     *   0.1209870E+2,0.3806320E+1,0.1272890E+1,0.3196460E+1,0.7478130E+
     *   00,0.2199660E+0,0.8230990E-1 /
      DATA (A(I,5),I=1,10) / 0.2068880E+4,0.3106500E+3,0.7068300E+2,
     *   0.1986110E+2,0.6299300E+1,0.2127030E+1,0.4727970E+1,0.1190340E+
     *   01,0.3594120E+0,0.1267510E+0 /
      DATA (A(I,6),I=1,10) / 0.3047520E+4,0.4573700E+3,0.1039490E+3,
     *   0.2921020E+2,0.9286660E+1,0.3163930E+1,0.7868270E+1,0.1881290E+
     *   01,0.5442490E+0,0.1687140E+0 /
      DATA (A(I,7),I=1,10) / 0.4173510E+4,0.6274580E+3,0.1429020E+3,
     *   0.4023430E+2,0.1282020E+2,0.4390440E+1,0.1162640E+2,0.2716280E+
     *   01,0.7722180E+0,0.2120310E+0 /
      DATA (A(I,8),I=1,10) / 0.5484670E+4,0.8252350E+3,0.1880470E+3,
     *   0.5296450E+2,0.1689760E+2,0.5799640E+1,0.1553960E+2,0.3599930E+
     *   01,0.1013760E+1,0.2700060E+0 /
      DATA (A(I,9),I=1,10) / 0.7001710E+4,0.1051370E+4,0.2392860E+3,
     *   0.6739740E+2,0.2152000E+2,0.7403100E+1,0.2084800E+2,0.4808310E+
     *   01,0.1344070E+1,0.3581510E+0 /
      DATA (A(I,10),I=1,10) / 0.8425850E+4,0.1268520E+4,0.2896210E+3,
     *   0.8185900E+2,0.2625150E+2,0.9094720E+1,0.2653210E+2,0.6101760E+
     *   01,0.1696270E+1,0.4458190E+0 /
      DATA (A(I,11),I=1,16) / 0.9993200E+4,0.1499890E+4,0.3419510E+3,
     *   0.9467960E+2,0.2973450E+2,0.1000630E+2,0.1509630E+3,0.3558780E+
     *   02,0.1116830E+2,0.3902010E+1,0.1381770E+1,0.4663820E+0,
     *   0.4979660E+0,0.8435290E-1,0.6663500E-1,0.2595440E-1 /
      DATA (A(I,12),I=1,16) / 0.1172280E+5,0.1759930E+4,0.4008460E+3,
     *   0.1128070E+3,0.3599970E+2,0.1218280E+2,0.1891800E+3,0.4521190E+
     *   02,0.1435630E+2,0.5138860E+1,0.1906520E+1,0.7058870E+0,
     *   0.9293400E+0,0.2690350E+0,0.1173790E+0,0.4210610E-1 /
      DATA (A(I,13),I=1,16) / 0.1398310E+5,0.2098750E+4,0.4777050E+3,
     *   0.1343600E+3,0.4287090E+2,0.1451890E+2,0.2396680E+3,0.5744190E+
     *   02,0.1828590E+2,0.6599140E+1,0.2490490E+1,0.9445450E+0,
     *   0.1277900E+1,0.3975900E+0,0.1600950E+0,0.5565770E-1 /
      DATA (A(I,14),I=1,16) / 0.1611590E+5,0.2425580E+4,0.5538670E+3,
     *   0.1563400E+3,0.5006830E+2,0.1701780E+2,0.2927180E+3,0.6987310E+
     *   02,0.2233630E+2,0.8150390E+1,0.3134580E+1,0.1225430E+1,
     *   0.1727380E+1,0.5729220E+0,0.2221920E+0,0.7783690E-1 /
      DATA (A(I,15),I=1,16) / 0.1941330E+5,0.2909420E+4,0.6613640E+3,
     *   0.1857590E+3,0.5919430E+2,0.2003100E+2,0.3394780E+3,0.8101010E+
     *   02,0.2587800E+2,0.9452210E+1,0.3665660E+1,0.1467460E+1,
     *   0.2156230E+1,0.7489970E+0,0.2831450E+0,0.9983170E-1 /
      DATA (A(I,16),I=1,16) / 0.2191710E+5,0.3301490E+4,0.7541460E+3,
     *   0.2127110E+3,0.6798960E+2,0.2305150E+2,0.4237350E+3,0.1007100E+
     *   03,0.3215990E+2,0.1180790E+2,0.4631100E+1,0.1870250E+1,
     *   0.2615840E+1,0.9221670E+0,0.3412870E+0,0.1171670E+0 /
      DATA (A(I,17),I=1,16) / 0.2518010E+5,0.3780350E+4,0.8604740E+3,
     *   0.2421450E+3,0.7733490E+2,0.2624700E+2,0.4917650E+3,0.1169840E+
     *   03,0.3741530E+2,0.1378340E+2,0.5452150E+1,0.2225880E+1,
     *   0.3186490E+1,0.1144270E+1,0.4203770E+0,0.1426570E+0 /
      DATA (A(I,18),I=1,16) / 0.2834830E+5,0.4257620E+4,0.9698570E+3,
     *   0.2732630E+3,0.8736950E+2,0.2968670E+2,0.5758910E+3,0.1368160E+
     *   03,0.4380980E+2,0.1620940E+2,0.6460840E+1,0.2651140E+1,
     *   0.3860280E+1,0.1413730E+1,0.5166460E+0,0.1738880E+0 /
      DATA (S(I,1),I=1,4) / 0.3349460E-1,0.2347270E+0,0.8137570E+0,
     *   0.1000000E+1 /
      DATA (S(I,2),I=1,4) / 0.2376600E-1,0.1546790E+0,0.4696300E+0,
     *   0.1000000E+1 /
      DATA (S(I,3),I=1,10) / 0.2142610E-2,0.1620890E-1,0.7731560E-1,
     *   0.2457860E+0,0.4701890E+0,0.3454710E+0,-.3509170E-1,-.1912330E+
     *   00,0.1083990E+1,0.1000000E+1 /
      DATA (S(I,4),I=1,10) / 0.1944760E-2,0.1483510E-1,0.7209050E-1,
     *   0.2371540E+0,0.4691990E+0,0.3565200E+0,-.1126490E+0,-.2295060E+
     *   00,0.1186920E+1,0.1000000E+1 /
      DATA (S(I,5),I=1,10) / 0.1866270E-2,0.1425150E-1,0.6955160E-1,
     *   0.2325730E+0,0.4670790E+0,0.3634310E+0,-.1303940E+0,-.1307890E+
     *   00,0.1130940E+1,0.1000000E+1 /
      DATA (S(I,6),I=1,10) / 0.1834740E-2,0.1403730E-1,0.6884260E-1,
     *   0.2321840E+0,0.4679410E+0,0.3623120E+0,-.1193320E+0,-.1608540E+
     *   00,0.1143460E+1,0.1000000E+1 /
      DATA (S(I,7),I=1,10) / 0.1834770E-2,0.1399460E-1,0.6858660E-1,
     *   0.2322410E+0,0.4690700E+0,0.3604550E+0,-.1149610E+0,-.1691170E+
     *   00,0.1145850E+1,0.1000000E+1 /
      DATA (S(I,8),I=1,10) / 0.1831070E-2,0.1395020E-1,0.6844510E-1,
     *   0.2327140E+0,0.4701930E+0,0.3585210E+0,-.1107780E+0,-.1480260E+
     *   00,0.1130770E+1,0.1000000E+1 /
      DATA (S(I,9),I=1,10) / 0.1819620E-2,0.1391610E-1,0.6840530E-1,
     *   0.2331860E+0,0.4712670E+0,0.3566190E+0,-.1085070E+0,-.1464520E+
     *   00,0.1128690E+1,0.1000000E+1 /
      DATA (S(I,10),I=1,10) / 0.1884350E-2,0.1433690E-1,0.7010960E-1,
     *   0.2373730E+0,0.4730070E+0,0.3484010E+0,-.1071180E+0,-.1461640E+
     *   00,0.1127770E+1,0.1000000E+1 /
      DATA (S(I,11),I=1,16) / 0.1937660E-2,0.1480700E-1,0.7270550E-1,
     *   0.2526290E+0,0.4932420E+0,0.3131690E+0,-.3542080E-2,-.4395880E-
     *   01,-.1097520E+0,0.1873980E+0,0.6466990E+0,0.3060580E+0,-.
     *   2485030E+0,-.1317040E+0,0.1233520E+1,0.1000000E+1 /
      DATA (S(I,12),I=1,16) / 0.1977830E-2,0.1511400E-1,0.7391080E-1,
     *   0.2491910E+0,0.4879280E+0,0.3196620E+0,-.3237170E-2,-.4100790E-
     *   01,-.1126000E+0,0.1486330E+0,0.6164970E+0,0.3648290E+0,-.
     *   2122900E+0,-.1079850E+0,0.1175840E+1,0.1000000E+1 /
      DATA (S(I,13),I=1,16) / 0.1942670E-2,0.1485990E-1,0.7284940E-1,
     *   0.2468300E+0,0.4872580E+0,0.3234960E+0,-.2926190E-2,-.3740830E-
     *   01,-.1144870E+0,0.1156350E+0,0.6125950E+0,0.3937990E+0,-.
     *   2276060E+0,0.1445830E-2,0.1092790E+1,0.1000000E+1 /
      DATA (S(I,14),I=1,16) / 0.1959480E-2,0.1492880E-1,0.7284780E-1,
     *   0.2461300E+0,0.4859140E+0,0.3250020E+0,-.2780940E-2,-.3571460E-
     *   01,-.1149850E+0,0.9356340E-1,0.6030170E+0,0.4189590E+0,-.
     *   2446300E+0,0.4315720E-2,0.1098180E+1,0.1000000E+1 /
      DATA (S(I,15),I=1,16) / 0.1851600E-2,0.1420620E-1,0.6999950E-1,
     *   0.2400790E+0,0.4847620E+0,0.3352000E+0,-.2782170E-2,-.3604990E-
     *   01,-.1166310E+0,0.9683280E-1,0.6144180E+0,0.4037980E+0,-.
     *   2529230E+0,0.3285170E-1,0.1081250E+1,0.1000000E+1 /
      DATA (S(I,16),I=1,16) / 0.1869240E-2,0.1423030E-1,0.6969620E-1,
     *   0.2384870E+0,0.4833070E+0,0.3380740E+0,-.2376770E-2,-.3169300E-
     *   01,-.1133170E+0,0.5609000E-1,0.5922550E+0,0.4550060E+0,-.
     *   2503740E+0,0.6695700E-1,0.1054510E+1,0.1000000E+1 /
      DATA (S(I,17),I=1,16) / 0.1832960E-2,0.1403420E-1,0.6909740E-1,
     *   0.2374520E+0,0.4830340E+0,0.3398560E+0,-.2297390E-2,-.3071370E-
     *   01,-.1125280E+0,0.4501630E-1,0.5893530E+0,0.4652060E+0,-.
     *   2518270E+0,0.6158900E-1,0.1060180E+1,0.1000000E+1 /
      DATA (S(I,18),I=1,16) / 0.1825260E-2,0.1396860E-1,0.6870730E-1,
     *   0.2362040E+0,0.4822140E+0,0.3420430E+0,-.2159720E-2,-.2907750E-
     *   01,-.1108270E+0,0.2769990E-1,0.5776130E+0,0.4886880E+0,-.
     *   2555920E+0,0.3780660E-1,0.1080560E+1,0.1000000E+1 /
      DATA (P(I,3),I=1,4) / 0.8941510E-2,0.1410090E+0,0.9453640E+0,
     *   0.1000000E+1 /
      DATA (P(I,4),I=1,4) / 0.5598020E-1,0.2615510E+0,0.7939720E+0,
     *   0.1000000E+1 /
      DATA (P(I,5),I=1,4) / 0.7459760E-1,0.3078470E+0,0.7434570E+0,
     *   0.1000000E+1 /
      DATA (P(I,6),I=1,4) / 0.6899910E-1,0.3164240E+0,0.7443080E+0,
     *   0.1000000E+1 /
      DATA (P(I,7),I=1,4) / 0.6757970E-1,0.3239070E+0,0.7408950E+0,
     *   0.1000000E+1 /
      DATA (P(I,8),I=1,4) / 0.7087430E-1,0.3397530E+0,0.7271590E+0,
     *   0.1000000E+1 /
      DATA (P(I,9),I=1,4) / 0.7162870E-1,0.3459120E+0,0.7224700E+0,
     *   0.1000000E+1 /
      DATA (P(I,10),I=1,4) / 0.7190960E-1,0.3495130E+0,0.7199400E+0,
     *   0.1000000E+1 /
      DATA (P(I,11),I=1,10) / 0.5001660E-2,0.3551090E-1,0.1428250E+0,
     *   0.3386200E+0,0.4515790E+0,0.2732710E+0,-.2302250E-1,0.9503590E+
     *   00,0.5985790E-1,0.1000000E+1 /
      DATA (P(I,12),I=1,10) / 0.4928130E-2,0.3498880E-1,0.1407250E+0,
     *   0.3336420E+0,0.4449400E+0,0.2692540E+0,-.2241920E-1,0.1922710E+
     *   00,0.8461810E+0,0.1000000E+1 /
      DATA (P(I,13),I=1,10) / 0.4602850E-2,0.3319900E-1,0.1362820E+0,
     *   0.3304760E+0,0.4491460E+0,0.2657040E+0,-.1751260E-1,0.2445330E+
     *   00,0.8049340E+0,0.1000000E+1 /
      DATA (P(I,14),I=1,10) / 0.4438260E-2,0.3266790E-1,0.1347210E+0,
     *   0.3286780E+0,0.4496400E+0,0.2613720E+0,-.1779510E-1,0.2535390E+
     *   00,0.8006690E+0,0.1000000E+1 /
      DATA (P(I,15),I=1,10) / 0.4564620E-2,0.3369360E-1,0.1397550E+0,
     *   0.3393620E+0,0.4509210E+0,0.2385860E+0,-.1776530E-1,0.2740580E+
     *   00,0.7854210E+0,0.1000000E+1 /
      DATA (P(I,16),I=1,10) / 0.4061010E-2,0.3068130E-1,0.1304520E+0,
     *   0.3272050E+0,0.4528510E+0,0.2560420E+0,-.1451050E-1,0.3102630E+
     *   00,0.7544830E+0,0.1000000E+1 /
      DATA (P(I,17),I=1,10) / 0.3989400E-2,0.3031770E-1,0.1298800E+0,
     *   0.3279510E+0,0.4535270E+0,0.2521540E+0,-.1429930E-1,0.3235720E+
     *   00,0.7435070E+0,0.1000000E+1 /
      DATA (P(I,18),I=1,10) / 0.3806650E-2,0.2923050E-1,0.1264670E+0,
     *   0.3235100E+0,0.4548960E+0,0.2566300E+0,-.1591970E-1,0.3246460E+
     *   00,0.7439900E+0,0.1000000E+1 /
C
C     DIFFUSE EXPONENTS FOR H ---> AR. (+ FOR LI-->AR, ++ FOR H,HE
C     (A_PL(HE)=0.0,A_PL(NE)=0.0,A_PL(AR)=0.0)
C
      DATA APL / 0.0360E+0,0.0000E+0,0.00740E+0,0.02070E+0,0.03150E+0,
     *   0.04380E+0,0.06390E+0,0.08450E+0,0.10760E+0,0.0000E+0,0.00760E+
     *   0,0.01460E+0,0.03180E+0,0.03310E+0,0.03480E+0,0.04050E+0,
     *   0.04830E+0,0.0000E+0 /
C
C       READ IN STAR PARAMETERS FOR H-AR PLACED IN ASTAR ARRAY
C       DLS 12/24/87
C
      DATA ASTAR(1) / 0.1100000E+1 /
      DATA ASTAR(2) / 0.1100000E+1 /
      DATA ASTAR(3) / 0.2000000E+0 /
      DATA ASTAR(4) / 0.4000000E+0 /
      DATA ASTAR(5) / 0.6000000E+0 /
      DATA ASTAR(6) / 0.8000000E+0 /
      DATA ASTAR(7) / 0.8000000E+0 /
      DATA ASTAR(8) / 0.8000000E+0 /
      DATA ASTAR(9) / 0.8000000E+0 /
      DATA ASTAR(10) / 0.8000000E+0 /
      DATA ASTAR(11) / 0.1750000E+0 /
      DATA ASTAR(12) / 0.1750000E+0 /
      DATA ASTAR(13) / 0.3250000E+0 /
      DATA ASTAR(14) / 0.4500000E+0 /
      DATA ASTAR(15) / 0.5500000E+0 /
      DATA ASTAR(16) / 0.6500000E+0 /
      DATA ASTAR(17) / 0.7500000E+0 /
      DATA ASTAR(18) / 0.8500000E+0 /
C
C       NOW INITIALIZE ALL OF THE NON-"R" DEPENDANT VALUES RATHER THAN
C       RECOMPUTING THEM MXPTS**3 TIMES IN THE Z,Y,X LOOP OVER THE
C       ORBITAL VALUE MATRIX (DENSIT)
C
C       THESE VALUES ARE INDEPENDANT OF ATOM TYPE, BEING DEPENDANT ONLY
C       ON THE ROW OF THE PERIODIC TABLE AND WHETHER IT IS "S" OR "P"
C       INITIALIZE THEM HERE AND ACCESS THEM WITHIN THE NAT LOOP, BEFORE
C       ENTERING THE LOOP OVER THE GRID (CUBE) POINTS.
C
C       CALC HAS BEEN CONVERTED TO UPPERCASE IN THE DRIVER ROUTINE
C
      NBASIS = 0
      DO 10 J = 1, NAT
         IF (IAN(J).LE.2) THEN
            NBASIS = NBASIS+2
C
C       CHECK FOR '6-31++G...'
C
            IF (INDEX(CALC,'++').NE.0) NBASIS = NBASIS+1
C
C       CATCH 6-31G**, 6-31+G**, AND 6-31++G** (P'S ON H'S)
C
            IF (INDEX(CALC,'**').NE.0.OR.INDEX(CALC,'P').NE.0) NBASIS =
     *         NBASIS+3
         ELSE
            NBASIS = NBASIS+9
            IF (IAN(J).GE.11) NBASIS = NBASIS+4
            IF (INDEX(CALC,'+').NE.0) NBASIS = NBASIS+4
            IF (INDEX(CALC,'*').NE.0.OR.INDEX(CALC,'D').NE.0) NBASIS =
     *          NBASIS+6
         ENDIF
   10 CONTINUE
      WRITE (ILST,*) ' BASIS SET IS ',CALC,' NUMBER OF AOS IS ',NBASIS
C
C       READ IN EIGENVECTORS
C
      IF (IONEMO.NE.0) THEN
         READ (IRD,20) (V(I,1),I=1,NBASIS)
      ELSE
         READ (IRD,20,END=30) ((V(I,J),I=1,NBASIS),J=1,NBASIS)
   20    FORMAT (8F10.6)
      ENDIF
   30 CONTINUE
      DO 40 J = 1, NAT
         XYZ(1,J) = C(J,1)
         XYZ(2,J) = C(J,2)
         XYZ(3,J) = C(J,3)
   40 CONTINUE
C
C       COMPUTE X,Y, AND Z VALUES FOR THE GRID (CUBE) POINTS.
C
      X(1) = XMI
      Y(1) = YMI
      Z(1) = ZMI
      DO 50 I = 2, MXPTS
         X(I) = XINC+X(I-1)
         Y(I) = YINC+Y(I-1)
         Z(I) = ZINC+Z(I-1)
   50 CONTINUE
C
C       ZERO THE ORBITAL VALUE ARRAY
C
      DO 80 IZ = 1, MXPTS
         DO 70 IY = 1, MXPTS
            DO 60 IX = 1, MXPTS
               DENSIT(IX,IY,IZ) = 0.0E+0
   60       CONTINUE
   70    CONTINUE
   80 CONTINUE
C
C       INITIALIZE THE AO COUNTER AND LOOP OVER ATOMS, IAT IS THE ATOMIC
C       NUMBER
C
C       FOR THE PRESENT, WILL ONLY PLOT THE MO SPECIFIED BY MONE
C
      MO = MONE
      KNTBAS = 1
      DO 940 I = 1, NAT
         IAT = IAN(I)
C
C       COMPUTE XDEL,YDEL,AND ZDEL (I.E. DELTA X,Y, AND Z FROM THE ATOM
C       TO EACH POINT ON THE GRID.  ONLY MXPTS VALUES FOR EACH SINCE,
C       FOR INSTANCE, EVERY POINT ON A PARTICULAR XY PLANE IS THE SAME
C       DELTA Z VALUE FROM THE POINT.  THEREFORE YOU HAVE ONLY ONE VALUE
C       FOR THE ENTIRE PLANE FOR DELTA Z, INSTEAD OF (FOR MXPTS=51)
C       2601.  AGAIN, BY COMPUTING THIS HERE, RATHER THAN INSIDE THE
C       LOOP WE CUT DOWN THESE SUBTRACTIONS AND MULTIPLICATIONS BY
C       A FACTOR OF 2601 TO 1. THIS HAS A SUBSTANTIAL EFFECT ON THE
C       SPEED OF THE COMPUTATIONS.
C
         DO 90 IXYZ = 1, MXPTS
            XDEL(IXYZ) = X(IXYZ)-XYZ(1,I)
            XDELSQ(IXYZ) = XDEL(IXYZ)*XDEL(IXYZ)
   90    CONTINUE
         DO 100 IXYZ = 1, MXPTS
            YDEL(IXYZ) = Y(IXYZ)-XYZ(2,I)
            YDELSQ(IXYZ) = YDEL(IXYZ)*YDEL(IXYZ)
  100    CONTINUE
         DO 110 IXYZ = 1, MXPTS
            ZDEL(IXYZ) = Z(IXYZ)-XYZ(3,I)
            ZDELSQ(IXYZ) = ZDEL(IXYZ)*ZDEL(IXYZ)
  110    CONTINUE
         WRITE (ILST,'(2(A,I5))') 
     *       'PROCESSING ATOM NUMBER ',I,' ATOMIC NUMBER ',IAT
C
C
C       **********************
C       ***                ***
C       ***    H TO HE     ***
C       ***                ***
C       **********************
C
C       THE EXPONENTIATIONS RELATED TO XDEL,YDEL,AND ZDEL
C       WILL BE MULTIPLIED IN THE LOOP, RATHER THAN
C       MAXPTS**3 SEPARATE EXPONENTIATIONS OVER R.  THESE
C       LM*3*MXPTS VALUES ARE THE UNIQUE ONES FOR H-HE.
C       (L+M GAUSSIANS*3 (X,Y AND Z) * MXPTS PLANES )
C       L,M REFER TO GAUSSIANS IN A K-LMG BASIS
C
         IF (IAT.LE.2) THEN
C
C     INNER 1S
C
            DO 120 IG = 1, LD
               AO1S(IG) = S(IG,IAT)*(A(IG,IAT)**0.750E+0)*0.822240980E+0
     *            *V(KNTBAS,MO)
  120       CONTINUE
C
C     OUTER - LM IS L+M FROM K-LMG,THIS WILL HELP ADDING NEW BASIS SETS
C
            AO1S(LM) = S(LM,IAT)*(A(LM,IAT)**0.750E+0)*0.822240980E+0*V(
     *         KNTBAS+1,MO)
C
C     PR0E-NEGATE EXPONENT FOR EXPONENTIAL FUNCTION
C
            DO 130 IG = 1, LM
               ANEG(IG) = -A(IG,IAT)
  130       CONTINUE
C
C       IG IS THE COUNTER OVER THE L+M GAUSSIANS
C
            DO 150 IG = 1, LM
               DO 140 IXYZ = 1, MXPTS
C
C       X, Y AND Z DEPENDANT PARTS OF THE EXPONENTIAL TERM, THE R
C       INDEPENDANT PART (AO1S) IS MULTIPLIED IN HERE, RATHER THAN
C       INSIDE THE LOOP.  NOTHING ELSE GETS MULTIPLIED BY THE EXP.
C       PART, SO NOTHING IS BEING CORRUPTED.  THE GIST OF WHAT
C       IS EQUIVALENTLY BEING DONE INSIDE THE INNER LOOP IS
C       E^(-ALPHA(IG)*R^2)*AO1S(IG) - WHERE THE EXPONENTIAL TERM
C       IS FACTORED INTO E^(-ALPHA*X^2)*E^(-ALPHA*Y^2)*E^(-ALPHA*Z^2)
C
                  E1SX(IXYZ,IG) = GEXP(XDELSQ(IXYZ)*ANEG(IG))*AO1S(IG)
                  E1SY(IXYZ,IG) = GEXP(YDELSQ(IXYZ)*ANEG(IG))
                  E1SZ(IXYZ,IG) = GEXP(ZDELSQ(IXYZ)*ANEG(IG))
  140          CONTINUE
  150       CONTINUE
            DO 190 IZ = 1, MXPTS
               DO 180 IY = 1, MXPTS
                  DO 170 IG = 1, LM
                     DO 160 IX = 1, MXPTS
C
C       CONTR IS THE SUM OF CONTRIBUTIONS OVER THIS YZ PLANE FOR THIS
C       ATOM.   WHEN FINISHED, SUM INTO THE ORBITAL VALUE ARRAY.
C
                        DENSIT(IX,IY,IZ) = DENSIT(IX,IY,IZ)+E1SY(IY,IG)*
     *                     E1SZ(IZ,IG)*E1SX(IX,IG)
  160                CONTINUE
  170             CONTINUE
  180          CONTINUE
  190       CONTINUE
            KNTBAS = KNTBAS+2
C
C       CHECK FOR ++ ANION DISFFUSE FUNCTION, DIFFUSE S ORBITAL ON H
C       USE SAME VARIABLES AS ABOVE FOR CONVENIENCE
C
            IF (INDEX(CALC,'++').NE.0) THEN
               AO1S(1) = (APL(IAT)**0.750E+0)*0.822240980E+0*V(KNTBAS,MO
     *            )
               APNEG = -APL(IAT)
               DO 200 IXYZ = 1, MXPTS
                  E1SX(IXYZ,1) = GEXP(XDELSQ(IXYZ)*APNEG)*AO1S(1)
                  E1SY(IXYZ,1) = GEXP(YDELSQ(IXYZ)*APNEG)
                  E1SZ(IXYZ,1) = GEXP(ZDELSQ(IXYZ)*APNEG)
  200          CONTINUE
C
C       LOOP OVER THE "CUBE" Z,Y, AND X
C
               DO 230 IZ = 1, MXPTS
                  DO 220 IY = 1, MXPTS
                     DO 210 IX = 1, MXPTS
C
C    SUM THE DIFFUSE ORBITAL CONTRIBUTION INTO THE ORBITAL VALUE ARRAY
C    (ONLY ONE PRIMITIVE )
C
                        DENSIT(IX,IY,IZ) = DENSIT(IX,IY,IZ)+E1SX(IX,1)*
     *                     E1SY(IY,1)*E1SZ(IZ,1)
  210                CONTINUE
  220             CONTINUE
  230          CONTINUE
               KNTBAS = KNTBAS+1
C
C       DONE WITH ++ CODE
C
            ENDIF
C
C       CODE FOR ** FUNCTIONS, P SET ON H,HE
C
            IF (INDEX(CALC,'**').NE.0.OR.INDEX(CALC,'P').NE.0) THEN
               ASNEG = -ASTAR(IAT)
               DO 240 IXYZ = 1, MXPTS
                  E2SPX(IXYZ,1) = GEXP(XDELSQ(IXYZ)*ASNEG)
                  E2SPY(IXYZ,1) = GEXP(YDELSQ(IXYZ)*ASNEG)
                  E2SPZ(IXYZ,1) = GEXP(ZDELSQ(IXYZ)*ASNEG)
  240          CONTINUE
               AO2P(1) = (ASTAR(IAT)**1.250E+0)*1.425410940E+0
               AO2PX(1) = AO2P(1)*V(KNTBAS,MO)
               AO2PY(1) = AO2P(1)*V(KNTBAS+1,MO)
               AO2PZ(1) = AO2P(1)*V(KNTBAS+2,MO)
               DO 250 IXYZ = 1, MXPTS
                  CNS2PX(IXYZ,1) = AO2PX(1)*XDEL(IXYZ)
                  CNS2PY(IXYZ,1) = AO2PY(1)*YDEL(IXYZ)
                  CNS2PZ(IXYZ,1) = AO2PZ(1)*ZDEL(IXYZ)
  250          CONTINUE
C
C       LOOP OVER THE "CUBE" Z,Y, AND X
C
               DO 280 IZ = 1, MXPTS
                  DO 270 IY = 1, MXPTS
                     DO 260 IX = 1, MXPTS
C
C       SUM THE DIFFUSE ORBITAL CONTRIBUTIONS INTO THE ORBITAL VALUE ARR
C       (ONLY ONE PRIMITIVE )
C
                        DENSIT(IX,IY,IZ) = DENSIT(IX,IY,IZ)+(CNS2PY(IY,1
     *                     )+CNS2PZ(IZ,1)+CNS2PX(IX,1))*E2SPX(IX,1)*
     *                     E2SPY(IY,1)*E2SPZ(IZ,1)
  260                CONTINUE
  270             CONTINUE
  280          CONTINUE
               KNTBAS = KNTBAS+3
C
C       DONE WITH ** CODE
C
            ENDIF
C
C       **********************
C       ***                ***
C       ***    LI TO NE    ***
C       ***                ***
C       **********************
C
         ELSEIF (IAT.LE.10) THEN
C
C       AO1S(1-KD) ARE THE KD "CONSTANTS" FOR THE 1S
C       ORBITAL, EVERYTHING THAT IS INDEPENDANT OF R IS IN THERE.
C       (1 PER GAUSSIAN PRIMITIVE)
C
            DO 290 IG = 1, KD
               AO1S(IG) = S(IG,IAT)*(A(IG,IAT)**0.750E+0)*0.712705470E+0
     *            *V(KNTBAS,MO)
  290       CONTINUE
C
C       INNER S FUNCTION
C
            DO 300 IG = 1, LD
               AO2S(IG) = S(KD+IG,IAT)*(A(KD+IG,IAT)**0.750E+0)*
     *            0.712705470E+0*V(KNTBAS+1,MO)
  300       CONTINUE
C
C       OUTER S FUNCTION
C
            AO2S(LM) = S(KD+LM,IAT)*(A(KD+LM,IAT)**0.750E+0)*
     *         0.712705470E+0*V(KNTBAS+5,MO)
C
C       INNER P FUNCTION
C
            DO 310 IG = 1, LD
               AO2P(IG) = P(IG,IAT)*(A(KD+IG,IAT)**1.250E+0)*
     *            1.425410940E+0
  310       CONTINUE
C
C       OUTER P FUNCTION
C
            AO2P(LM) = P(LM,IAT)*(A(KD+LM,IAT)**1.250E+0)*1.425410940E+0
C
C       INNER P * COEFFICIENT (PX)
C
            DO 320 IG = 1, LD
               AO2PX(IG) = AO2P(IG)*V(KNTBAS+2,MO)
  320       CONTINUE
C
C       OUTER P * COEFFICIENT (PX)
C
            AO2PX(LM) = AO2P(LM)*V(KNTBAS+6,MO)
C
C       INNER P * COEFFICIENT (PY)
C
            DO 330 IG = 1, LD
               AO2PY(IG) = AO2P(IG)*V(KNTBAS+3,MO)
  330       CONTINUE
C
C       OUTER P * COEFFICIENT (PY)
C
            AO2PY(LM) = AO2P(LM)*V(KNTBAS+7,MO)
C
C       INNER P * COEFFICIENT (PZ)
C
            DO 340 IG = 1, LD
               AO2PZ(IG) = AO2P(IG)*V(KNTBAS+4,MO)
  340       CONTINUE
C
C       OUTER P * COEFFICIENT (PZ)
C
            AO2PZ(LM) = AO2P(LM)*V(KNTBAS+8,MO)
C
C       AO2S(1-LM) CORRESPONDS TO THE LM GAUSSIANS FOR THE 2S ORBITAL
C       AO2PX(1-LM) "                            " FOR THE 2PX ""
C       ETC.
C       YOU NEED DIFFERENT "CONSTANTS" FOR EACH DUE TO THE DIFFERENT
C       VALUE OF THE WAVEFUNCTION FOR EACH ORBITAL. AT THE END WE WILL
C       HAVE A 3 SETS OF NUMBERS ALL MULTIPLIED BY THE EXPONENTIAL,
C       WE CAN ADD THEM FIRST, THEN CALC. AND MULT. BY THE EXP.
C
            DO 360 IG = 1, LM
               DO 350 IXYZ = 1, MXPTS
C
C       X - DEPENDANT PART ( 2PX )
C
                  CNS2PX(IXYZ,IG) = AO2PX(IG)*XDEL(IXYZ)
C
C       Y - DEPENDANT PART ( 2PY )
C
                  CNS2PY(IXYZ,IG) = AO2PY(IG)*YDEL(IXYZ)
C
C       AO2PZ(7->10) ARE THE 2 S MULTIPLIERS.  THEY CAN BE ADDED TO ANY
C       ONE OF THE 3 MULTIPLIERS OUTSIDE THE LOOP, SINCE THEY HAVE NO
C       X,Y, OR Z DEPENDANCE. THESE 4 ALL HAVE THE SAME EXP TO MULTIPLY
C       BY AND BE ADDED TO THE DENSITY MATRIX, SO THEY CAN BE ADDED
C       BEFORE, AND THEN MULTIPLIED BY THE EXP ALL AT ONCE.
C
                  CNS2PZ(IXYZ,IG) = AO2PZ(IG)*ZDEL(IXYZ)+AO2S(IG)
  350          CONTINUE
  360       CONTINUE
C
C       EVALUATE 2S, 2P, AND 1S
C
C       MINUS ALPHA FOR 1S:
C
            DO 370 IG = 1, KD
               ANEG(IG) = -A(IG,IAT)
  370       CONTINUE
C
C       MINUS ALPHA FOR 2SP:
C
            DO 380 IG = KD+1, KD+LM
               ANEG(IG) = -A(IG,IAT)
  380       CONTINUE
C
C       PRECOMPUTE EXP(-A*DELO**2) WHERE DELO**2 IS
C       DELTA-X**2, DELTA-Y**2, AND DELTA-Z**2
C       THIS WILL BE K*3(X,Y, AND Z)*MXPTS FOR K-LMG CALCULATIONS
C       ON THE 1S CORE.
C
            DO 400 IG = 1, KD
               DO 390 IXYZ = 1, MXPTS
                  E1SX(IXYZ,IG) = GEXP(XDELSQ(IXYZ)*ANEG(IG))*AO1S(IG)
                  E1SY(IXYZ,IG) = GEXP(YDELSQ(IXYZ)*ANEG(IG))
                  E1SZ(IXYZ,IG) = GEXP(ZDELSQ(IXYZ)*ANEG(IG))
  390          CONTINUE
  400       CONTINUE
C
C       NOW THE VALENCE EXPONENTS, THIS WILL BE (L+M)*3(X,Y, AND Z)*
C       MXPTS COMPUTATIONS FOR K-LMG CALCULATIONS. THIS SAVES DOING
C       (L+M)*1(R)*MXPTS**3 COMPUTATIONS INSIDE THE LOOP.  THIS CODE
C       WAS REPONSIBLE FOR A TEST STO-3G CASE (F2) GOING FROM 1:45
C       TO 0:13 (MIN:SEC).
C
            DO 420 IG = 1, LM
               DO 410 IXYZ = 1, MXPTS
C
C       X^2 PART OF EXPONENTIAL
C
                  E2SPX(IXYZ,IG) = GEXP(XDELSQ(IXYZ)*ANEG(KD+IG))
C
C       Y^2 PART OF EXPONENTIAL
C
                  E2SPY(IXYZ,IG) = GEXP(YDELSQ(IXYZ)*ANEG(KD+IG))
C
C       Z^2 PART OF EXPONENTIAL
C
                  E2SPZ(IXYZ,IG) = GEXP(ZDELSQ(IXYZ)*ANEG(KD+IG))
  410          CONTINUE
  420       CONTINUE
C
C       LOOP OVER THE "CUBE" Z,Y, AND X
C
            DO 480 IZ = 1, MXPTS
C
C       NOW, LOOP OVER X, INCLUDE X-DEPENDANT CONTRIBUTIONS, COMPUTE THE
C       ORBITAL VALUE AND SUM INTO THE DENSITY ARRAY.
C
               DO 470 IY = 1, MXPTS
                  DO 440 IG = 1, KD
                     DO 430 IX = 1, MXPTS
C
C       COMPUTE AND SUM THE ORBITAL CONTRIBUTIONS INTO THE ORBITAL
C       VALUE ARRAY
C
C       FIRST THE KD GAUSSIANS FOR THE 1S:
C
                        DENSIT(IX,IY,IZ) = DENSIT(IX,IY,IZ)+E1SY(IY,IG)*
     *                     E1SZ(IZ,IG)*E1SX(IX,IG)
  430                CONTINUE
  440             CONTINUE
C
C       NEXT, THE LM (L-INNER-M-OUTER) FOR THE 2SP:
C
                  DO 460 IG = 1, LM
C
C       LOOP OVER X FOR THE 2SP SHELL
C
C       SUM THE NON-EXPONENTIAL PARTS OF THE 2SP SHELL.  THE 2S HAS
C       ALREADY BEEN SUMMED IN, AS IT HAS NO DEPENDANCE ON X,Y, OR Z.
C
C       MULT EXP(Z**2)*EXP(Y**2)*EXP(Z**2)
C
                     DO 450 IX = 1, MXPTS
                        DENSIT(IX,IY,IZ) = DENSIT(IX,IY,IZ)+(CNS2PY(IY,
     *                     IG)+CNS2PZ(IZ,IG)+CNS2PX(IX,IG))*E2SPX(IX,IG)
     *                     *E2SPY(IY,IG)*E2SPZ(IZ,IG)
  450                CONTINUE
  460             CONTINUE
  470          CONTINUE
  480       CONTINUE
            KNTBAS = KNTBAS+9
C
C       THIS WILL BE DONE FOR BOTH 6-31+G AND 6-31++G, AS WELL AS EITHER
C       CASE SUPPLEMENTED WITH '*' , I.E. 6-31+G*, 6-31++G**, ETC....
C
            IF (INDEX(CALC,'+').NE.0) THEN
               APNEG = -APL(IAT)
               DO 490 IXYZ = 1, MXPTS
                  E2SPX(IXYZ,1) = GEXP(XDELSQ(IXYZ)*APNEG)
                  E2SPY(IXYZ,1) = GEXP(YDELSQ(IXYZ)*APNEG)
                  E2SPZ(IXYZ,1) = GEXP(ZDELSQ(IXYZ)*APNEG)
  490          CONTINUE
               AO2S(1) = V(KNTBAS,MO)*(APL(IAT)**0.750E+0)*0.712705470E+
     *            0
               AO2P(1) = (APL(IAT)**1.250E+0)*1.425410940E+0
               AO2PX(1) = AO2P(1)*V(KNTBAS+1,MO)
               AO2PY(1) = AO2P(1)*V(KNTBAS+2,MO)
               AO2PZ(1) = AO2P(1)*V(KNTBAS+3,MO)
               DO 500 IXYZ = 1, MXPTS
                  CNS2PX(IXYZ,1) = AO2PX(1)*XDEL(IXYZ)
                  CNS2PY(IXYZ,1) = AO2PY(1)*YDEL(IXYZ)
C
C       AO2S(1->4) ARE THE 2 S MULTIPLIERS.  THEY CAN BE ADDED TO ANY
C       ONE OF THE 3 MULTIPLIERS OUTSIDE THE LOOP, SINCE THEY HAVE NO
C       X,Y, OR Z DEPENDANCE. THESE 4 ALL HAVE THE SAME EXP TO MULTIPLY
C       BY AND BE ADDED TO THE DENSITY MATRIX, SO THEY CAN BE ADDED
C       BEFORE, AND THEN MULTIPLIED BY THE EXP ALL AT ONCE.
C
                  CNS2PZ(IXYZ,1) = AO2PZ(1)*ZDEL(IXYZ)+AO2S(1)
  500          CONTINUE
C
C       LOOP OVER THE "CUBE" Z,Y, AND X
C
               DO 530 IZ = 1, MXPTS
                  DO 520 IY = 1, MXPTS
                     DO 510 IX = 1, MXPTS
C
C       SUM THE DIFFUSE ORBITAL CONTRIBUTIONS INTO THE ORBITAL VALUE ARR
C       (ONLY ONE PRIMITIVE )
C
                        DENSIT(IX,IY,IZ) = DENSIT(IX,IY,IZ)+(CNS2PY(IY,1
     *                     )+CNS2PZ(IZ,1)+CNS2PX(IX,1))*E2SPX(IX,1)*
     *                     E2SPY(IY,1)*E2SPZ(IZ,1)
  510                CONTINUE
  520             CONTINUE
  530          CONTINUE
               KNTBAS = KNTBAS+4
            ENDIF
C
C       THIS WILL BE DONE FOR EITHER *, OR ** WAVEFUNCTIONS, WITH ANY
C       COMBINATION OF (OR LACK OF) + FUNCTIONS
C
            IF (INDEX(CALC,'*').NE.0.OR.INDEX(CALC,'D').NE.0) THEN
               ASNEG = -ASTAR(IAT)
               DO 540 IXYZ = 1, MXPTS
                  EXPDX(IXYZ) = GEXP(XDELSQ(IXYZ)*ASNEG)
                  EXPDY(IXYZ) = GEXP(YDELSQ(IXYZ)*ASNEG)
                  EXPDZ(IXYZ) = GEXP(ZDELSQ(IXYZ)*ASNEG)
  540          CONTINUE
               TMP = (ASTAR(IAT)**1.750E+0)*1.645922780E+0
               AODXX = V(KNTBAS,MO)*TMP
               AODYY = V(KNTBAS+1,MO)*TMP
               AODZZ = V(KNTBAS+2,MO)*TMP
               AODXY = V(KNTBAS+3,MO)*TMP
               AODXZ = V(KNTBAS+4,MO)*TMP
               AODYZ = V(KNTBAS+5,MO)*TMP
C
C       CNSXY AND CNSXZ WILL BE ADDED OUTSIDE THE INNER LOOP, THEN
C       THE SUM WILL BE MULTIPLIED BY XDEL(IX) IN THE LOOP.  THAT SHOULD
C       SAVE ONE ADDITION IN THE INNER LOOP.
C
               DO 550 IXYZ = 1, MXPTS
                  CNSXX(IXYZ) = AODXX*XDELSQ(IXYZ)
                  CNSYY(IXYZ) = AODYY*YDELSQ(IXYZ)
                  CNSZZ(IXYZ) = AODZZ*ZDELSQ(IXYZ)
                  CNSXY(IXYZ) = AODXY*YDEL(IXYZ)
                  CNSXZ(IXYZ) = AODXZ*ZDEL(IXYZ)
                  CNSYZ(IXYZ) = AODYZ*YDEL(IXYZ)
  550          CONTINUE
C
C       LOOP OVER THE "CUBE" Z,Y, AND X
C
               DO 580 IZ = 1, MXPTS
                  DO 570 IY = 1, MXPTS
                     DO 560 IX = 1, MXPTS
C
C       SUM THE DIFFUSE ORBITAL CONTRIBUTIONS INTO THE ORBITAL VALUE ARR
C       (ONLY ONE PRIMITIVE )
C
                        DENSIT(IX,IY,IZ) = DENSIT(IX,IY,IZ)+(CNSZZ(IZ)+
     *                     CNSYZ(IY)*ZDEL(IZ)+CNSYY(IY)+CNSXX(IX)+
     *                     XDEL(IX)*CNSXZ(IZ)+CNSXY(IY))*EXPDX(IX)*EXPDY
     *                     (IY)*EXPDZ(IZ)
  560                CONTINUE
  570             CONTINUE
  580          CONTINUE
               KNTBAS = KNTBAS+6
            ENDIF
C
C       **********************
C       ***                ***
C       ***    NA TO AR    ***
C       ***                ***
C       **********************
C
         ELSEIF (IAT.LE.18) THEN
            DO 590 IG = 1, KD
               IG2 = KD+IG
C
C       AO1S(1-KD) ARE THE KD GAUSSION "CONSTANTS" FOR THE 1S
C       ORBITAL, EVERYTHING THAT IS INDEPENDANT OF R IS IN THERE.
C
               AO1S(IG) = S(IG,IAT)*(A(IG,IAT)**0.750E+0)*0.712705470E+0
     *            *V(KNTBAS,MO)
C
C       CORE 2SP SHELL.... 2S:
C
               AO2S(IG) = S(IG2,IAT)*(A(IG2,IAT)**0.750E+0)*0.712705470E
     *            +0*V(KNTBAS+1,MO)
C
C       CORE 2SP SHELL.... 2P (GENERAL PARTS):
C
               AO2P(IG) = P(IG,IAT)*(A(IG2,IAT)**1.250E+0)*1.425410940E+
     *            0
C
C       CORE 2SP SHELL.... 2PX:
C
               AO2PX(IG) = AO2P(IG)*V(KNTBAS+2,MO)
C
C       CORE 2SP SHELL.... 2PY:
C
               AO2PY(IG) = AO2P(IG)*V(KNTBAS+3,MO)
C
C       CORE 2SP SHELL.... 2PZ:
C
               AO2PZ(IG) = AO2P(IG)*V(KNTBAS+4,MO)
  590       CONTINUE
C
C       INNER S FUNCTION
C
            DO 600 IG = 1, LD
               AO3S(IG) = S(K2+IG,IAT)*(A(K2+IG,IAT)**0.750E+0)*
     *            0.712705470E+0*V(KNTBAS+5,MO)
  600       CONTINUE
C
C       OUTER S FUNCTION
C
            AO3S(LM) = S(K2+LM,IAT)*(A(K2+LM,IAT)**0.750E+0)*
     *         0.712705470E+0*V(KNTBAS+9,MO)
C
C       INNER P FUNCTION
C
            DO 610 IG = 1, LD
               AO3P(IG) = P(KD+IG,IAT)*(A(K2+IG,IAT)**1.250E+0)*
     *            1.425410940E+0
  610       CONTINUE
C
C       OUTER P FUNCTION
C
            AO3P(LM) = P(KD+LM,IAT)*(A(K2+LM,IAT)**1.250E+0)*
     *         1.425410940E+0
C
C       INNER P * COEFFICIENT (PX)
C
            DO 620 IG = 1, LD
               AO3PX(IG) = AO3P(IG)*V(KNTBAS+6,MO)
  620       CONTINUE
C
C       OUTER P * COEFFICIENT (PX)
C
            AO3PX(LM) = AO3P(LM)*V(KNTBAS+10,MO)
C
C       INNER P * COEFFICIENT (PY)
C
            DO 630 IG = 1, LD
               AO3PY(IG) = AO3P(IG)*V(KNTBAS+7,MO)
  630       CONTINUE
C
C       OUTER P * COEFFICIENT (PY)
C
            AO3PY(LM) = AO3P(LM)*V(KNTBAS+11,MO)
C
C       INNER P * COEFFICIENT (PZ)
C
            DO 640 IG = 1, LD
               AO3PZ(IG) = AO3P(IG)*V(KNTBAS+8,MO)
  640       CONTINUE
C
C       OUTER P * COEFFICIENT (PZ)
C
            AO3PZ(LM) = AO3P(LM)*V(KNTBAS+12,MO)
C
C       CORE 2P X,Y,AND Z DEPENDANT PARTS
C
            DO 660 IG = 1, KD
               DO 650 IXYZ = 1, MXPTS
C
C       CORE 2PX - X DEPENDANT PART
C
                  CNS2PX(IXYZ,IG) = AO2PX(IG)*XDEL(IXYZ)
C
C       CORE 2PY - Y DEPENDANT PART
C
                  CNS2PY(IXYZ,IG) = AO2PY(IG)*YDEL(IXYZ)
C
C       AO2S(1->6) ARE THE 2S MULTIPLIERS.  THEY CAN BE ADDED TO ANY
C       ONE OF THE 3 MULTIPLIERS OUTSIDE THE LOOP, SINCE THEY HAVE NO
C       X,Y, OR Z DEPENDANCE. THESE 4 ALL HAVE THE SAME EXP TO MULTIPLY
C       BY AND BE ADDED TO THE DENSITY MATRIX, SO THEY CAN BE ADDED
C       BEFORE, AND THEN MULTIPLIED BY THE EXP ALL AT ONCE.
C
C       CORE 2PZ - Z DEPENDANT PART - S ADDED IN HERE
C
                  CNS2PZ(IXYZ,IG) = AO2PZ(IG)*ZDEL(IXYZ)+AO2S(IG)
  650          CONTINUE
  660       CONTINUE
C
C       AO3S  CORRESPONDS TO THE 4 GAUSSIANS FOR THE 3S ORBITAL
C       AO3PX "                            " FOR THE 3PX ""
C       ETC.
C       YOU NEED DIFFERENT "CONSTANTS" FOR EACH DUE TO THE DIFFERENT
C       VALUE OF THE WAVEFUNCTION FOR EACH ORBITAL. AT THE END WE WILL
C       HAVE A 3 SETS OF NUMBERS ALL MULTIPLIED BY THE EXPONENTIAL,
C       WE CAN ADD THEM FIRST, THEN CALC. AND MULT. BY THE EXP.
C
            DO 680 IG = 1, LM
               DO 670 IXYZ = 1, MXPTS
                  CNS3PX(IXYZ,IG) = AO3PX(IG)*XDEL(IXYZ)
                  CNS3PY(IXYZ,IG) = AO3PY(IG)*YDEL(IXYZ)
C
C       AO3S(1->4) ARE THE 3S MULTIPLIERS.  THEY CAN BE ADDED TO ANY
C       ONE OF THE 3 MULTIPLIERS OUTSIDE THE LOOP, SINCE THEY HAVE NO
C       X,Y, OR Z DEPENDANCE. THESE 4 ALL HAVE THE SAME EXP TO MULTIPLY
C       BY AND BE ADDED TO THE DENSITY MATRIX, SO THEY CAN BE ADDED
C       BEFORE, AND THEN MULTIPLIED BY THE EXP ALL AT ONCE.
C
                  CNS3PZ(IXYZ,IG) = AO3PZ(IG)*ZDEL(IXYZ)+AO3S(IG)
  670          CONTINUE
  680       CONTINUE
C
C       EVALUATE 1S, 2SP, AND 3SP EXPONENTIALS
C
C       MINUS ALPHA FOR 1S:
C
            DO 690 IG = 1, KD
               ANEG(IG) = -A(IG,IAT)
  690       CONTINUE
C
C       MINUS ALPHA FOR 2SP:
C
            DO 700 IG = 1, KD
               ANEG(KD+IG) = -A(KD+IG,IAT)
  700       CONTINUE
C
C       MINUS ALPHA FOR 3SP:
C
            DO 710 IG = 1, LM
               ANEG(K2+IG) = -A(K2+IG,IAT)
  710       CONTINUE
C
C       PRECOMPUTE EXP(-A*DELO**2) WHERE DELO**2 IS
C       DELTA-X**2, DELTA-Y**2, AND DELTA-Z**2
C       THIS WILL BE K*3(X,Y, AND Z)*MXPTS FOR K-LMG CALCULATIONS
C       ON THE 1S CORE.
C
C       THE AO CONTRIBUTION CAN BE MULTIPLIED IN RIGHT HERE, SINCE
C       THERE IS NO X,Y, OR Z DEPENDANCE ON IT, AND NOTHING ELSE
C       NEEDS TO BE MULTIPLIED INTO THE EXPONENTIAL TERM.
C
            DO 730 IG = 1, KD
               DO 720 IXYZ = 1, MXPTS
C
C       E^(-ALPHA(1S)*X^2)
C
                  E1SX(IXYZ,IG) = GEXP(XDELSQ(IXYZ)*ANEG(IG))*AO1S(IG)
C
C       E^(-ALPHA(1S)*Y^2)
C
                  E1SY(IXYZ,IG) = GEXP(YDELSQ(IXYZ)*ANEG(IG))
C
C       E^(-ALPHA(1S)*Z^2)
C
                  E1SZ(IXYZ,IG) = GEXP(ZDELSQ(IXYZ)*ANEG(IG))
  720          CONTINUE
  730       CONTINUE
            DO 750 IG = 1, KD
               DO 740 IXYZ = 1, MXPTS
C
C       E^(-ALPHA(2SP)*X^2)
C
                  E2SPX(IXYZ,IG) = GEXP(XDELSQ(IXYZ)*ANEG(KD+IG))
C
C       E^(-ALPHA(2SP)*Y^2)
C
                  E2SPY(IXYZ,IG) = GEXP(YDELSQ(IXYZ)*ANEG(KD+IG))
C
C       E^(-ALPHA(2SP)*Z^2)
C
                  E2SPZ(IXYZ,IG) = GEXP(ZDELSQ(IXYZ)*ANEG(KD+IG))
  740          CONTINUE
  750       CONTINUE
C
C       NOW THE VALENCE EXPONENTS, THIS WILL BE (L+M)*3(X,Y, AND Z)*
C       MXPTS COMPUTATIONS FOR K-LMG CALCULATIONS. THIS SAVES DOING
C       (L+M)*1(R)*MXPTS**3 COMPUTATIONS INSIDE THE LOOP.  THIS CODE
C       WAS REPONSIBLE FOR A TEST STO-3G CASE (F2) GOING FROM 1:45
C       TO 0:13 (MIN:SEC).
C
            DO 770 IG = 1, LM
               DO 760 IXYZ = 1, MXPTS
C
C       E^(-ALPHA(3SP)*X^2)
C
                  E3SPX(IXYZ,IG) = GEXP(XDELSQ(IXYZ)*ANEG(K2+IG))
C
C       E^(-ALPHA(3SP)*Y^2)
C
                  E3SPY(IXYZ,IG) = GEXP(YDELSQ(IXYZ)*ANEG(K2+IG))
C
C       E^(-ALPHA(3SP)*Z^2)
C
                  E3SPZ(IXYZ,IG) = GEXP(ZDELSQ(IXYZ)*ANEG(K2+IG))
  760          CONTINUE
  770       CONTINUE
C
C
C       LOOP OVER THE "CUBE" Z,Y, AND X
C
            DO 830 IZ = 1, MXPTS
C
C       NOW, LOOP OVER X, INCLUDE X-DEPENDANT CONTRIBUTIONS, COMPUTE THE
C       ORBITAL VALUE AND SUM INTO THE DENSITY ARRAY.
C
               DO 820 IY = 1, MXPTS
                  DO 790 IG = 1, KD
C
C       SUM THE NON-EXPONENTIAL PARTS OF THE 2SP SHELL.  THE 2S HAS
C       ALREADY BEEN SUMMED IN, AS IT HAS NO DEPENDANCE ON X,Y, OR Z.
C
                     DO 780 IX = 1, MXPTS
C
C       COMPUTE AND SUM THE ORBITAL CONTRIBUTIONS INTO THE ORBITAL
C       VALUE ARRAY
C
C       FIRST THE KD GAUSSIANS FOR THE 1S:
C
                        DENSIT(IX,IY,IZ) = DENSIT(IX,IY,IZ)+E1SY(IY,IG)*
     *                     E1SZ(IZ,IG)*E1SX(IX,IG)
C
C       NOW THE KD FOR THE CORE 2SP SHELL
C
                        DENSIT(IX,IY,IZ) = DENSIT(IX,IY,IZ)+(CNS2PY(IY,
     *                     IG)+CNS2PZ(IZ,IG)+CNS2PX(IX,IG))*E2SPX(IX,IG)
     *                     *E2SPY(IY,IG)*E2SPZ(IZ,IG)
  780                CONTINUE
  790             CONTINUE
C
C       NOW THE LM (L-INNER-M-OUTER) FOR THE VALENCE 3SP SHELL
C
                  DO 810 IG = 1, LM
C
C       NOW DO THE SAME FOR THE 3SP
C
                     DO 800 IX = 1, MXPTS
                        DENSIT(IX,IY,IZ) = DENSIT(IX,IY,IZ)+(CNS3PY(IY,
     *                     IG)+CNS3PZ(IZ,IG)+CNS3PX(IX,IG))*E3SPX(IX,IG)
     *                     *E3SPY(IY,IG)*E3SPZ(IZ,IG)
  800                CONTINUE
  810             CONTINUE
  820          CONTINUE
  830       CONTINUE
            KNTBAS = KNTBAS+13
C
C       THIS WILL BE DONE FOR BOTH 6-31+G AND 6-31++G, AS WELL AS EITHER
C       CASE SUPPLEMENTED WITH '*' , I.E. 6-31+G*, 6-31++G**, ETC....
C
            IF (INDEX(CALC,'+').NE.0) THEN
               APNEG = -APL(IAT)
               DO 840 IXYZ = 1, MXPTS
                  E3SPX(IXYZ,1) = GEXP(XDELSQ(IXYZ)*APNEG)
                  E3SPY(IXYZ,1) = GEXP(YDELSQ(IXYZ)*APNEG)
                  E3SPZ(IXYZ,1) = GEXP(ZDELSQ(IXYZ)*APNEG)
  840          CONTINUE
C
C       USE SAME VARIABLES AS USED ABOVE, JUST FOR CONVENIENCE
C       THEY DON'T NEED ARRAYS NOW, BUT SOMEONE MAY WANT A + SHELL
C       OR D SHELL WITH MORE THAN ONE PRIMITIVE
C
               AO3S(1) = V(KNTBAS,MO)*(APL(IAT)**0.750E+0)*0.712705470E+
     *            0
               AO3P(1) = (APL(IAT)**1.250E+0)*1.425410940E+0
               AO3PX(1) = AO3P(1)*V(KNTBAS+1,MO)
               AO3PY(1) = AO3P(1)*V(KNTBAS+2,MO)
               AO3PZ(1) = AO3P(1)*V(KNTBAS+3,MO)
               DO 850 IXYZ = 1, MXPTS
                  CNS3PX(IXYZ,1) = AO3PX(1)*XDEL(IXYZ)
                  CNS3PY(IXYZ,1) = AO3PY(1)*YDEL(IXYZ)
C
C       AO3S(1) IS THE S MULTIPLIER.  IT CAN BE ADDED TO ANY
C       ONE OF THE 3 MULTIPLIERS OUTSIDE THE LOOP, SINCE THEY HAVE NO
C       X,Y, OR Z DEPENDANCE. THESE 4 ALL HAVE THE SAME EXP TO MULTIPLY
C       BY AND BE ADDED TO THE DENSITY MATRIX, SO THEY CAN BE ADDED
C       BEFORE, AND THEN MULTIPLIED BY THE EXP ALL AT ONCE.
C
                  CNS3PZ(IXYZ,1) = AO3PZ(1)*ZDEL(IXYZ)+AO3S(1)
  850          CONTINUE
C
C       LOOP OVER THE "CUBE" Z,Y, AND X
C
               DO 880 IZ = 1, MXPTS
                  DO 870 IY = 1, MXPTS
                     DO 860 IX = 1, MXPTS
C
C       SUM THE DIFFUSE ORBITAL CONTRIBUTIONS INTO THE ORBITAL VALUE ARR
C       (ONLY ONE PRIMITIVE )
C
                        DENSIT(IX,IY,IZ) = DENSIT(IX,IY,IZ)+(CNS3PY(IY,1
     *                     )+CNS3PZ(IZ,1)+CNS3PX(IX,1))*E3SPX(IX,1)*
     *                     E3SPY(IY,1)*E3SPZ(IZ,1)
  860                CONTINUE
  870             CONTINUE
  880          CONTINUE
               KNTBAS = KNTBAS+4
            ENDIF
C
C       THIS WILL BE DONE FOR EITHER *, OR ** WAVEFUNCTIONS, WITH ANY
C       COMBINATION OF (OR LACK OF) + FUNCTIONS
C
            IF (INDEX(CALC,'*').NE.0.OR.INDEX(CALC,'D').NE.0) THEN
               ASNEG = -ASTAR(IAT)
               DO 890 IXYZ = 1, MXPTS
                  EXPDX(IXYZ) = GEXP(XDELSQ(IXYZ)*ASNEG)
                  EXPDY(IXYZ) = GEXP(YDELSQ(IXYZ)*ASNEG)
                  EXPDZ(IXYZ) = GEXP(ZDELSQ(IXYZ)*ASNEG)
  890          CONTINUE
               TMP = (ASTAR(IAT)**1.750E+0)*1.645922780E+0
               AODXX = V(KNTBAS,MO)*TMP
               AODYY = V(KNTBAS+1,MO)*TMP
               AODZZ = V(KNTBAS+2,MO)*TMP
               AODXY = V(KNTBAS+3,MO)*TMP
               AODXZ = V(KNTBAS+4,MO)*TMP
               AODYZ = V(KNTBAS+5,MO)*TMP
C
C       CNSXY AND CNSXZ WILL BE ADDED OUTSIDE THE INNER LOOP, THEN
C       THE SUM WILL BE MULTIPLIED BY XDEL(IX) IN THE LOOP.  THAT SHOULD
C       SAVE ONE ADDITION IN THE INNER LOOP.
C
               DO 900 IXYZ = 1, MXPTS
                  CNSXX(IXYZ) = AODXX*XDELSQ(IXYZ)
                  CNSYY(IXYZ) = AODYY*YDELSQ(IXYZ)
                  CNSZZ(IXYZ) = AODZZ*ZDELSQ(IXYZ)
                  CNSXY(IXYZ) = AODXY*YDEL(IXYZ)
                  CNSXZ(IXYZ) = AODXZ*ZDEL(IXYZ)
                  CNSYZ(IXYZ) = AODYZ*YDEL(IXYZ)
  900          CONTINUE
C
C       LOOP OVER THE "CUBE" Z,Y, AND X
C
               DO 930 IZ = 1, MXPTS
C
C       THE Z^2,YZ, AND Y^2 ARE ALL "CONSTANT" WITHING THE X LOOP, THEY
C       CAN BE ADDED NOW
C
                  DO 920 IY = 1, MXPTS
                     DO 910 IX = 1, MXPTS
C
C       SUM THE DIFFUSE ORBITAL CONTRIBUTIONS INTO THE ORBITAL VALUE ARR
C       (ONLY ONE PRIMITIVE )
C
                        DENSIT(IX,IY,IZ) = DENSIT(IX,IY,IZ)+(CNSZZ(IZ)+
     *                     CNSYZ(IY)*ZDEL(IZ)+CNSYY(IY)+CNSXX(IX)+
     *                     XDEL(IX)*CNSXZ(IZ)+CNSXY(IY))*EXPDX(IX)*EXPDY
     *                     (IY)*EXPDZ(IZ)
  910                CONTINUE
  920             CONTINUE
  930          CONTINUE
               KNTBAS = KNTBAS+6
            ENDIF
C
         ENDIF
  940 CONTINUE
      RETURN
      END
C
C
      SUBROUTINE MO321G
C
      IMPLICIT REAL (A-H,O-Z)
C
      PARAMETER (MXPTS=51)
      PARAMETER (KD=3,LD=2,MD=1,K3=KD*3,K2=KD*2)
      PARAMETER (MAXATM=50)
      PARAMETER (MAXORB=200)
C
C       WHEN GENERATING A NEW WAVEFUNCTION ROUTINE, DEFINE THE
C       PARAMETERS KDIM,LDIM,AND MDIM FOR THE K-LMG WAVEFUNCTION
C       BEING USED - 3-21G WILL BE 6,3, AND 1 RESPECTIVELY.  REPLACE
C       THE DATA STATEMENTS FOR A,S,P, AND APLUS AS NECESSARY,
C       AND CHANGE ALL OCCURANCES OF 3-21 TO WHATEVER BASIS SET
C       YOU ARE DEFINING.   YOU WILL ALSO WANT TO REMOVE ANY CODE
C       NOT APPLICABLE TO THE BASIS (** FOR 3-21, FOR EXAMPLE), OR
C       NOT IMPLIMENTED AT THE PRESENT TIME.
C
C       DAN SEVERANCE 12/24/87
C
      PARAMETER (KLM=(KD+LD+MD),LM=(MD+LD),LM3=LM*3,LM2=LM*2)
C
C       WRITTEN AS A STAND-ALONE SUBROUTINE FOR CALCULATING ORBITAL
C       VALUES. WORKED TO REDUCE REDUNDANT COMPUTATION OF POWERS,
C       SQUARE ROOTS, AND EXPONENTIALS.  ALSO WRITTEN TO ALLOW
C       EASY VECTORIZATION BY COMPILERS ON VECTOR MACHINES.
C
C       DAN SEVERANCE 12/22/87-->1/6/88 MANY MODIFICATIONS
C
      COMMON /SPLIT/ XMI,YMI,ZMI,XINC,YINC,ZINC,MONE,MOLAST,NAT,
     *   IONEMO
      COMMON /DENS/ DENSIT(MXPTS,MXPTS,MXPTS),V(MAXORB,MAXORB),C(MAXATM,
     *   3),OCMO(MAXORB),IAN(MAXATM)
      COMMON /IO/ IRD,ILST
      DIMENSION XDEL(MXPTS),YDEL(MXPTS),ZDEL(MXPTS),AO1S(KD)
      DIMENSION XDELSQ(MXPTS),YDELSQ(MXPTS),ZDELSQ(MXPTS)
      DIMENSION X(MXPTS),Y(MXPTS),Z(MXPTS),XYZ(3,50),ANEG(50)
      DIMENSION AO2S(KD),AO2P(KD),AO2PX(KD),AO2PY(KD),AO2PZ(KD)
      DIMENSION AO3S(KD),AO3P(KD),AO3PX(KD),AO3PY(KD),AO3PZ(KD)
      DIMENSION CNS2PX(MXPTS,KD),CNS2PY(MXPTS,KD),CNS2PZ(MXPTS,KD)
      DIMENSION CNS3PX(MXPTS,KD),CNS3PY(MXPTS,KD),CNS3PZ(MXPTS,KD)
      DIMENSION E1SX(MXPTS,KD),E1SY(MXPTS,KD),E1SZ(MXPTS,KD)
      DIMENSION E2SPX(MXPTS,KD),E2SPY(MXPTS,KD),E2SPZ(MXPTS,KD)
      DIMENSION E3SPX(MXPTS,KD),E3SPY(MXPTS,KD),E3SPZ(MXPTS,KD)
      DIMENSION EXPDX(MXPTS),EXPDY(MXPTS),EXPDZ(MXPTS)
      DIMENSION CNSXX(MXPTS),CNSYY(MXPTS),CNSZZ(MXPTS)
      DIMENSION CNSYZ(MXPTS),CNSXY(MXPTS),CNSXZ(MXPTS)
C
C
C     FORM OF S,P, AND D GAUSSIAN FUNCTIONS:
C       J.COMP.CHEM.,7,359,(1986)
C
C       (^ DENOTES "RAISED TO THE POWER" I.E. R^2 IS R SQUARED)
C       S    TYPE : (   2 * ALPHA   /PI   )^(3/4)       * E^(-ALPHA*R^2)
C       PX   TYPE : ( 128 * ALPHA^5 /PI^3 )^(1/4) * X   * E^(-ALPHA*R^2)
C       DYZ  TYPE : (2048 * ALPHA^7 /PI^3 )^(1/4) * Y*Z * E^(-ALPHA*R^2)
C
C       FOR PY, REPLACE X WITH Y
C       FOR DX^2 REPLACE Y AND Z WITH X AND X, I.E. Y*Z-->X*X=X^2
C
C       THE FOLLOWING COMMENT LINES (EDITED) COURTESY AOSV:
C       WRITTEN BY JMB  NOV. 1986
C
C----------------------------------------------------------------------
C     -------------------------------------------------------
C     | FOR THE FORM OF THE SPLIT VALENCE AO FUNCTIONS SEE: |
C     -------------------------------------------------------
C     JACS,102,939,(1980)  - INNER SHELL H-NE. | ALSO CONTAINS THE FORM
C     JACS,104,2798,(1982) - INNER SHELL NA-AR | OF THE FUNCTIONS
C     J. CHEM. PHYS.,51,2657(1962) - GAUSSIANS | INCLUDING GAUSSIANS.
C     J. CHEM. PHYS.,80,3265(1984) - | GENERAL DESCRIPTION OF DIFFUSE
C                                    | FUNCTIONS  AND THEIR COEFS.
C
C     THE DATA STATEMENTS CONTAINING THE EXPONENTS, S AND P COEFFS
C     WERE READ FROM THE FILE 321.GBS IN THE G86 DISTRIBUTION.
C
C     NUMERICAL VALUES FOR THE DIFFUSE EXPONENTS WERE TAKEN FROM:
C     "AB INITIO MOLECULAR ORBITAL THEORY",HEHRE,RADOM,SCHLEYER,
C     POPLE, JOHN WILEY & SONS, 1986. PG. 87.
C
C     NORMALIZATION - FOR HYDROGEN AND HELIUM (ZETA = 1.0 FOR OTHERS)
C
C     AO1S = AO1S*(ZETA**1.50E+0)*((2.0E+0/PI)**0.750E+0)
C     ((2.0E+0/PI)**0.750E+0) = 0.712705470E+0
C     ZETA = 1.100E+0
C     ZETA*SQRT(ZETA) = 1.153689730E+0
C     1.153689730E+0*0.712705470E+0 = 0.822240980E+0
C
C----------------------------------------------------------------------
C
      COMMON /BASIS/ CALC
      CHARACTER CALC*20
C
C     THE DIMENSIONS OF THE FOLLOWING ARRAYS ARE:
C     9 IS (NROW-1)*KD+LM, WHERE H,HE IS FIRST, ETC. NROW=3
C     18 IS THE MAXIMUM ATOMIC NUMBER IMPLEMENTED
C     6 IS (NROW-2)*KD+LM
C
      DIMENSION A(9,18),S(9,18),P(6,18),APLUS(18),ASTAR(18)
      DATA (A(I,1),I=1,3) / 0.5447180E+1,0.8245470E+0,0.1831920E+0 /
      DATA (A(I,2),I=1,3) / 0.1362670E+2,0.1999350E+1,0.3829930E+0 /
      DATA (A(I,3),I=1,6) / 0.3683820E+2,0.5481720E+1,0.1113270E+1,
     *   0.5402050E+0,0.1022550E+0,0.2856450E-1 /
      DATA (A(I,4),I=1,6) / 0.7188760E+2,0.1072890E+2,0.2222050E+1,
     *   0.1295480E+1,0.2688810E+0,0.7735010E-1 /
      DATA (A(I,5),I=1,6) / 0.1164340E+3,0.1743140E+2,0.3680160E+1,
     *   0.2281870E+1,0.4652480E+0,0.1243280E+0 /
      DATA (A(I,6),I=1,6) / 0.1722560E+3,0.2591090E+2,0.5533350E+1,
     *   0.3664980E+1,0.7705450E+0,0.1958570E+0 /
      DATA (A(I,7),I=1,6) / 0.2427660E+3,0.3648510E+2,0.7814490E+1,
     *   0.5425220E+1,0.1149150E+1,0.2832050E+0 /
      DATA (A(I,8),I=1,6) / 0.3220370E+3,0.4843080E+2,0.1042060E+2,
     *   0.7402940E+1,0.1576200E+1,0.3736840E+0 /
      DATA (A(I,9),I=1,6) / 0.4138010E+3,0.6224460E+2,0.1343400E+2,
     *   0.9777590E+1,0.2086170E+1,0.4823830E+0 /
      DATA (A(I,10),I=1,6) / 0.5157240E+3,0.7765380E+2,0.1681360E+2,
     *   0.1248300E+2,0.2664510E+1,0.6062500E+0 /
      DATA (A(I,11),I=1,9) / 0.5476130E+3,0.8206780E+2,0.1769170E+2,
     *   0.1754070E+2,0.3793980E+1,0.9064410E+0,0.5018240E+0,0.6094580E-
     *   1,0.2443490E-1 /
      DATA (A(I,12),I=1,9) / 0.6528410E+3,0.9838050E+2,0.2129960E+2,
     *   0.2337270E+2,0.5199530E+1,0.1315080E+1,0.6113490E+0,0.1418410E+
     *   0,0.4640110E-1 /
      DATA (A(I,13),I=1,9) / 0.7757370E+3,0.1169520E+3,0.2533260E+2,
     *   0.2947960E+2,0.6633140E+1,0.1726750E+1,0.9461600E+0,0.2025060E+
     *   0,0.6390880E-1 /
      DATA (A(I,14),I=1,9) / 0.9106550E+3,0.1373360E+3,0.2976010E+2,
     *   0.3667160E+2,0.8317290E+1,0.2216450E+1,0.1079130E+1,0.3024220E+
     *   0,0.9333920E-1 /
      DATA (A(I,15),I=1,9) / 0.1054900E+4,0.1591950E+3,0.3453040E+2,
     *   0.4428660E+2,0.1010190E+2,0.2739970E+1,0.1218650E+1,0.3955460E+
     *   0,0.1228110E+0 /
      DATA (A(I,16),I=1,9) / 0.1210620E+4,0.1827470E+3,0.3966730E+2,
     *   0.5222360E+2,0.1196290E+2,0.3289110E+1,0.1223840E+1,0.4573030E+
     *   0,0.1422690E+0 /
      DATA (A(I,17),I=1,9) / 0.1376400E+4,0.2078570E+3,0.4515540E+2,
     *   0.6080140E+2,0.1397650E+2,0.3887100E+1,0.1352990E+1,0.5269550E+
     *   0,0.1667140E+0 /
      DATA (A(I,18),I=1,9) / 0.1553710E+4,0.2346780E+3,0.5101210E+2,
     *   0.7004530E+2,0.1614730E+2,0.4534920E+1,0.1542090E+1,0.6072670E+
     *   0,0.1953730E+0 /
      DATA (S(I,1),I=1,3) / 0.1562850E+0,0.9046910E+0,0.1000000E+1 /
      DATA (S(I,2),I=1,3) / 0.1752300E+0,0.8934830E+0,0.1000000E+1 /
      DATA (S(I,3),I=1,6) / 0.6966860E-1,0.3813460E+0,0.6817020E+0,-.
     *   2321270E+0,0.1143390E+1,0.1000000E+1 /
      DATA (S(I,4),I=1,6) / 0.6442630E-1,0.3660960E+0,0.6959340E+0,-.
     *   4210640E+0,0.1224070E+1,0.1000000E+1 /
      DATA (S(I,5),I=1,6) / 0.6296050E-1,0.3633040E+0,0.6972550E+0,-.
     *   3686620E+0,0.1199440E+1,0.1000000E+1 /
      DATA (S(I,6),I=1,6) / 0.6176690E-1,0.3587940E+0,0.7007130E+0,-.
     *   3958970E+0,0.1215840E+1,0.1000000E+1 /
      DATA (S(I,7),I=1,6) / 0.5986570E-1,0.3529550E+0,0.7065130E+0,-.
     *   4133010E+0,0.1224420E+1,0.1000000E+1 /
      DATA (S(I,8),I=1,6) / 0.5923940E-1,0.3515000E+0,0.7076580E+0,-.
     *   4044530E+0,0.1221560E+1,0.1000000E+1 /
      DATA (S(I,9),I=1,6) / 0.5854830E-1,0.3493080E+0,0.7096320E+0,-.
     *   4073270E+0,0.1223140E+1,0.1000000E+1 /
      DATA (S(I,10),I=1,6) / 0.5814300E-1,0.3479510E+0,0.7107140E+0,-.
     *   4099220E+0,0.1224310E+1,0.1000000E+1 /
      DATA (S(I,11),I=1,9) / 0.6749110E-1,0.3935050E+0,0.6656050E+0,-.
     *   1119370E+0,0.2546540E+0,0.8444170E+0,-.2196600E+0,0.1089120E+1,
     *   0.1000000E+1 /
      DATA (S(I,12),I=1,9) / 0.6759820E-1,0.3917780E+0,0.6666610E+0,-.
     *   1102460E+0,0.1841190E+0,0.8963990E+0,-.3611010E+0,0.1215050E+1,
     *   0.1000000E+1 /
      DATA (S(I,13),I=1,9) / 0.6683470E-1,0.3890610E+0,0.6694680E+0,-.
     *   1079020E+0,0.1462450E+0,0.9237300E+0,-.3203270E+0,0.1184120E+1,
     *   0.1000000E+1 /
      DATA (S(I,14),I=1,9) / 0.6608230E-1,0.3862290E+0,0.6723800E+0,-.
     *   1045110E+0,0.1074100E+0,0.9514460E+0,-.3761080E+0,0.1251650E+1,
     *   0.1000000E+1 /
      DATA (S(I,15),I=1,9) / 0.6554070E-1,0.3840360E+0,0.6745410E+0,-.
     *   1021300E+0,0.8159220E-1,0.9697880E+0,-.3714950E+0,0.1270990E+1,
     *   0.1000000E+1 /
      DATA (S(I,16),I=1,9) / 0.6500710E-1,0.3820400E+0,0.6765450E+0,-.
     *   1003100E+0,0.6508770E-1,0.9814550E+0,-.2860890E+0,0.1228060E+1,
     *   0.1000000E+1 /
      DATA (S(I,17),I=1,9) / 0.6458270E-1,0.3803630E+0,0.6781900E+0,-.
     *   9876390E-1,0.5113380E-1,0.9913370E+0,-.2224010E+0,0.1182520E+1,
     *   0.1000000E+1 /
      DATA (S(I,18),I=1,9) / 0.6417070E-1,0.3787970E+0,0.6797520E+0,-.
     *   9746610E-1,0.3905690E-1,0.9999160E+0,-.1768660E+0,0.1146900E+1,
     *   0.1000000E+1 /
      DATA (P(I,3),I=1,3) / 0.1615460E+0,0.9156630E+0,0.1000000E+1 /
      DATA (P(I,4),I=1,3) / 0.2051320E+0,0.8825280E+0,0.1000000E+1 /
      DATA (P(I,5),I=1,3) / 0.2311520E+0,0.8667640E+0,0.1000000E+1 /
      DATA (P(I,6),I=1,3) / 0.2364600E+0,0.8606190E+0,0.1000000E+1 /
      DATA (P(I,7),I=1,3) / 0.2379720E+0,0.8589530E+0,0.1000000E+1 /
      DATA (P(I,8),I=1,3) / 0.2445860E+0,0.8539550E+0,0.1000000E+1 /
      DATA (P(I,9),I=1,3) / 0.2466800E+0,0.8523210E+0,0.1000000E+1 /
      DATA (P(I,10),I=1,3) / 0.2474600E+0,0.8517430E+0,0.1000000E+1 /
      DATA (P(I,11),I=1,6) / 0.1282330E+0,0.4715330E+0,0.6042730E+0,
     *   0.9066490E-2,0.9972020E+0,0.1000000E+1 /
      DATA (P(I,12),I=1,6) / 0.1210140E+0,0.4628100E+0,0.6069070E+0,
     *   0.2426330E-1,0.9866730E+0,0.1000000E+1 /
      DATA (P(I,13),I=1,6) / 0.1175740E+0,0.4611740E+0,0.6055350E+0,
     *   0.5193830E-1,0.9726600E+0,0.1000000E+1 /
      DATA (P(I,14),I=1,6) / 0.1133550E+0,0.4575780E+0,0.6074270E+0,
     *   0.6710300E-1,0.9568830E+0,0.1000000E+1 /
      DATA (P(I,15),I=1,6) / 0.1108510E+0,0.4564950E+0,0.6069360E+0,
     *   0.9158230E-1,0.9349240E+0,0.1000000E+1 /
      DATA (P(I,16),I=1,6) / 0.1096460E+0,0.4576490E+0,0.6042610E+0,
     *   0.1647770E+0,0.8708550E+0,0.1000000E+1 /
      DATA (P(I,17),I=1,6) / 0.1085980E+0,0.4586820E+0,0.6019620E+0,
     *   0.2192160E+0,0.8223210E+0,0.1000000E+1 /
      DATA (P(I,18),I=1,6) / 0.1076190E+0,0.4595760E+0,0.6000410E+0,
     *   0.2556870E+0,0.7898420E+0,0.1000000E+1 /
C
C       PARAMETERS FOR H-AR PLACED IN ASTAR ARRAY
C       DLS 3/29/89 3-21G(*) HAS 0.0 EXPONENTS FOR H->NE
C
      DATA ASTAR(1) / 0.0000000E+0 /
      DATA ASTAR(2) / 0.0000000E+0 /
      DATA ASTAR(3) / 0.0000000E+0 /
      DATA ASTAR(4) / 0.0000000E+0 /
      DATA ASTAR(5) / 0.0000000E+0 /
      DATA ASTAR(6) / 0.0000000E+0 /
      DATA ASTAR(7) / 0.0000000E+0 /
      DATA ASTAR(8) / 0.0000000E+0 /
      DATA ASTAR(9) / 0.0000000E+0 /
      DATA ASTAR(10) / 0.0000000E+0 /
      DATA ASTAR(11) / 0.1750000E+0 /
      DATA ASTAR(12) / 0.1750000E+0 /
      DATA ASTAR(13) / 0.3250000E+0 /
      DATA ASTAR(14) / 0.4500000E+0 /
      DATA ASTAR(15) / 0.5500000E+0 /
      DATA ASTAR(16) / 0.6500000E+0 /
      DATA ASTAR(17) / 0.7500000E+0 /
      DATA ASTAR(18) / 0.8500000E+0 /
C
C     DIFFUSE EXPONENTS FOR H ---> AR. (+ FOR LI-->AR, ++ FOR H,HE
C     (A_PLUS(HE)=0.0E+0,A_PLUS(NE)=0.0E+0,A_PLUS(AR)=0.0E+0)
C
      DATA APLUS / 0.0360E+0,0.00E+0,0.00740E+0,0.02070E+0,0.03150E+0,
     *   0.04380E+0,0.06390E+0,0.08450E+0,0.10760E+0,0.00E+0,0.00760E+0,
     *   0.01460E+0,0.03180E+0,0.03310E+0,0.03480E+0,0.04050E+0,0.04830E
     *   +0,0.0000E+0 /
C
C       NOW INITIALIZE ALL OF THE NON-"R" DEPENDANT VALUES RATHER THAN
C       RECOMPUTING THEM MXPTS**3 TIMES IN THE Z,Y,X LOOP OVER THE
C       ORBITAL VALUE MATRIX (DENSIT)
C
C       THESE VALUES ARE ALSO INDEPENDANT OF ATOM TYPE, AND ONLY DEPENDA
C       ON THE ROW OF THE PERIODIC TABLE AND WHETHER IT IS "S" OR "P"
C       INITIALIZE THEM HERE AND ACCESS THEM WITHIN THE NAT LOOP, BEFORE
C       ENTERING THE LOOP OVER THE GRID (CUBE) POINTS.
C
C       DETERMINE THE NUMBER OF AOS AND READ IN THE WAVEFUNCTION.
C
C       WHEN THE (*) FUNCTIONS ARE IMPLEMENTED, NECCESSARY CHECKS WILL
C       HAVE TO BE ADDED HERE.....  1/88 DLS
C
      N = 0
      DO 10 J = 1, NAT
         IF (IAN(J).LE.2) THEN
            N = N+2
C
C       CHECK FOR '3-21++G...' OR 3-21G* VARIATIONS
C
            IF (INDEX(CALC,'++').NE.0) N = N+1
         ELSE
            N = N+9
            IF (INDEX(CALC,'+').NE.0) N = N+4
            IF (IAN(J).GE.11) THEN
               N = N+4
               IF (INDEX(CALC,'*').NE.0) N = N+6
            ENDIF
         ENDIF
   10 CONTINUE
C
C       READ IN EIGENVECTORS
C       IT IS ASSUMED THAT THE EIGENVECTORS HAVE BEEN NORMALIZED TO 1
C       ELECTRON WITH THE OVERLAP MATRIX INCLUDED.
C
      IF (IONEMO.NE.0) THEN
         READ (IRD,20) (V(I,1),I=1,N)
      ELSE
         READ (IRD,20,END=30) ((V(I,J),I=1,N),J=1,N)
   20    FORMAT (8F10.6)
      ENDIF
   30 CONTINUE
      DO 40 J = 1, NAT
         XYZ(1,J) = C(J,1)
         XYZ(2,J) = C(J,2)
         XYZ(3,J) = C(J,3)
   40 CONTINUE
C
C       COMPUTE X,Y, AND Z VALUES FOR THE GRID (CUBE) POINTS.
C
      X(1) = XMI
      Y(1) = YMI
      Z(1) = ZMI
      DO 50 I = 2, MXPTS
         X(I) = XINC+X(I-1)
         Y(I) = YINC+Y(I-1)
         Z(I) = ZINC+Z(I-1)
   50 CONTINUE
C
C       ZERO THE ORBITAL VALUE ARRAY
C
      DO 80 IZ = 1, MXPTS
         DO 70 IY = 1, MXPTS
            DO 60 IX = 1, MXPTS
               DENSIT(IX,IY,IZ) = 0.0E+0
   60       CONTINUE
   70    CONTINUE
   80 CONTINUE
C
C       INITIALIZE THE AO COUNTER AND LOOP OVER ATOMS, IAT IS THE ATOMIC
C       NUMBER OF THE I'TH ATOM
C
C       FOR THE PRESENT, ONLY PLOTTING THE MO SPECIFIED BY MONE
C
      MO = MONE
      M = 1
      DO 840 I = 1, NAT
         IAT = IAN(I)
C
C       COMPUTE XDEL,YDEL,AND ZDEL (I.E. DELTA X,Y, AND Z FROM THE ATOM
C       TO EACH POINT ON THE GRID.  ONLY MXPTS VALUES FOR EACH SINCE,
C       FOR INSTANCE, EVERY POINT ON A PARTICULAR XY PLANE IS THE SAME
C       DELTA Z VALUE FROM THE POINT.  THEREFORE YOU HAVE ONLY ONE VALUE
C       FOR THE ENTIRE PLANE FOR DELTA Z, INSTEAD OF (FOR MXPTS=51)
C       2601.  AGAIN, BY COMPUTING THIS HERE, RATHER THAN INSIDE THE
C       LOOP WE CUT DOWN THESE SUBTRACTIONS AND MULTIPLICATIONS BY
C       A FACTOR OF 2601 TO 1. THIS HAS A SUBSTANTIAL EFFECT ON THE
C       SPEED OF THE COMPUTATIONS.
C
         DO 90 IXYZ = 1, MXPTS
            XDEL(IXYZ) = X(IXYZ)-XYZ(1,I)
            XDELSQ(IXYZ) = XDEL(IXYZ)*XDEL(IXYZ)
   90    CONTINUE
         DO 100 IXYZ = 1, MXPTS
            YDEL(IXYZ) = Y(IXYZ)-XYZ(2,I)
            YDELSQ(IXYZ) = YDEL(IXYZ)*YDEL(IXYZ)
  100    CONTINUE
         DO 110 IXYZ = 1, MXPTS
            ZDEL(IXYZ) = Z(IXYZ)-XYZ(3,I)
            ZDELSQ(IXYZ) = ZDEL(IXYZ)*ZDEL(IXYZ)
  110    CONTINUE
         WRITE (ILST,'(2(A,I5))') 
     *       'PROCESSING ATOM NUMBER ',I,' ATOMIC NUMBER ',IAT
C
C
C       **********************
C       ***                ***
C       ***     H TO HE    ***
C       ***                ***
C       **********************
C
C       THE INDIVIDUAL EXPONENTIATIONS RELATED TO XDEL,YDEL,
C       AND ZDEL WILL BE PRECOMPUTED, STORED IN ARRAYS, AND
C       THE VALUES MULTIPLIED IN THE LOOP, RATHER THAN
C       MAXPTS**3 SEPARATE EXPONENTIATIONS OVER R.  THESE
C       LM*3*MXPTS VALUES ARE THE UNIQUE ONES FOR H-HE.
C       (L+M GAUSSIANS*3 (X,Y AND Z) * MXPTS PLANES )
C       L,M REFER TO GAUSSIANS IN A K-LMG BASIS
C
         IF (IAT.LE.2) THEN
C
C     INNER 1S
C
            DO 120 IG = 1, LD
               AO1S(IG) = S(IG,IAT)*(A(IG,IAT)**0.750E+0)*0.822240980E+0
     *            *V(M,MO)
  120       CONTINUE
C
C     OUTER - LM IS L+M FROM K-LMG,THIS WILL HELP ADDING NEW BASIS SETS
C
            AO1S(LM) = S(LM,IAT)*(A(LM,IAT)**0.750E+0)*0.822240980E+0*V(
     *         M+1,MO)
C
C     PRE-NEGATE EXPONENT FOR EXPONENTIAL FUNCTION
C
            DO 130 IG = 1, LM
               ANEG(IG) = -A(IG,IAT)
  130       CONTINUE
C
C       IG IS THE COUNTER OVER THE L+M GAUSSIANS
C
            DO 150 IG = 1, LM
               DO 140 IXYZ = 1, MXPTS
C
C       X, Y AND Z DEPENDANT PARTS OF THE EXPONENTIAL TERM, THE R
C       INDEPENDANT PART (AO1S) IS MULTIPLIED IN HERE, RATHER THAN
C       INSIDE THE LOOP.  NOTHING ELSE GETS MULTIPLIED BY THE EXP.
C       PART, SO NOTHING IS BEING CORRUPTED.  THE GIST OF WHAT
C       IS EQUIVALENTLY BEING DONE INSIDE THE INNER LOOP IS
C       E^(-ALPHA(IG)*R^2)*AO1S(IG) - WHERE THE EXPONENTIAL TERM
C       IS FACTORED INTO E^(-ALPHA*X^2)*E^(-ALPHA*Y^2)*E^(-ALPHA*Z^2)
C
                  E1SX(IXYZ,IG) = GEXP(XDELSQ(IXYZ)*ANEG(IG))*AO1S(IG)
                  E1SY(IXYZ,IG) = GEXP(YDELSQ(IXYZ)*ANEG(IG))
                  E1SZ(IXYZ,IG) = GEXP(ZDELSQ(IXYZ)*ANEG(IG))
  140          CONTINUE
  150       CONTINUE
            DO 190 IZ = 1, MXPTS
               DO 180 IY = 1, MXPTS
                  DO 170 IG = 1, LM
                     DO 160 IX = 1, MXPTS
C
C       CONTR IS THE SUM OF CONTRIBUTIONS OVER THIS YZ PLANE FOR THIS
C       ATOM.   WHEN FINISHED, SUM INTO THE ORBITAL VALUE ARRAY.
C
                        DENSIT(IX,IY,IZ) = DENSIT(IX,IY,IZ)+E1SY(IY,IG)*
     *                     E1SZ(IZ,IG)*E1SX(IX,IG)
  160                CONTINUE
  170             CONTINUE
  180          CONTINUE
  190       CONTINUE
            M = M+2
C
C       CHECK FOR ++ ANION DISFFUSE FUNCTION, DIFFUSE S ORBITAL ON H
C       USE SAME VARIABLES AS ABOVE FOR CONVENIENCE
C
            IF (INDEX(CALC,'++').NE.0) THEN
               AO1S(1) = (APLUS(IAT)**0.750E+0)*0.822240980E+0*V(M,MO)
               APNEG = -APLUS(IAT)
               DO 200 IXYZ = 1, MXPTS
                  E1SX(IXYZ,1) = GEXP(XDELSQ(IXYZ)*APNEG)*AO1S(1)
                  E1SY(IXYZ,1) = GEXP(YDELSQ(IXYZ)*APNEG)
                  E1SZ(IXYZ,1) = GEXP(ZDELSQ(IXYZ)*APNEG)
  200          CONTINUE
C
C       LOOP OVER THE "CUBE" Z,Y, AND X
C
               DO 230 IZ = 1, MXPTS
                  DO 220 IY = 1, MXPTS
                     DO 210 IX = 1, MXPTS
C
C    SUM THE DIFFUSE ORBITAL CONTRIBUTION INTO THE ORBITAL VALUE ARRAY
C    (ONLY ONE PRIMITIVE )
C
                        DENSIT(IX,IY,IZ) = DENSIT(IX,IY,IZ)+E1SX(IX,1)*
     *                     E1SY(IY,1)*E1SZ(IZ,1)
  210                CONTINUE
  220             CONTINUE
  230          CONTINUE
               M = M+1
C
C       DONE WITH ++ CODE
C
            ENDIF
C
C       **********************
C       ***                ***
C       ***    LI TO NE    ***
C       ***                ***
C       **********************
C
         ELSEIF (IAT.LE.10) THEN
C
C       AO1S(1-KD) ARE THE KD "CONSTANTS" FOR THE 1S
C       ORBITAL, EVERYTHING THAT IS INDEPENDANT OF R IS IN THERE.
C       (1 PER GAUSSIAN PRIMITIVE)
C
            DO 240 IG = 1, KD
               AO1S(IG) = S(IG,IAT)*(A(IG,IAT)**0.750E+0)*0.712705470E+0
     *            *V(M,MO)
  240       CONTINUE
C
C       INNER S FUNCTION
C
            DO 250 IG = 1, LD
               AO2S(IG) = S(KD+IG,IAT)*(A(KD+IG,IAT)**0.750E+0)*
     *            0.712705470E+0*V(M+1,MO)
  250       CONTINUE
C
C       OUTER S FUNCTION
C
            AO2S(LM) = S(KD+LM,IAT)*(A(KD+LM,IAT)**0.750E+0)*
     *         0.712705470E+0*V(M+5,MO)
C
C       INNER P FUNCTION
C
            DO 260 IG = 1, LD
               AO2P(IG) = P(IG,IAT)*(A(KD+IG,IAT)**1.250E+0)*
     *            1.425410940E+0
  260       CONTINUE
C
C       OUTER P FUNCTION
C
            AO2P(LM) = P(LM,IAT)*(A(KD+LM,IAT)**1.250E+0)*1.425410940E+0
C
C       INNER P * COEFFICIENT (PX)
C
            DO 270 IG = 1, LD
               AO2PX(IG) = AO2P(IG)*V(M+2,MO)
  270       CONTINUE
C
C       OUTER P * COEFFICIENT (PX)
C
            AO2PX(LM) = AO2P(LM)*V(M+6,MO)
C
C       INNER P * COEFFICIENT (PY)
C
            DO 280 IG = 1, LD
               AO2PY(IG) = AO2P(IG)*V(M+3,MO)
  280       CONTINUE
C
C       OUTER P * COEFFICIENT (PY)
C
            AO2PY(LM) = AO2P(LM)*V(M+7,MO)
C
C       INNER P * COEFFICIENT (PZ)
C
            DO 290 IG = 1, LD
               AO2PZ(IG) = AO2P(IG)*V(M+4,MO)
  290       CONTINUE
C
C       OUTER P * COEFFICIENT (PZ)
C
            AO2PZ(LM) = AO2P(LM)*V(M+8,MO)
C
C       AO2S(1-LM) CORRESPONDS TO THE LM GAUSSIANS FOR THE 2S ORBITAL
C       AO2PX(1-LM) "                            " FOR THE 2PX ""
C       ETC. (LM IS L+M FROM THE K-LMG BASIS USED )
C       YOU NEED DIFFERENT "CONSTANTS" FOR EACH DUE TO THE DIFFERENT
C       VALUE OF THE WAVEFUNCTION FOR EACH ORBITAL. AT THE END WE WILL
C       HAVE A 3 SETS OF NUMBERS ALL MULTIPLIED BY THE EXPONENTIAL,
C       WE CAN ADD THEM FIRST, THEN CALC. AND MULT. BY THE EXP.
C
            DO 310 IG = 1, LM
               DO 300 IXYZ = 1, MXPTS
C
C       X - DEPENDANT PART ( 2PX )
C
                  CNS2PX(IXYZ,IG) = AO2PX(IG)*XDEL(IXYZ)
C
C       Y - DEPENDANT PART ( 2PY )
C
                  CNS2PY(IXYZ,IG) = AO2PY(IG)*YDEL(IXYZ)
C
C       AO2PZ(7->10) ARE THE 2S MULTIPLIERS.  THEY CAN BE ADDED TO ANY
C       ONE OF THE 3 MULTIPLIERS OUTSIDE THE LOOP, SINCE THEY HAVE NO
C       X,Y, OR Z DEPENDANCE. THESE 4 ALL HAVE THE SAME EXP TO MULTIPLY
C       BY AND BE ADDED TO THE DENSITY MATRIX, SO THEY CAN BE ADDED
C       BEFORE, AND THEN MULTIPLIED BY THE EXP ALL AT ONCE.
C
                  CNS2PZ(IXYZ,IG) = AO2PZ(IG)*ZDEL(IXYZ)+AO2S(IG)
  300          CONTINUE
  310       CONTINUE
C
C       EVALUATE 2S, 2P, AND 1S
C
C       MINUS ALPHA FOR 1S:
C
            DO 320 IG = 1, KD
               ANEG(IG) = -A(IG,IAT)
  320       CONTINUE
C
C       MINUS ALPHA FOR 2SP:
C
            DO 330 IG = KD+1, KD+LM
               ANEG(IG) = -A(IG,IAT)
  330       CONTINUE
C
C       PRECOMPUTE EXP(-A*DELO**2) WHERE DELO**2 IS
C       DELTA-X**2, DELTA-Y**2, AND DELTA-Z**2
C       THIS WILL BE K*3(X,Y, AND Z)*MXPTS FOR K-LMG CALCULATIONS
C       ON THE 1S CORE.
C
            DO 350 IG = 1, KD
               DO 340 IXYZ = 1, MXPTS
                  E1SX(IXYZ,IG) = GEXP(XDELSQ(IXYZ)*ANEG(IG))*AO1S(IG)
                  E1SY(IXYZ,IG) = GEXP(YDELSQ(IXYZ)*ANEG(IG))
                  E1SZ(IXYZ,IG) = GEXP(ZDELSQ(IXYZ)*ANEG(IG))
  340          CONTINUE
  350       CONTINUE
C
C       NOW THE VALENCE EXPONENTS, THIS WILL BE (L+M)*3(X,Y, AND Z)*
C       MXPTS COMPUTATIONS FOR K-LMG CALCULATIONS. THIS SAVES DOING
C       (L+M)*1(R)*MXPTS**3 COMPUTATIONS INSIDE THE LOOP.  THIS CODE
C       WAS REPONSIBLE FOR A TEST STO-3G CASE (F2) GOING FROM 1:45
C       TO 0:13 (MIN:SEC).
C
            DO 370 IG = 1, LM
               DO 360 IXYZ = 1, MXPTS
C
C       X^2 PART OF EXPONENTIAL
C
                  E2SPX(IXYZ,IG) = GEXP(XDELSQ(IXYZ)*ANEG(KD+IG))
C
C       Y^2 PART OF EXPONENTIAL
C
                  E2SPY(IXYZ,IG) = GEXP(YDELSQ(IXYZ)*ANEG(KD+IG))
C
C       Z^2 PART OF EXPONENTIAL
C
                  E2SPZ(IXYZ,IG) = GEXP(ZDELSQ(IXYZ)*ANEG(KD+IG))
  360          CONTINUE
  370       CONTINUE
C
C       LOOP OVER THE "CUBE" Z,Y, AND X
C
            DO 430 IZ = 1, MXPTS
C
C       SUM THE NON-EXPONENTIAL PARTS OF THE 2SP SHELL.  THE 2S HAS
C       ALREADY BEEN SUMMED IN, AS IT HAS NO DEPENDANCE ON X,Y, OR Z.
C
               DO 420 IY = 1, MXPTS
C
C       NOW, LOOP OVER X, INCLUDE X-DEPENDANT CONTRIBUTIONS, COMPUTE THE
C       ORBITAL VALUE AND SUM INTO THE DENSITY ARRAY.
C
                  DO 390 IG = 1, KD
                     DO 380 IX = 1, MXPTS
C
C       COMPUTE AND SUM THE ORBITAL CONTRIBUTIONS INTO THE ORBITAL
C       VALUE ARRAY
C
C       FIRST THE KD GAUSSIANS FOR THE 1S:
C
                        DENSIT(IX,IY,IZ) = DENSIT(IX,IY,IZ)+E1SY(IY,IG)*
     *                     E1SZ(IZ,IG)*E1SX(IX,IG)
  380                CONTINUE
  390             CONTINUE
C
C       NEXT, THE LM (L-INNER-M-OUTER) FOR THE 2SP:
C
                  DO 410 IG = 1, LM
                     DO 400 IX = 1, MXPTS
                        DENSIT(IX,IY,IZ) = DENSIT(IX,IY,IZ)+(CNS2PY(IY,
     *                     IG)+CNS2PZ(IZ,IG)+CNS2PX(IX,IG))*E2SPX(IX,IG)
     *                     *E2SPY(IY,IG)*E2SPZ(IZ,IG)
  400                CONTINUE
  410             CONTINUE
  420          CONTINUE
  430       CONTINUE
            M = M+9
C
C       THIS WILL BE DONE FOR BOTH 3-21+G AND 3-21++G
C
            IF (INDEX(CALC,'+').NE.0) THEN
               APNEG = -APLUS(IAT)
               DO 440 IXYZ = 1, MXPTS
                  E2SPX(IXYZ,1) = GEXP(XDELSQ(IXYZ)*APNEG)
                  E2SPY(IXYZ,1) = GEXP(YDELSQ(IXYZ)*APNEG)
                  E2SPZ(IXYZ,1) = GEXP(ZDELSQ(IXYZ)*APNEG)
  440          CONTINUE
               AO2S(1) = V(M,MO)*(APLUS(IAT)**0.750E+0)*0.712705470E+0
               AO2P(1) = (APLUS(IAT)**1.250E+0)*1.425410940E+0
               AO2PX(1) = AO2P(1)*V(M+1,MO)
               AO2PY(1) = AO2P(1)*V(M+2,MO)
               AO2PZ(1) = AO2P(1)*V(M+3,MO)
               DO 450 IXYZ = 1, MXPTS
                  CNS2PX(IXYZ,1) = AO2PX(1)*XDEL(IXYZ)
                  CNS2PY(IXYZ,1) = AO2PY(1)*YDEL(IXYZ)
C
C       AO2S(1->4) ARE THE 2 S MULTIPLIERS.  THEY CAN BE ADDED TO ANY
C       ONE OF THE 3 MULTIPLIERS OUTSIDE THE LOOP, SINCE THEY HAVE NO
C       X,Y, OR Z DEPENDANCE. THESE 4 ALL HAVE THE SAME EXP TO MULTIPLY
C       BY AND BE ADDED TO THE DENSITY MATRIX, SO THEY CAN BE ADDED
C       BEFORE, AND THEN MULTIPLIED BY THE EXP ALL AT ONCE.
C
                  CNS2PZ(IXYZ,1) = AO2PZ(1)*ZDEL(IXYZ)+AO2S(1)
  450          CONTINUE
C
C       LOOP OVER THE "CUBE" Z,Y, AND X
C
               DO 480 IZ = 1, MXPTS
                  DO 470 IY = 1, MXPTS
                     DO 460 IX = 1, MXPTS
C
C       SUM THE DIFFUSE ORBITAL CONTRIBUTIONS INTO THE ORBITAL VALUE
C       ARRAY (ONLY ONE PRIMITIVE )
C
                        DENSIT(IX,IY,IZ) = DENSIT(IX,IY,IZ)+(CNS2PY(IY,1
     *                     )+CNS2PZ(IZ,1)+CNS2PX(IX,1))*E2SPX(IX,1)*
     *                     E2SPY(IY,1)*E2SPZ(IZ,1)
  460                CONTINUE
  470             CONTINUE
  480          CONTINUE
               M = M+4
            ENDIF
C
C       **********************
C       ***                ***
C       ***    NA TO AR    ***
C       ***                ***
C       **********************
C
         ELSEIF (IAT.LE.18) THEN
            DO 490 IG = 1, KD
               IG2 = KD+IG
C
C       AO1S(1-KD) ARE THE KD GAUSSION "CONSTANTS" FOR THE 1S
C       ORBITAL, EVERYTHING THAT IS INDEPENDANT OF R IS IN THERE.
C
               AO1S(IG) = S(IG,IAT)*(A(IG,IAT)**0.750E+0)*0.712705470E+0
     *            *V(M,MO)
C
C       CORE 2SP SHELL.... 2S:
C
               AO2S(IG) = S(IG2,IAT)*(A(IG2,IAT)**0.750E+0)*0.712705470E
     *            +0*V(M+1,MO)
C
C       CORE 2SP SHELL.... 2P (GENERAL PARTS):
C
               AO2P(IG) = P(IG,IAT)*(A(IG2,IAT)**1.250E+0)*1.425410940E+
     *            0
C
C       CORE 2SP SHELL.... 2PX:
C
               AO2PX(IG) = AO2P(IG)*V(M+2,MO)
C
C       CORE 2SP SHELL.... 2PY:
C
               AO2PY(IG) = AO2P(IG)*V(M+3,MO)
C
C       CORE 2SP SHELL.... 2PZ:
C
               AO2PZ(IG) = AO2P(IG)*V(M+4,MO)
  490       CONTINUE
C
C       INNER S FUNCTION
C
            DO 500 IG = 1, LD
               AO3S(IG) = S(K2+IG,IAT)*(A(K2+IG,IAT)**0.750E+0)*
     *            0.712705470E+0*V(M+5,MO)
  500       CONTINUE
C
C       OUTER S FUNCTION
C
            AO3S(LM) = S(K2+LM,IAT)*(A(K2+LM,IAT)**0.750E+0)*
     *         0.712705470E+0*V(M+9,MO)
C
C       INNER P FUNCTION
C
            DO 510 IG = 1, LD
               AO3P(IG) = P(KD+IG,IAT)*(A(K2+IG,IAT)**1.250E+0)*
     *            1.425410940E+0
  510       CONTINUE
C
C       OUTER P FUNCTION
C
            AO3P(LM) = P(KD+LM,IAT)*(A(K2+LM,IAT)**1.250E+0)*
     *         1.425410940E+0
C
C       INNER P * COEFFICIENT (PX)
C
            DO 520 IG = 1, LD
               AO3PX(IG) = AO3P(IG)*V(M+6,MO)
  520       CONTINUE
C
C       OUTER P * COEFFICIENT (PX)
C
            AO3PX(LM) = AO3P(LM)*V(M+10,MO)
C
C       INNER P * COEFFICIENT (PY)
C
            DO 530 IG = 1, LD
               AO3PY(IG) = AO3P(IG)*V(M+7,MO)
  530       CONTINUE
C
C       OUTER P * COEFFICIENT (PY)
C
            AO3PY(LM) = AO3P(LM)*V(M+11,MO)
C
C       INNER P * COEFFICIENT (PZ)
C
            DO 540 IG = 1, LD
               AO3PZ(IG) = AO3P(IG)*V(M+8,MO)
  540       CONTINUE
C
C       OUTER P * COEFFICIENT (PZ)
C
            AO3PZ(LM) = AO3P(LM)*V(M+12,MO)
C
C       CORE 2P X,Y,AND Z DEPENDANT PARTS
C
            DO 560 IG = 1, KD
               DO 550 IXYZ = 1, MXPTS
C
C       CORE 2PX - X DEPENDANT PART
C
                  CNS2PX(IXYZ,IG) = AO2PX(IG)*XDEL(IXYZ)
C
C       CORE 2PY - Y DEPENDANT PART
C
                  CNS2PY(IXYZ,IG) = AO2PY(IG)*YDEL(IXYZ)
C
C       AO2S(1->6) ARE THE 2S MULTIPLIERS.  THEY CAN BE ADDED TO ANY
C       ONE OF THE 3 MULTIPLIERS OUTSIDE THE LOOP, SINCE THEY HAVE NO
C       X,Y, OR Z DEPENDANCE. THESE 4 ALL HAVE THE SAME EXP TO MULTIPLY
C       BY AND BE ADDED TO THE DENSITY MATRIX, SO THEY CAN BE ADDED
C       BEFORE, AND THEN MULTIPLIED BY THE EXP ALL AT ONCE.
C
C       CORE 2PZ - Z DEPENDANT PART - S ADDED IN HERE
C
                  CNS2PZ(IXYZ,IG) = AO2PZ(IG)*ZDEL(IXYZ)+AO2S(IG)
  550          CONTINUE
  560       CONTINUE
C
C       AO3S  CORRESPONDS TO THE 4 GAUSSIANS FOR THE 3S ORBITAL
C       AO3PX "                            " FOR THE 3PX ""
C       ETC.
C       YOU NEED DIFFERENT "CONSTANTS" FOR EACH DUE TO THE DIFFERENT
C       VALUE OF THE WAVEFUNCTION FOR EACH ORBITAL. AT THE END WE WILL
C       HAVE A 3 SETS OF NUMBERS ALL MULTIPLIED BY THE EXPONENTIAL,
C       WE CAN ADD THEM FIRST, THEN CALC. AND MULT. BY THE EXP.
C
            DO 580 IG = 1, LM
               DO 570 IXYZ = 1, MXPTS
                  CNS3PX(IXYZ,IG) = AO3PX(IG)*XDEL(IXYZ)
                  CNS3PY(IXYZ,IG) = AO3PY(IG)*YDEL(IXYZ)
C
C       AO3S(1->4) ARE THE 3S MULTIPLIERS.  THEY CAN BE ADDED TO ANY
C       ONE OF THE 3 MULTIPLIERS OUTSIDE THE LOOP, SINCE THEY HAVE NO
C       X,Y, OR Z DEPENDANCE. THESE 4 ALL HAVE THE SAME EXP TO MULTIPLY
C       BY AND BE ADDED TO THE DENSITY MATRIX, SO THEY CAN BE ADDED
C       BEFORE, AND THEN MULTIPLIED BY THE EXP ALL AT ONCE.
C
                  CNS3PZ(IXYZ,IG) = AO3PZ(IG)*ZDEL(IXYZ)+AO3S(IG)
  570          CONTINUE
  580       CONTINUE
C
C       EVALUATE 1S, 2SP, AND 3SP EXPONENTIALS
C
C       MINUS ALPHA FOR 1S:
C
            DO 590 IG = 1, KD
               ANEG(IG) = -A(IG,IAT)
  590       CONTINUE
C
C       MINUS ALPHA FOR 2SP:
C
            DO 600 IG = 1, KD
               ANEG(KD+IG) = -A(KD+IG,IAT)
  600       CONTINUE
C
C       MINUS ALPHA FOR 3SP:
C
            DO 610 IG = 1, LM
               ANEG(K2+IG) = -A(K2+IG,IAT)
  610       CONTINUE
C
C       PRECOMPUTE EXP(-A*DELO**2) WHERE DELO**2 IS
C       DELTA-X**2, DELTA-Y**2, AND DELTA-Z**2
C       THIS WILL BE K*3(X,Y, AND Z)*MXPTS FOR K-LMG CALCULATIONS
C       ON THE 1S CORE.
C
C       THE AO CONTRIBUTION CAN BE MULTIPLIED IN RIGHT HERE, SINCE
C       THERE IS NO X,Y, OR Z DEPENDANCE ON IT, AND NOTHING ELSE
C       NEEDS TO BE MULTIPLIED INTO THE EXPONENTIAL TERM.
C
            DO 630 IG = 1, KD
               DO 620 IXYZ = 1, MXPTS
C
C       E^(-ALPHA(1S)*X^2)
C
                  E1SX(IXYZ,IG) = GEXP(XDELSQ(IXYZ)*ANEG(IG))*AO1S(IG)
C
C       E^(-ALPHA(1S)*Y^2)
C
                  E1SY(IXYZ,IG) = GEXP(YDELSQ(IXYZ)*ANEG(IG))
C
C       E^(-ALPHA(1S)*Z^2)
C
                  E1SZ(IXYZ,IG) = GEXP(ZDELSQ(IXYZ)*ANEG(IG))
  620          CONTINUE
  630       CONTINUE
            DO 650 IG = 1, KD
               DO 640 IXYZ = 1, MXPTS
C
C       E^(-ALPHA(2SP)*X^2)
C
                  E2SPX(IXYZ,IG) = GEXP(XDELSQ(IXYZ)*ANEG(KD+IG))
C
C       E^(-ALPHA(2SP)*Y^2)
C
                  E2SPY(IXYZ,IG) = GEXP(YDELSQ(IXYZ)*ANEG(KD+IG))
C
C       E^(-ALPHA(2SP)*Z^2)
C
                  E2SPZ(IXYZ,IG) = GEXP(ZDELSQ(IXYZ)*ANEG(KD+IG))
  640          CONTINUE
  650       CONTINUE
C
C       NOW THE VALENCE EXPONENTS, THIS WILL BE (L+M)*3(X,Y, AND Z)*
C       MXPTS COMPUTATIONS FOR K-LMG CALCULATIONS. THIS SAVES DOING
C       (L+M)*1(R)*MXPTS**3 COMPUTATIONS INSIDE THE LOOP.  THIS CODE
C       WAS REPONSIBLE FOR A TEST STO-3G CASE (F2) GOING FROM 1:45
C       TO 0:13 (MIN:SEC).
C
            DO 670 IG = 1, LM
               DO 660 IXYZ = 1, MXPTS
C
C       E^(-ALPHA(3SP)*X^2)
C
                  E3SPX(IXYZ,IG) = GEXP(XDELSQ(IXYZ)*ANEG(K2+IG))
C
C       E^(-ALPHA(3SP)*Y^2)
C
                  E3SPY(IXYZ,IG) = GEXP(YDELSQ(IXYZ)*ANEG(K2+IG))
C
C       E^(-ALPHA(3SP)*Z^2)
C
                  E3SPZ(IXYZ,IG) = GEXP(ZDELSQ(IXYZ)*ANEG(K2+IG))
  660          CONTINUE
  670       CONTINUE
C
C       LOOP OVER THE "CUBE" Z,Y, AND X
C
            DO 730 IZ = 1, MXPTS
C
C       SUM THE NON-EXPONENTIAL PARTS OF THE 2SP SHELL.  THE 2S HAS
C       ALREADY BEEN SUMMED IN, AS IT HAS NO DEPENDANCE ON X,Y, OR Z.
C
               DO 720 IY = 1, MXPTS
                  DO 690 IG = 1, KD
C
C       NOW, LOOP OVER X, INCLUDE X-DEPENDANT CONTRIBUTIONS, COMPUTE THE
C       ORBITAL VALUE AND SUM INTO THE DENSITY ARRAY.
C
                     DO 680 IX = 1, MXPTS
C
C       COMPUTE AND SUM THE ORBITAL CONTRIBUTIONS INTO THE ORBITAL
C       VALUE ARRAY
C
C       FIRST THE KD GAUSSIANS FOR THE 1S:
C
                        DENSIT(IX,IY,IZ) = DENSIT(IX,IY,IZ)+E1SY(IY,IG)*
     *                     E1SZ(IZ,IG)*E1SX(IX,IG)
C
C       NOW THE KD FOR THE CORE 2SP SHELL
C
                        DENSIT(IX,IY,IZ) = DENSIT(IX,IY,IZ)+(CNS2PY(IY,
     *                     IG)+CNS2PZ(IZ,IG)+CNS2PX(IX,IG))*E2SPX(IX,IG)
     *                     *E2SPY(IY,IG)*E2SPZ(IZ,IG)
  680                CONTINUE
  690             CONTINUE
C
C       NOW THE LM (L-INNER-M-OUTER) FOR THE VALENCE 3SP SHELL
C
                  DO 710 IG = 1, LM
                     DO 700 IX = 1, MXPTS
                        DENSIT(IX,IY,IZ) = DENSIT(IX,IY,IZ)+(CNS3PY(IY,
     *                     IG)+CNS3PZ(IZ,IG)+CNS3PX(IX,IG))*E3SPX(IX,IG)
     *                     *E3SPY(IY,IG)*E3SPZ(IZ,IG)
  700                CONTINUE
  710             CONTINUE
  720          CONTINUE
  730       CONTINUE
            M = M+13
C
C       THIS WILL BE DONE FOR (*) WAVEFUNCTIONS, WITH ANY
C       COMBINATION OF (OR LACK OF) + FUNCTIONS
C
            IF (INDEX(CALC,'*').NE.0.OR.INDEX(CALC,'D').NE.0) THEN
               ASNEG = -ASTAR(IAT)
               DO 740 IXYZ = 1, MXPTS
                  EXPDX(IXYZ) = GEXP(XDELSQ(IXYZ)*ASNEG)
                  EXPDY(IXYZ) = GEXP(YDELSQ(IXYZ)*ASNEG)
                  EXPDZ(IXYZ) = GEXP(ZDELSQ(IXYZ)*ASNEG)
  740          CONTINUE
               TMP = (ASTAR(IAT)**1.750E+0)*1.645922780E+0
               AODXX = V(M,MO)*TMP
               AODYY = V(M+1,MO)*TMP
               AODZZ = V(M+2,MO)*TMP
               AODXY = V(M+3,MO)*TMP
               AODXZ = V(M+4,MO)*TMP
               AODYZ = V(M+5,MO)*TMP
C
C       CNSXY AND CNSXZ WILL BE ADDED OUTSIDE THE INNER LOOP, THEN
C       THE SUM WILL BE MULTIPLIED BY XDEL(IX) IN THE LOOP.  THAT SHOULD
C       SAVE ONE ADDITION IN THE INNER LOOP.
C
               DO 750 IXYZ = 1, MXPTS
                  CNSXX(IXYZ) = AODXX*XDELSQ(IXYZ)
                  CNSYY(IXYZ) = AODYY*YDELSQ(IXYZ)
                  CNSZZ(IXYZ) = AODZZ*ZDELSQ(IXYZ)
                  CNSXY(IXYZ) = AODXY*YDEL(IXYZ)
                  CNSXZ(IXYZ) = AODXZ*ZDEL(IXYZ)
                  CNSYZ(IXYZ) = AODYZ*YDEL(IXYZ)
  750          CONTINUE
C
C       LOOP OVER THE "CUBE" Z,Y, AND X
C
               DO 780 IZ = 1, MXPTS
C
C       THE Z^2,YZ, AND Y^2 ARE ALL "CONSTANT" WITHING THE X LOOP, THEY
C       CAN BE ADDED NOW
C
                  DO 770 IY = 1, MXPTS
                     DO 760 IX = 1, MXPTS
C
C       SUM THE DIFFUSE ORBITAL CONTRIBUTIONS INTO THE ORBITAL VALUE ARR
C       (ONLY ONE PRIMITIVE )
C
                        DENSIT(IX,IY,IZ) = DENSIT(IX,IY,IZ)+(CNSZZ(IZ)+
     *                     CNSYZ(IY)*ZDEL(IZ)+CNSYY(IY)+CNSXX(IX)+
     *                     XDEL(IX)*CNSXZ(IZ)+CNSXY(IY))*EXPDX(IX)*EXPDY
     *                     (IY)*EXPDZ(IZ)
  760                CONTINUE
  770             CONTINUE
  780          CONTINUE
               M = M+6
            ENDIF
C
C       THIS WILL BE DONE FOR BOTH 3-21+G AND 3-21++G, AS WELL AS EITHER
C       CASE SUPPLEMENTED WITH '*' , I.E. 3-21+G*, 3-21++G*, ETC....
C
            IF (INDEX(CALC,'+').NE.0) THEN
               APNEG = -APLUS(IAT)
               DO 790 IXYZ = 1, MXPTS
                  E3SPX(IXYZ,1) = GEXP(XDELSQ(IXYZ)*APNEG)
                  E3SPY(IXYZ,1) = GEXP(YDELSQ(IXYZ)*APNEG)
                  E3SPZ(IXYZ,1) = GEXP(ZDELSQ(IXYZ)*APNEG)
  790          CONTINUE
C
C       USE SAME VARIABLES AS USED ABOVE, JUST FOR CONVENIENCE
C       THEY DON'T NEED ARRAYS NOW, BUT SOMEONE MAY WANT A + SHELL
C       OR D SHELL WITH MORE THAN ONE PRIMITIVE
C
               AO3S(1) = V(M,MO)*(APLUS(IAT)**0.750E+0)*0.712705470E+0
               AO3P(1) = (APLUS(IAT)**1.250E+0)*1.425410940E+0
               AO3PX(1) = AO3P(1)*V(M+1,MO)
               AO3PY(1) = AO3P(1)*V(M+2,MO)
               AO3PZ(1) = AO3P(1)*V(M+3,MO)
               DO 800 IXYZ = 1, MXPTS
                  CNS3PX(IXYZ,1) = AO3PX(1)*XDEL(IXYZ)
                  CNS3PY(IXYZ,1) = AO3PY(1)*YDEL(IXYZ)
C
C       AO3S(1) IS THE S MULTIPLIER.  IT CAN BE ADDED TO ANY
C       ONE OF THE 3 MULTIPLIERS OUTSIDE THE LOOP, SINCE THEY HAVE NO
C       X,Y, OR Z DEPENDANCE. THESE 4 ALL HAVE THE SAME EXP TO MULTIPLY
C       BY AND BE ADDED TO THE DENSITY MATRIX, SO THEY CAN BE ADDED
C       BEFORE, AND THEN MULTIPLIED BY THE EXP ALL AT ONCE.
C
                  CNS3PZ(IXYZ,1) = AO3PZ(1)*ZDEL(IXYZ)+AO3S(1)
  800          CONTINUE
C
C       LOOP OVER THE "CUBE" Z,Y, AND X
C
               DO 830 IZ = 1, MXPTS
                  DO 820 IY = 1, MXPTS
                     DO 810 IX = 1, MXPTS
C
C       SUM THE DIFFUSE ORBITAL CONTRIBUTIONS INTO THE ORBITAL VALUE ARR
C       (ONLY ONE PRIMITIVE )
C
                        DENSIT(IX,IY,IZ) = DENSIT(IX,IY,IZ)+(CNS3PY(IY,1
     *                     )+CNS3PZ(IZ,1)+CNS3PX(IX,1))*E3SPX(IX,1)*
     *                     E3SPY(IY,1)*E3SPZ(IZ,1)
  810                CONTINUE
  820             CONTINUE
  830          CONTINUE
               M = M+4
            ENDIF
         ENDIF
  840 CONTINUE
      RETURN
      END
C
C
      SUBROUTINE STOMO
      IMPLICIT REAL (A-H,O-Z)
      PARAMETER (MXPTS=51)
      PARAMETER (MAXATM=50)
      PARAMETER (MAXORB=200)
C
C       REEPLACES ORIGINAL ROUTINE TO EVALUATE RADIAL PARTS OF STO-3G
C       WAVEFUNCTIONS AS WELL AS THE MAINLINE CODE FOR THE REST.
C       REF FROM ORIGINAL CODE (FUNCTION AO):
C
C       ATOMS H TO AR ARE HANDLED ACCORDING TO J CHEM PHYS 51, 2657
C       (1969), 52, 2769 (1970).
C       W.L. JORGENSEN - MARCH,1976, JULY, 1977.
C
C       REWRITTEN INTO A STAND-ALONE SUBROUTINE FOR CALCULATING ORBITAL
C       VALUES. REWORKED TO REDUCE REDUNDANT COMPUTATION OF POWERS AND
C       SQUARE ROOTS, AS WELL AS PUTTING INTO AN EASILY VECTORIZABLE
C       FORM FOR VECTOR MACHINES ( WE HAVE A CYBER 205 AT PURDUE ).
C
C       DAN SEVERANCE 12/10/87
C
C       AFTER DISCUSSION WITH JIM BRIGGS WHO HAS BEEN USING THIS PROGRAM
C       ON DR. JORGENSEN'S GOULD FOR THE LAST FEW WEEKS; THE AO
C       COMPUTATION AND WAVEFUNCTION READ WAS MOVED FROM THE MAIN LINE
C       TO THE RESPECTIVE SUBROUTINES.  WHEN A WAVEFUNCTION IS ADDED,
C       ONLY THE SUBROUTINE SHOULD NEED SIGNIFICANT MODIFICATION, THE
C       MAINLINE SHOULD ONLY NEED TO HAVE THE CALL ADDED.
C       DAN SEVERANCE 1/17/88
C
      COMMON /SPLIT/ XMI,YMI,ZMI,XINC,YINC,ZINC,MONE,MOLAST,NAT,
     *   IONEMO
      COMMON /DENS/ DENSIT(MXPTS,MXPTS,MXPTS),V(MAXORB,MAXORB),C(MAXATM,
     *   3),OCMO(MAXORB),IAN(MAXATM)
      COMMON /IO/ IRD,ILST
      COMMON /BASIS/ CALC
      CHARACTER CALC*20
      DIMENSION XDEL(MXPTS),YDEL(MXPTS),ZDEL(MXPTS),EXP2SP(MXPTS,9)
      DIMENSION XDELSQ(MXPTS),YDELSQ(MXPTS),ZDELSQ(MXPTS)
      DIMENSION XYZ(3,50),ZSQ(3),EXP1S(MXPTS,9)
      DIMENSION X(MXPTS),Y(MXPTS),Z(MXPTS),CNSTX(MXPTS,6),CNSTY(MXPTS,6)
      DIMENSION Z1(18),Z2(18),Z3(18),A(3,3),D(3,3),D2P(3),D3P(3)
      DIMENSION CNSTNS(3,3),CNSTNP(3,3),VNORM(27),CNSTZ(MXPTS,6)
      DIMENSION ANEG(9),EXP3SP(MXPTS,9)
      DATA Z1 / 1.240E+0,1.690E+0,2.690E+0,3.680E+0,4.680E+0,5.670E+0,
     *   6.670E+0,7.660E+0,8.650E+0,9.640E+0,10.610E+0,11.590E+0,12.560E
     *   +0,13.530E+0,14.50E+0,15.470E+0,16.430E+0,17.40E+0 /
      DATA Z2 / 0.0E+0,0.0E+0,0.80E+0,1.150E+0,1.50E+0,1.720E+0,1.950E+0
     *   ,2.250E+0,2.550E+0,2.880E+0,3.480E+0,3.90E+0,4.360E+0,4.830E+0,
     *   5.310E+0,5.790E+0,6.260E+0,6.740E+0 /
      DATA Z3 / 10*0.0E+0,1.750E+0,1.70E+0,1.70E+0,1.750E+0,1.90E+0,
     *   2.050E+0,2.10E+0,2.330E+0 /
C
C       12/10/87 DAN SEVERANCE
C
      DATA A(1,1),A(1,2),A(1,3) / 1.098180E-1,4.057710E-1,2.227660E+0 /
      DATA A(2,1),A(2,2),A(2,3) / 7.513860E-2,2.310310E-1,9.942030E-1 /
      DATA A(3,1),A(3,2),A(3,3) / 5.272660E-2,1.347150E-1,4.828540E-1 /
      DATA D(1,1),D(1,2),D(1,3) / 4.446350E-1,5.353280E-1,1.543290E-1 /
      DATA D(2,1),D(2,2),D(2,3) / 7.001150E-1,3.995130E-1,-9.996720E-2 /
      DATA D(3,1),D(3,2),D(3,3) / 9.003980E-1,2.255950E-1,-2.196200E-1 /
      DATA D2P / 3.919570E-1,6.076840E-1,1.559160E-1 /
      DATA D3P / 4.620010E-1,5.951670E-1,1.058760E-2 /
C
C       NOW INITIALIZE ALL OF THE NON-"R" DEPENDANT VALUES RATHER THAN
C       RECOMPUTING THEM MXPTS**3 TIMES IN THE Z,Y,X LOOP OVER THE
C       ORBITAL VALUE MATRIX (DENSIT)
C
C       THESE VALUES ARE ALSO INDEPENDANT OF ATOM TYPE, ONLY DEPENDANT
C       ON THE ROW OF THE PERIODIC TABLE AND WHETHER IT IS "S" OR "P"
C       INITIALIZE THEM HERE AND ACCESS THEM WITHIN THE NAT LOOP, BEFORE
C       ENTERING THE LOOP OVER THE GRID (CUBE) POINTS.
C
      N = 0
      DO 10 I = 1, NAT
         N = N+1
         IF (IAN(I).GT.2) THEN
            N = N+4
            IF (IAN(I).GT.10) THEN
               N = N+4
               IF (IAN(I).GT.18) THEN
                  WRITE (ILST,*)
     *               'ATOMIC NUMBERS > 18 NOT YET IMPLEMENTED'
               ENDIF
            ENDIF
         ENDIF
   10 CONTINUE
C
C       READ IN EIGENVECTORS
C       IT IS ASSUMED THAT THE EIGENVECTORS HAVE BEEN NORMALIZED TO 1
C       ELECTRON WITH THE OVERLAP MATRIX INCLUDED.
C
      IF (IONEMO.NE.0) THEN
         READ (IRD,20) (V(I,1),I=1,N)
      ELSE
         READ (IRD,20,END=30) ((V(I,J),I=1,N),J=1,N)
   20    FORMAT (8F10.6)
      ENDIF
   30 CONTINUE
      MO = MONE
      WRITE (ILST,*) ' BASIS SET IS ',CALC,' NUMBER OF AOS IS ',N
C     Write(ILST,*)' plotting MO number ',MO
      DO 40 J = 1, NAT
         XYZ(1,J) = C(J,1)
         XYZ(2,J) = C(J,2)
         XYZ(3,J) = C(J,3)
   40 CONTINUE
      X(1) = XMI
      Y(1) = YMI
      Z(1) = ZMI
      DO 50 I = 2, MXPTS
         X(I) = XINC+X(I-1)
         Y(I) = YINC+Y(I-1)
         Z(I) = ZINC+Z(I-1)
   50 CONTINUE
C
C       THESE ARE THE FIRST PART OF THE EQUATIONS FOR 1S->3P, CALCULATE
C       THEM ONCE AND ONCE ONLY, REAL POWERS (**X.XX) ARE VERY SLOW
C       COMPUTATIONS.  THESE "CONTANTS" NS AND NP WILL BE MULTIPLIED
C       BY THE ATOM DEPENDANT - R INDEPENDANT VALUES TO FORM ONE SINGLE
C       CONSTANT FOR MULTIPLICATION WITHIN THE CUBE LOOP.  THIS WILL
C       SPEED COMPUTATION CONSIDERABLY RATHER THAN DOING ALL OF THIS IN
C       THE LOOP.
C
C       FIRST THE "NS" ORBITAL "CONSTANTS"
C
      CNSTNS(1,1) = (A(1,1)**0.750E+0)*D(1,1)
      CNSTNS(2,1) = (A(1,2)**0.750E+0)*D(1,2)
      CNSTNS(3,1) = (A(1,3)**0.750E+0)*D(1,3)
      CNSTNS(1,2) = (A(2,1)**0.750E+0)*D(2,1)
      CNSTNS(2,2) = (A(2,2)**0.750E+0)*D(2,2)
      CNSTNS(3,2) = (A(2,3)**0.750E+0)*D(2,3)
      CNSTNS(1,3) = (A(3,1)**0.750E+0)*D(3,1)
      CNSTNS(2,3) = (A(3,2)**0.750E+0)*D(3,2)
      CNSTNS(3,3) = (A(3,3)**0.750E+0)*D(3,3)
C
C       NOW FOR THE "NP" ORBITALS (THE SECOND ARG IS THE QUANTUM NUMBER)
C       THE QUANTUM NUMBER RANGES FROM 2,3 SINCE THERE IS NO 1P ORBITAL
C
      CNSTNP(1,2) = (A(2,1)**1.250E+0)*D2P(1)
      CNSTNP(2,2) = (A(2,2)**1.250E+0)*D2P(2)
      CNSTNP(3,2) = (A(2,3)**1.250E+0)*D2P(3)
      CNSTNP(1,3) = (A(3,1)**1.250E+0)*D3P(1)
      CNSTNP(2,3) = (A(3,2)**1.250E+0)*D3P(2)
      CNSTNP(3,3) = (A(3,3)**1.250E+0)*D3P(3)
C
C       ZERO THE ORBITAL VALUE ARRAY
C
      DO 80 IZ = 1, MXPTS
         DO 70 IY = 1, MXPTS
            DO 60 IX = 1, MXPTS
               DENSIT(IX,IY,IZ) = 0.0E+0
   60       CONTINUE
   70    CONTINUE
   80 CONTINUE
C
C       INITIALIZE THE AO COUNTER AND LOOP OVER ATOMS, IAT IS THE ATOMIC
C       NUMBER
C
      M = 1
      DO 310 I = 1, NAT
         IAT = IAN(I)
C
C       COMPUTE XDEL,YDEL,AND ZDEL (I.E. DELTA X,Y, AND Z FROM THE ATOM
C       TO EACH POINT ON THE GRID.  ONLY MXPTS VALUES FOR EACH SINCE,
C       FOR INSTANCE, EVERY POINT ON A PARTICULAR XY PLANE IS THE SAME
C       DELTA Z VALUE FROM THE POINT.  THEREFORE YOU HAVE ONLY ONE VALUE
C       FOR THE ENTIRE PLANE FOR DELTA Z, INSTEAD OF (FOR MXPTS=51)
C       2601.  AGAIN, BY COMPUTING THIS HERE, RATHER THAN INSIDE THE
C       LOOP WE CUT DOWN THESE SUBTRACTIONS AND MULTIPLICATIONS BY
C       A FACTOR OF 2601 TO 1. THIS HAS A SUBSTANTIAL EFFECT ON THE
C       SPEED OF THE COMPUTATIONS.
C
         DO 90 IXYZ = 1, MXPTS
            XDEL(IXYZ) = X(IXYZ)-XYZ(1,I)
            XDELSQ(IXYZ) = XDEL(IXYZ)*XDEL(IXYZ)
   90    CONTINUE
         DO 100 IXYZ = 1, MXPTS
            YDEL(IXYZ) = Y(IXYZ)-XYZ(2,I)
            YDELSQ(IXYZ) = YDEL(IXYZ)*YDEL(IXYZ)
  100    CONTINUE
         DO 110 IXYZ = 1, MXPTS
            ZDEL(IXYZ) = Z(IXYZ)-XYZ(3,I)
            ZDELSQ(IXYZ) = ZDEL(IXYZ)*ZDEL(IXYZ)
  110    CONTINUE
C
C       FIRST THE H, HE ATOMS
C
         WRITE (ILST,'(2(A,I5))') 
     *       'PROCESSING ATOM NUMBER ',I,' ATOMIC NUMBER ',IAT
         IF (IAT.LE.2) THEN
C
C       NOW CALCULATE THE NORMALIZATION FACTORS WHICH ARE ATOM TYPE AND
C       QUANTUM NUMBER DEPENDANT.  MULTIPLY BY THE "CONSTANTS" FOR THE
C       PARTICULAR A.O. AND THE EIGENVECTOR FOR THAT A.O. IN THE M.O.
C       SINCE IT IS ALSO POSITION INDEPENDANT.  THE R*COS(THETA) (XDEL,
C       YDEL,AND ZDEL) WILL HAVE TO BE DONE INSIDE THE X,Y,AND Z LOOPS
C       RESPECTIVELY FOR ATOMS WHICH HAVE "P" ORBITALS, HERE WE DON'T
C       NEED TO WORRY.
C
            ZN = Z1(IAT)
            ZSQRT = SQRT(ZN)
            ZSQ(1) = ZN*ZN
C
C       ZN * SQRT(ZN) * (2.0E+0/PI)**0.750E+0
C
            RNORM = ZN*ZSQRT*0.712705470E+0
            VNORM(1) = RNORM*V(M,MO)*CNSTNS(1,1)
            VNORM(2) = RNORM*V(M,MO)*CNSTNS(2,1)
            VNORM(3) = RNORM*V(M,MO)*CNSTNS(3,1)
C
C       THERE IS ONLY THE 1S TO EVALUATE
C
            ANEG(1) = -A(1,1)*ZSQ(1)
            ANEG(2) = -A(1,2)*ZSQ(1)
            ANEG(3) = -A(1,3)*ZSQ(1)
C
C       THE EXPONENTIATIONS RELATED TO XDEL,YDEL,AND ZDEL
C       THEY WILL BE MULTIPLIED IN THE LOOP, RATHER THAN
C       MAXPTS**3 SEPARATE EXPONENTIATIONS OVER R.  THESE
C       9*MXPTS MAKE THE UNIQUE ONES FOR 1S. (3 GAUSSIANS*
C       3 CARTESIAN COORDS * MXPTS PLANES )
C
            DO 120 IXYZ = 1, MXPTS
               EXP1S(IXYZ,1) = GEXP(XDELSQ(IXYZ)*ANEG(1))*VNORM(1)
               EXP1S(IXYZ,2) = GEXP(XDELSQ(IXYZ)*ANEG(2))*VNORM(2)
               EXP1S(IXYZ,3) = GEXP(XDELSQ(IXYZ)*ANEG(3))*VNORM(3)
               EXP1S(IXYZ,4) = GEXP(YDELSQ(IXYZ)*ANEG(1))
               EXP1S(IXYZ,5) = GEXP(YDELSQ(IXYZ)*ANEG(2))
               EXP1S(IXYZ,6) = GEXP(YDELSQ(IXYZ)*ANEG(3))
               EXP1S(IXYZ,7) = GEXP(ZDELSQ(IXYZ)*ANEG(1))
               EXP1S(IXYZ,8) = GEXP(ZDELSQ(IXYZ)*ANEG(2))
               EXP1S(IXYZ,9) = GEXP(ZDELSQ(IXYZ)*ANEG(3))
  120       CONTINUE
            DO 160 IZ = 1, MXPTS
               DO 150 IY = 1, MXPTS
                  DO 140 IG = 1, 3
                     DO 130 IX = 1, MXPTS
C
C       CONTR IS THE SUM OF CONTRIBUTIONS OVER THIS YZ PLANE FOR THIS
C       ATOM.   WHEN FINISHED, SUM INTO THE ORBITAL VALUE ARRAY.
C
                        DENSIT(IX,IY,IZ) = DENSIT(IX,IY,IZ)+EXP1S(IY,IG+
     *                     3)*EXP1S(IZ,IG+6)*EXP1S(IX,IG)
  130                CONTINUE
  140             CONTINUE
  150          CONTINUE
  160       CONTINUE
            M = M+1
         ELSEIF (IAT.LE.10) THEN
            ZN = Z1(IAT)
            ZSQRT = SQRT(ZN)
            ZSQ(1) = ZN*ZN
C
C       ZN * SQRT(ZN) * (2.0E+0/PI)**0.750E+0
C
C       VNORM(1-3) THE THE THREE GAUSSION "CONSTANTS" FOR THE 1S
C       ORBITAL, EVERYTHING THAT IS INDEPENDANT OF R IS IN THERE.
C
            RNORM = ZN*ZSQRT*0.712705470E+0
            VNORM(1) = RNORM*V(M,MO)*CNSTNS(1,1)
            VNORM(2) = RNORM*V(M,MO)*CNSTNS(2,1)
            VNORM(3) = RNORM*V(M,MO)*CNSTNS(3,1)
C
C       CALC ZETA(2SP) SQUARED AND SQRT FOR CONSTANTS.
C
            ZN = Z2(IAT)
            ZSQRT = SQRT(ZN)
            ZSQ(2) = ZN*ZN
C
C       ZN*ZN * SQRT(ZN) * ((128.0E+0/PI**3)**0.250E+0)
C
C       VNORM(4-6) CORRESPONDS TO THE 3 GAUSSIANS FOR THE 2S ORBITAL
C       VNORM(7-9) "                            " FOR THE 2PX ""
C       VNORM(10-12) "" FOR THE 2PY ""
C       VNORM(13-15) "" FOR THE 2PZ ""
C       YOU NEED DIFFERENT "CONSTANTS" FOR EACH DUE TO THE DIFFERENT
C       VALUE OF THE WAVEFUNCTION FOR EACH ORBITAL. AT THE END WE WILL
C       HAVE A 3 SETS OF NUMBERS ALL MULTIPLIED BY THE EXPONENTIAL,
C       WE CAN ADD THEM FIRST, THEN CALC. AND MULT. BY THE EXP.
C
            RNORM = ZN*ZSQRT*0.712705470E+0
            VNORM(4) = RNORM*CNSTNS(1,2)*V(M+1,MO)
            VNORM(5) = RNORM*CNSTNS(2,2)*V(M+1,MO)
            VNORM(6) = RNORM*CNSTNS(3,2)*V(M+1,MO)
            RNORM = ZSQ(2)*ZSQRT*1.425410940E+0
            VNORM(7) = RNORM*CNSTNP(1,2)
            VNORM(8) = RNORM*CNSTNP(2,2)
            VNORM(9) = RNORM*CNSTNP(3,2)
            VNORM(10) = VNORM(7)*V(M+3,MO)
            VNORM(11) = VNORM(8)*V(M+3,MO)
            VNORM(12) = VNORM(9)*V(M+3,MO)
            VNORM(13) = VNORM(7)*V(M+4,MO)
            VNORM(14) = VNORM(8)*V(M+4,MO)
            VNORM(15) = VNORM(9)*V(M+4,MO)
            VNORM(7) = VNORM(7)*V(M+2,MO)
            VNORM(8) = VNORM(8)*V(M+2,MO)
            VNORM(9) = VNORM(9)*V(M+2,MO)
            DO 170 IXYZ = 1, MXPTS
               CNSTX(IXYZ,1) = VNORM(7)*XDEL(IXYZ)
               CNSTX(IXYZ,2) = VNORM(8)*XDEL(IXYZ)
               CNSTX(IXYZ,3) = VNORM(9)*XDEL(IXYZ)
               CNSTY(IXYZ,1) = VNORM(10)*YDEL(IXYZ)
               CNSTY(IXYZ,2) = VNORM(11)*YDEL(IXYZ)
               CNSTY(IXYZ,3) = VNORM(12)*YDEL(IXYZ)
               CNSTZ(IXYZ,1) = VNORM(13)*ZDEL(IXYZ)+VNORM(4)
               CNSTZ(IXYZ,2) = VNORM(14)*ZDEL(IXYZ)+VNORM(5)
               CNSTZ(IXYZ,3) = VNORM(15)*ZDEL(IXYZ)+VNORM(6)
  170       CONTINUE
C
C       EVALUATE 2S, 2P, AND 1S
C
C       MINUS ALPHA FOR 1S:
C
            ANEG(1) = -A(1,1)*ZSQ(1)
            ANEG(2) = -A(1,2)*ZSQ(1)
            ANEG(3) = -A(1,3)*ZSQ(1)
C
C       MINUS ALPHA FOR 2SP:
C
            ANEG(4) = -A(2,1)*ZSQ(2)
            ANEG(5) = -A(2,2)*ZSQ(2)
            ANEG(6) = -A(2,3)*ZSQ(2)
C
C       PRECOMPUTE EXP(-A*Z**2*DELO**2) WHERE DELO**2 IS
C       DELTA-X**2, DELTA-Y**2, AND DELTA-Z**2
C
            DO 180 IXYZ = 1, MXPTS
               EXP1S(IXYZ,1) = GEXP(XDELSQ(IXYZ)*ANEG(1))*VNORM(1)
               EXP1S(IXYZ,2) = GEXP(XDELSQ(IXYZ)*ANEG(2))*VNORM(2)
               EXP1S(IXYZ,3) = GEXP(XDELSQ(IXYZ)*ANEG(3))*VNORM(3)
               EXP1S(IXYZ,4) = GEXP(YDELSQ(IXYZ)*ANEG(1))
               EXP1S(IXYZ,5) = GEXP(YDELSQ(IXYZ)*ANEG(2))
               EXP1S(IXYZ,6) = GEXP(YDELSQ(IXYZ)*ANEG(3))
               EXP1S(IXYZ,7) = GEXP(ZDELSQ(IXYZ)*ANEG(1))
               EXP1S(IXYZ,8) = GEXP(ZDELSQ(IXYZ)*ANEG(2))
               EXP1S(IXYZ,9) = GEXP(ZDELSQ(IXYZ)*ANEG(3))
  180       CONTINUE
            DO 190 IXYZ = 1, MXPTS
               EXP2SP(IXYZ,1) = GEXP(XDELSQ(IXYZ)*ANEG(4))
               EXP2SP(IXYZ,2) = GEXP(XDELSQ(IXYZ)*ANEG(5))
               EXP2SP(IXYZ,3) = GEXP(XDELSQ(IXYZ)*ANEG(6))
               EXP2SP(IXYZ,4) = GEXP(YDELSQ(IXYZ)*ANEG(4))
               EXP2SP(IXYZ,5) = GEXP(YDELSQ(IXYZ)*ANEG(5))
               EXP2SP(IXYZ,6) = GEXP(YDELSQ(IXYZ)*ANEG(6))
               EXP2SP(IXYZ,7) = GEXP(ZDELSQ(IXYZ)*ANEG(4))
               EXP2SP(IXYZ,8) = GEXP(ZDELSQ(IXYZ)*ANEG(5))
               EXP2SP(IXYZ,9) = GEXP(ZDELSQ(IXYZ)*ANEG(6))
  190       CONTINUE
C
C       LOOP OVER THE "CUBE" Z,Y, AND X
C
            DO 230 IZ = 1, MXPTS
               DO 220 IY = 1, MXPTS
                  DO 210 IG = 1, 3
                     DO 200 IX = 1, MXPTS
C
C       SUM THE ORBITAL CONTRIBUTIONS INTO THE ORBITAL VALUE ARRAY
C
C       FIRST THE 3 GAUSSIANS FOR THE 1S:
C
                        DENSIT(IX,IY,IZ) = DENSIT(IX,IY,IZ)+EXP1S(IX,IG)
     *                     *EXP1S(IY,IG+3)*EXP1S(IZ,IG+6)
C
C       NEXT, THE 3 FOR THE 2SP:
C
                        DENSIT(IX,IY,IZ) = DENSIT(IX,IY,IZ)+(CNSTY(IY,IG
     *                     )+CNSTZ(IZ,IG)+CNSTX(IX,IG))*EXP2SP(IX,IG)*
     *                     EXP2SP(IY,IG+3)*EXP2SP(IZ,IG+6)
  200                CONTINUE
  210             CONTINUE
  220          CONTINUE
  230       CONTINUE
            M = M+5
         ELSEIF (IAT.LE.18) THEN
            ZN = Z1(IAT)
            ZSQRT = SQRT(ZN)
            ZSQ(1) = ZN*ZN
C
C       ZN * SQRT(ZN) * (2.0E+0/PI)**0.750E+0
C
C       VNORM(1-3) THE THE THREE GAUSSION "CONSTANTS" FOR THE 1S
C       ORBITAL, EVERYTHING THAT IS INDEPENDANT OF R IS IN THERE.
C
            RNORM = ZN*ZSQRT*0.712705470E+0
            VNORM(1) = RNORM*V(M,MO)*CNSTNS(1,1)
            VNORM(2) = RNORM*V(M,MO)*CNSTNS(2,1)
            VNORM(3) = RNORM*V(M,MO)*CNSTNS(3,1)
C
C       CALC ZETA(2SP) SQUARED AND SQRT FOR CONSTANTS.
C
            ZN = Z2(IAT)
            ZSQRT = SQRT(ZN)
            ZSQ(2) = ZN*ZN
C
C       ZN*ZN * SQRT(ZN) * ((128.0E+0/PI**3)**0.250E+0)
C
C       VNORM(4-6) CORRESPONDS TO THE 3 GAUSSIANS FOR THE 2S ORBITAL
C       VNORM(7-9) "                            " FOR THE 2PX ""
C       VNORM(10-12) "" FOR THE 2PY ""
C       VNORM(13-15) "" FOR THE 2PZ ""
C       YOU NEED DIFFERENT "CONSTANTS" FOR EACH DUE TO THE DIFFERENT
C       VALUE OF THE WAVEFUNCTION FOR EACH ORBITAL. AT THE END WE WILL
C       HAVE A 3 SETS OF NUMBERS ALL MULTIPLIED BY THE EXPONENTIAL,
C       WE CAN ADD THEM FIRST, THEN CALC. AND MULT. BY THE EXP.
C
            RNORM = ZN*ZSQRT*0.712705470E+0
            VNORM(4) = RNORM*CNSTNS(1,2)*V(M+1,MO)
            VNORM(5) = RNORM*CNSTNS(2,2)*V(M+1,MO)
            VNORM(6) = RNORM*CNSTNS(3,2)*V(M+1,MO)
            RNORM = ZSQ(2)*ZSQRT*1.425410940E+0
            VNORM(7) = RNORM*CNSTNP(1,2)
            VNORM(8) = RNORM*CNSTNP(2,2)
            VNORM(9) = RNORM*CNSTNP(3,2)
            VNORM(10) = VNORM(7)*V(M+3,MO)
            VNORM(11) = VNORM(8)*V(M+3,MO)
            VNORM(12) = VNORM(9)*V(M+3,MO)
            VNORM(13) = VNORM(7)*V(M+4,MO)
            VNORM(14) = VNORM(8)*V(M+4,MO)
            VNORM(15) = VNORM(9)*V(M+4,MO)
            VNORM(7) = VNORM(7)*V(M+2,MO)
            VNORM(8) = VNORM(8)*V(M+2,MO)
            VNORM(9) = VNORM(9)*V(M+2,MO)
C
C       NOW THE 3RD ROW STUFF
C
            ZN = Z3(IAT)
            ZSQRT = SQRT(ZN)
            ZSQ(3) = ZN*ZN
C
C       ZN*ZN * SQRT(ZN) * ((128.0E+0/PI**3)**0.250E+0)
C
C       VNORM(16-18) CORRESPONDS TO THE 3 GAUSSIANS FOR THE 2S ORBITAL
C       VNORM(19-21) "                            " FOR THE 2PX ""
C       VNORM(22-24) "" FOR THE 2PY ""
C       VNORM(25-27) "" FOR THE 2PZ ""
C       YOU NEED DIFFERENT "CONSTANTS" FOR EACH DUE TO THE DIFFERENT
C       VALUE OF THE WAVEFUNCTION FOR EACH ORBITAL. AT THE END WE WILL
C       HAVE A 3 SETS OF NUMBERS ALL MULTIPLIED BY THE EXPONENTIAL,
C       WE CAN ADD THEM FIRST, THEN CALC. AND MULT. BY THE EXP.
C
            RNORM = ZN*ZSQRT*0.712705470E+0
            VNORM(16) = RNORM*CNSTNS(1,3)*V(M+5,MO)
            VNORM(17) = RNORM*CNSTNS(2,3)*V(M+5,MO)
            VNORM(18) = RNORM*CNSTNS(3,3)*V(M+5,MO)
            RNORM = ZSQ(3)*ZSQRT*1.425410940E+0
            VNORM(19) = RNORM*CNSTNP(1,3)
            VNORM(20) = RNORM*CNSTNP(2,3)
            VNORM(21) = RNORM*CNSTNP(3,3)
            VNORM(22) = VNORM(19)*V(M+7,MO)
            VNORM(23) = VNORM(20)*V(M+7,MO)
            VNORM(24) = VNORM(21)*V(M+7,MO)
            VNORM(25) = VNORM(19)*V(M+8,MO)
            VNORM(26) = VNORM(20)*V(M+8,MO)
            VNORM(27) = VNORM(21)*V(M+8,MO)
            VNORM(19) = VNORM(19)*V(M+6,MO)
            VNORM(20) = VNORM(20)*V(M+6,MO)
            VNORM(21) = VNORM(21)*V(M+6,MO)
            DO 240 IXYZ = 1, MXPTS
               CNSTX(IXYZ,1) = VNORM(7)*XDEL(IXYZ)
               CNSTX(IXYZ,2) = VNORM(8)*XDEL(IXYZ)
               CNSTX(IXYZ,3) = VNORM(9)*XDEL(IXYZ)
               CNSTY(IXYZ,1) = VNORM(10)*YDEL(IXYZ)
               CNSTY(IXYZ,2) = VNORM(11)*YDEL(IXYZ)
               CNSTY(IXYZ,3) = VNORM(12)*YDEL(IXYZ)
               CNSTZ(IXYZ,1) = VNORM(13)*ZDEL(IXYZ)+VNORM(4)
               CNSTZ(IXYZ,2) = VNORM(14)*ZDEL(IXYZ)+VNORM(5)
               CNSTZ(IXYZ,3) = VNORM(15)*ZDEL(IXYZ)+VNORM(6)
               CNSTX(IXYZ,4) = VNORM(19)*XDEL(IXYZ)
               CNSTX(IXYZ,5) = VNORM(20)*XDEL(IXYZ)
               CNSTX(IXYZ,6) = VNORM(21)*XDEL(IXYZ)
               CNSTY(IXYZ,4) = VNORM(22)*YDEL(IXYZ)
               CNSTY(IXYZ,5) = VNORM(23)*YDEL(IXYZ)
               CNSTY(IXYZ,6) = VNORM(24)*YDEL(IXYZ)
               CNSTZ(IXYZ,4) = VNORM(25)*ZDEL(IXYZ)+VNORM(16)
               CNSTZ(IXYZ,5) = VNORM(26)*ZDEL(IXYZ)+VNORM(17)
               CNSTZ(IXYZ,6) = VNORM(27)*ZDEL(IXYZ)+VNORM(18)
  240       CONTINUE
C
C       EVALUATE 2S, 2P, AND 1S
C
C       MINUS ALPHA FOR 1S:
C
            ANEG(1) = -A(1,1)*ZSQ(1)
            ANEG(2) = -A(1,2)*ZSQ(1)
            ANEG(3) = -A(1,3)*ZSQ(1)
C
C       MINUS ALPHA FOR 2SP:
C
            ANEG(4) = -A(2,1)*ZSQ(2)
            ANEG(5) = -A(2,2)*ZSQ(2)
            ANEG(6) = -A(2,3)*ZSQ(2)
C
C       MINUS ALPHA FOR 3SP:
C
            ANEG(7) = -A(3,1)*ZSQ(3)
            ANEG(8) = -A(3,2)*ZSQ(3)
            ANEG(9) = -A(3,3)*ZSQ(3)
C
C       PRECOMPUTE EXP(-A*Z**2*DELO**2) WHERE DELO**2 IS
C       DELTA-X**2, DELTA-Y**2, AND DELTA-Z**2
C
            DO 250 IXYZ = 1, MXPTS
               EXP1S(IXYZ,1) = GEXP(XDELSQ(IXYZ)*ANEG(1))*VNORM(1)
               EXP1S(IXYZ,2) = GEXP(XDELSQ(IXYZ)*ANEG(2))*VNORM(2)
               EXP1S(IXYZ,3) = GEXP(XDELSQ(IXYZ)*ANEG(3))*VNORM(3)
               EXP1S(IXYZ,4) = GEXP(YDELSQ(IXYZ)*ANEG(1))
               EXP1S(IXYZ,5) = GEXP(YDELSQ(IXYZ)*ANEG(2))
               EXP1S(IXYZ,6) = GEXP(YDELSQ(IXYZ)*ANEG(3))
               EXP1S(IXYZ,7) = GEXP(ZDELSQ(IXYZ)*ANEG(1))
               EXP1S(IXYZ,8) = GEXP(ZDELSQ(IXYZ)*ANEG(2))
               EXP1S(IXYZ,9) = GEXP(ZDELSQ(IXYZ)*ANEG(3))
  250       CONTINUE
            DO 260 IXYZ = 1, MXPTS
               EXP2SP(IXYZ,1) = GEXP(XDELSQ(IXYZ)*ANEG(4))
               EXP2SP(IXYZ,2) = GEXP(XDELSQ(IXYZ)*ANEG(5))
               EXP2SP(IXYZ,3) = GEXP(XDELSQ(IXYZ)*ANEG(6))
               EXP2SP(IXYZ,4) = GEXP(YDELSQ(IXYZ)*ANEG(4))
               EXP2SP(IXYZ,5) = GEXP(YDELSQ(IXYZ)*ANEG(5))
               EXP2SP(IXYZ,6) = GEXP(YDELSQ(IXYZ)*ANEG(6))
               EXP2SP(IXYZ,7) = GEXP(ZDELSQ(IXYZ)*ANEG(4))
               EXP2SP(IXYZ,8) = GEXP(ZDELSQ(IXYZ)*ANEG(5))
               EXP2SP(IXYZ,9) = GEXP(ZDELSQ(IXYZ)*ANEG(6))
               EXP3SP(IXYZ,1) = GEXP(XDELSQ(IXYZ)*ANEG(7))
               EXP3SP(IXYZ,2) = GEXP(XDELSQ(IXYZ)*ANEG(8))
               EXP3SP(IXYZ,3) = GEXP(XDELSQ(IXYZ)*ANEG(9))
               EXP3SP(IXYZ,4) = GEXP(YDELSQ(IXYZ)*ANEG(7))
               EXP3SP(IXYZ,5) = GEXP(YDELSQ(IXYZ)*ANEG(8))
               EXP3SP(IXYZ,6) = GEXP(YDELSQ(IXYZ)*ANEG(9))
               EXP3SP(IXYZ,7) = GEXP(ZDELSQ(IXYZ)*ANEG(7))
               EXP3SP(IXYZ,8) = GEXP(ZDELSQ(IXYZ)*ANEG(8))
               EXP3SP(IXYZ,9) = GEXP(ZDELSQ(IXYZ)*ANEG(9))
  260       CONTINUE
C
C       LOOP OVER THE "CUBE" Z,Y, AND X
C
            DO 300 IZ = 1, MXPTS
               DO 290 IY = 1, MXPTS
                  DO 280 IG = 1, 3
                     DO 270 IX = 1, MXPTS
C
C       SUM THE ORBITAL CONTRIBUTIONS INTO THE ORBITAL VALUE ARRAY
C
C       FIRST THE 3 GAUSSIANS FOR THE 1S:
C
                        DENSIT(IX,IY,IZ) = DENSIT(IX,IY,IZ)+EXP1S(IX,IG)
     *                     *EXP1S(IY,IG+3)*EXP1S(IZ,IG+6)
C
C       NEXT, THE 3 FOR THE 2SP:
C
                        DENSIT(IX,IY,IZ) = DENSIT(IX,IY,IZ)+(CNSTX(IX,IG
     *                     )+CNSTY(IY,IG)+CNSTZ(IZ,IG))*EXP2SP(IX,IG)*
     *                     EXP2SP(IY,IG+3)*EXP2SP(IZ,IG+6)
C
C       NEXT, THE 3 FOR THE 3SP:
C
                        DENSIT(IX,IY,IZ) = DENSIT(IX,IY,IZ)+(CNSTY(IY,IG
     *                     +3)+CNSTZ(IZ,IG+3)+CNSTX(IX,IG+3))*EXP3SP(IX,
     *                     IG)*EXP3SP(IY,IG+3)*EXP3SP(IZ,IG+6)
  270                CONTINUE
  280             CONTINUE
  290          CONTINUE
  300       CONTINUE
            M = M+9
         ENDIF
  310 CONTINUE
      RETURN
      END
C
C
      SUBROUTINE MOSEMI
      IMPLICIT REAL (A-H,O-Z)
      PARAMETER (MXPTS=51)
      PARAMETER (MAXATM=50)
      PARAMETER (MAXORB=200)
C
C       REEPLACES ORIGINAL ROUTINE TO EVALUATE RADIAL PARTS OF STO-3G
C       WAVEFUNCTIONS AS WELL AS THE MAINLINE CODE FOR THE REST.
C       REF FROM ORIGINAL CODE (FUNCTION AO):
C
C       ATOMS H TO AR ARE HANDLED ACCORDING TO J CHEM PHYS 51, 2657
C       (1969), 52, 2769 (1970).
C       W.L. JORGENSEN - MARCH,1976, JULY, 1977.
C
C       REWRITTEN INTO A STAND-ALONE SUBROUTINE FOR CALCULATING ORBITAL
C       VALUES. REWORKED TO REDUCE REDUNDANT COMPUTATION OF POWERS AND
C       SQUARE ROOTS, AS WELL AS PUTTING INTO AN EASILY VECTORIZABLE
C       FORM FOR VECTOR MACHINES ( WE HAVE A CYBER 205 AT PURDUE, IT
C       WILL ALSO WORK AS WELL ON OTHER MACHINES).
C
C       DAN SEVERANCE 12/10/87
C
C       AFTER DISCUSSION WITH JIM BRIGGS WHO HAS BEEN USING THIS PROGRAM
C       ON DR. JORGENSEN'S GOULD FOR THE LAST FEW WEEKS; THE AO
C       COMPUTATION AND WAVEFUNCTION READ WAS MOVED FROM THE MAIN LINE
C       TO THE RESPECTIVE SUBROUTINES.  WHEN A WAVEFUNCTION IS ADDED,
C       ONLY THE SUBROUTINE SHOULD NEED SIGNIFICANT MODIFICATION, THE
C       MAINLINE SHOULD ONLY NEED TO HAVE THE CALL ADDED.
C       DAN SEVERANCE 1/17/88
C
      COMMON /SPLIT/ XMI,YMI,ZMI,XINC,YINC,ZINC,MONE,MOLAST,NAT,
     *   IONEMO
      COMMON /DENS/ DENSIT(MXPTS,MXPTS,MXPTS),V(MAXORB,MAXORB),C(MAXATM,
     *   3),OCMO(MAXORB),IAN(MAXATM)
      COMMON /IO/ IRD,ILST
      DIMENSION XDEL(MXPTS),YDEL(MXPTS),ZDEL(MXPTS),EXPP(MXPTS,9)
      DIMENSION XDELSQ(MXPTS),YDELSQ(MXPTS),ZDELSQ(MXPTS)
      DIMENSION XYZ(3,50),ZSQ(3),EXPS(MXPTS,9),ANEG(6)
      DIMENSION X(MXPTS),Y(MXPTS),Z(MXPTS),CNSTX(MXPTS,3),CNSTY(MXPTS,3)
      DIMENSION Z1(50),Z2(50),Z3(50),A(3,3),D(3,3),D2P(3),D3P(3)
      DIMENSION CNSTNS(3,3),CNSTNP(3,3),VNORM(12),CNSTZ(MXPTS,3)
C
C       12/10/87 DAN SEVERANCE
C
      DATA A(1,1),A(1,2),A(1,3) / 1.098180E-1,4.057710E-1,2.227660E+0 /
      DATA A(2,1),A(2,2),A(2,3) / 7.513860E-2,2.310310E-1,9.942030E-1 /
      DATA A(3,1),A(3,2),A(3,3) / 5.272660E-2,1.347150E-1,4.828540E-1 /
      DATA D(1,1),D(1,2),D(1,3) / 4.446350E-1,5.353280E-1,1.543290E-1 /
      DATA D(2,1),D(2,2),D(2,3) / 7.001150E-1,3.995130E-1,-9.996720E-2 /
      DATA D(3,1),D(3,2),D(3,3) / 9.003980E-1,2.255950E-1,-2.196200E-1 /
      DATA D2P / 3.919570E-1,6.076840E-1,1.559160E-1 /
      DATA D3P / 4.620010E-1,5.951670E-1,1.058760E-2 /
C
C       NOW INITIALIZE ALL OF THE NON-"R" DEPENDANT VALUES RATHER THAN
C       RECOMPUTING THEM MXPTS**3 TIMES IN THE Z,Y,X LOOP OVER THE
C       ORBITAL VALUE MATRIX (DENSIT)
C
C       THESE VALUES ARE ALSO INDEPENDANT OF ATOM TYPE, ONLY DEPENDANT
C       ON THE ROW OF THE PERIODIC TABLE AND WHETHER IT IS "S" OR "P"
C       INITIALIZE THEM HERE AND ACCESS THEM WITHIN THE NAT LOOP, BEFORE
C       ENTERING THE LOOP OVER THE GRID (CUBE) POINTS.
C
      WRITE (ILST,*) ' EVALUATING THE SEMIEMPIRICAL WAVEFUNCTION'
      NAO = 0
      DO 10 I = 1, NAT
         NAO = NAO+1
         IF (IAN(I).GT.2) THEN
            NAO = NAO+3
            IF (IAN(I).GT.18) THEN
               WRITE (ILST,*) ' ATOMIC NUMBERS > 18 NOT YET IMPLEMENTED'
            ENDIF
         ENDIF
   10 CONTINUE
      WRITE (ILST,*) NAO,' WAVEFUNCTIONS TO BE PROCESSED '
C
C       READ IN EIGENVECTORS
C       IT IS ASSUMED THAT THE EIGENVECTORS HAVE BEEN NORMALIZED TO 1
C       ELECTRON WITH THE OVERLAP MATRIX INCLUDED. FOR SEMI-EMPIRICAL
C       WAVEFUNCTIONS, THIS REQUIRES THE LOWDIN TRANSFORMATION AS
C       IMPLEMENTED IN THE MOPAC ROUTINE MULT.
C
C       THE ZETA VALUES ARE ASSUMED TO BE AT THE END OF THE DATA, ONE
C       FOR EACH ATOM NUMBER. (8F10.6)
C
      IF (IONEMO.NE.0) THEN
         READ (IRD,40) (V(I,1),I=1,NAO)
         DO 20 I = 1, NAT
            READ (IRD,40) Z1(I),Z2(I),Z3(I)
   20    CONTINUE
      ELSE
         READ (IRD,40,END=50) ((V(I,J),I=1,NAO),J=1,NAO)
         DO 30 I = 1, NAT
            READ (IRD,40) Z1(I),Z2(I),Z3(I)
   30    CONTINUE
   40    FORMAT (8F10.6)
      ENDIF
   50 CONTINUE
      MO = MONE
      DO 60 J = 1, NAT
         XYZ(1,J) = C(J,1)
         XYZ(2,J) = C(J,2)
         XYZ(3,J) = C(J,3)
   60 CONTINUE
      X(1) = XMI
      Y(1) = YMI
      Z(1) = ZMI
      DO 70 I = 2, MXPTS
         X(I) = XINC+X(I-1)
         Y(I) = YINC+Y(I-1)
         Z(I) = ZINC+Z(I-1)
   70 CONTINUE
C
C       THESE ARE THE FIRST PART OF THE EQUATIONS FOR 1S->3P, CALCULATE
C       THEM ONCE AND ONCE ONLY, REAL POWERS (**X.XX) ARE VERY SLOW
C       COMPUTATIONS.  THESE "CONTANTS" NS AND NP WILL BE MULTIPLIED
C       BY THE ATOM DEPENDANT - R INDEPENDANT VALUES TO FORM ONE SINGLE
C       CONSTANT FOR MULTIPLICATION WITHIN THE CUBE LOOP.  THIS WILL
C       SPEED COMPUTATION CONSIDERABLY RATHER THAN DOING ALL OF THIS IN
C       THE LOOP.
C
C       FIRST THE "NS" ORBITAL "CONSTANTS"
C
      CNSTNS(1,1) = (A(1,1)**0.750E+0)*D(1,1)
      CNSTNS(2,1) = (A(1,2)**0.750E+0)*D(1,2)
      CNSTNS(3,1) = (A(1,3)**0.750E+0)*D(1,3)
      CNSTNS(1,2) = (A(2,1)**0.750E+0)*D(2,1)
      CNSTNS(2,2) = (A(2,2)**0.750E+0)*D(2,2)
      CNSTNS(3,2) = (A(2,3)**0.750E+0)*D(2,3)
      CNSTNS(1,3) = (A(3,1)**0.750E+0)*D(3,1)
      CNSTNS(2,3) = (A(3,2)**0.750E+0)*D(3,2)
      CNSTNS(3,3) = (A(3,3)**0.750E+0)*D(3,3)
C
C       NOW FOR THE "NP" ORBITALS (THE SECOND ARG IS THE QUANTUM NUMBER)
C       THE QUANTUM NUMBER RANGES FROM 2,3 SINCE THERE IS NO 1P ORBITAL
C
      CNSTNP(1,2) = (A(2,1)**1.250E+0)*D2P(1)
      CNSTNP(2,2) = (A(2,2)**1.250E+0)*D2P(2)
      CNSTNP(3,2) = (A(2,3)**1.250E+0)*D2P(3)
      CNSTNP(1,3) = (A(3,1)**1.250E+0)*D3P(1)
      CNSTNP(2,3) = (A(3,2)**1.250E+0)*D3P(2)
      CNSTNP(3,3) = (A(3,3)**1.250E+0)*D3P(3)
C
C       ZERO THE ORBITAL VALUE ARRAY
C
      DO 100 IZ = 1, MXPTS
         DO 90 IY = 1, MXPTS
            DO 80 IX = 1, MXPTS
               DENSIT(IX,IY,IZ) = 0.0E+0
   80       CONTINUE
   90    CONTINUE
  100 CONTINUE
C
C       INITIALIZE THE AO COUNTER AND LOOP OVER ATOMS, IAT IS THE ATOMIC
C       NUMBER
C
      M = 1
      DO 330 I = 1, NAT
         IAT = IAN(I)
C
C       COMPUTE XDEL,YDEL,AND ZDEL (I.E. DELTA X,Y, AND Z FROM THE ATOM
C       TO EACH POINT ON THE GRID.  ONLY MXPTS VALUES FOR EACH SINCE,
C       FOR INSTANCE, EVERY POINT ON A PARTICULAR XY PLANE IS THE SAME
C       DELTA Z VALUE FROM THE POINT.  THEREFORE YOU HAVE ONLY ONE VALUE
C       FOR THE ENTIRE PLANE FOR DELTA Z, INSTEAD OF (FOR MXPTS=51)
C       2601.  AGAIN, BY COMPUTING THIS HERE, RATHER THAN INSIDE THE
C       LOOP WE CUT DOWN THESE SUBTRACTIONS AND MULTIPLICATIONS BY
C       A FACTOR OF 2601 TO 1. THIS HAS A SUBSTANTIAL EFFECT ON THE
C       SPEED OF THE COMPUTATIONS.
C
         DO 110 IXYZ = 1, MXPTS
            XDEL(IXYZ) = X(IXYZ)-XYZ(1,I)
            XDELSQ(IXYZ) = XDEL(IXYZ)*XDEL(IXYZ)
  110    CONTINUE
         DO 120 IXYZ = 1, MXPTS
            YDEL(IXYZ) = Y(IXYZ)-XYZ(2,I)
            YDELSQ(IXYZ) = YDEL(IXYZ)*YDEL(IXYZ)
  120    CONTINUE
         DO 130 IXYZ = 1, MXPTS
            ZDEL(IXYZ) = Z(IXYZ)-XYZ(3,I)
            ZDELSQ(IXYZ) = ZDEL(IXYZ)*ZDEL(IXYZ)
  130    CONTINUE
C
C       FIRST THE H, HE ATOMS
C
         WRITE (ILST,'(2(A,I5))') 
     *       'PROCESSING ATOM NUMBER ',I,' ATOMIC NUMBER ',IAT
         IF (IAT.LE.2) THEN
C
C       NOW CALCULATE THE NORMALIZATION FACTORS WHICH ARE ATOM TYPE AND
C       QUANTUM NUMBER DEPENDANT.  MULTIPLY BY THE "CONSTANTS" FOR THE
C       PARTICULAR A.O. AND THE EIGENVECTOR FOR THAT A.O. IN THE M.O.
C       SINCE IT IS ALSO POSITION INDEPENDANT.  THE R*COS(THETA) (XDEL,
C       YDEL,AND ZDEL) WILL HAVE TO BE DONE INSIDE THE X,Y,AND Z LOOPS
C       RESPECTIVELY FOR ATOMS WHICH HAVE "P" ORBITALS, HERE WE DON'T
C       NEED TO WORRY.
C
            ZN = Z1(I)
            ZSQRT = SQRT(ZN)
            ZSQ(1) = ZN*ZN
C
C       ZN * SQRT(ZN) * (2.0E+0/PI)**0.750E+0
C
            RNORM = ZN*ZSQRT*0.712705470E+0
            VNORM(1) = RNORM*V(M,MO)*CNSTNS(1,1)
            VNORM(2) = RNORM*V(M,MO)*CNSTNS(2,1)
            VNORM(3) = RNORM*V(M,MO)*CNSTNS(3,1)
C
C       THERE IS ONLY THE 1S TO EVALUATE
C
            ANEG(1) = -A(1,1)*ZSQ(1)
            ANEG(2) = -A(1,2)*ZSQ(1)
            ANEG(3) = -A(1,3)*ZSQ(1)
C
C       THE EXPONENTIATIONS RELATED TO XDEL,YDEL,AND ZDEL
C       THEY WILL BE MULTIPLIED IN THE LOOP, RATHER THAN
C       MAXPTS**3 SEPARATE EXPONENTIATIONS OVER R.  THESE
C       9*MXPTS MAKE THE UNIQUE ONES FOR 1S. (3 GAUSSIANS*
C       3 CARTESIAN COORDS * MXPTS PLANES )
C
            DO 140 IXYZ = 1, MXPTS
               EXPS(IXYZ,1) = GEXP(XDELSQ(IXYZ)*ANEG(1))*VNORM(1)
               EXPS(IXYZ,2) = GEXP(XDELSQ(IXYZ)*ANEG(2))*VNORM(2)
               EXPS(IXYZ,3) = GEXP(XDELSQ(IXYZ)*ANEG(3))*VNORM(3)
               EXPS(IXYZ,4) = GEXP(YDELSQ(IXYZ)*ANEG(1))
               EXPS(IXYZ,5) = GEXP(YDELSQ(IXYZ)*ANEG(2))
               EXPS(IXYZ,6) = GEXP(YDELSQ(IXYZ)*ANEG(3))
               EXPS(IXYZ,7) = GEXP(ZDELSQ(IXYZ)*ANEG(1))
               EXPS(IXYZ,8) = GEXP(ZDELSQ(IXYZ)*ANEG(2))
               EXPS(IXYZ,9) = GEXP(ZDELSQ(IXYZ)*ANEG(3))
  140       CONTINUE
            DO 180 IZ = 1, MXPTS
               DO 170 IY = 1, MXPTS
                  DO 160 IG = 1, 3
                     DO 150 IX = 1, MXPTS
C
C       CONTR IS THE SUM OF CONTRIBUTIONS OVER THIS YZ PLANE FOR THIS
C       ATOM.   WHEN FINISHED, SUM INTO THE ORBITAL VALUE ARRAY.
C
                        DENSIT(IX,IY,IZ) = DENSIT(IX,IY,IZ)+EXPS(IY,IG+3
     *                     )*EXPS(IZ,IG+6)*EXPS(IX,IG)
  150                CONTINUE
  160             CONTINUE
  170          CONTINUE
  180       CONTINUE
            M = M+1
         ELSEIF (IAT.LE.10) THEN
C
C       CALC ZETA(2S) SQUARED AND SQRT FOR CONSTANTS.
C
            ZN = Z1(I)
            ZSQRT = SQRT(ZN)
            ZSQ(1) = ZN*ZN
C
C       ZN * SQRT(ZN) * (2.0E+0/PI)**0.750E+0
C
C       VNORM(1-3) THE THE THREE GAUSSION "CONSTANTS" FOR THE 2S
C       ORBITAL, EVERYTHING THAT IS INDEPENDANT OF R IS IN THERE.
C
            RNORM = ZN*ZSQRT*0.712705470E+0
            VNORM(1) = RNORM*V(M,MO)*CNSTNS(1,2)
            VNORM(2) = RNORM*V(M,MO)*CNSTNS(2,2)
            VNORM(3) = RNORM*V(M,MO)*CNSTNS(3,2)
C
C       CALC ZETA(2P) SQUARED AND SQRT FOR CONSTANTS.
C
            ZN = Z2(I)
            ZSQRT = SQRT(ZN)
            ZSQ(2) = ZN*ZN
C
C       ZN*ZN * SQRT(ZN) * ((128.0E+0/PI**3)**0.250E+0)
C
C       VNORM(1-3) CORRESPONDS TO THE 3 GAUSSIANS FOR THE 2S ORBITAL
C       VNORM(4-6) "                            " FOR THE 2PX ""
C       VNORM(7-9) "" FOR THE 2PY ""
C       VNORM(10-12) "" FOR THE 2PZ ""
C       YOU NEED DIFFERENT "CONSTANTS" FOR EACH DUE TO THE DIFFERENT
C       VALUE OF THE WAVEFUNCTION FOR EACH ORBITAL. AT THE END WE WILL
C       HAVE A 3 SETS OF NUMBERS ALL MULTIPLIED BY THE EXPONENTIAL,
C       WE CAN ADD THEM FIRST, THEN CALC. AND MULT. BY THE EXP.
C
            RNORM = ZSQ(2)*ZSQRT*1.425410940E+0
            AOP1 = RNORM*CNSTNP(1,2)
            AOP2 = RNORM*CNSTNP(2,2)
            AOP3 = RNORM*CNSTNP(3,2)
            VNORM(4) = AOP1*V(M+1,MO)
            VNORM(5) = AOP2*V(M+1,MO)
            VNORM(6) = AOP3*V(M+1,MO)
            VNORM(7) = AOP1*V(M+2,MO)
            VNORM(8) = AOP2*V(M+2,MO)
            VNORM(9) = AOP3*V(M+2,MO)
            VNORM(10) = AOP1*V(M+3,MO)
            VNORM(11) = AOP2*V(M+3,MO)
            VNORM(12) = AOP3*V(M+3,MO)
            DO 190 IXYZ = 1, MXPTS
               CNSTX(IXYZ,1) = VNORM(4)*XDEL(IXYZ)
               CNSTX(IXYZ,2) = VNORM(5)*XDEL(IXYZ)
               CNSTX(IXYZ,3) = VNORM(6)*XDEL(IXYZ)
               CNSTY(IXYZ,1) = VNORM(7)*YDEL(IXYZ)
               CNSTY(IXYZ,2) = VNORM(8)*YDEL(IXYZ)
               CNSTY(IXYZ,3) = VNORM(9)*YDEL(IXYZ)
               CNSTZ(IXYZ,1) = VNORM(10)*ZDEL(IXYZ)
               CNSTZ(IXYZ,2) = VNORM(11)*ZDEL(IXYZ)
               CNSTZ(IXYZ,3) = VNORM(12)*ZDEL(IXYZ)
  190       CONTINUE
C
C       EVALUATE 2S AND 2P
C
C       MINUS ALPHA FOR 2S:
C
            ANEG(1) = -A(2,1)*ZSQ(1)
            ANEG(2) = -A(2,2)*ZSQ(1)
            ANEG(3) = -A(2,3)*ZSQ(1)
C
C       MINUS ALPHA FOR 2P:
C
            ANEG(4) = -A(2,1)*ZSQ(2)
            ANEG(5) = -A(2,2)*ZSQ(2)
            ANEG(6) = -A(2,3)*ZSQ(2)
C
C       PRECOMPUTE EXP(-A*Z**2*DELO**2) WHERE DELO**2 IS
C       DELTA-X**2, DELTA-Y**2, AND DELTA-Z**2
C
            DO 200 IXYZ = 1, MXPTS
               EXPS(IXYZ,1) = GEXP(XDELSQ(IXYZ)*ANEG(1))*VNORM(1)
               EXPS(IXYZ,2) = GEXP(XDELSQ(IXYZ)*ANEG(2))*VNORM(2)
               EXPS(IXYZ,3) = GEXP(XDELSQ(IXYZ)*ANEG(3))*VNORM(3)
               EXPS(IXYZ,4) = GEXP(YDELSQ(IXYZ)*ANEG(1))
               EXPS(IXYZ,5) = GEXP(YDELSQ(IXYZ)*ANEG(2))
               EXPS(IXYZ,6) = GEXP(YDELSQ(IXYZ)*ANEG(3))
               EXPS(IXYZ,7) = GEXP(ZDELSQ(IXYZ)*ANEG(1))
               EXPS(IXYZ,8) = GEXP(ZDELSQ(IXYZ)*ANEG(2))
               EXPS(IXYZ,9) = GEXP(ZDELSQ(IXYZ)*ANEG(3))
  200       CONTINUE
            DO 210 IXYZ = 1, MXPTS
               EXPP(IXYZ,1) = GEXP(XDELSQ(IXYZ)*ANEG(4))
               EXPP(IXYZ,2) = GEXP(XDELSQ(IXYZ)*ANEG(5))
               EXPP(IXYZ,3) = GEXP(XDELSQ(IXYZ)*ANEG(6))
               EXPP(IXYZ,4) = GEXP(YDELSQ(IXYZ)*ANEG(4))
               EXPP(IXYZ,5) = GEXP(YDELSQ(IXYZ)*ANEG(5))
               EXPP(IXYZ,6) = GEXP(YDELSQ(IXYZ)*ANEG(6))
               EXPP(IXYZ,7) = GEXP(ZDELSQ(IXYZ)*ANEG(4))
               EXPP(IXYZ,8) = GEXP(ZDELSQ(IXYZ)*ANEG(5))
               EXPP(IXYZ,9) = GEXP(ZDELSQ(IXYZ)*ANEG(6))
  210       CONTINUE
C
C       LOOP OVER THE "CUBE" Z,Y, AND X
C
            DO 250 IZ = 1, MXPTS
               DO 240 IY = 1, MXPTS
                  DO 230 IG = 1, 3
                     DO 220 IX = 1, MXPTS
C
C       SUM THE ORBITAL CONTRIBUTIONS INTO THE ORBITAL VALUE ARRAY
C
C       FIRST THE 3 GAUSSIANS FOR THE 2S:
C
                        DENSIT(IX,IY,IZ) = DENSIT(IX,IY,IZ)+EXPS(IY,IG+3
     *                     )*EXPS(IZ,IG+6)*EXPS(IX,IG)
C
C       NEXT, THE 3 FOR THE 2P:
C
                        DENSIT(IX,IY,IZ) = DENSIT(IX,IY,IZ)+(CNSTY(IY,IG
     *                     )+CNSTZ(IZ,IG)+CNSTX(IX,IG))*EXPP(IX,IG)*EXPP
     *                     (IY,IG+3)*EXPP(IZ,IG+6)
  220                CONTINUE
  230             CONTINUE
  240          CONTINUE
  250       CONTINUE
            M = M+4
         ELSEIF (IAT.LE.18) THEN
C
C       CALC ZETA(3S) SQUARED AND SQRT FOR CONSTANTS.
C
            ZN = Z1(I)
            ZSQRT = SQRT(ZN)
            ZSQ(1) = ZN*ZN
C
C       ZN * SQRT(ZN) * (2.0E+0/PI)**0.750E+0
C
C       VNORM(1-3) THE THE THREE GAUSSION "CONSTANTS" FOR THE 3S
C       ORBITAL, EVERYTHING THAT IS INDEPENDANT OF R IS IN THERE.
C
            RNORM = ZN*ZSQRT*0.712705470E+0
            VNORM(1) = RNORM*V(M,MO)*CNSTNS(1,3)
            VNORM(2) = RNORM*V(M,MO)*CNSTNS(2,3)
            VNORM(3) = RNORM*V(M,MO)*CNSTNS(3,3)
C
C       CALC ZETA(2P) SQUARED AND SQRT FOR CONSTANTS.
C
            ZN = Z2(I)
            ZSQRT = SQRT(ZN)
            ZSQ(2) = ZN*ZN
C
C       ZN*ZN * SQRT(ZN) * ((128.0E+0/PI**3)**0.250E+0)
C
C       VNORM(1-3) CORRESPONDS TO THE 3 GAUSSIANS FOR THE 3S ORBITAL
C       VNORM(4-6) "                            " FOR THE 3PX ""
C       VNORM(7-9) "" FOR THE 3PY ""
C       VNORM(10-12) "" FOR THE 3PZ ""
C       YOU NEED DIFFERENT "CONSTANTS" FOR EACH DUE TO THE DIFFERENT
C       VALUE OF THE WAVEFUNCTION FOR EACH ORBITAL. AT THE END WE WILL
C       HAVE A 3 SETS OF NUMBERS ALL MULTIPLIED BY THE EXPONENTIAL,
C       WE CAN ADD THEM FIRST, THEN CALC. AND MULT. BY THE EXP.
C
            RNORM = ZSQ(2)*ZSQRT*1.425410940E+0
            AOP1 = RNORM*CNSTNP(1,3)
            AOP2 = RNORM*CNSTNP(2,3)
            AOP3 = RNORM*CNSTNP(3,3)
            VNORM(4) = AOP1*V(M+1,MO)
            VNORM(5) = AOP2*V(M+1,MO)
            VNORM(6) = AOP3*V(M+1,MO)
            VNORM(7) = AOP1*V(M+2,MO)
            VNORM(8) = AOP2*V(M+2,MO)
            VNORM(9) = AOP3*V(M+2,MO)
            VNORM(10) = AOP1*V(M+3,MO)
            VNORM(11) = AOP2*V(M+3,MO)
            VNORM(12) = AOP3*V(M+3,MO)
            DO 260 IXYZ = 1, MXPTS
               CNSTX(IXYZ,1) = VNORM(4)*XDEL(IXYZ)
               CNSTX(IXYZ,2) = VNORM(5)*XDEL(IXYZ)
               CNSTX(IXYZ,3) = VNORM(6)*XDEL(IXYZ)
               CNSTY(IXYZ,1) = VNORM(7)*YDEL(IXYZ)
               CNSTY(IXYZ,2) = VNORM(8)*YDEL(IXYZ)
               CNSTY(IXYZ,3) = VNORM(9)*YDEL(IXYZ)
               CNSTZ(IXYZ,1) = VNORM(10)*ZDEL(IXYZ)
               CNSTZ(IXYZ,2) = VNORM(11)*ZDEL(IXYZ)
               CNSTZ(IXYZ,3) = VNORM(12)*ZDEL(IXYZ)
  260       CONTINUE
C
C       EVALUATE 3S AND 3P
C
C       MINUS ALPHA FOR 3S:
C
            ANEG(1) = -A(3,1)*ZSQ(1)
            ANEG(2) = -A(3,2)*ZSQ(1)
            ANEG(3) = -A(3,3)*ZSQ(1)
C
C       MINUS ALPHA FOR 3P:
C
            ANEG(4) = -A(3,1)*ZSQ(2)
            ANEG(5) = -A(3,2)*ZSQ(2)
            ANEG(6) = -A(3,3)*ZSQ(2)
C
C       PRECOMPUTE EXP(-A*Z**2*DELO**2) WHERE DELO**2 IS
C       DELTA-X**2, DELTA-Y**2, AND DELTA-Z**2
C
            DO 270 IXYZ = 1, MXPTS
               EXPS(IXYZ,1) = GEXP(XDELSQ(IXYZ)*ANEG(1))*VNORM(1)
               EXPS(IXYZ,2) = GEXP(XDELSQ(IXYZ)*ANEG(2))*VNORM(2)
               EXPS(IXYZ,3) = GEXP(XDELSQ(IXYZ)*ANEG(3))*VNORM(3)
               EXPS(IXYZ,4) = GEXP(YDELSQ(IXYZ)*ANEG(1))
               EXPS(IXYZ,5) = GEXP(YDELSQ(IXYZ)*ANEG(2))
               EXPS(IXYZ,6) = GEXP(YDELSQ(IXYZ)*ANEG(3))
               EXPS(IXYZ,7) = GEXP(ZDELSQ(IXYZ)*ANEG(1))
               EXPS(IXYZ,8) = GEXP(ZDELSQ(IXYZ)*ANEG(2))
               EXPS(IXYZ,9) = GEXP(ZDELSQ(IXYZ)*ANEG(3))
  270       CONTINUE
            DO 280 IXYZ = 1, MXPTS
               EXPP(IXYZ,1) = GEXP(XDELSQ(IXYZ)*ANEG(4))
               EXPP(IXYZ,2) = GEXP(XDELSQ(IXYZ)*ANEG(5))
               EXPP(IXYZ,3) = GEXP(XDELSQ(IXYZ)*ANEG(6))
               EXPP(IXYZ,4) = GEXP(YDELSQ(IXYZ)*ANEG(4))
               EXPP(IXYZ,5) = GEXP(YDELSQ(IXYZ)*ANEG(5))
               EXPP(IXYZ,6) = GEXP(YDELSQ(IXYZ)*ANEG(6))
               EXPP(IXYZ,7) = GEXP(ZDELSQ(IXYZ)*ANEG(4))
               EXPP(IXYZ,8) = GEXP(ZDELSQ(IXYZ)*ANEG(5))
               EXPP(IXYZ,9) = GEXP(ZDELSQ(IXYZ)*ANEG(6))
  280       CONTINUE
C
C       LOOP OVER THE "CUBE" Z,Y, AND X
C
            DO 320 IZ = 1, MXPTS
               DO 310 IY = 1, MXPTS
                  DO 300 IG = 1, 3
                     DO 290 IX = 1, MXPTS
C
C       SUM THE ORBITAL CONTRIBUTIONS INTO THE ORBITAL VALUE ARRAY
C
C       FIRST THE 3 GAUSSIANS FOR THE 3S:
C
                        DENSIT(IX,IY,IZ) = DENSIT(IX,IY,IZ)+EXPS(IY,IG+3
     *                     )*EXPS(IZ,IG+6)*EXPS(IX,IG)
C
C       NEXT, THE 3 FOR THE 3P:
C
                        DENSIT(IX,IY,IZ) = DENSIT(IX,IY,IZ)+(CNSTY(IY,IG
     *                     )+CNSTZ(IZ,IG)+CNSTX(IX,IG))*EXPP(IX,IG)*EXPP
     *                     (IY,IG+3)*EXPP(IZ,IG+6)
  290                CONTINUE
  300             CONTINUE
  310          CONTINUE
  320       CONTINUE
            M = M+4
         ENDIF
  330 CONTINUE
      RETURN
      END
      function GEXP(X)
      IF(X.GE.-19.0) THEN
         GEXP = EXP(X)
      ELSE
         GEXP = 0.0E+0
      ENDIF
      RETURN
      END
