C *********************************************************** C ******************** PSI/77 - PART 1 ******************** C *********************************************************** C C WILLIAM L. JORGENSEN C PURDUE UNIVERSITY, DEPARTMENT OF CHEMISTRY C WEST LAFAYETTE, INDIANA 47907, USA C PHONE 317-494-8824 C C PROGRAM HAS BEEN MODIFIED TO PRODUCE D TYPE ORBITALS. C EIGENVECTORS FROM EHT CALCULATION ARE REQUIRED. C ORBITALS ARE IN THE ORDER OF D(Z**2,XZ,YZ,X**2-Y**2,XY) C C J. GAO, JAN., 1986 C COMMON /RPLOT/ N,CO(3),CM,THE,GAM,PHI,X(850),Y(850),Z(850),IDASH, * SCALE,PERZ,KA,KO,ZO COMMON /HLCOM/ NAC(850),NC,NCURV,NCDASH(850) COMMON /IO/ IRD,ILST,IDSK1,IDSK2,IDSK3 C C THIS PROGRAM ACCEPTS A WAVEFUNCTION AND DETERMINES C ELECTRON DENSITY OR ORBITAL VALUES IN THE PLANES C WHICH ARE SPECIFIED, THE EIGENVECTORS AND MOLECULAR C COORDINATES ARE READ FROM DISK FILE 21 OR MAY BE CARD INPUT. C CLOSED CONTOURS AT THE SPECIFIED LEVEL(S) C ARE DRAWN THROUGH THE DENSITY OR ORBITAL VALUE C ARRAYS AND THE POINTS(IN ANGSTROMS) FOR EACH CONTOUR ARE STORED C AS OUTPUT IN DISK FILE 22. THE TOTAL C NUMBER OF POINTS AND THE NUMBER OF CURVES ALONG C WITH THE NUMBER OF POINTS IN EACH CURVE AND AN C INDICATOR FOR WHETHER THE CONTOUR LEVEL IS POSITIVE C OR NEGATIVE ARE STORED IN THE OUTPUT DISK FILE C 23. FOR PROPER HIDDEN LINE ELIMINATION C IN THE 3-D DRAWING PROGRAM, IT IS ESSENTIAL THAT ALL CONTOURS C WHICH ARE CONSTRUCTED BE CLOSED. IF THIS CONDITION IS NOT MET, C AN ERROR MESSAGE IS PRINTED. C C CARD INPUT IS DESCRIBED IN SUBROUTINE THREED. C C NOT ALL OF THE VARIABLES IN COMMON/RPLOT ARE C ESSENTIAL TO THIS PROGRAM, BUT HAVE BEEN RETAINED C FOR COMPATIBILITY WITH THE HIDDEN LINE PLOTTING C PROGRAM IN CASE THE TWO ARE MERGED C IRD = 5 ILST = 6 IDSK1 = 22 IDSK2 = 23 IDSK3 = 21 THE = 0. GAM = 0. PHI = 0. NC = 0 NCURV = 0 C C MAKE 3-D CONTOUR MAP C CALL THREED C C AFTER THE CURVES HAVE BEEN CREATED THE POINT AND C CURVE TOTAL INFORMATION IS STORED IN THE FILE, POINT C STOP END SUBROUTINE THREED C C PROGRAM TO PLOT VALENCE SHELL C OR CHARGE DENSITIES IN 2 OR 3 DIMENSIONS FOR ATOMS H - AR C W.L. JORGENSEN - 12/19/71, MAJOR REVISION 5/15/77 C C THIS VERSION ONLY WORKS FOR EHT AND EHT(N3D) TYPES C JOHN NASH - 4/28/90 C parameter (MXPTS=41) CHARACTER CALC*6,AUTO*4 INTEGER WFTYP COMMON /RPLOT/ NNN,CO(3),CM,THE,GAM,PHI,XX,YY,ZZ,IDASH,SCALE,PERZ, * KA,KO,ZZZ COMMON /HLCOM/ NAC(850),NC,NCURV,NCDASH(850) COMMON /IO/ IRD,ILST,IDSK1,IDSK2,IDSK3 COMMON /DENS/ DENS(MXPTS,MXPTS,MXPTS),CONTR(100),RAD(100) * ,RNORM(100),V(100,100),C(50,3),IAN(49),RCUT(18),OCMO(100) COMMON /MINMAX/ XMIN,XMAX,YMIN,YMAX,ZMIN,ZMAX COMMON /SPLIT/ XMI,YMI,ZMI,XINC,YINC,ZINC,MONE,MOLAST * ,NORB,NOCUT,NAT COMMON /EHT/ Z1S(18),Z2SP(18),Z3S(8),Z3P(8),Z3D(8) DIMENSION ZO(31),Z1(31),XX(850),YY(850),ZZ(850),CMIN(3),CMAX(3) DIMENSION CTR(15),SAV(MXPTS,MXPTS) C C SLATER EXPONENTS ARE FROM HOFFMANN + SUMMERVILLE, JACS, 1976. C DATA Z1S / 1.3,1.7,2.7,3.7,4.7,5.7,6.7,7.7,8.7,9.7,10.7,11.7,12.7, * 13.7,14.7,15.7,16.7,17.7 / DATA Z2SP / 0.0,0.0,0.65,0.975,1.3,1.625,1.95,2.275,2.425,2.925, * 3.425,3.925,4.425,4.925,5.425,5.925,6.425,6.925 / DATA Z3S / 0.733,0.95,1.167,1.383,1.75,2.122,2.183,2.25 / DATA Z3P / 0.733,0.95,1.167,1.383,1.30,1.827,1.733,2.25 / DATA Z3D / 0.0,0.0,0.0,1.383,1.40,1.50,2.033,0.0 / DATA RCUT / 5.8,4.7,13.2,9.5,7.5,6.2,5.4,4.7,4.5,3.8,14.0,11.3,9.5 * ,8.4,8.3,7.8,6.1,5.6 / DATA AU,PI,RT3 / 0.529167,3.1415926536,1.7320508 / C C V IS THE EIGENVECTOR MATRIX. DENS IS THE CHARGE DENSITY OR C ORBITAL VALUE MATRIX. IAN IS THE ARRAY OF ATOMIC NUMBERS C FOR THE ATOMS WHOSE COORDINATES, READ IN IN ANGSTROMS, ARE IN C C. CTR IS THE ARRAY OF CONTOUR LEVELS IN AU. C C CARD INPUT IS AS FOLLOWS - C CARD 1 - TYPE OF WAVEFUNCTION = EHT, EHT3D.(1A6) C EHT3D INCLUDES 3D ORBITALS ON SI, P, S, CL. C CC 6-9 = AUTO FOR AUTOMATIC DETERMINATION OF PLOTTING C PLANES AND LIMITS, ETC. IF USED, NOW SKIP TO C AFTER CARD 9. THE OTHER PARAMETERS ON THIS CARD ARE C STILL ACTIVE. C CC 11-15 IRDPOP = INCREMENT BETWEEN PLANES IN ANGSTROMS C WHEN AUTO IS USED. C CC 16 IONEMO = 1 IF ONE MO ONLY IS TO BE CARD INPUT. C SEE BELOW FOR AO ORDER. C CC 17 NOCUT = 1 IF A CUTOFF OF 0.0005 AU IS NOT TO BE C USED IN EVALUATING AN AO. THIS SLOWS EXECUTION. C CARD 2 - X,Y AND NO. OF FIRST SET OF PLANES (=1 FOR 2-D) - 2A1, C CC 5-6 = +1 OR -1 IF ORBITAL IS SYMMETRIC IN THE YZ PLANE C - +1 FOR SYMMETRIC, -1 FOR ANTISYMMETRIC (OPTIONAL). C CARD 3 - XMIN (=0 FOR AUTO SCALING) XMAX YMIN YMAX - 4F10.6 C CARD 4 - ZMIN ZMAX - 2F10.6 C CARD 5 - X',Y'AND NO. OF SECOND SET OF PLANES (=0 FOR 2-D) - 2A1 C CC 5-6 = +1 OR -1 IF ORBITAL IS SYMMETRIC IN THE YZ PLANE C - +1 FOR SYMMETRIC, -1 FOR ANTISYMMETRIC (OPTIONAL). C CARD 6 - X'MIN X'MAX Y'MIN Y'MAX (OMIT FOR 2-D) - 4F10.6 C CARD 7 - Z'MIN Z'MAX (OMIT FOR 2-D) 2F10.6 C CARD 8 - NO. OF CONTOURS, ICONN, ICONZ, NORB -4I2 C ICONN = 1 FOR NEGATIVE CONTOURS, TOO. C ICONZ=1 FOR ZERO CONTOUR AND NORB=2 FOR CHARGE DENSITY PLOT. C CARD 9 - CONTOUR LEVELS (ONLY POSITIVE NEED BE SPECIFIED WHEN C ICONN AND ICONZ ARE USED - 8F10.6). C ***IF IONEMO=1, THE GEOMETRY CARDS AND A 99 CARD FOLLOWED BY C THE ONE MO TO BE PLOTTED (8F10.6) ARE INSERTED HERE *** C LAST CARD - FIRST MO TO BE SUMMED, LAST MO, SCALE (FOR AUTO SCALE) C - 2I2, F10.6 - MONE=MOLAST=0 FOR ALL OCCUPIED MO'S. C FOR AN ORBITAL VALUE PLOT, MONE=MOLAST=THE NUMBER OF C THE DESIRED MO. IF IONEMO=1, MONE=MOLAST=01. C C----------------------------------------------------------------------- C WHEN TRANSITION ELEMENTS ARE INCLUDED IN YOUR MO PLOTS, C USE ATOMIC NUMBER RANGING FROM 19-53 FOR EACH ATOM. THE AO C COEFFICIENTS ARE TO BE READ IN FROM YOUR INFILE. FOR EACH C TS METAL, ONE ADDITIONAL CARD FOLLOWING THE XYZ COORDINATE C FOR THE TRANSITION METAL IS REQUIRED. THE FORMAT IS AS FOLLOW: C A2 ATOMIC SYMBOL C I3 NUMBER OF VALENCE ELECTRONS C I3 PRINCIPAL QUANTUM NUMBER OF S SHELL C F6.3 S SHELL EXPONENT C I3 PRINCIPAL QUANTUM NUMBER FOR THE P SHELL C F6.3 P SHELL EXPONENT C I3 PRINCIPAL QUANTUM NUMBER FOR THE D SHELL C F6.3 D SHELL EXPONENT 1 C F6.4 COEFFICIENT FOR THE FIRST D EXPONENT C F6.3 D SHELL EXPONENT 2 C F6.4 COEFFICIENT FOR THE SECOND D EXPONENT C----------------------------------------------------------------------- C C INITIALIZE FOR AUTOMATIC PROCESSING C NCT = 1 CTR(1) = 0.075 NORB = 1 ICONN = 1 ICONZ = 0 XMIN = 0.0 X1MIN = 0.0 NZ = 0 NZ1 = 0 ISYM = 0 ISYM1 = 0 ISAV = 0 N3D = 0 CTRFACT = 1.0 READ (IRD,10) CALC,AUTO,IRDPOP,IONEMO,NOCUT WRITE(*,*)' CALC,ETC ',CALC,AUTO,IRDPOP,IONEMO,NOCUT 10 FORMAT (1A6,1A4,1X,I1,4X,2I1) IF (CALC.EQ.'EHT3D ') N3D = 1 IF (CALC.NE.'EHT '.AND.CALC.NE.'EHT3D ') THEN WRITE (ILST,*) 'PSI1 IS NOT SET UP TO DO THE CALCULATIONS ON' WRITE (ILST,*) 'THIS BASIS SET. YOU CAN USE EHT AND EHT3D ONLY' STOP ENDIF IF (AUTO.NE.'AUTO') THEN C C READ ABSCISSA AND ORDINATE MIN AND MAX - 0 = AUTO. C READ (IRD,30) XMIN,XMAX READ (IRD,30) YMIN,YMAX READ (IRD,30) ZMIN,ZMAX ENDIF 30 FORMAT (2F10.6) C C OBTAIN PLOTTING DATA FROM DISK C IRD1 = IDSK3 IF (IONEMO.EQ.1) IRD1 = IRD READ (IRD1,200) ICHG WRITE(*,*)' CHARGE = ',ICHG NOEL = -ICHG N = 0 NBASIS = 0 NAT = 0 NDCALS = 0 C C READ COORDINATES AND DETERMINE NORMALIZING FACTORS FOR AOS. C THE FORM OF THE SLATER ORBITALS MAY BE FOUND IN I.G. CSIZMADIA, C THEORY AND PRACTICE OF MO CALCULATIONS,ELSEVIER,1976, P 313. C 40 NAT = NAT+1 READ (IRD1,50) IN,(C(NAT,J),J=1,3) WRITE(*,*)' ATOM ',IN,(C(NAT,J),J=1,3) 50 FORMAT (I2,8X,3F10.6) IF (IN.NE.99) THEN IAN(NAT) = IN N3FLG = 0 IF (IN.LE.2) THEN N = N+1 NBASIS = NBASIS+2 NOEL = NOEL+IN RNORM(N) = SQRT((Z1S(IN)**3)/PI) GO TO 40 ENDIF IF (WFTYP.EQ.1) THEN NOEL = NOEL+2 N = N+1 RNORM(N) = SQRT((Z1S(IN)**3)/PI) ENDIF IF (IN.GE.11) THEN IF (WFTYP.EQ.0) GO TO 80 ENDIF NOEL = NOEL+8 IF (IN.LT.11) NOEL = NOEL-10+IN RN = SQRT((Z2SP(IN)**5)/PI) N = N+1 RNORM(N) = RN/RT3 60 DO 70 I = 1, 3 N = N+1 RNORM(N) = RN 70 CONTINUE IF (N3FLG.EQ.1) GO TO 90 IF (IN.LT.11) GO TO 40 80 IF (IN.GT.18) GO TO 110 IN = IN-10 NOEL = NOEL+IN N = N+1 RNORM(N) = SQRT(2.*(Z3S(IN)**7)/(45.*PI)) RN = SQRT(2.*(Z3P(IN)**7)/(15.*PI)) N3FLG = 1 GO TO 60 90 IF (IN.LT.4) GO TO 40 IF (N3D.EQ.0) GO TO 40 RN = SQRT(2.*(Z3D(IN)**7)/(3.*PI)) DO 100 I = 1, 5 N = N+1 RNORM(N) = RN 100 CONTINUE I = N-4 RNORM(I) = RN*0.5/RT3 I = N-1 RNORM(I) = RN*0.5 GO TO 40 110 CALL DORBIT (IN,N,NOEL,NAT,IRD1) GO TO 40 ENDIF NAT = NAT-1 CALL DRAMNP (C,NAT,CO,ICM,CMIN,CMAX) C C CONVERT COORDINATES TO ATOMIC UNITS C C IT'S FASTER TO MULTIPLY THAN TO DIVIDE ON MOST COMPUTERS C AUINV = 1.0/AU DO 130 I = 1, NAT DO 120 J = 1, 3 C(I,J) = C(I,J)*AUINV 120 CONTINUE 130 CONTINUE WRITE(*,*)' ATOM ',IN,(C(NAT,J),J=1,3) C C NAT=NO. OF ATOMS% NOEL=NO. OF ELECTRONS% N= NO. C OF BASIS FUNCTIONS. C NORB=1,2 FOR ORBITAL VALUES OR CHARGE DENSITIES C NMO = (NOEL+1)/2 C C READ IN EIGENVECTORS C IT IS ASSUMED THAT THE EIGENVECTORS HAVE BEEN NORMALIZED TO 1 C ELECTRON WITH THE OVERLAP MATRIX INCLUDED. C C ALTHOUGH THIS DOES NOT APPEAR NECESSARY - JOHNN J. NASH 4/29/90 C WRITE (ILST,140) N,NAT,NOEL 140 FORMAT (' NBASIS = ',I3,' NATOMS = ',I3,' NOEL = ',I3) IF (IONEMO.NE.0) THEN READ (IRD,150) (V(I,1),I=1,N) WRITE(*,*)' EIGENVECTORS ',(V(I,1),I=1,N) ELSE READ (IDSK3,150) ((V(I,J),I=1,N),J=1,N) 150 FORMAT (8f10.6) ENDIF C C ESTABLISH MO OCCUPANCIES. THESE MAY ALSO BE READ C IN FOR OPEN SHELLS C C C READ IN MONE AND MOLAST C READ (IRD,190) MONE,MOLAST,SCALE WRITE(*,*)' MONE,MOLAST,SCALE ',MONE,MOLAST,SCALE 190 FORMAT (2I2,F10.6) C C READ THE NO. OF CONTOURS AND THEN THE LEVELS. C READ (IRD,200) NCT,ICONN,ICONZ,NORB WRITE(*,*)' NCT,ICONN,ICONZ,NORB ',NCT,ICONN,ICONZ,NORB 200 FORMAT (4I2) READ (IRD,210) (CTR(I),I=1,NCT) WRITE(*,*)' CTR ',CTR(1) 210 FORMAT (8F10.6) C C DEFAULT VALUES C IF (MONE.EQ.0) MONE = 1 IF (MOLAST.EQ.0) MOLAST = NMO IF (IONEMO.EQ.1) MOLAST = 1 IF (SCALE.EQ.0.0) SCALE = 1.0 NUM = MOLAST-MONE+1 C C DEFAULT VALUES C IF (NCT.EQ.0) CTR(1) = 0.075 IF (NCT.EQ.0) NCT = 1 IF (NORB.EQ.0) NORB = 1 IF (NORB.EQ.1) ICONN = 1 WRITE (ILST,20) CALC,IONEMO,NORB,IRDPOP 20 FORMAT (' CALC = ',1A6,' IONEMO = ',I2,' NORB = ',I2,' IRDPOP=', * I2) C C NORMALIZE MO S TO PROPER OCCUPANCIES C C THIS WILL NOW BE DONE AT THE END SO MO'S CAN BE SIMULTANEOUSLY C BE PLOTTED, IMPOSSIBLE IF IT'S DONE HERE C C DO 180 I = 1, N C FACT = SQRT(OCMO(I)) C DO 170 J = 1, N C V(J,I) = FACT*V(J,I) C 170 CONTINUE C 180 CONTINUE C ENDIF C ENDIF C C GREATER RESOLUTION IS OBTAINABLE BY INCREASING THE C SIZE OF MXPTS WHICH IS THE DIMENSIONS OF THE DENSITY C OR ORBITAL VALUE ARRAY. OF COURSE, THE DIMENSION C STATEMENTS WILL ALSO HAVE TO BE APPROPRIATELY FIXED. C MXPTS MUST BE ODD. C C MXPTS = 51 C see the parameter statement at top for this parameter C SPACES = FLOAT(MXPTS-1)*AU C C DETERMINE DEFAULT PLANES C C C DETERMINE DEFAULT RANGES, IF XMIN = 0. INCREASING SCALE C INCREASES THE RANGE OF THE PLOT. C WRITE(*,*)' AUTO, SCALE ',AUTO,SCALE IF (AUTO.EQ.'AUTO') THEN R = 1.3*SCALE XMIN = CMIN(1)-R YMIN = CMIN(2)-R XMAX = CMAX(1)+R YMAX = CMAX(2)+R ZMIN = CMIN(3)-R ZMAX = CMAX(3)+R ENDIF XMI = XMIN/AU YMI = YMIN/AU ZMI = ZMIN/AU XINC = (XMAX-XMIN)/SPACES YINC = (YMAX-YMIN)/SPACES ZINC = (ZMAX-ZMIN)/SPACES WRITE (ILST,220) 'X',XMIN,XMAX,'Y',YMIN,YMAX,'Z',ZMIN,ZMAX 220 FORMAT ('0',3(2X,A1,2F10.4)) C C SYMMETRY SET-UP C ISFLAG = 0 C MXXPTS = MXPTS C C MAKE A MAP FOR EACH SLICE C IF (CALC.EQ.'EHT '.OR.CALC.EQ.'EHT3D ') THEN CALL EHTMO (N3D) ELSE WRITE (*,*) ' NO COMPUTATIONS DONE' ENDIF RETURN END C C SUBROUTINE DRAMNP (C,NAT,CO,ICM,CMIN,CMAX) DIMENSION C(50,3),CO(3),CMIN(3),CMAX(3) C C THIS ROUTINE DETERMINES CO AND CM WHICH ARE USED C FOR AUTOMATIC SCALING OF PLOTTING AREAS C C THE ROUTINE WAS ADAPTED FROM THE ROUTINE CALLED C DRAMOL WHICH IS USED IN THE HIDDEN LINE PART OF C THE PROGRAM C CM = -100. DO 20 I = 1, 3 CMI = 100. CMA = -100. DO 10 J = 1, NAT P = C(J,I) CMI = MIN(CMI,P) CMA = MAX(CMA,P) 10 CONTINUE CO(I) = (CMA+CMI)/2. CMIN(I) = CMI CMAX(I) = CMA P = CMA-CMI IF (P.GT.CM) THEN ICM = I CM = P ENDIF 20 CONTINUE RETURN END C C SUBROUTINE DORBIT (IN,N,NOEL,NAT,IRD) PARAMETER (MXPTS=41) C C SETS UP D ORBITALS C J. GAO, JAN., 1986 C C COEFFECIENTS OF ORBITALS : C NS = SQRT( (2ZS)^(2N+1)/(4(2N)!*PI) ) * R^(N-1)EXP(-ZSR) C C NP = SQRT( 3(2ZP)^(2N+1)/(4(2N)!*PI) ) * R^(N-2)EXP(-ZPR)X C C ND = SQRT( 15(2ZD)^(2N+1)/(4(2N)!*PI) ) --- AND ADJUSTMENTS C FOR EACH D'S C COMMON /DENS/ DENS(MXPTS,MXPTS,MXPTS),CONTR(100),RAD(100) * ,RNORM(100),V(100,100),C(50,3),IAN(49),RCUT(18),OCMO(100) COMMON /DORB/ NSHELL(25),ZTS(25),ZTP(25),ZTD(25,2),DCONST(25,2), * RRNORM(100) DIMENSION ZD(2),CD(2) DATA PI,RT3 / 3.1415926536,1.7320508 / C READ (IRD,10) SYM,NVLE,NS,ZS,NP,ZP,ND,ZD(1),CD(1),ZD(2),CD(2) 10 FORMAT (A2,I3,I3,F6.3,I3,F6.3,I3,F6.3,F6.4,F6.3,F6.4) NSHELL(NAT) = ND ZTS(NAT) = ZS ZTP(NAT) = ZP ZTD(NAT,1) = ZD(1) ZTD(NAT,2) = ZD(2) DCONST(NAT,1) = CD(1) DCONST(NAT,2) = CD(2) NOEL = NOEL+NVLE N = N+1 IF (NS.LT.3.OR.NP.LT.3.OR.ND.LT.3) THEN WRITE (*,*) ' PROGRAM STOPPED IN SUBROUTINE DORBIT' STOP ENDIF GO TO (20,20,20,30,40,50), NS 20 RNORM(N) = SQRT(2.*ZS**7/(45.*PI)) GO TO 60 30 RNORM(N) = SQRT(ZS**9/(315.*PI)) GO TO 60 40 RNORM(N) = SQRT(2.*ZS**11/(14175.*PI)) GO TO 60 50 RNORM(N) = SQRT(2.*ZS**13/(467775.*PI)) 60 GO TO (70,70,70,90,110,130), NP 70 RN = SQRT(2.*ZP**7/(15.*PI)) DO 80 I = 1, 3 N = N+1 RNORM(N) = RN 80 CONTINUE GO TO 150 90 RN = SQRT(ZP**9/(105.*PI)) DO 100 I = 1, 3 N = N+1 RNORM(N) = RN 100 CONTINUE GO TO 150 110 RN = SQRT(2.*ZP**11/(4725.*PI)) DO 120 I = 1, 3 N = N+1 RNORM(N) = RN 120 CONTINUE GO TO 150 130 RN = SQRT(2.*ZP**13/(155925.*PI)) DO 140 I = 1, 3 N = N+1 RNORM(N) = RN 140 CONTINUE 150 GO TO (160,160,160,180,200), ND 160 RN = SQRT(2.*ZD(1)**7/(3.*PI)) RRN = SQRT(2.*ZD(2)**7/(3.*PI)) DO 170 I = 1, 5 N = N+1 RNORM(N) = RN RRNORM(N) = RRN 170 CONTINUE I = N-4 RNORM(I) = RNORM(I)*0.5/RT3 RRNORM(I) = RRNORM(I)*0.5/RT3 I = N-1 RNORM(I) = 0.5*RNORM(I) RRNORM(I) = 0.5*RRNORM(I) GO TO 220 180 RN = SQRT(ZD(1)**9/(21.*PI)) RRN = SQRT(ZD(2)**9/(21.*PI)) DO 190 I = 1, 5 N = N+1 RNORM(N) = RN RRNORM(N) = RRN 190 CONTINUE I = N-4 RNORM(I) = 0.5*RNORM(I)/RT3 RRNORM(I) = 0.5*RNORM(I)/RT3 I = I-1 RNORM(I) = 0.5*RNORM(I) RRNORM(I) = 0.5*RRNORM(I) GO TO 220 200 RN = SQRT(2.*ZD(1)**11/(945.*PI)) RRN = SQRT(2.*ZD(2)**11/(945.*PI)) DO 210 I = 1, 5 N = N+1 RNORM(N) = RN RRNORM(N) = RRN 210 CONTINUE I = N-4 RNORM(I) = 0.5*RNORM(I)/RT3 RRNORM(I) = 0.5*RRNORM(I)/RT3 I = I-1 RNORM(I) = RNORM(I)*0.5 RRNORM(I) = RRNORM(I)*0.5 220 CONTINUE RETURN END C SUBROUTINE EHTMO (N3D) C PARAMETER (MXPTS=41) COMMON /DORB/ NSHELL(25),ZTS(25),ZTP(25),ZTD(25,2),DCONST(25,2), * RRNORM(100) COMMON /SPLIT/ XMI,YMI,ZMI,XINC,YINC,ZINC,MONE,MOLAST,NORB,NOCUT, * NAT COMMON /DENS/ DENS(MXPTS,MXPTS,MXPTS),CONTR(100),RAD(100), * RNORM(100),V2(100,100),C(50,3),IAN(49),RCUT(18),OCMO(100) COMMON /MINMAX/ XMIN,XMAX,YMIN,YMAX,ZMIN,ZMAX COMMON /EHT/ Z1S(18),Z2SP(18),Z3S(8),Z3P(8),Z3D(8) COMMON /IO/ IRD,ILST,IDSK1,IDSK2,IDSK3 DIMENSION XDEL(MXPTS),YDEL(MXPTS),ZDEL(MXPTS) DIMENSION XDELSQ(MXPTS),YDELSQ(MXPTS),ZDELSQ(MXPTS) DIMENSION PTMP(MXPTS,3),V(100,100) DIMENSION YDELRSQ(MXPTS),IND(100),X(MXPTS),Y(MXPTS),Z(MXPTS) DIMENSION TMPYZ(MXPTS) C C THIS SUBROUTINE IS A SIMPLIFIED VERSION OF THE DENSITY C COMPUTATION SECTION OF PSI77. C THIS ROUTINE ONLY CALCULATES EHT(3D) TYPE MO'S, THUS REMOVING C THE LARGE NUMBER OF IF - GOTO'S IN THE ORIGINAL. C C THE PURPOSE FOR SPLITTING IT INTO SEPARATE MODULES IS SIMPLE, C THERE IS NO REASON TO CHECK EACH TIME THROUGH THE LOOP WHICH C TYPE OF WAVE FUNCTION IS BEING USED, ONCE YOU KNOW IT'S EHT(3D) C IT'S NEVER GOING TO CHANGE DURING THE PROGRAM RUN. C DAN SEVERANCE 9/3/87 C NOCUT = 0 DO 10 NZZ = 1, MXPTS DO 10 NY = 1, MXPTS DO 10 NX = 1, MXPTS DENS(NX,NY,NZZ) = 0.0 10 CONTINUE WRITE (*,*) ' MXPTS ',MXPTS X(1) = XMI Y(1) = YMI Z(1) = ZMI DO 20 I = 2, MXPTS X(I) = XINC+X(I-1) Y(I) = YINC+Y(I-1) Z(I) = ZINC+Z(I-1) 20 CONTINUE NAO = 0 Do 501 NS = 1,nat in = ian(ns) IF( IN .LT. 2 ) THEN NAO = NAO + 1 ELSEIF(IN .LT. 18) THEN NAO = NAO + 4 ELSEIF(IN .GT. 18) THEN NAO = NAO + 9 ENDIF 501 CONTINUE do 589 m=1,na0 write(*,*)' rnorm(',m,') = ',rnorm(m) 589 continue DO 500 MO=MONE,MOLAST DO 500 M = 1, NAO V(M,MO) = V2(M,MO)*RNORM(M) 500 CONTINUE M = 1 MO = MONE DO 250 NS = 1, NAT VTMP = V(M,MO) IN = IAN(NS) WRITE (*,*) ' NS, IAN(NS) ',NS,IN C C RCUT CONTAINS THE DISTANCES BEYOND WHICH THE RADIAL PART C OF THE MOST DIFFUSE AO FOR AN ATOM IS .LT. .0005 AU. FOR C NORMAL CONTOUR LEVELS (.GT. .01 AU) SUCH CONTRIBUTIONS ARE C NEGLIGIBLE. THIS FEATURE IS DISABLED BY NOCUT. C C IF (NOCUT.EQ.1.OR.R(NS).LE.RCUT(IN)) THEN C DISABLED FOR NOW C c3 = c(ns,3) c2 = c(ns,2) c1 = c(ns,1) DO 30 NXYZ = 1, MXPTS ZDEL(NXYZ) = Z(NXYZ)-C3 ZDELSQ(NXYZ) = ZDEL(NXYZ)*ZDEL(NXYZ) 30 CONTINUE DO 32 NXYZ = 1, MXPTS YDEL(NXYZ) = Y(NXYZ)-C2 YDELSQ(NXYZ) = YDEL(NXYZ)*YDEL(NXYZ) 32 CONTINUE DO 33 NXYZ = 1, MXPTS XDEL(NXYZ) = X(NXYZ)-C1 XDELSQ(NXYZ) = XDEL(NXYZ)*XDEL(NXYZ) 33 CONTINUE IF (IN.LE.2) THEN Z1SIN = -Z1S(IN) DO 90 NZ = 1, MXPTS ZDSQ = ZDELSQ(NZ) DO 80 NY = 1, MXPTS YDRSQ = YDELSQ(NY)+ZDSQ DO 50 NX = 1, MXPTS R = XDELSQ(NX)+YDRSQ R = SQRT(R) R = EXP(Z1SIN*R) DENS(NX,NY,NZ) = DENS(NX,NY,NZ)+R*VTMP 50 CONTINUE 80 CONTINUE 90 CONTINUE M = M+1 ELSEIF (IN.LE.10) THEN C C 2S AND 2P ORBITALS C Z2SPIN = -Z2SP(IN) VTMPX = V(M+1,MO) VTMPY = V(M+2,MO) VTMPZ = V(M+3,MO) DO 100 NXYZ = 1, MXPTS PTMP(NXYZ,1) = VTMPX*XDEL(NXYZ) PTMP(NXYZ,2) = VTMPY*YDEL(NXYZ) PTMP(NXYZ,3) = VTMPZ*ZDEL(NXYZ) 100 CONTINUE DO 170 NZ = 1, MXPTS ZDSQ = ZDELSQ(NZ) TMPP = PTMP(NZ,3) DO 160 NY = 1, MXPTS YDRSQ = YDELSQ(NY)+ZDSQ tmp = PTMP(NY,2) + TMPP DO 120 NX = 1, MXPTS RSQ2 = XDELSQ(NX)+YDRSQ R = SQRT(RSQ2) RADPS = EXP(Z2SPIN*R) DENS(NX,NY,NZ) = DENS(NX,NY,NZ)+RADPS* * (PTMP(NX,1)+TMP ) RADPS = RADPS*R DENS(NX,NY,NZ) = DENS(NX,NY,NZ)+RADPS*VTMP 120 CONTINUE C C The next and previous code was changed such that RADP(NX) has the C exponent only to start, then rad = radp*r, instead of : C Rad = exp C radp = rad * (stuff moved into PTMP) C rad = rad *(stuff moved into Vtmp) C C Now V is premultiplied by the Rnorm factor, only PTMP need calc. to C Include Xdel,Ydel, and Zdel C 160 CONTINUE 170 CONTINUE M = M+4 ELSEIF (IN.LE.18) THEN IN = IN - 10 C C 3S AND 3P ORBITALS C Z3SIN = -Z3S(IN) Z3PIN = -Z3P(IN) VTMPX = V(M+1,MO) VTMPY = V(M+2,MO) VTMPZ = V(M+3,MO) DO 1001 NXYZ = 1, MXPTS PTMP(NXYZ,1) = VTMPX*XDEL(NXYZ) PTMP(NXYZ,2) = VTMPY*YDEL(NXYZ) PTMP(NXYZ,3) = VTMPZ*ZDEL(NXYZ) 1001 CONTINUE DO 1701 NZ = 1, MXPTS ZDSQ = ZDELSQ(NZ) TMPP = PTMP(NZ,3) DO 1601 NY = 1, MXPTS YDRSQ = YDELSQ(NY)+ZDSQ tmp = PTMP(NY,2) + TMPP DO 1201 NX = 1, MXPTS RSQ2 = XDELSQ(NX)+YDRSQ R = SQRT(RSQ2) RADS = EXP(Z3SIN*R) RADS = RADS * RSQ2 DENS(NX,NY,NZ) = DENS(NX,NY,NZ)+RADS*VTMP RADP = EXP(Z3PIN*R) RADP = RADP * R DENS(NX,NY,NZ) = DENS(NX,NY,NZ)+RADP* * (PTMP(NX,1)+TMP ) 1201 CONTINUE C C The next and previous code was changed such that RADP(NX) has the C exponent only to start, then rad = radp*r, instead of : C Rad = exp C radp = rad * (stuff moved into PTMP) C rad = rad *(stuff moved into Vtmp) C C Now V is premultiplied by the Rnorm factor, only PTMP need calc. to C Include Xdel,Ydel, and Zdel C 1601 CONTINUE 1701 CONTINUE M = M+4 ELSE DO 1702 NZ = 1, MXPTS ZDSQ = ZDELSQ(NZ) DO 1602 NY = 1, MXPTS YDRSQ = YDELSQ(NY)+ZDSQ DO 1202 NX = 1, MXPTS RSQ2 = XDELSQ(NX)+YDRSQ R = SQRT(RSQ2) CALL MOCALC(IN,M,NS,XDEL(NX),YDEL(NY),ZDEL(NZ) * ,DENVAR,R) DENS(NX,NY,NZ) = DENS(NX,NY,NZ)+DENVAR 1202 CONTINUE 1602 CONTINUE 1702 CONTINUE M = M + 9 ENDIF 250 CONTINUE OPEN(17,FORM='UNFORMATTED',STATUS='NEW') WRITE(17)MXPTS,NAT WRITE(17)(IAN(I),I=1,NAT) WRITE(17)((C(I,J),J=1,3),I=1,NAT) WRITE(17)(((DENS(NX,NY,NZ),NX=1,MXPTS),NY=1,MXPTS),NZ=1,MXPTS) WRITE(17)XMIN,XMAX,YMIN,YMAX,ZMIN,ZMAX CLOSE(17) RETURN END C C SUBROUTINE MOCALC (IN,M,NS,XDEL,YDEL,ZDEL,CONTR,R) C C EVALUATES SLATER TYPE ORBITALS C J. GAO, JAN., 1986 C C VALENCE ORBITALS ARE: (3D,4S,4P), (4D,5S,5P), & (5D,6S,6P) C PARAMETER(MXPTS=41) COMMON /SPLIT/ XMI,YMI,ZMI,XINC,YINC,ZINC,MONE,MOLAST,NORB,NOCUT, * NAT COMMON /DENS/ DENS(MXPTS,MXPTS,MXPTS),JUNK(100),RAD(100), * RNORM(100),V(100,100),C(50,3),IAN(49),RCUT(18),OCMO(100) COMMON /RPLOT/ NNN,CO(3),CM,THE,GAM,PHI,XX,YY,ZZ,IDASH,SCALE,PERZ, * KA,KO,ZZZ COMMON /DORB/ NSHELL(25),ZTS(25),ZTP(25),ZTD(25,2),DCONST(25,2), * RRNORM(100) DIMENSION XX(850),YY(850),ZZ(850) C ND = NSHELL(NS) GO TO (10,10,10,20,30) ND 10 RR = R*R RADS = RR*R*RNORM(M)*EXP(-ZTS(NS)*R) RADP = RR*RNORM(M+1)*EXP(-ZTP(NS)*R) RADD1 = EXP(-ZTD(NS,1)*R) RADD2 = EXP(-ZTD(NS,2)*R) GO TO 40 20 RR = R*R RADS = RR*RR*RNORM(M)*EXP(-ZTS(NS)*R) RADP = RR*R*RNORM(M+1)*EXP(-ZTP(NS)*R) RADD1 = R*EXP(-ZTD(NS,1)*R) RADD2 = R*EXP(-ZTD(NS,2)*R) GO TO 40 30 RR = R*R RADS = RR*RR*R*RNORM(M)*EXP(-ZTS(NS)*R) RADP = RR*RR*RNORM(M+1)*EXP(-ZTP(NS)*R) RADD1 = RR*EXP(-ZTD(NS,1)*R) RADD2 = RR*EXP(-ZTD(NS,2)*R) 40 CONTR = V(M,MONE)*RADS CONTR = CONTR+V(M+1,MONE)*XDEL*RADP CONTR = CONTR+V(M+2,MONE)*YDEL*RADP CONTR = CONTR+V(M+3,MONE)*ZDEL*RADP MJ = M+4 CONTR = CONTR+V(MJ,MONE)*(XDEL**2-YDEL**2)*(DCONST(NS, * 1)*RNORM(MJ)*RADD1+DCONST(NS,2)*RRNORM(MJ)*RADD2) MJ = M+5 CONTR = CONTR+V(MJ,MONE)*(3.*ZDEL**2-R*R)*(DCONST(NS,1)* * RNORM(MJ)*RADD1+DCONST(NS,2)*RRNORM(MJ)*RADD2) MJ = M+6 CONTR = CONTR+V(MJ,MONE)*XDEL*YDEL*(DCONST(NS,1)*RNORM * (MJ)*RADD1+DCONST(NS,2)*RRNORM(MJ)*RADD2) MJ = M+7 CONTR = CONTR+V(MJ,MONE)*XDEL*ZDEL*(DCONST(NS,1)*RNORM * (MJ)*RADD1+DCONST(NS,2)*RRNORM(MJ)*RADD2) MJ = M+8 CONTR = CONTR+V(MJ,MONE)*YDEL*ZDEL*(DCONST(NS,1)*RNORM * (MJ)*RADD1+DCONST(NS,2)*RRNORM(MJ)*RADD2) RETURN END