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 STOMO (calc,natoms,norbs,MAXATOMS,MAXPTS,mo,ian, & eigvec,xyz,x,y,z,density,numgpx,numgpy,numgpz) IMPLICIT REAL (A-H,O-Z) PARAMETER (MXPTS=100) PARAMETER (MAXATM=200) PARAMETER (MAXORB=400) 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 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 COMMON /BASIS/ CALC CHARACTER CALC*20 DIMENSION XDEL(MXPTS),YDEL(MXPTS),ZDEL(MXPTS),EXP2SP(MXPTS,9) DIMENSION XDELSQ(MXPTS),YDELSQ(MXPTS),ZDELSQ(MXPTS) c DIMENSION XYZ(3,50),ZSQ(3),EXP1S(MXPTS,9) DIMENSION ZSQ(3),EXP1S(MXPTS,9) c DIMENSION X(MXPTS),Y(MXPTS),Z(MXPTS),CNSTX(MXPTS,6),CNSTY(MXPTS,6) DIMENSION 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) 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/ 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 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 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, natoms 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 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) (V(I,1),I=1,N) c ELSE c READ (IRD,20,END=30) ((V(I,J),I=1,N),J=1,N) c 20 FORMAT (8F10.6) c ENDIF c 30 CONTINUE c MO = MONE c WRITE (ILST,*) ' BASIS SET IS ',CALC,' NUMBER OF AOS IS ',N cC Write(ILST,*)' plotting MO number ',MO 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 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 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, 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 M = 1 DO 310 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 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*Evectors(M,MO)*CNSTNS(1,1) VNORM(2) = RNORM*Evectors(M,MO)*CNSTNS(2,1) VNORM(3) = RNORM*Evectors(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, numgpx 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) 120 CONTINUE DO 122 IXYZ = 1, numgpy EXP1S(IXYZ,4) = GEXP(YDELSQ(IXYZ)*ANEG(1)) EXP1S(IXYZ,5) = GEXP(YDELSQ(IXYZ)*ANEG(2)) EXP1S(IXYZ,6) = GEXP(YDELSQ(IXYZ)*ANEG(3)) 122 CONTINUE DO 124 IXYZ = 1, numgpz EXP1S(IXYZ,7) = GEXP(ZDELSQ(IXYZ)*ANEG(1)) EXP1S(IXYZ,8) = GEXP(ZDELSQ(IXYZ)*ANEG(2)) EXP1S(IXYZ,9) = GEXP(ZDELSQ(IXYZ)*ANEG(3)) 124 CONTINUE DO 160 IZ = 1, numgpz DO 150 IY = 1, numgpy DO 140 IG = 1, 3 DO 130 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)+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*Evectors(M,MO)*CNSTNS(1,1) VNORM(2) = RNORM*Evectors(M,MO)*CNSTNS(2,1) VNORM(3) = RNORM*Evectors(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)*Evectors(M+1,MO) VNORM(5) = RNORM*CNSTNS(2,2)*Evectors(M+1,MO) VNORM(6) = RNORM*CNSTNS(3,2)*Evectors(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)*Evectors(M+3,MO) VNORM(11) = VNORM(8)*Evectors(M+3,MO) VNORM(12) = VNORM(9)*Evectors(M+3,MO) VNORM(13) = VNORM(7)*Evectors(M+4,MO) VNORM(14) = VNORM(8)*Evectors(M+4,MO) VNORM(15) = VNORM(9)*Evectors(M+4,MO) VNORM(7) = VNORM(7)*Evectors(M+2,MO) VNORM(8) = VNORM(8)*Evectors(M+2,MO) VNORM(9) = VNORM(9)*Evectors(M+2,MO) DO 170 IXYZ = 1, numgpx CNSTX(IXYZ,1) = VNORM(7)*XDEL(IXYZ) CNSTX(IXYZ,2) = VNORM(8)*XDEL(IXYZ) CNSTX(IXYZ,3) = VNORM(9)*XDEL(IXYZ) 170 CONTINUE DO 172 IXYZ = 1, numgpy CNSTY(IXYZ,1) = VNORM(10)*YDEL(IXYZ) CNSTY(IXYZ,2) = VNORM(11)*YDEL(IXYZ) CNSTY(IXYZ,3) = VNORM(12)*YDEL(IXYZ) 172 CONTINUE DO 174 IXYZ = 1, numgpz 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) 174 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, numgpx 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) 180 CONTINUE DO 182 IXYZ = 1, numgpy EXP1S(IXYZ,4) = GEXP(YDELSQ(IXYZ)*ANEG(1)) EXP1S(IXYZ,5) = GEXP(YDELSQ(IXYZ)*ANEG(2)) EXP1S(IXYZ,6) = GEXP(YDELSQ(IXYZ)*ANEG(3)) 182 CONTINUE DO 184 IXYZ = 1, numgpz EXP1S(IXYZ,7) = GEXP(ZDELSQ(IXYZ)*ANEG(1)) EXP1S(IXYZ,8) = GEXP(ZDELSQ(IXYZ)*ANEG(2)) EXP1S(IXYZ,9) = GEXP(ZDELSQ(IXYZ)*ANEG(3)) 184 CONTINUE DO 190 IXYZ = 1, numgpx EXP2SP(IXYZ,1) = GEXP(XDELSQ(IXYZ)*ANEG(4)) EXP2SP(IXYZ,2) = GEXP(XDELSQ(IXYZ)*ANEG(5)) EXP2SP(IXYZ,3) = GEXP(XDELSQ(IXYZ)*ANEG(6)) 190 CONTINUE DO 192 IXYZ = 1, numgpy EXP2SP(IXYZ,4) = GEXP(YDELSQ(IXYZ)*ANEG(4)) EXP2SP(IXYZ,5) = GEXP(YDELSQ(IXYZ)*ANEG(5)) EXP2SP(IXYZ,6) = GEXP(YDELSQ(IXYZ)*ANEG(6)) 192 CONTINUE DO 194 IXYZ = 1, numgpz EXP2SP(IXYZ,7) = GEXP(ZDELSQ(IXYZ)*ANEG(4)) EXP2SP(IXYZ,8) = GEXP(ZDELSQ(IXYZ)*ANEG(5)) EXP2SP(IXYZ,9) = GEXP(ZDELSQ(IXYZ)*ANEG(6)) 194 CONTINUE C C LOOP OVER THE "CUBE" Z,Y, AND X C DO 230 IZ = 1, numgpz DO 220 IY = 1, numgpy DO 210 IG = 1, 3 DO 200 IX = 1, numgpx 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*Evectors(M,MO)*CNSTNS(1,1) VNORM(2) = RNORM*Evectors(M,MO)*CNSTNS(2,1) VNORM(3) = RNORM*Evectors(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)*Evectors(M+1,MO) VNORM(5) = RNORM*CNSTNS(2,2)*Evectors(M+1,MO) VNORM(6) = RNORM*CNSTNS(3,2)*Evectors(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)*Evectors(M+3,MO) VNORM(11) = VNORM(8)*Evectors(M+3,MO) VNORM(12) = VNORM(9)*Evectors(M+3,MO) VNORM(13) = VNORM(7)*Evectors(M+4,MO) VNORM(14) = VNORM(8)*Evectors(M+4,MO) VNORM(15) = VNORM(9)*Evectors(M+4,MO) VNORM(7) = VNORM(7)*Evectors(M+2,MO) VNORM(8) = VNORM(8)*Evectors(M+2,MO) VNORM(9) = VNORM(9)*Evectors(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)*Evectors(M+5,MO) VNORM(17) = RNORM*CNSTNS(2,3)*Evectors(M+5,MO) VNORM(18) = RNORM*CNSTNS(3,3)*Evectors(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)*Evectors(M+7,MO) VNORM(23) = VNORM(20)*Evectors(M+7,MO) VNORM(24) = VNORM(21)*Evectors(M+7,MO) VNORM(25) = VNORM(19)*Evectors(M+8,MO) VNORM(26) = VNORM(20)*Evectors(M+8,MO) VNORM(27) = VNORM(21)*Evectors(M+8,MO) VNORM(19) = VNORM(19)*Evectors(M+6,MO) VNORM(20) = VNORM(20)*Evectors(M+6,MO) VNORM(21) = VNORM(21)*Evectors(M+6,MO) DO 240 IXYZ = 1, numgpx CNSTX(IXYZ,1) = VNORM(7)*XDEL(IXYZ) CNSTX(IXYZ,2) = VNORM(8)*XDEL(IXYZ) CNSTX(IXYZ,3) = VNORM(9)*XDEL(IXYZ) CNSTX(IXYZ,4) = VNORM(19)*XDEL(IXYZ) CNSTX(IXYZ,5) = VNORM(20)*XDEL(IXYZ) CNSTX(IXYZ,6) = VNORM(21)*XDEL(IXYZ) 240 CONTINUE DO 242 IXYZ = 1, numgpy CNSTY(IXYZ,1) = VNORM(10)*YDEL(IXYZ) CNSTY(IXYZ,2) = VNORM(11)*YDEL(IXYZ) CNSTY(IXYZ,3) = VNORM(12)*YDEL(IXYZ) CNSTY(IXYZ,4) = VNORM(22)*YDEL(IXYZ) CNSTY(IXYZ,5) = VNORM(23)*YDEL(IXYZ) CNSTY(IXYZ,6) = VNORM(24)*YDEL(IXYZ) 242 CONTINUE DO 244 IXYZ = 1, numgpz 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) 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) 244 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, numgpx 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) 250 CONTINUE DO 252 IXYZ = 1, numgpy EXP1S(IXYZ,4) = GEXP(YDELSQ(IXYZ)*ANEG(1)) EXP1S(IXYZ,5) = GEXP(YDELSQ(IXYZ)*ANEG(2)) EXP1S(IXYZ,6) = GEXP(YDELSQ(IXYZ)*ANEG(3)) 252 CONTINUE DO 254 IXYZ = 1, numgpz EXP1S(IXYZ,7) = GEXP(ZDELSQ(IXYZ)*ANEG(1)) EXP1S(IXYZ,8) = GEXP(ZDELSQ(IXYZ)*ANEG(2)) EXP1S(IXYZ,9) = GEXP(ZDELSQ(IXYZ)*ANEG(3)) 254 CONTINUE DO 260 IXYZ = 1, numgpx EXP2SP(IXYZ,1) = GEXP(XDELSQ(IXYZ)*ANEG(4)) EXP2SP(IXYZ,2) = GEXP(XDELSQ(IXYZ)*ANEG(5)) EXP2SP(IXYZ,3) = GEXP(XDELSQ(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)) 260 CONTINUE DO 262 IXYZ = 1, numgpy EXP2SP(IXYZ,4) = GEXP(YDELSQ(IXYZ)*ANEG(4)) EXP2SP(IXYZ,5) = GEXP(YDELSQ(IXYZ)*ANEG(5)) EXP2SP(IXYZ,6) = GEXP(YDELSQ(IXYZ)*ANEG(6)) EXP3SP(IXYZ,4) = GEXP(YDELSQ(IXYZ)*ANEG(7)) EXP3SP(IXYZ,5) = GEXP(YDELSQ(IXYZ)*ANEG(8)) EXP3SP(IXYZ,6) = GEXP(YDELSQ(IXYZ)*ANEG(9)) 262 CONTINUE DO 264 IXYZ = 1, numgpz 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,7) = GEXP(ZDELSQ(IXYZ)*ANEG(7)) EXP3SP(IXYZ,8) = GEXP(ZDELSQ(IXYZ)*ANEG(8)) EXP3SP(IXYZ,9) = GEXP(ZDELSQ(IXYZ)*ANEG(9)) 264 CONTINUE C C LOOP OVER THE "CUBE" Z,Y, AND X C DO 300 IZ = 1, numgpz DO 290 IY = 1, numgpy DO 280 IG = 1, 3 DO 270 IX = 1, numgpx 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 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