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