Fritz
S. Allen
and Edward A. Walters
University of New Mexico Albuquerque, 87131
I
I
TWO Computer Problems Dealing with Graphic Display in Chemistry
For the past several years we have taught a course introducing upper division chemistry students to the Fortran language and the application of numerical techniques to chemical prohlems. One phase of this course has centered around the use of graphic display and we have developed two exercises which the students have enjoyed solving. A number of problems have appeared (1-3) which utilize the peripheral plotting facilities of the computer. We have avoided this because of the costs involved and also because plotter operations often differ between computer centers. By employing the line printer as our graphic display device, we gain the added advantage that the routines developed can be used on any machine with a Fortran compiler and line printer. In the course, the line printer displays have a dual role. We introduce them after the class has been exposed to the various types of formats as an exercise in extensive and critical use of formats. At the completion of a problem the students have usually mastered the concepts of formatting. The second use of the displays is to demonstrate that even with the calculating ability of a computer "a picture is worth a thousand words." We try to emphasize that plots are easy to make, useful, and much easier to interpret than tables of data. One difficulty often encountered in a programming class with an exercise as complex as the generation of a line printer plot, is that the problem needs to be stated in such a way that i t seems solvable and yet the solution is not merely a simple Iine by line transliteration of the problem statement into Fortran. We have circumvented this difficulty by giving the problem as a schematic flow diagram which specifies the appropriate sequence of logical operations while leaving the mechanics of each step to the student. This technique seems to work quite well. In a class with a typical spectrum of ability, all of the students understand the problem, while only a few could have drawn the flow diagram properly. Slower students can be assisted over rough sections. If this is done by the instructor i t hecomes clear very quickly how well students are grasping the intricacies of programming. In any case, the problem remains challenging to the best students and is completed by nearly everyone. Two different types of problems have been used. The first is a procedure to produce a contour plot of a hydrogen-like atomic orbital in a plane normal to the X, Y, or Z axis. The second is a vector display routine which has wide application in a variety of chemical prohlems, and may consequently serve as a convenient subroutine for later use. The flow diagram which is handed out to the class for the contour plot of the hydrogen-like orbitals is shown in Figure 1. In this problem, the students first determine the symbols which will comprise the plot. The symbols are used to denote different intensities on the contour plot. The students commonly divide the range of the plot into 20 increments and indicate each increment with a specific symhol. By beginning with and then alternating the symbol "blank" in the symbol array the final plot is easy to interpret. The second, fourth, and subsequent even numbered members of the array are usually letters or numbers in increasing order. The region of space which is to be the domain of the plot is the cube contained between +25 on the X, Y, and 596
/
Journal of Chemical Education
Figure 1. Flow Diagram for Contour Plots of Atomic Orbitals. Students who elect to produce contour plots of the hydrogen-like orbitals are given the above flow chart as a statement of the problem
Z axes. The line printer is basically a Cartesian coordinate device. We have used this coordinate system to make the generation of the plot as simple as possible. We convert to polar coordinates only to compute the value of the wave function. For the time being it is convenient to assume that all plots are made for a plane Z = constant, i.e., a plane parallel to the XY plane. The wave function of the chosen orhital will be scaled such that the orbital nearly fills the cube. The printing spaces on a typical line printer are higher than they are wide; as a result using equal increments in the X and Y direction results in a rectangular plot which gives a somewhat distorted view. To prevent this we divide the Y direction into 51 printer lines separated by one unit and each line which is the X direction into 81 elements separated by 0.625 units. The resulting 51 X 81 character array fits on one page of standard computer paper and is nearly square; so a good representation of an orbital is produced. The plot is generated by printing an 81 element array P with an A format 51 times. The next operation specified by the flow diagram is to read the value of Z which determines the location of the plot plane. For purposes of termination Z is compared to 50 and if i t is less a plot will be generated. The first element of the plot is a title specifying the nature of the orbital and any other information. The title must be formatted such that the plot and title start a t the top of a new page. After printing the title the next eight instmctions set u p the loop structure which counts the lines and
number of elements in each line to produce a plot. Most of the students formulate the loops with DO statements instead of the IF statements implied in the flow chart. Once the counting apparatus is set up, the various values of X , Y, and Z are used to evaluate the appropriate wave function. Tabulations of the wave functions are available in several sources (4-6). These functions are usually presented in spherical polar cwrdinates. Hence a coordinate change is necessary. Exceptional students have written their own conversion subprogramsbut we have usually given the class a subprogram .for this purpose. The subroutine we provide, which also computes the probability and weiehts i t bv the volume element.. is eiven in the ADpendix. > t a t e m i t 60 in this subroutine is where the scaling occurs. Multi~licationbv a FACTOR neater than one increases R and shrinks thk plot. ~ecre'sing R enlarges the plot. From the calculated value of the probability, the student next must determine in which interval between the minimum and maximum the current value lies and then store the symbol for that interval in the element of P which corresponds to the current value of X. The 81 elements of P represent all of the possible values of X a t a given value of Y and Z. Each location in P(X) denotes a point on the plane Z = constant. At each of these points the probability density is calculated and the result is compared to the interval from 0 to PSIMAX. If we assume that we have chosen 20 different symbols in the symbol array, each symbol then corresponds to a 5% increment. To choose the appropriate symbol for the current value of the probability density we determine which 5% increment from O-PSIMAX contains the value, say the Kth and put the Kth element of the symbol array into P(X). This process is repeated for all 81 values of X . When the calculation of P is completed, the vector P is printed as a single line of type. Subsequently Y is incremented by one and the process begins anew. When Y has been varied from -25 to +25 the plot is completed and control transfers to statement 1. Values for FACTOR and PSIMAX are given in the table for several orbitals. Thus far we have discussed plots made in a plane Z = constant. It is no more effort to produce plots in planes normal to the X or Y axis. All that is required is permutation of X and Z or of Y and Z in the call to subroutine CALC. Thus, one clear advantage of using a subroutine to make the coordinate change is that X , Y, and Z can be permuted a t will. The students are generally very intrigued by this problem and many produced plots as shown in Figure 2. Many of them write programs which are much more elegant than the basic problem. On their own initiative students have expanded the plot from an 81 x 51 array of characters to 130 X 82. This character array is nearly the maximum size for the line printer and provides for increased resolution. T o implement this, the students must overide the line printer spacing and provide the required carriage control. Other students have written routines which when supplied with the quantum numbers generate the wave function, while some have taken linear combinations of atomic orbitals to make hybrid orbital plots. In any case the students who have done this problem have been very gratified to find that the shapes of the hydrogen-like orbitals do indeed come from the wave functions of quantum mechanics and confirm the diagrams in texts whose origins are reasonably nebulous. A second problem we have developed involving graphic nresentation is a subroutine for the disolav of a curve. ?he routine allows for simultaneous display i f several sets of data of identical but arbitraw leneth. The flow diaeram we provide for the students is &ven:n Figure 3. We enter the subroutine and provide the following information in the data list or in common: the data to be displayed, X (multiple data sets are stored sequentially in X);the total number of elements in X , NEL; a minimum value for the
Problem 1 Program Parameters.
Orbital
FACTOR
PSIMAX
Is
0.3
0 043
1.3 1.4
0.024 0.031
-'~
3~ 3d
Appropriate values for FACTOR and PSIMAX.
Figure 2. A Computer Plot of the 4 d ( z 2 ) Orbital. The characters in the above plot are scaled such that the asterisks represent the maximum electron density. The plot has bean made in the YZ plane at the nucleus With the Z axis horizontal.
elements of X , XMIN; a similar maximum, XMAX; and the number of vectors to be displayed, NPL. This subroutine generates a plot by repeatedly printing a vector A which contains 101 elements; thus the abscissa runs in the direction the paper comes out of the printer. This technique obviates the necessity of scaling both the X and Y cwrdinates. The first element is always a which generates the abscissa. After entry, i t is necessary to create an array, A, an array, INDEX, and a symbol array. The latter should be filled with appropriate characters so that different graphs will have different symbols. Next we must compute the number of points in one graph, NPTS, and then use XMAX and XMIN to determine the step size on the ordinate. We check to make certain that the plot interval is not zero nor negative before proceeding. If it is, an error message must be generated before return. Assuming the interval is positive a heading should be printed which should include the number of points and a scale ranging from XMIN to XMAX. The heading should be printed a t the top of a new page. Next A should be filled with "+" and printed with a lOlAl format to generate the ordinate. At this point A(2) through A(101) should be filled with " "or "blank" to clear the field of the plot. The l's, J's and K's in the flow chart form the indices of the loops which ensure that A is printed NPTS times and contains an element for each of the NPL plots every time. As in the previous problem the students most often formulate the loop structure with DO statements. The next statement of consequence generates a nearly square grid of "+"s in the field of the plot. We have chosen to insert a "+" every six rows and every ten columns. This divides the ordinate into tenths and the grid produced looks approximately square. The heart of the plotting routine is contained in the next series of instructions. First, the Ith element of X is examined to determine the proper position on the ordinate for plotting. By comparison of X(1) with the interval XMAX-XMIN, a number L is generated such that A(L) is then set equal to SYMBOL (J). This process is repeated
"+"
Volume 57, Number 9, September 1974 / 597
-
mwmc nd. r o r n r r l ~ ~ m
ms
IRL/"?L
-
3. Flow Diaaram for a Vector Display Subroutine. Students who elect to produce a vector display subroutine are given the above flow chart as a Statement of the problem.
for all NPL vectors to be plotted. The L determined in each case is stored in the array INDEX for future reference. When A has been completed it is printed and then it is necessary to set all of the symbols back to "blank." The array INDEX contains the locations of all of the symbols and it is referenced to reinitialize A. A check is performed to see if a grid was produced for this particular value of I and if so these locations too are filled with "blank." At this point, A is ready to be used to print the next set of points. The subroutine returns when A has been printed NPTS times. After writing this program many of the students modify it to give more esthetically satisfying plots. They commonly find XMAX and XMIN internally, print index values on the abscissa, and generally introduce refinements. The output of a typical student program is given in Figure 4. Toward the end of the semester each student is assigned a special project. These projects are mostly determined by the students themselves, perhaps after some discussion with us. Many of the projects are centered around the graphic display problems. For example, we have had students use the vector display routines to present simultaneously the concentrations of the species involved in multiple equilibria, to show plots of frictional coefficients for molecules of differing shapes (7). and to plot titration curves with the first and second derivatives. By using the contour plots of the hydrogen-like orbitals, students have produced comprehensive pictures for all of the orbitals up to n = 4 and by stacking transparencies made from plots in parallel planes a very good three dimensional view can be obtained. The extent to which the graphic display IOUtines are used in special projects indicates to us that they are very successful and highly motivating problems. Acknowledgment
We would like to thank Mr. Wayne Trott and Mr. Mark Martin for allowing us to reproduce their solution to the above problems. 598
/
Journal of Chemical Education
Figure 4. Multiple Data Set Display. Typical student output tor simultaneous display of several vectors. The curves shown are Sin(A). CoslA) 1). Note that in this plot the student has shifted the axis to and Sin(A the center of the field to better accommodate the data to be plotted and has also dispensed with the grid points to avoid a cluttered appearence.
+
Appendix
SUBROUTINE CALC (X,Y,Z,PSI) FACTOR = 1.4 PI= 3.14159 PI2 = PI/2. R = SQRT(Z**2 + X**2 Y**2)
+
-. ..
~~
20 THETA = ARCOS(Z/R) 30 IF(ABS(X1-0.01)40,40,50 40 PHI = PI-SIGN(PI2,Y) r.n m fin 50 PHI= ATAN(Y/X) + PD-SIGN(PI2,X) 60 R=ReFACTOR PSI = (1/(3l.*SQRT(2*PI)))*R**2*EXP(-R/3)*(SW(THETA))**2*COS(2*PHI) -v
. V " "
PSI = (R*PSI)**Z*SIN(THETA) RETURN END Literature Cited (I) w e d m a n , H. G. J,., Ch~l,~i".R. 0 ,and Feuerhacker, D. G..J . CHEM. EDUC.. 41.354(19641.
(21 Onyzlo, E. A., J. CHEM. EDUC., 42,150 (1965). (31 Bordsss, W.T., and Linnett. J. W . J . CHEM. EDUC.. 47,672 (19701. (41 Eving, H., Walter. J . . and Kimhlo, G. E., "Quantum Chemistry," John Wlley & Sons, hc.,Now York. 1944. (51 Psulinp, L., and Wilsn, E. B.. "lntmdudlon to Quantum Mechanics," McGmwHlllElmkCo.. NewYork. 1335. (61 Leuine, I. N.."Quantum Chemirfry,"Val. 1.AllynandBaconlne., Baton. 1910. (71 Van Holde, K. E.. "Physical Bioehemintry," Prenfiee~Hall,Hc.. Englewmd Cliffs. NcuJerpey, 1971.
A-
Volume 5 1 , Number 9, September 1974
/
599
