C C *********************************************************** C ******************** PSI/88 - PART 3 ******************** C *********************************************************** C C Version 1.0 Any questions to the author should specify C the version being used. C PROGRAM PSI2 C C PROGRAM FOR PLOTTING CHARGE DENSITIES OR ORBITAL VALUES C IN TWO OR THREE DIMENSIONS WITH OR WITHOUT HIDDEN LINE C ELIMINATION. C C WILLIAM L. JORGENSEN C Daniel L. Severance C DEPARTMENT OF CHEMISTRY, Yale University C New Haven, CT 06511, USA C C MODIFIED FOR PLOTTING INPUT SYMBOLS C JAN., 1986 C JIALI GAO C C Heavy Modifications of the Hidden Line Section to Increase C Speed by a Factor of 10 or More 8-87 to 3-88. Introduction C of Code to Properly Handle Indentations and Doughnuts as Well. C Dan Severance, Purdue. 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 William Jorgensen at Purdue University C The name of the University or William Jorgensen may not be used to endorse C or promote products derived from this software without specific prior C written permission. The author is 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 AC IS THE ARRAY OF ALL POINTS IN THE CURVES C WHICH WERE GENERATED IN THE FIRST PART OF C THE PROGRAM. IT IS NECESSARY TO HAVE ALL OF C THEM IN CORE AT ONCE FOR PROPER HIDDEN LINE C ELIMINATION. NAC CONTAINS THE NUMBER OF POINTS C IN EACH CURVE. NBEG CONTAINS THE STARTING C SUBSCRIPT FOR EACH CURVE IN AC. ACMIN AND ACMAX C CONTAIN THE EXTREMA FOR EACH CURVE AFTER IT C HAS BEEN TRANSFORMED. C ** THE DIMENSION OF AC MUST BE CHANGED IF NC .GT.75000 ** C ** ARRAYS DIMENSIONED 1024 MUST BE CHANGED IF NCURV .GT. 1024 C THE, GAM AND PHI ARE THE ROTATION ANGLES C IN THE XY, YZ AND XZ PLANES RESPECTIVELY. C ONCE THE CURVES HAVE BEEN GENERATED IN PART ONE C MANY VALUES OF THE, GAM AND PHI MAY BE TRIED C UNTIL THE PLOT ADOPTS THE VIEW YOU DESIRE. C THE SIZE OF THE PLOT IS CONTROLLED BY SCALE C AS BEFORE. PERZ PROVIDES THE PLOT WITH PERSPECTIVE. C A LITTLE PERSPECTIVE IS ALMOST ALWAYS NEEDED C TO RELIEVE ANY AMBIGUITY IN THE STEREOCHEMISTRY. C THE DEFAULT VALUE OF 0.0 FOR PERZ WILL CAUSE C AN ADEQUATE DEGREE OF PERSPECTIVE TO BE PROVIDED. C PERSPECTIVE IS ELIMINATED BY TAKING PERZ VERY LARGE. C C CARD INPUT - C CARD 1 - TITLE (A120) C CARD 2 - SUBTITLE (CENTERED A40) C CARD 3 - IRDXYZ = 01, IF COORDINATES ARE TO BE READ FROM CARDS C THIS OPTION NO LONGER DOES ANYTHING, COORDINATES ARE C ALWAYS READ FROM THE INPUT FILE (IRDXYZ ALWAYS SET TO 01 C MOL = 01, FOR DRAWING OF STRUCTURE ONLY C NHL = 01 FOR NO HIDDEN LINE ELIMINATION (3I2) C CARD 4 - NO. OF BONDS (I2) =0 FOR AUTOMATIC C CARD 5 - BOND PAIRS (2I2) IF NO. OF BONDS SPECIFIED C CARD 6 - MOLECULAR TITLE AND COORDINATES, FOLLOWED BY 99 CARD C C LAST CARD - THETA, GAMMA AND PHI ROTATION ANGLES, SCALE 4F10.6 C C LAST CARD + 1 - IF ANY ATOM HAS ATOMIC NUMBER GREATER THAN 18, C THE ATOMIC SYMBOL IS READ IN. (A2,I1, # OF LETTER C EVERY ATOM NEEDS ONE CARD! C LAST CARD + 2 - AN OPTION CARD IS NEEDED TO CONTROL COLOR PLOT C OR DASHLINE PLOT. 0 FOR DASHLINE, 1 FOR COLOR, C AND 2 FOR BOTH. C DEFAULT -1- COLOR. C COMMON /POINTS/ AC(75000,3) COMMON /HLCOM/ NAC(1024),NBEG(1024),ACMIN(1024,3),ACMAX(1024,3) COMMON /RPLOT/ N,CO(3),CM,THE,GAM,PHI,X(1024),Y(1024),Z(1024), * IDASH,SCALE,PERZ COMMON /BONDS/ NBND,IB(35),IC(35) C DIMENSION IAN(60),C(1024,3),IATSMN(60) CHARACTER TITLE*120,SUBT*40,IATSYM(60)*2 EQUIVALENCE (C(1,1),X(1)) DATA IRD,ILST/ 5,6 / C C INITIALIZE PLOTTING C CALL PLOTS (0) C READ (IRD,10) TITLE,SUBT 10 FORMAT (A/A) READ (IRD,20) IRDXYZ,MOL,NHL,FC IRDXYZ=1 C C THE DEFAULT FACTOR IS TO OPEN A 6.5"X6.5" VIEWPORT C IF (FC.EQ.0.0E+0) FC = 0.650E+0 CALL FACTOR (FC) READ (IRD,20) NBND 20 FORMAT (3I2,F4.3) IF (NBND.NE.0) THEN DO 30 I = 1, NBND READ (IRD,20) IB(I),IC(I) 30 CONTINUE ENDIF C C THE COORDINATES ARE READ IN FROM THE PLOTTING FILE C THIS IS THE ONLY INFORMATION REQUIRED FROM THE PLOT C FILE FOR THIS PART OF THE PROGRAM. C N = 1 READ (IRD,20) JUNK 40 READ (IRD,50) IAN(N),(C(N,J),J=1,3) 50 FORMAT (I2,8X,3F10.6) IF (IAN(N).NE.99) THEN N = N+1 GO TO 40 ENDIF N = N-1 READ (IRD,60) THE,GAM,PHI,SCALE 60 FORMAT (4F10.4) IF (SCALE.EQ.0.0E+0) SCALE = 1.0E+0 INIT = 1 C C READ IN ATOMIC LABELS FOR ATOMS WITH IAN > 18 C DO 70 I = 1, N IF (IAN(I).GT.18) READ (IRD,80) IATSYM(I),IATSMN(I) 70 CONTINUE C C The field of IATSMN(I) was changed from 1 to 2 (I2) to allow C more than 9 atoms to be user-defined. Although it is unlikely C that this would happen with a typical organic molecule, the EHMO C program will undoubtedly need more at some point. C C JJN 11-21-90 C 80 FORMAT (A2,I2) C READ (IRD,'(I2)',ERR=90) JLINES GO TO 100 90 JLINES = 1 100 CONTINUE C CALL DRAMOL (C,N,IAN,TITLE,CO,THE,GAM,PHI,CM,PERZ,SCALE,INIT, * IATSYM,IATSMN) IF (MOL.NE.1) CALL HIDPLT (SUBT,NHL,JLINES) C C IF PLOTTING TO A TERMINAL SCREEN, YOU MAY NEED A ROUTINE HERE C TO MAKE THE TERMINAL WAIT FOR YOU TO FINISH BEFORE CLEARING THE C SCREEN. (WAIT FOR CARRIAGE-RETURN, OR MOUSE-CLICK, ETC. BEFORE C EXECUTING THE LAST CALL TO PLOT. C CALL PLOT (0.0E+0,0.0E+0,999) WRITE (ILST,110) TITLE,SUBT 110 FORMAT (1X,A/) WRITE (ILST,120) THE,GAM,PHI,SCALE 120 FORMAT ('0','THETA = ',F6.1,2X,'GAMMA = ',F6.1,2X,'PHI = ',F6.1/1X * ,'SCALE = ',F6.3/) STOP END C C SUBROUTINE DRAMOL (C,NAT,IAN,TITLE,CO,THE,GAM,PHI,CM,PERZ,SCALE, * INIT,IATSYM,IATSMN) COMMON /MIN/ YMOLMN COMMON /BONDS/ NBND,IB,IC DIMENSION C(1024,3),IAN(60),IATSMN(60),CO(3),NCHAR1(18) DIMENSION IB(35),IC(35),BX(4),BY(4),COVRAD(18) CHARACTER*(*) TITLE CHARACTER*2 ATS(18),IATSYM(60) C C ROUTINE TO DRAW MOLECULAR FRAMEWORK C INIT=1 TO INITIALIZE - CO,CM,PERZ RETURNED. C C IS DESTROYED. NAT = NO. OF ATOMS. C DATA ATS / 'H ','HE','LI','BE','B ','C ','N ','O ','F ','NE','NA', * 'MG','AL','SI','P ','S ','CL','AR'/ C C THIS ARRAY CONTAINS APPROXIMATE COVALENT RADIUS DATA FOR C AUTOMATIC DETERMINATION OF BONDING. C DATA COVRAD / 0.350E+0,0.0E+0,1.40E+0,1.060E+0,0.840E+0,0.770E+0, * 0.740E+0,0.740E+0,0.640E+0,0.0E+0,1.570E+0,1.450E+0,1.30E+0, * 1.220E+0,1.20E+0,1.140E+0,1.10E+0,0.0E+0 / C C THIS ARRAY CONTAINS THE NUMBER OF LETTERS IN THE ATOMIC SYMBOL C FOR EACH ATOM IN ARRAY IATS. C DATA NCHAR1 / 1,2,2,2,1,1,1,1,1,2,2,2,2,2,1,1,2,2 / IF (INIT.EQ.1) THEN CM = -100.0E+0 DO 20 I = 1, 3 CMI = 100.0E+0 CMA = -100.0E+0 DO 10 J = 1, NAT P = C(J,I) CMI = MIN(CMI,P) CMA = MAX(CMA,P) 10 CONTINUE CO(I) = (CMA+CMI)*0.50E+0 P = CMA-CMI CM = MAX(CM,P) 20 CONTINUE IF (CM.LT.0.10E+0) CM = 2.50E+0 PERZ = 10.0E+0*CM ENDIF C C DETERMINE BONDS C IF (NBND.EQ.0) THEN NMI = NAT-1 IF (NMI.NE.0) THEN DO 40 I = 1, NMI IPI = I+1 DO 30 J = IPI, NAT P = SQRT((C(I,1)-C(J,1))**2+(C(I,2)-C(J,2))**2+(C(I,3) * -C(J,3))**2) C C CHECK TO SEE IF THE ATOMS ARE WITHIN 10 PERCENT OF THE COVALENT C BONDING DISTANCE, IF SO, DRAW THE BOND C COVBND = (COVRAD(IAN(I))+COVRAD(IAN(J)))*1.10E+0 IF (P.LE.COVBND) THEN NBND = NBND+1 IB(NBND) = I IC(NBND) = J ENDIF 30 CONTINUE 40 CONTINUE ENDIF ENDIF C C TRANSFORM C C CALL ROTCOR (C,NAT,THE,GAM,PHI,PERZ,CO) CL = CM*SCALE ADMI = -CL PII = CL*0.20E+0 C C MAKE TITLE C CALL SYMBOL (0.10E+0,0.10E+0,0.300E+0,TITLE,0.00E+0,120) C C THE NEXT LINE MAY NEED MODIFICATION - SHIFT THE ORIGIN TO CENTER T C PLOTS MORE, AS WELL AS MOVE AWAY FROM THE TITLE C CALL PLOT (1.250E+0,1.250E+0,-3) C C DETERMINE LOWEST POINT IN MOLECULE C YMOLMN = 100.0E+0 DO 50 I = 1, NAT YMOLMN = MIN(YMOLMN,C(I,2)) 50 CONTINUE C C MAKE PLOTTING BOX C C TO PUT A BOX AROUND THE PLOT C COMMENT OUT THE CALL PLOT STATMENTS C C CALL PLOT (0.00E+0,0.00E+0,3) C CALL PLOT (10.00E+0,0.00E+0,2) C CALL PLOT (10.00E+0,10.00E+0,2) C CALL PLOT (0.00E+0,10.00E+0,2) C CALL PLOT (0.00E+0,0.00E+0,-2) C C MAKE ATOM SYMBOLS C DO 100 I = 1, NAT C C ******************************************* C C USE THIS SECTION TO PICK OUT THE ATOMS TO PLOT C FOR 2-D PLOTS C C IF (I.EQ.1) GOTO 50 C IF (I.EQ.2) GOTO 50 C GOTO 11 C 50 CONTINUE C C *************************************** C IANI = IAN(I) IF (IANI.LE.18) THEN NP = NCHAR1(IANI) ELSE NP = IATSMN(I) ENDIF X = ((C(I,1)+CL)/CL)*5.00E+0-0.0710E+0 Y = ((C(I,2)+CL)/CL)*5.00E+0-0.0710E+0+0.01250E+0 C C X = ((C(I,1)+CL)/CL)*5.00E+0-0.0710E+0-0.0700E+0 C Y = ((C(I,2)+CL)/CL)*5.00E+0-0.0710E+0-0.0500E+0 C IF (NP.EQ.2) X = X-0.0500E+0 C C DRAW ATOMS 4 TIMES TO DARKEN (FOR BAD PENS, ETC.) C DO 100 J = 1, 4 GO TO (90,60,70,80), J 60 Y = Y+0.0050E+0 X = X+0.0050E+0 GO TO 90 70 X = X+0.0050E+0 Y = Y-0.0050E+0 GO TO 90 80 Y = Y-0.0050E+0 X = X-0.0050E+0 90 IF (IANI.GT.18) THEN CALL SYMBOL (X,Y,0.250E+0,IATSYM(I),0.00E+0,NP) ELSE CALL SYMBOL (X,Y,0.250E+0,ATS(IANI),0.00E+0,NP) ENDIF 100 CONTINUE C C *********************************************** C C UNCOMMENTING THIS LINE WILL SKIP THE PLOTTING OF BONDS C C GOTO 120 C C *********************************************** C C DRAW BONDS C IF (NBND.NE.0) THEN BX(3) = ADMI BY(3) = ADMI BX(4) = PII BY(4) = PII DO 110 I = 1, NBND IBI = IB(I) ICI = IC(I) BX(1) = C(IBI,1) BX(2) = C(ICI,1) BY(1) = C(IBI,2) BY(2) = C(ICI,2) CALL LVGAP (BX(1),BY(1),BX(2),BY(2),PII,IER) C C PLOT BONDS 4 TIMES SLIGHTLY OFFSET, THIS WILL DARKEN THEM IN C ON PEN-PLOTTERS AS WELL AS LASER PRINTERS. C IF (IER.NE.1) THEN CALL DSHLIN (BX,BY,2,NULL,NULL,0) BY(1) = BY(1)+.0050E+0 BY(2) = BY(2)+.0050E+0 CALL DSHLIN (BX,BY,2,NULL,NULL,0) BX(1) = BX(1)+.0050E+0 BX(2) = BX(2)+.0050E+0 CALL DSHLIN (BX,BY,2,NULL,NULL,0) BY(1) = BY(1)-0.0050E+0 BY(2) = BY(2)-0.0050E+0 CALL DSHLIN (BX,BY,2,NULL,NULL,0) ENDIF 110 CONTINUE ENDIF RETURN END C C SUBROUTINE LVGAP (X1,Y1,X2,Y2,PII,IER) DATA R / 0.280E+0 / C C CALCULATE GAPS AROUND ATOM SYMBOLS C DX = X2-X1 DY = Y2-Y1 DXY2 = DX**2+DY**2 DIST = SQRT(DXY2) TH = ATAN(ABS(DY/DX)) XP = R*PII*COS(TH) YP = R*PII*SIN(TH) IF (X1.GT.X2) XP = -XP IF (Y1.GT.Y2) YP = -YP X1 = X1+XP Y1 = Y1+YP X2 = X2-XP Y2 = Y2-YP PDIST = XP**2+YP**2 PDIST = SQRT(PDIST) IF (DIST.LT.PDIST) THEN IER = 1 ELSE IER = 0 ENDIF RETURN END C C SUBROUTINE ROTCOR (A,N,THE,GAM,PHI,PERZ,CO) DIMENSION A(1024,3),CO(3) DATA PI / 3.141592653610E+0 / RPI = PI/180.0E+0 C C TRANSFORM TO CENTERED SYSTEM C DO 20 I = 1, 3 DO 10 J = 1, N A(J,I) = A(J,I)-CO(I) 10 CONTINUE 20 CONTINUE C C XY ROTATION C RT = THE*RPI CT = COS(RT) ST = SIN(RT) DO 30 I = 1, N RT = CT*A(I,1)-ST*A(I,2) A(I,2) = ST*A(I,1)+CT*A(I,2) A(I,1) = RT 30 CONTINUE C C ZX ROTATION C RT = PHI*RPI CT = COS(RT) ST = SIN(RT) DO 40 I = 1, N RT = CT*A(I,3)-ST*A(I,1) A(I,1) = ST*A(I,3)+CT*A(I,1) A(I,3) = RT 40 CONTINUE C C YZ ROTATION C RT = GAM*RPI CT = COS(RT) ST = SIN(RT) DO 50 I = 1, N RT = CT*A(I,2)-ST*A(I,3) A(I,3) = ST*A(I,2)+CT*A(I,3) A(I,2) = RT 50 CONTINUE C C GIVE PERSPECTIVE C DO 60 I = 1, N RT = PERZ-A(I,3) IF (RT.NE.0.0E+0) THEN RT = PERZ/RT A(I,1) = RT*A(I,1) A(I,2) = RT*A(I,2) ENDIF 60 CONTINUE RETURN END C C SUBROUTINE HIDPLT (SUBT,NHL,JLINES) IMPLICIT REAL (A-H,O-Z) C C SUBROUTINE TO COMPUTE HIDDEN LINE ELIMINATION ON CONTOURS C COMMON /DLT/ YINT(75000),SLPINV(75000),PLANES(1024,3),NCHID COMMON /POINTS/ AC(75000,3) COMMON /HLCOM/ NAC,NBEG,ACMIN,ACMAX COMMON /MIN/ YMOLMN COMMON /RPLOT/ N,CO(3),CM,THE,GAM,PHI,X(1024),Y(1024),Z(1024), * IDASH,SCALE,PERZ CHARACTER*(*) SUBT C C LOGICAL PARLEL - DESIGNATES IF TWO PLANES ARE PARALLEL C C 1/88 MODIFIED HIDDEN LINE ALGORITHM FOR SPEED AND CORRECTION C OF THE TREATMENT OF SOME SPECIAL CASES - DAN SEVERANCE, PURDUE C DIMENSION ACMIN(1024,3),ACMAX(1024,3),C(1024,3),AMA(3),AMI(3) DIMENSION XPL(952),YPL(952),ACNMIN(1024,3),ACNMAX(1024,3) INTEGER NAC(1024),NCDASH(1024),PENCOD INTEGER IPOS(1024),NBEG(1024),NCHID(1024),IPARA(1024) LOGICAL HIDDEN,PARLEL,FIRST EQUIVALENCE (C(1,1),X(1)) DATA ILST,IDSK1,IDSK2,IDSK3 / 6,22,23,24 /,PENCOD / 1 / C C READ CURVE DATA FROM DISK FILES, TRANSFORM CURVES AND C DETERMINE EXTREMA FOR EACH CURVE C OPEN (IDSK1,FILE='FOR022',STATUS='OLD') OPEN (IDSK2,FILE='FOR023',STATUS='OLD') OPEN (IDSK3,FILE='FOR024',STATUS='OLD') READ (IDSK2,40) NC,NCURV,NC1,NCURV1,NC2,NCURV2,NC3,NCURV3 C C SET POINTER FOR PARALLEL PLANES, I.E. IF IPARA(I)= C IPARA(J), THEN PLANES I AND J ARE PARALLEL. THIS IS USED C BY THE HIDDEN LINE ROUTINE. C DO 10 I = 1, NCURV1 IPARA(I) = 1 10 CONTINUE DO 20 I = (NCURV1+1), NCURV2 IPARA(I) = 2 20 CONTINUE DO 30 I = (NCURV2+1), NCURV3 IPARA(I) = 3 30 CONTINUE 40 FORMAT (10I6) C C CHECK TO SEE IF ENOUGH ARRAY SPACE IS ALLOCATED FOR THE C CONTOURS.. C IF (NC.GT.75000) THEN WRITE (ILST,50) NC 50 FORMAT (1X,'%%%%%% ARRAY AC MUST BE REDIMENSIONED %%%%%%'/,1X * ,'%%%%%% AC NEEDS',I7,' BUT HAS ONLY 75000 %%%%%%') CALL NEWPEN (1) CALL PLOT (0.0E+0,0.0E+0,999) STOP ENDIF C C CHECK THAT THE DIMENSIONS OF INDIVIDUAL CONTOUR ARRAYS ARE LARGE C ENOUGH.. C IF (NCURV.GT.1024) THEN CALL NEWPEN (1) CALL PLOT (0.0E+0,0.0E+0,999) WRITE (ILST,60) NCURV 60 FORMAT (1X, * '%%%%%% ARRAYS DIMENSIONED AT 1024 NEED TO BE %%%%%%'/,1X, * '%%%%%% AS',I6,' %%%%%%') STOP ENDIF C C READ IN THE INDEXES TO THE CONTOUR FILE, I.E. NAC(I) IS THE C NUMBER OF POINTS IN CONTOUR (I), AND NCDASH(I) IS THE "COLOR" C OF CONTOUR (I).... C READ (IDSK2,100) (NAC(I),NCDASH(I),I=1,NCURV) DO 80 I = 1, NCURV C C PERHAPS THIS SHOULD BE DONE SOME WAY OTHER THAN HITTING EOF C THERE ARE AT MOST NCURV ELEMENTS IN FILE 24 C READ (IDSK3,70,END=90) NCHID(I) 70 FORMAT (I6) 80 CONTINUE READ (IDSK3,70,END=90) NCHID(NCURV+1) WRITE (ILST,*) 'ERROR -- ATTEMPTING TO READ MORE POINTS THAN' WRITE (ILST,*) * 'POSSIBLE FROM FILE 24, PLEASE CHECK FILE INTEGRITY' STOP 90 NCHID2 = I-1 100 FORMAT (20I4) C C PLOTTING PARAMETERS C DSH = 0.120E+0 GAP = 0.100E+0 HITE = 0.250E+0 XSBT = 5.00E+0-15.0E+0*HITE KI = 0 C DO 150 I = 1, NCURV K = NAC(I) IF (K.GT.1024) THEN WRITE (ILST,110) K 110 FORMAT (1X,'%%%%%% ARRAYS X,Y,Z DIMENSIONED AT 1024. %%%%%%' * /,1X,'%%%%%% NEED ',I4,'... ALL 1024 AND 952 %%%%%%'/, * 1X,'%%%%%% ARRAYS NEED TO BE CHANGED %%%%%%') CALL NEWPEN (1) CALL PLOT (0.00E+0,0.00E+0,999) STOP ENDIF C C READ IN A CONTOUR C READ (IDSK1,120) (X(M),Y(M),Z(M),M=1,K) 120 FORMAT (8F10.6) AMA(1) = -100.0E+0 AMA(2) = -100.0E+0 AMA(3) = -100.0E+0 AMI(1) = 100.0E+0 AMI(2) = 100.0E+0 AMI(3) = 100.0E+0 C C ROTATE THE CONTOUR USING THETA,PHI,GAMMA..... C CALL ROTCOR (C,K,THE,GAM,PHI,PERZ,CO) C C DETERMINE THE PLANE EQUATION FOR EACH CONTOUR TO USE IN C THE HIDDEN LINE ROUTINE. C CALL PLANEQ (K,I) NBEG(I) = KI+1 C C DETERMINE X,Y, AND Z MIN AND MAX FOR EACH CONTOUR C DO 130 J = 1, K KI = KI+1 P = C(J,1) AMA(1) = MAX(AMA(1),P) AMI(1) = MIN(AMI(1),P) AC(KI,1) = P P = C(J,2) AMA(2) = MAX(AMA(2),P) AMI(2) = MIN(AMI(2),P) AC(KI,2) = P P = C(J,3) AMA(3) = MAX(AMA(3),P) AMI(3) = MIN(AMI(3),P) AC(KI,3) = P 130 CONTINUE DO 140 II = 1, 3 ACMIN(I,II) = AMI(II) ACMAX(I,II) = AMA(II) 140 CONTINUE 150 CONTINUE C C COMPUTE THE INVERSE OF THE SLOPE AND THE INTERCEPT C FOR EACH LINE SEGMENT FORMING THE CONTOURS. C DO 160 IK = 2, NC DENOM = AC(IK,2)-AC((IK-1),2) IF (ABS(DENOM).GT.0.00010E+0) THEN SLPINV(IK) = (AC(IK,1)-AC((IK-1),1))/DENOM ELSE SLPINV(IK) = 1000.0E+0 ENDIF IF (ABS(SLPINV(IK)).GT.0.00010E+0) THEN YINT(IK) = AC(IK,2)-(AC(IK,1)/SLPINV(IK)) ELSE YINT(IK) = -1000.00E+0 ENDIF 160 CONTINUE C C GENERATE ARRAY OF POINTERS TO THE START OF INDIVIDUAL CONTOURS. C THESE ARE USED TO KEEP TRACK OF ALL CONTOURS CONTAINED IN A C PLANE FOR THE HIDDEN LINE ELIMINATION. C KI = NCURV+1 NBEG(KI) = NBEG(NCURV)+NAC(NCURV) CL = CM*SCALE ADMI = -CL PII = 2.00E+0*CL*0.10E+0 C C HIDDEN LINE ELIMINATION C NCCNT = 1 DO 180 NK = 1, NCHID2 C C FIND TOTAL MIN AND MAX FOR EACH PLANE C ACNMAX(NK,1) = -1000.00E+0 ACNMAX(NK,2) = -1000.00E+0 ACNMAX(NK,3) = -1000.00E+0 ACNMIN(NK,1) = 1000.00E+0 ACNMIN(NK,2) = 1000.00E+0 ACNMIN(NK,3) = 1000.00E+0 C C REDUCE TO ONE PLANE EQUATION PER PLANE... C PLANES(NK,1) = PLANES(NCCNT,1) PLANES(NK,2) = PLANES(NCCNT,2) PLANES(NK,3) = PLANES(NCCNT,3) IPARA(NK) = IPARA(NCHID(NK)) DO 170 J = NCCNT, NCHID(NK) ACNMAX(NK,1) = MAX(ACNMAX(NK,1),ACMAX(J,1)) ACNMAX(NK,2) = MAX(ACNMAX(NK,2),ACMAX(J,2)) ACNMAX(NK,3) = MAX(ACNMAX(NK,3),ACMAX(J,3)) ACNMIN(NK,1) = MIN(ACNMIN(NK,1),ACMIN(J,1)) ACNMIN(NK,2) = MIN(ACNMIN(NK,2),ACMIN(J,2)) ACNMIN(NK,3) = MIN(ACNMIN(NK,3),ACMIN(J,3)) 170 CONTINUE NCCNT = NCHID(NK)+1 180 CONTINUE C MCCNT = 1 DO 350 MK = 1, NCHID2 C C COMPUTE THE Z VALUE OF THE PLANE AT X=1, Y=1 TO COMPARE PARALLEL C PLANES ( IF THIS ONE'S Z VALUE IS IN FRONT OF THE Z VALUE OF A C PLANE PARALLEL TO IT, ALL POINTS ON THE PLANE ARE IN FRONT OF C THE PARLEL PLANE, AND THUS CANNOT BE HIDDEN BY ANY SURFACE ON C THE PLANE. PLANES(MK,1)=(A/C), PLANES(MK,2)=(B/C), C PLANES(MK,3)=(D/C), Z= - ( (A/C)X + (B/C)Y + (D/C) ), SO FOR X=1 C AND Y=1, Z = - ( PLANES(MK,1)+PLANES(MK,2)+PLANES(MK,3) ) C ZTST = PLANES(MK,1)+PLANES(MK,2)+PLANES(MK,3) IPARA1 = IPARA(MK) DO 340 I = MCCNT, NCHID(MK) K = NAC(I) NUM = 0 IS = 1 IF = K IF (NHL.NE.1) THEN NCCNT = 1 DO 200 NK = 1, NCHID2 ZNK = PLANES(NK,1)+PLANES(NK,2)+PLANES(NK,3) PARLEL = IPARA1.EQ.IPARA(NK) IF ((ZTST.LT.ZNK).OR.(.NOT.PARLEL)) THEN FIRST = .TRUE. DO 190 J = NCCNT, NCHID(NK) IF (FIRST) THEN IF (ACMIN(I,1).LT.ACMAX(J,1)) THEN IF (ACMAX(I,1).GT.ACMIN(J,1)) THEN IF (ACMIN(I,2).LT.ACMAX(J,2)) THEN IF (ACMAX(I,2).GT.ACMIN(J,2)) THEN IF (ACMIN(I,3).LT.ACMAX(J,3)) THEN NUM = NUM+1 IPOS(NUM) = NK FIRST = .FALSE. ENDIF ENDIF ENDIF ENDIF ENDIF ENDIF 190 CONTINUE ENDIF NCCNT = NCHID(NK)+1 200 CONTINUE ENDIF C C IPOS CONTAINS THE POSSIBLE COVERING CURVES, PLACE CURVE IN C C NBE = NBEG(I)-1 DO 210 L = 1, K NBE = NBE+1 C(L,1) = AC(NBE,1) C(L,2) = AC(NBE,2) C(L,3) = AC(NBE,3) 210 CONTINUE C C TAKE EACH POINT AND SEE IF IT IS BEHIND A IPOS C IF IT IS, MARK IT WITH 1000 FOR ITS X. C IF (NUM.NE.0) THEN DO 240 L = 1, K X1 = X(L) Y1 = Y(L) Z1 = Z(L) DO 220 J = 1, NUM NCRV = IPOS(J) Z2 = X1*PLANES(NCRV,1)+Y1*PLANES(NCRV,2)+ * PLANES(NCRV,3) IF (Z2.GT.(Z1+0.20)) THEN C C IF (Z1.LT.ACNMAX(NCRV,3)) THEN C IF (Z2.LT.ACNMAX(NCRV,3)) THEN IF (Z2.GT.ACNMIN(NCRV,3)) THEN IF (X1.LT.ACNMAX(NCRV,1)) THEN IF (X1.GT.ACNMIN(NCRV,1)) THEN IF (Y1.LT.ACNMAX(NCRV,2)) THEN IF (Y1.GT.ACNMIN(NCRV,2)) THEN IF (HIDDEN(X1,Y1,NCRV)) GO TO 230 ENDIF ENDIF ENDIF ENDIF ENDIF ENDIF C C ENDIF C ENDIF 220 CONTINUE GO TO 240 230 X(L) = 1000.0E+0 240 CONTINUE ENDIF C C PLOT WHAT S LEFT OF THE CURVE C NSEC = 1 IF (NCDASH(I).EQ.0) LINCOL = 2 IF (NCDASH(I).NE.0) LINCOL = 3 IF (JLINES.GE.1.AND.LINCOL.NE.PENCOD) CALL NEWPEN (LINCOL) PENCOD = LINCOL NSEC = 1 IF (NCDASH(I).NE.0) NSEC = 0 IF (JLINES.EQ.1) NSEC = 0 IF (NUM.EQ.0) GO TO 310 K1 = 1 250 DO 260 J = K1, K IF (X(J).NE.1000.0E+0) GO TO 270 260 CONTINUE GO TO 340 270 IS = J DO 280 J = IS, K IF (X(J).EQ.1000.0E+0) GO TO 290 280 CONTINUE IF = K GO TO 300 290 IF = J-1 300 IF (IF.EQ.IS) GO TO 330 310 JJ = 0 DO 320 II = IS, IF JJ = JJ+1 XPL(JJ) = X(II) YPL(JJ) = Y(II) 320 CONTINUE XPL(JJ+1) = ADMI YPL(JJ+1) = ADMI XPL(JJ+2) = PII YPL(JJ+2) = PII CALL DSHLIN (XPL,YPL,JJ,DSH,GAP,NSEC) 330 IF (IF.NE.K) THEN K1 = IF+1 GO TO 250 ENDIF 340 CONTINUE MCCNT = NCHID(MK)+1 350 CONTINUE CALL NEWPEN (1) C C PLOT SUBTITLE C THE SUBTITLE IS PLACED SLIGHTLY BELOW THE MOLECULE C AND IS GENERALLY USED TO REPORT THE ENERGY OF THE C ORBITAL C YMIN = 100.0E+0 DO 360 I = 1, NCURV YMIN = MIN(YMIN,ACMIN(I,2)) 360 CONTINUE YMIN = MIN(YMIN,YMOLMN) YMIN = 5.00E+0*(YMIN+CL)/CL-0.800E+0 CALL SYMBOL (XSBT,YMIN,HITE,SUBT,0.00E+0,40) RETURN END C C LOGICAL FUNCTION HIDDEN (X,Y,NCRV) C C GIVEN THE POINT X,Y WHICH IS BEHIND THE PLANE CONTAINING C THE SURFACE NCRV, DETERMINE IF THE AREA OF THE SURFACE C COVERS THE POINT. C C THIS ROUTINE IS THE GUTS OF THE HIDDEN LINE ELIMINATION. C THE ALGORITHM IT USES WAS DEVELOPED BY C W.L. JORGENSEN AND D.E. BARTH IN APRIL 1972. C C MAJOR MODIFICATIONS TO THE ALGORITHM WERE MADE IN OCTOBER 1987 C BY DANIEL L. SEVERANCE TO INCREASE THE EFFICIENCY. THE HEIGHT C OF THE CURVE AT POINT X,Y IS CALCULATED BY SUBSTITUTING INTO C THE PLANE EQN. FOR THE CURVE. THIS IS DONE OUTSIDE OF THE C ROUTINE TO SCREEN OUT CURVES WHICH WILL NOT BE HIDING THE POINT C ANYWAY. THIS ROUTINE THEN DETERMINES FOR THOSE CURVES WITH C LARGER 'Z' VALUES, IF THE POINT IS ACTUALLY BEHIND A FILLED C PORTION OF THE CURVE, OR IN AN OPEN AREA. THIS HAS NOW BEEN C MADE INTO A LOGICAL FUNCTION (TRUE/FALSE RETURNED). C D.L. SEVERANCE C C JAN. '88 MODIFIED TO IMPROVE PERFORMANCE, CLEANED UP "CLIPPING" C TO USE XOR (.NEQV.) LOGIC RATHER THAN IF STATEMENTS. D.L.S. C LOGICAL LT1,LT2,EVEN COMMON /DLT/ YINT(75000),SLPINV(75000),PLANES(1024,3),NCHID(1024) COMMON /POINTS/ AC(75000,3) COMMON /HLCOM/ NAC(1024),NNBEG(1024),ACMIN(1024,3),ACMAX(1024,3) IF (NCRV.EQ.1) THEN NBEG = 1 ELSE NBEG = NCHID(NCRV-1)+1 ENDIF EVEN = .TRUE. HIDDEN = .FALSE. DO 20 IK = NBEG, NCHID(NCRV) NBEG2 = NNBEG(IK) K = NAC(IK)-1 C C THIS SECTION DETERMINES WHICH LINES MAKE WINDOWS ABOUT THE POIN C USE A SINGLE LOGICAL ASSIGNMENT FOR EFFICIENCY. C LT2 = (Y.LT.AC(NBEG2,2)) DO 10 N2 = NBEG2+1, NBEG2+K LT1 = LT2 LT2 = (Y.LT.AC(N2,2)) IF (LT1.NEQV.LT2) THEN EVEN = .NOT.EVEN XP = (Y-YINT(N2))*SLPINV(N2) C C XP IS THE X VALUE FOR THE C INTERSECTION POINT OF A LINE AT Y=Y WITH THE C CURVE IN QUESTION, NCRV. (THERE MAY BE ANY C NUMBER OF INTERSECTIONS.) C IF (X.GE.XP) THEN HIDDEN = .NOT.HIDDEN ENDIF ENDIF 10 CONTINUE 20 CONTINUE HIDDEN = HIDDEN.AND.EVEN RETURN END C C SUBROUTINE PLANEQ (IDIM,NCTR) COMMON /DLT/ YINT(75000),SLPINV(75000),PLANES(1024,3),NCH(1024) COMMON /RPLOT/ N,CO(3),CM,THE,GAM,PHI,X(1024),Y(1024),Z(1024), * IDASH,SCALE,PERZ DIMENSION V(2,3),D(4) C C GENERATE THE COEFFICIENTS OF THE PLANE EQUATION GIVEN 3 INPUT C POINTS KNOWN TO LIE IN THE PLANE C I2 = IDIM/3 I3 = I2*2+1 I2 = I2+1 V(1,1) = X(I2)-X(1) V(1,2) = Y(I2)-Y(1) V(1,3) = Z(I2)-Z(1) V(2,1) = X(I3)-X(1) V(2,2) = Y(I3)-Y(1) V(2,3) = Z(I3)-Z(1) C C NOW WE HAVE THE TWO LINES CALC. CROSS PRODUCT C D(1) = V(1,2)*V(2,3)-V(1,3)*V(2,2) D(2) = V(1,3)*V(2,1)-V(1,1)*V(2,3) D(3) = V(1,1)*V(2,2)-V(1,2)*V(2,1) C C NOW THE VALUE OF D BY SUBSTITUTING VALUES FROM C(1) C D(4) = -(D(1)*X(1)+D(2)*Y(1)+D(3)*Z(1)) C C TO FIND Z GIVEN X AND Y, THE EQUATION REARRANGES TO C Z = ( D(1)*X + D(2)*Y + D(4) ) / D(3) C REWRITING AS Z = (D(1)/D(3))*X + (D(2)/D(3))*Y + D(4)/D(3) C THIS NEXT CODE CALCS THE MULTIPLIERS ABOVE IN PARENTHESIS C D3INV = -1.00E+0/D(3) PLANES(NCTR,1) = D(1)*D3INV PLANES(NCTR,2) = D(2)*D3INV PLANES(NCTR,3) = D(4)*D3INV RETURN END