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 MO321G (calc,natoms,norbs,MAXATOMS,MAXPTS,mo,ian, & eigvec,xyz,x,y,z,density,numgpx,numgpy,numgpz) C IMPLICIT REAL (A-H,O-Z) C PARAMETER (MXPTS=100) PARAMETER (MAXATM=200) PARAMETER (MAXORB=400) PARAMETER (KD=3,LD=2,MD=1,K3=KD*3,K2=KD*2) 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 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 DIMENSION XDEL(MXPTS),YDEL(MXPTS),ZDEL(MXPTS),AO1S(KD) DIMENSION XDELSQ(MXPTS),YDELSQ(MXPTS),ZDELSQ(MXPTS) c DIMENSION X(MXPTS),Y(MXPTS),Z(MXPTS),XYZ(3,50),ANEG(50) DIMENSION ANEG(maxatm) 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 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/ 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 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 put the eigvec array into the 2D evectors form used bu psi88 c the original v array from psi88 is now the evectors array 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 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, natoms 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 cC READ IN EIGENVECTORS cC IT IS ASSUMED THAT THE EIGENVECTORS HAVE BEEN NORMALIZED TO 1 cC ELECTRON WITH THE OVERLAP MATRIX INCLUDED. cC c IF (IONEMO.NE.0) THEN c READ (IRD,20) (Evectors(I,1),I=1,N) c ELSE c READ (IRD,20,END=30) ((Evectors(I,J),I=1,N),J=1,N) c 20 FORMAT (8F10.6) c ENDIF c 30 CONTINUE c DO 40 J = 1, NAT 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 OF THE I'TH ATOM C C FOR THE PRESENT, ONLY PLOTTING THE MO SPECIFIED BY MONE C c MO = MONE M = 1 DO 840 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 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 * *Evectors(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*Evectors(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, numgpx 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) 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 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*Evectors(M,MO) APNEG = -APLUS(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 2 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 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 * *Evectors(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*Evectors(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*Evectors(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)*Evectors(M+2,MO) 270 CONTINUE C C OUTER P * COEFFICIENT (PX) C AO2PX(LM) = AO2P(LM)*Evectors(M+6,MO) C C INNER P * COEFFICIENT (PY) C DO 280 IG = 1, LD AO2PY(IG) = AO2P(IG)*Evectors(M+3,MO) 280 CONTINUE C C OUTER P * COEFFICIENT (PY) C AO2PY(LM) = AO2P(LM)*Evectors(M+7,MO) C C INNER P * COEFFICIENT (PZ) C DO 290 IG = 1, LD AO2PZ(IG) = AO2P(IG)*Evectors(M+4,MO) 290 CONTINUE C C OUTER P * COEFFICIENT (PZ) C AO2PZ(LM) = AO2P(LM)*Evectors(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, numgpx C C X - DEPENDANT PART ( 2PX ) C CNS2PX(IXYZ,IG) = AO2PX(IG)*XDEL(IXYZ) 300 CONTINUE DO 302 IXYZ = 1, numgpy C C Y - DEPENDANT PART ( 2PY ) C CNS2PY(IXYZ,IG) = AO2PY(IG)*YDEL(IXYZ) 302 CONTINUE DO 304 IXYZ = 1, numgpz 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) 304 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, numgpx E1SX(IXYZ,IG) = GEXP(XDELSQ(IXYZ)*ANEG(IG))*AO1S(IG) 340 CONTINUE DO 342 IXYZ = 1, numgpy E1SY(IXYZ,IG) = GEXP(YDELSQ(IXYZ)*ANEG(IG)) 342 CONTINUE DO 344 IXYZ = 1, numgpz E1SZ(IXYZ,IG) = GEXP(ZDELSQ(IXYZ)*ANEG(IG)) 344 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, numgpx C C X^2 PART OF EXPONENTIAL C E2SPX(IXYZ,IG) = GEXP(XDELSQ(IXYZ)*ANEG(KD+IG)) 360 CONTINUE DO 362 IXYZ = 1, numgpy C C Y^2 PART OF EXPONENTIAL C E2SPY(IXYZ,IG) = GEXP(YDELSQ(IXYZ)*ANEG(KD+IG)) 362 CONTINUE DO 364 IXYZ = 1, numgpz C C Z^2 PART OF EXPONENTIAL C E2SPZ(IXYZ,IG) = GEXP(ZDELSQ(IXYZ)*ANEG(KD+IG)) 364 CONTINUE 370 CONTINUE C C LOOP OVER THE "CUBE" Z,Y, AND X C DO 430 IZ = 1, numgpz 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, numgpy 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, 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) 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, 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) 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, numgpx E2SPX(IXYZ,1) = GEXP(XDELSQ(IXYZ)*APNEG) 440 CONTINUE DO 442 IXYZ = 1, numgpy E2SPY(IXYZ,1) = GEXP(YDELSQ(IXYZ)*APNEG) 442 CONTINUE DO 444 IXYZ = 1, numgpz E2SPZ(IXYZ,1) = GEXP(ZDELSQ(IXYZ)*APNEG) 444 CONTINUE AO2S(1) = Evectors(M,MO)*(APLUS(IAT)**0.750E+0) & *0.712705470E+0 AO2P(1) = (APLUS(IAT)**1.250E+0)*1.425410940E+0 AO2PX(1) = AO2P(1)*Evectors(M+1,MO) AO2PY(1) = AO2P(1)*Evectors(M+2,MO) AO2PZ(1) = AO2P(1)*Evectors(M+3,MO) DO 450 IXYZ = 1, numgpx CNS2PX(IXYZ,1) = AO2PX(1)*XDEL(IXYZ) 450 CONTINUE DO 452 IXYZ = 1, numgpy CNS2PY(IXYZ,1) = AO2PY(1)*YDEL(IXYZ) 452 CONTINUE DO 454 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) 454 CONTINUE C C LOOP OVER THE "CUBE" Z,Y, AND X C DO 480 IZ = 1, numgpz DO 470 IY = 1, numgpy DO 460 IX = 1, numgpx 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 * *Evectors(M,MO) C C CORE 2SP SHELL.... 2S: C AO2S(IG) = S(IG2,IAT)*(A(IG2,IAT)**0.750E+0)*0.712705470E * +0*Evectors(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)*Evectors(M+2,MO) C C CORE 2SP SHELL.... 2PY: C AO2PY(IG) = AO2P(IG)*Evectors(M+3,MO) C C CORE 2SP SHELL.... 2PZ: C AO2PZ(IG) = AO2P(IG)*Evectors(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*Evectors(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*Evectors(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)*Evectors(M+6,MO) 520 CONTINUE C C OUTER P * COEFFICIENT (PX) C AO3PX(LM) = AO3P(LM)*Evectors(M+10,MO) C C INNER P * COEFFICIENT (PY) C DO 530 IG = 1, LD AO3PY(IG) = AO3P(IG)*Evectors(M+7,MO) 530 CONTINUE C C OUTER P * COEFFICIENT (PY) C AO3PY(LM) = AO3P(LM)*Evectors(M+11,MO) C C INNER P * COEFFICIENT (PZ) C DO 540 IG = 1, LD AO3PZ(IG) = AO3P(IG)*Evectors(M+8,MO) 540 CONTINUE C C OUTER P * COEFFICIENT (PZ) C AO3PZ(LM) = AO3P(LM)*Evectors(M+12,MO) C C CORE 2P X,Y,AND Z DEPENDANT PARTS C DO 560 IG = 1, KD DO 550 IXYZ = 1, numgpx C C CORE 2PX - X DEPENDANT PART C CNS2PX(IXYZ,IG) = AO2PX(IG)*XDEL(IXYZ) 550 CONTINUE DO 552 IXYZ = 1, numgpy C C CORE 2PY - Y DEPENDANT PART C CNS2PY(IXYZ,IG) = AO2PY(IG)*YDEL(IXYZ) 552 CONTINUE 4 DO 554 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) 554 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, numgpx CNS3PX(IXYZ,IG) = AO3PX(IG)*XDEL(IXYZ) 570 CONTINUE DO 572 IXYZ = 1, numgpy CNS3PY(IXYZ,IG) = AO3PY(IG)*YDEL(IXYZ) 572 CONTINUE DO 574 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) 574 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, numgpx C C E^(-ALPHA(1S)*X^2) C E1SX(IXYZ,IG) = GEXP(XDELSQ(IXYZ)*ANEG(IG))*AO1S(IG) 620 CONTINUE DO 622 IXYZ = 1, numgpy C C E^(-ALPHA(1S)*Y^2) C E1SY(IXYZ,IG) = GEXP(YDELSQ(IXYZ)*ANEG(IG)) 622 CONTINUE DO 624 IXYZ = 1, numgpz C C E^(-ALPHA(1S)*Z^2) C E1SZ(IXYZ,IG) = GEXP(ZDELSQ(IXYZ)*ANEG(IG)) 624 CONTINUE 630 CONTINUE DO 650 IG = 1, KD DO 640 IXYZ = 1, numgpx C C E^(-ALPHA(2SP)*X^2) C E2SPX(IXYZ,IG) = GEXP(XDELSQ(IXYZ)*ANEG(KD+IG)) 640 CONTINUE DO 642 IXYZ = 1, numgpy C C E^(-ALPHA(2SP)*Y^2) C E2SPY(IXYZ,IG) = GEXP(YDELSQ(IXYZ)*ANEG(KD+IG)) 642 CONTINUE DO 644 IXYZ = 1, numgpz C C E^(-ALPHA(2SP)*Z^2) C E2SPZ(IXYZ,IG) = GEXP(ZDELSQ(IXYZ)*ANEG(KD+IG)) 644 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, numgpx C C E^(-ALPHA(3SP)*X^2) C E3SPX(IXYZ,IG) = GEXP(XDELSQ(IXYZ)*ANEG(K2+IG)) 660 CONTINUE DO 662 IXYZ = 1, numgpy C C E^(-ALPHA(3SP)*Y^2) C E3SPY(IXYZ,IG) = GEXP(YDELSQ(IXYZ)*ANEG(K2+IG)) 662 CONTINUE DO 664 IXYZ = 1, numgpz C C E^(-ALPHA(3SP)*Z^2) C E3SPZ(IXYZ,IG) = GEXP(ZDELSQ(IXYZ)*ANEG(K2+IG)) 664 CONTINUE 670 CONTINUE C C LOOP OVER THE "CUBE" Z,Y, AND X C DO 730 IZ = 1, numgpz 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, numgpy 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, 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) 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, 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) 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, numgpx EXPDX(IXYZ) = GEXP(XDELSQ(IXYZ)*ASNEG) 740 CONTINUE DO 742 IXYZ = 1, numgpy EXPDY(IXYZ) = GEXP(YDELSQ(IXYZ)*ASNEG) 742 CONTINUE DO 744 IXYZ = 1, numgpz EXPDZ(IXYZ) = GEXP(ZDELSQ(IXYZ)*ASNEG) 744 CONTINUE TMP = (ASTAR(IAT)**1.750E+0)*1.645922780E+0 AODXX = Evectors(M,MO)*TMP AODYY = Evectors(M+1,MO)*TMP AODZZ = Evectors(M+2,MO)*TMP AODXY = Evectors(M+3,MO)*TMP AODXZ = Evectors(M+4,MO)*TMP AODYZ = Evectors(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, numgpx CNSXX(IXYZ) = AODXX*XDELSQ(IXYZ) 750 CONTINUE DO 752 IXYZ = 1, numgpy CNSYY(IXYZ) = AODYY*YDELSQ(IXYZ) CNSXY(IXYZ) = AODXY*YDEL(IXYZ) CNSYZ(IXYZ) = AODYZ*YDEL(IXYZ) 752 CONTINUE DO 754 IXYZ = 1, numgpz CNSZZ(IXYZ) = AODZZ*ZDELSQ(IXYZ) CNSXZ(IXYZ) = AODXZ*ZDEL(IXYZ) 754 CONTINUE C C LOOP OVER THE "CUBE" Z,Y, AND X C DO 780 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 770 IY = 1, numgpy DO 760 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) 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, numgpx E3SPX(IXYZ,1) = GEXP(XDELSQ(IXYZ)*APNEG) 790 CONTINUE DO 792 IXYZ = 1, numgpy E3SPY(IXYZ,1) = GEXP(YDELSQ(IXYZ)*APNEG) 792 CONTINUE DO 794 IXYZ = 1, numgpz E3SPZ(IXYZ,1) = GEXP(ZDELSQ(IXYZ)*APNEG) 794 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(M,MO)*(APLUS(IAT)**0.750E+0)*0.712705470E+0 AO3P(1) = (APLUS(IAT)**1.250E+0)*1.425410940E+0 AO3PX(1) = AO3P(1)*Evectors(M+1,MO) AO3PY(1) = AO3P(1)*Evectors(M+2,MO) AO3PZ(1) = AO3P(1)*Evectors(M+3,MO) DO 800 IXYZ = 1, numgpx CNS3PX(IXYZ,1) = AO3PX(1)*XDEL(IXYZ) 800 CONTINUE DO 802 IXYZ = 1, numgpy CNS3PY(IXYZ,1) = AO3PY(1)*YDEL(IXYZ) 802 CONTINUE DO 804 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) 804 CONTINUE C C LOOP OVER THE "CUBE" Z,Y, AND X C DO 830 IZ = 1, numgpz DO 820 IY = 1, numgpy DO 810 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) 810 CONTINUE 820 CONTINUE 830 CONTINUE M = M+4 ENDIF ENDIF 840 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