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 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 subroutine is adapted from the psi1.f program (referred to as psi88) c that is available via anonymous ftp. c generally all modifications to this subroutine for the display_gauss c module are in lowercase text SUBROUTINE MO631G (calc,natoms,norbs,MAXATOMS,MAXPTS,mo,ian, & eigvec,xyz,x,y,z,density,numgpx,numgpy,numgpz) IMPLICIT REAL (A-H,O-Z) C C THIS DETERMINES THE DENSITY OF PLANES COMPUTED C PARAMETER (MXPTS=100) PARAMETER (MAXATM=200) PARAMETER (MAXORB=400) 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 c COMMON /SPLIT/ XMI,YMI,ZMI,XINC,YINC,ZINC,MONE,MOLAST,NAT, c * IONEMO c COMMON /DENS/ DENSIT(MXPTS,MXPTS,MXPTS),V(MAXORB,MAXORB),C(MAXATM, c * 3),OCMO(MAXORB),IAN(MAXATM) c COMMON /IO/ IRD,ILST c the 2d array c in this subroutine corresponds to xyz in the calling program c the 2d array xyz that was originally in this subroutine is the c transposed version of the c matrix, since going from the C (language) c to fortran will also cause the 2D matrices to be transposed, we will c pass the main xyz[maxatoms][3] to xyz(3,maxatoms) and omit the c matrix c to xyz conversion real*8 evectors(maxorb,maxorb),xyz(3,MAXATM),eigvec(maxorb*maxorb) real*8 x(MAXATM),y(MAXATM),z(MAXATM) real*8 densit(MXPTS,MXPTS,MXPTS),density(MXPTS,MXPTS,MXPTS) integer ian(MAXATM) data ilst /6/ DIMENSION XDEL(MXPTS),YDEL(MXPTS),ZDEL(MXPTS) DIMENSION XDELSQ(MXPTS),YDELSQ(MXPTS),ZDELSQ(MXPTS) c DIMENSION XYZ(3,50) c DIMENSION X(MXPTS),Y(MXPTS),Z(MXPTS) DIMENSION CNSXX(MXPTS),CNSYY(MXPTS),CNSZZ(MXPTS),AO1S(KD) DIMENSION ANEG(maxatm),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 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 put the eigvec array into the 2D evectors form used bu psi88 icnt = 1 do i=1,norbs do j=1,norbs evectors(j,i) = eigvec(icnt) icnt = icnt + 1 enddo enddo 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, natoms 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 cC cC READ IN EIGENVECTORS cC c IF (IONEMO.NE.0) THEN c READ (IRD,20) (V(I,1),I=1,NBASIS) c ELSE c READ (IRD,20,END=30) ((V(I,J),I=1,NBASIS),J=1,NBASIS) c 20 FORMAT (8F10.6) c ENDIF c 30 CONTINUE c DO 40 J = 1, natoms c XYZ(1,J) = C(J,1) c XYZ(2,J) = C(J,2) c XYZ(3,J) = C(J,3) c 40 CONTINUE cC cC COMPUTE X,Y, AND Z VALUES FOR THE GRID (CUBE) POINTS. cC c X(1) = XMI c Y(1) = YMI c Z(1) = ZMI c DO 50 I = 2, MXPTS c X(I) = XINC+X(I-1) c Y(I) = YINC+Y(I-1) c Z(I) = ZINC+Z(I-1) c 50 CONTINUE C C ZERO THE ORBITAL VALUE ARRAY C DO 80 IZ = 1, numgpz DO 70 IY = 1, numgpy DO 60 IX = 1, numgpx 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 c MO = MONE KNTBAS = 1 DO 940 I = 1, natoms 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, numgpx XDEL(IXYZ) = X(IXYZ)-XYZ(1,I) XDELSQ(IXYZ) = XDEL(IXYZ)*XDEL(IXYZ) 90 CONTINUE DO 100 IXYZ = 1, numgpy YDEL(IXYZ) = Y(IXYZ)-XYZ(2,I) YDELSQ(IXYZ) = YDEL(IXYZ)*YDEL(IXYZ) 100 CONTINUE DO 110 IXYZ = 1, numgpz 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 * *evectors(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*evectors(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 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 DO 150 IG = 1, LM DO 140 IXYZ = 1, numgpx E1SX(IXYZ,IG) = GEXP(XDELSQ(IXYZ)*ANEG(IG))*AO1S(IG) 140 CONTINUE DO 142 IXYZ = 1, numgpy E1SY(IXYZ,IG) = GEXP(YDELSQ(IXYZ)*ANEG(IG)) 142 CONTINUE DO 144 IXYZ = 1, numgpz E1SZ(IXYZ,IG) = GEXP(ZDELSQ(IXYZ)*ANEG(IG)) 144 CONTINUE 150 CONTINUE DO 190 IZ = 1, numgpz DO 180 IY = 1, numgpy DO 170 IG = 1, LM DO 160 IX = 1, numgpx 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* * Evectors(KNTBAS,MO) APNEG = -APL(IAT) DO 200 IXYZ = 1, numgpx E1SX(IXYZ,1) = GEXP(XDELSQ(IXYZ)*APNEG)*AO1S(1) 200 CONTINUE DO 202 IXYZ = 1, numgpy E1SY(IXYZ,1) = GEXP(YDELSQ(IXYZ)*APNEG) 202 CONTINUE DO 204 IXYZ = 1, numgpz E1SZ(IXYZ,1) = GEXP(ZDELSQ(IXYZ)*APNEG) 204 CONTINUE C C LOOP OVER THE "CUBE" Z,Y, AND X C DO 230 IZ = 1, numgpz DO 220 IY = 1, numgpy DO 210 IX = 1, numgpx 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, numgpx E2SPX(IXYZ,1) = GEXP(XDELSQ(IXYZ)*ASNEG) 240 CONTINUE DO 242 IXYZ = 1, numgpy E2SPY(IXYZ,1) = GEXP(YDELSQ(IXYZ)*ASNEG) 242 CONTINUE DO 244 IXYZ = 1, numgpz E2SPZ(IXYZ,1) = GEXP(ZDELSQ(IXYZ)*ASNEG) 244 CONTINUE AO2P(1) = (ASTAR(IAT)**1.250E+0)*1.425410940E+0 AO2PX(1) = AO2P(1)*Evectors(KNTBAS,MO) AO2PY(1) = AO2P(1)*Evectors(KNTBAS+1,MO) AO2PZ(1) = AO2P(1)*Evectors(KNTBAS+2,MO) DO 250 IXYZ = 1, numgpx CNS2PX(IXYZ,1) = AO2PX(1)*XDEL(IXYZ) 250 CONTINUE DO 252 IXYZ = 1, numgpy CNS2PY(IXYZ,1) = AO2PY(1)*YDEL(IXYZ) 252 CONTINUE DO 254 IXYZ = 1, numgpz CNS2PZ(IXYZ,1) = AO2PZ(1)*ZDEL(IXYZ) 254 CONTINUE C C LOOP OVER THE "CUBE" Z,Y, AND X C DO 280 IZ = 1, numgpz DO 270 IY = 1, numgpy DO 260 IX = 1, numgpx 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 * *Evectors(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*Evectors(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*Evectors(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)*Evectors(KNTBAS+2,MO) 320 CONTINUE C C OUTER P * COEFFICIENT (PX) C AO2PX(LM) = AO2P(LM)*Evectors(KNTBAS+6,MO) C C INNER P * COEFFICIENT (PY) C DO 330 IG = 1, LD AO2PY(IG) = AO2P(IG)*Evectors(KNTBAS+3,MO) 330 CONTINUE C C OUTER P * COEFFICIENT (PY) C AO2PY(LM) = AO2P(LM)*Evectors(KNTBAS+7,MO) C C INNER P * COEFFICIENT (PZ) C DO 340 IG = 1, LD AO2PZ(IG) = AO2P(IG)*Evectors(KNTBAS+4,MO) 340 CONTINUE C C OUTER P * COEFFICIENT (PZ) C AO2PZ(LM) = AO2P(LM)*Evectors(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, numgpx C C X - DEPENDANT PART ( 2PX ) C CNS2PX(IXYZ,IG) = AO2PX(IG)*XDEL(IXYZ) 350 CONTINUE DO 352 IXYZ = 1, numgpy C C Y - DEPENDANT PART ( 2PY ) C CNS2PY(IXYZ,IG) = AO2PY(IG)*YDEL(IXYZ) 352 CONTINUE DO 354 IXYZ = 1, numgpz 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) 354 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, numgpx E1SX(IXYZ,IG) = GEXP(XDELSQ(IXYZ)*ANEG(IG))*AO1S(IG) 390 CONTINUE DO 392 IXYZ = 1, numgpy E1SY(IXYZ,IG) = GEXP(YDELSQ(IXYZ)*ANEG(IG)) 392 CONTINUE DO 394 IXYZ = 1, numgpz E1SZ(IXYZ,IG) = GEXP(ZDELSQ(IXYZ)*ANEG(IG)) 394 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, numgpx C C X^2 PART OF EXPONENTIAL C E2SPX(IXYZ,IG) = GEXP(XDELSQ(IXYZ)*ANEG(KD+IG)) 410 CONTINUE DO 412 IXYZ = 1, numgpy C C Y^2 PART OF EXPONENTIAL C E2SPY(IXYZ,IG) = GEXP(YDELSQ(IXYZ)*ANEG(KD+IG)) 412 CONTINUE DO 414 IXYZ = 1, numgpz C C Z^2 PART OF EXPONENTIAL C E2SPZ(IXYZ,IG) = GEXP(ZDELSQ(IXYZ)*ANEG(KD+IG)) 414 CONTINUE 420 CONTINUE C C LOOP OVER THE "CUBE" Z,Y, AND X C DO 480 IZ = 1, numgpz 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, numgpy DO 440 IG = 1, KD DO 430 IX = 1, numgpx 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, numgpx 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, numgpx E2SPX(IXYZ,1) = GEXP(XDELSQ(IXYZ)*APNEG) 490 CONTINUE DO 492 IXYZ = 1, numgpy E2SPY(IXYZ,1) = GEXP(YDELSQ(IXYZ)*APNEG) 492 CONTINUE DO 494 IXYZ = 1, numgpz E2SPZ(IXYZ,1) = GEXP(ZDELSQ(IXYZ)*APNEG) 494 CONTINUE AO2S(1) = Evectors(KNTBAS,MO)*(APL(IAT)**0.750E+0)* * 0.712705470E+0 AO2P(1) = (APL(IAT)**1.250E+0)*1.425410940E+0 AO2PX(1) = AO2P(1)*Evectors(KNTBAS+1,MO) AO2PY(1) = AO2P(1)*Evectors(KNTBAS+2,MO) AO2PZ(1) = AO2P(1)*Evectors(KNTBAS+3,MO) DO 500 IXYZ = 1, numgpx CNS2PX(IXYZ,1) = AO2PX(1)*XDEL(IXYZ) 500 CONTINUE DO 502 IXYZ = 1, numgpy CNS2PY(IXYZ,1) = AO2PY(1)*YDEL(IXYZ) 502 CONTINUE DO 504 IXYZ = 1, numgpz 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) 504 CONTINUE C C LOOP OVER THE "CUBE" Z,Y, AND X C DO 530 IZ = 1, numgpz DO 520 IY = 1, numgpy DO 510 IX = 1, numgpx 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, numgpx EXPDX(IXYZ) = GEXP(XDELSQ(IXYZ)*ASNEG) 540 CONTINUE DO 542 IXYZ = 1, numgpy EXPDY(IXYZ) = GEXP(YDELSQ(IXYZ)*ASNEG) 542 CONTINUE DO 544 IXYZ = 1, numgpz EXPDZ(IXYZ) = GEXP(ZDELSQ(IXYZ)*ASNEG) 544 CONTINUE TMP = (ASTAR(IAT)**1.750E+0)*1.645922780E+0 AODXX = Evectors(KNTBAS,MO)*TMP AODYY = Evectors(KNTBAS+1,MO)*TMP AODZZ = Evectors(KNTBAS+2,MO)*TMP AODXY = Evectors(KNTBAS+3,MO)*TMP AODXZ = Evectors(KNTBAS+4,MO)*TMP AODYZ = Evectors(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, numgpx CNSXX(IXYZ) = AODXX*XDELSQ(IXYZ) 550 CONTINUE DO 552 IXYZ = 1, numgpy CNSYY(IXYZ) = AODYY*YDELSQ(IXYZ) CNSXY(IXYZ) = AODXY*YDEL(IXYZ) CNSYZ(IXYZ) = AODYZ*YDEL(IXYZ) 552 CONTINUE DO 554 IXYZ = 1, numgpz CNSZZ(IXYZ) = AODZZ*ZDELSQ(IXYZ) CNSXZ(IXYZ) = AODXZ*ZDEL(IXYZ) 554 CONTINUE C C LOOP OVER THE "CUBE" Z,Y, AND X C DO 580 IZ = 1, numgpz DO 570 IY = 1, numgpy DO 560 IX = 1, numgpx 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 * *Evectors(KNTBAS,MO) C C CORE 2SP SHELL.... 2S: C AO2S(IG) = S(IG2,IAT)*(A(IG2,IAT)**0.750E+0)*0.712705470E * +0*Evectors(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)*Evectors(KNTBAS+2,MO) C C CORE 2SP SHELL.... 2PY: C AO2PY(IG) = AO2P(IG)*Evectors(KNTBAS+3,MO) C C CORE 2SP SHELL.... 2PZ: C AO2PZ(IG) = AO2P(IG)*Evectors(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*Evectors(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*Evectors(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)*Evectors(KNTBAS+6,MO) 620 CONTINUE C C OUTER P * COEFFICIENT (PX) C AO3PX(LM) = AO3P(LM)*Evectors(KNTBAS+10,MO) C C INNER P * COEFFICIENT (PY) C DO 630 IG = 1, LD AO3PY(IG) = AO3P(IG)*Evectors(KNTBAS+7,MO) 630 CONTINUE C C OUTER P * COEFFICIENT (PY) C AO3PY(LM) = AO3P(LM)*Evectors(KNTBAS+11,MO) C C INNER P * COEFFICIENT (PZ) C DO 640 IG = 1, LD AO3PZ(IG) = AO3P(IG)*Evectors(KNTBAS+8,MO) 640 CONTINUE C C OUTER P * COEFFICIENT (PZ) C AO3PZ(LM) = AO3P(LM)*Evectors(KNTBAS+12,MO) C C CORE 2P X,Y,AND Z DEPENDANT PARTS C DO 660 IG = 1, KD DO 650 IXYZ = 1, numgpx C C CORE 2PX - X DEPENDANT PART C CNS2PX(IXYZ,IG) = AO2PX(IG)*XDEL(IXYZ) 650 CONTINUE DO 652 IXYZ = 1, numgpy C C CORE 2PY - Y DEPENDANT PART C CNS2PY(IXYZ,IG) = AO2PY(IG)*YDEL(IXYZ) 652 CONTINUE DO 654 IXYZ = 1, numgpz 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) 654 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, numgpx CNS3PX(IXYZ,IG) = AO3PX(IG)*XDEL(IXYZ) 670 CONTINUE DO 672 IXYZ = 1, numgpy CNS3PY(IXYZ,IG) = AO3PY(IG)*YDEL(IXYZ) 672 CONTINUE DO 674 IXYZ = 1, numgpz 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) 674 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, numgpx C C E^(-ALPHA(1S)*X^2) C E1SX(IXYZ,IG) = GEXP(XDELSQ(IXYZ)*ANEG(IG))*AO1S(IG) 720 CONTINUE DO 722 IXYZ = 1, numgpy C C E^(-ALPHA(1S)*Y^2) C E1SY(IXYZ,IG) = GEXP(YDELSQ(IXYZ)*ANEG(IG)) 722 CONTINUE DO 724 IXYZ = 1, numgpz C C E^(-ALPHA(1S)*Z^2) C E1SZ(IXYZ,IG) = GEXP(ZDELSQ(IXYZ)*ANEG(IG)) 724 CONTINUE 730 CONTINUE DO 750 IG = 1, KD DO 740 IXYZ = 1, numgpx C C E^(-ALPHA(2SP)*X^2) C E2SPX(IXYZ,IG) = GEXP(XDELSQ(IXYZ)*ANEG(KD+IG)) 740 CONTINUE DO 742 IXYZ = 1, numgpy C C E^(-ALPHA(2SP)*Y^2) C E2SPY(IXYZ,IG) = GEXP(YDELSQ(IXYZ)*ANEG(KD+IG)) 742 CONTINUE DO 744 IXYZ = 1, numgpz C C E^(-ALPHA(2SP)*Z^2) C E2SPZ(IXYZ,IG) = GEXP(ZDELSQ(IXYZ)*ANEG(KD+IG)) 744 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, numgpx C C E^(-ALPHA(3SP)*X^2) C E3SPX(IXYZ,IG) = GEXP(XDELSQ(IXYZ)*ANEG(K2+IG)) 760 CONTINUE DO 762 IXYZ = 1, numgpy C C E^(-ALPHA(3SP)*Y^2) C E3SPY(IXYZ,IG) = GEXP(YDELSQ(IXYZ)*ANEG(K2+IG)) 762 CONTINUE DO 764 IXYZ = 1, numgpz C C E^(-ALPHA(3SP)*Z^2) C E3SPZ(IXYZ,IG) = GEXP(ZDELSQ(IXYZ)*ANEG(K2+IG)) 764 CONTINUE 770 CONTINUE C C C LOOP OVER THE "CUBE" Z,Y, AND X C DO 830 IZ = 1, numgpz 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, numgpy 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, numgpx 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, numgpx 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, numgpx E3SPX(IXYZ,1) = GEXP(XDELSQ(IXYZ)*APNEG) 840 CONTINUE DO 842 IXYZ = 1, numgpy E3SPY(IXYZ,1) = GEXP(YDELSQ(IXYZ)*APNEG) 842 CONTINUE DO 844 IXYZ = 1, numgpz E3SPZ(IXYZ,1) = GEXP(ZDELSQ(IXYZ)*APNEG) 844 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) = Evectors(KNTBAS,MO)*(APL(IAT)**0.750E+0)* * 0.712705470E+0 AO3P(1) = (APL(IAT)**1.250E+0)*1.425410940E+0 AO3PX(1) = AO3P(1)*Evectors(KNTBAS+1,MO) AO3PY(1) = AO3P(1)*Evectors(KNTBAS+2,MO) AO3PZ(1) = AO3P(1)*Evectors(KNTBAS+3,MO) DO 850 IXYZ = 1, numgpx CNS3PX(IXYZ,1) = AO3PX(1)*XDEL(IXYZ) 850 CONTINUE DO 852 IXYZ = 1, numgpy CNS3PY(IXYZ,1) = AO3PY(1)*YDEL(IXYZ) 852 CONTINUE DO 854 IXYZ = 1, numgpz 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) 854 CONTINUE C C LOOP OVER THE "CUBE" Z,Y, AND X C DO 880 IZ = 1, numgpz DO 870 IY = 1, numgpy DO 860 IX = 1, numgpx 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, numgpx EXPDX(IXYZ) = GEXP(XDELSQ(IXYZ)*ASNEG) 890 CONTINUE DO 892 IXYZ = 1, numgpy EXPDY(IXYZ) = GEXP(YDELSQ(IXYZ)*ASNEG) 892 CONTINUE DO 894 IXYZ = 1, numgpz EXPDZ(IXYZ) = GEXP(ZDELSQ(IXYZ)*ASNEG) 894 CONTINUE TMP = (ASTAR(IAT)**1.750E+0)*1.645922780E+0 AODXX = Evectors(KNTBAS,MO)*TMP AODYY = Evectors(KNTBAS+1,MO)*TMP AODZZ = Evectors(KNTBAS+2,MO)*TMP AODXY = Evectors(KNTBAS+3,MO)*TMP AODXZ = Evectors(KNTBAS+4,MO)*TMP AODYZ = Evectors(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, numgpx CNSXX(IXYZ) = AODXX*XDELSQ(IXYZ) 900 CONTINUE DO 902 IXYZ = 1, numgpy CNSYY(IXYZ) = AODYY*YDELSQ(IXYZ) CNSXY(IXYZ) = AODXY*YDEL(IXYZ) CNSYZ(IXYZ) = AODYZ*YDEL(IXYZ) 902 CONTINUE DO 904 IXYZ = 1, numgpz CNSZZ(IXYZ) = AODZZ*ZDELSQ(IXYZ) CNSXZ(IXYZ) = AODXZ*ZDEL(IXYZ) 904 CONTINUE C C LOOP OVER THE "CUBE" Z,Y, AND X C DO 930 IZ = 1, numgpz 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, numgpy DO 910 IX = 1, numgpx 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 c now transpose the densit matrix before passing it back as density do 980 ix=1,numgpx do 982 iy=1,numgpy do 984 iz=1,numgpz density(iz,iy,ix) = densit(ix,iy,iz) 984 continue 982 continue 980 continue RETURN END c ***************************************************************************** c ***************************************************************************** C C function GEXP(X) c real*8 x IF(X.GE.-19.0) THEN GEXP = EXP(X) ELSE GEXP = 0.0E+0 ENDIF RETURN END