|
src
|
|
chk2psi.f,
chk2psi92.f,
chk2psi92.f.txt,
ctplot.f,
gksplot.f,
hpplot.f,
makefile,
preplot.f,
psi1.f,
psi2.f,
psicon.f,
psplot.f
|
|
|
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
|
