In this section, we'll learn how to do numerical calculations. As an example, we will plot out the cubic form of the van der Waals equation of state as a function of volume, V. In fact, the code shown below is the first step in solving problem 16-22 in Physical Chemistry: A Molecular Approach by McQuarrie and Simon.
program main
implicit none
integer :: i
real :: V,f,fp
real :: c3,c2,c1,c0 ! coefficients of the cubic eqn
real, parameter :: P=35.00, T= 142.69 ! P in atm, T in K
real, parameter :: R = 0.082058, & ! dm^3.atm/mol.K (Tbl 16.2)
a=1.3307,b=0.031830, & ! Ar params, Tbl 16.3 (in dm and atm)
mm= 39.948 ! molar mass of Ar in g/mol
! use cubic form of the van der Waals eqn (eq 16.10)
c3 = 1.0
c2 = -(b+R*T/P)
c1 = a/P
c0 = -a*b/P
write (*,*) '#conditions = ',P,T
write (*,*) '#vdW eqn: (',c3,') v^3 + (',c2,') v^2 '
write (*,*) '# + (',c1,') v + (',c0,')'
! NOTE: everything above here is exactly the same in 16-22a.f90 and 16-22b.f90
write (*,*) '# v','f(v)','0'
do v=0.0,+0.3,0.01
f=c3*(V**3) + c2*(V**2) + c1*V + c0
write (*,*) v, f, 0.0
end do
stop
end program main
Let's take it line by line:
program main
implicit none
should be familiar from the previous lesson. If not, please go back and review. You'll notice there is a blank space here…we can insert these wherever we want, with the goal of inmproving the readibility of the source code.
defines a variable that stores an "integer" (i.e., a number without a decimal point) named "i". We can use this to store values during subsequent calculations. There are several types of variables in Fortran:
integer,
real,
complex,
logical,
character. The meanings of these should be fairly obvious (eg.,
real stores numbers with decimal points,
complex stores complex numbers,
logical stores true/false values,
character stores text, etc.) Variable names (e.g., "i") must start with a letter, and can contain letters, numbers, and underscores ("_"). Each variable should have a unique name that is declared only once. F90 is
case-insensitive, in other words, it doesn't matter whether you use uppercase or lowercase letters when referring to variables throughout your program. That being said, I like to try to be consistent.
just like above defines real variables to store volume, V, and the values of the cubic function, "f", and its derivative "fp" (think "f-prime"). The comma tells the compiler that we want to declare another variable next. We can do this as many times as we like.
real :: c3,c2,c1,c0 ! coefficients of the cubic eqn
again defines 4 real variables to store the coefficients of the cubic equation. But then, you already knew that from reading the
comments for the human; everything after the "!" is ignored by the computer. It is always a good idea to remind yourself what variables store what, unless it is obvious. You can use comments anywhere you like to clarify the meaning of the program
real, parameter :: P=35.00, T= 142.69 ! P in atm, T in K
Here we define variables as
parameters, that is as fixed constants whose values don't change during the course of our calculation. E.g., the values of the pressure (P) and temperature (T) might be used throughout the program. If we wanted to change our program for another P (or T) all we have to do is change them once up here, rather than having to go through the entire program. Also, it is much easier to understand "P" in an equation than some arbitrary number. As before, we can use commas to define as many parameters as we wish.
real, parameter :: R = 0.082058, & ! dm^3.atm/mol.K (Tbl 16.2)
a=1.3307,b=0.031830, & ! Ar params, Tbl 16.3 (in dm and atm)
mm= 39.948 ! molar mass of Ar in g/mol
illustrates a few things. Here we have defined "R" to hold the value of the
gas constant, rather than having to type this in multiple times in our program. As before, the comma following the definition of R tells the compiler that we want to define another "real" variable/parameter. Notice the "
&". This tells the compiler to continue reading from the next line of the program. In this case, the second line has the "a" and "b" parameters that are used in the van der Waals equation. We could just let the line run on and on, but that looks messy on the screen. Finally, it lets us comment each of our variables, to remind ourselves what units they are in and what table we obtained them from.
! use cubic form of the van der Waals eqn (eq 16.10)
c3 = 1.0
c2 = -(b+R*T/P)
c1 = a/P
c0 = -a*b/P
After reminding ourselves with a comment that we are calculating the coefficients for the cubic form of the vdW equation, we perform the calculation. In this case, we give "c3" the value of "1.0". In the other cases, we perform some arithmetic. E.g., we set "c1" equal to "a
divided by P" Remember that multiplication/division have a higher precedence than addition/subtraction. Add parentheses as necessary to make things clear for you and for the compiler. One thing to be careful of is division of two integers; this results in an integer value (e.g., "3.0/2.0 = 1.5, but 3/2 = 1").
write (*,*) '#conditions = ',P,T
write (*,*) '#vdW eqn: (',c3,') v^3 + (',c2,') v^2 '
write (*,*) '# + (',c1,') v + (',c0,')'
! NOTE: everything above here is exactly the same in 16-22a.f90 and 16-22b.f90
write (*,*) '# v','f(v)','0'
We've seen "write" statements before. The difference here is that we mix strings (e.g. "# conditions = ") with variable values (e.g., "P,T"). Just keep them separated with commas and it will be fine. Essentially all this is doing is either writing out in text whatever is in quotes (#conditions =) or writing the current values for variables/parameters (like P or T….if you put quotes around P or T, you'll just write out the letters P or T instead of the number values). Listing each of this in the same line separated by commas will write out each component in the same row. Having multiple write statements, as above, will result in the components being written underneath those written from above.
do v=0.0,+0.3,0.01
......
end do
Here's something new. We can tell the compiler to
loop over a variable, to perform a block of code multiple times. Here, we tell the compiler to make "V=0.0", and increase it to "+0.3" in steps of "0.01". As long as V is less than +0.3, it will perform the calculation in …. (which can be a single line or many lines). Finally, when it reaches "+0.3", the program will continue onwards from the "end do" statement. Next, let's take a closer look at the "…..":
f=c3*(V**3) + c2*(V**2) + c1*V + c0
write (*,*) v, f, 0.0
In the first line we are simply calculating "f" which is the cubic equation form of the vdW equation, using the constants we computed above. "
V**3" means "V-to-the-third-power" (i.e., "
**" is the
exponentiation operator). Finally, we write out the values of v, f, and the constant zero.
You already know what these do.
To compile and run this program try:
pgf90 -o 16-22a.r 16-22a.f90
./16-22a.r > vdw.dat
xmgrace -nxy vdw.dat
The first line compiles your program that you've written (ending in .f90) into a program that you can run (ending in .r)
The second line runs the program, and creates and saves another file (ending in .dat) that contains all the data points calculated from running your compiled program (the one that ended in .r)
The third line plots the data from the specified file. If your data exists in multiple columns, the first column will be the x-axis. Any additional columns will be plotted as individual series. this format is denoted by the "-nxy"
If you've got your X11 connection set up properly (see [http://h205.wikidot.com] if you don't), you should see something like the following. We will use this as an initial guess for our Newton-Raphson solver, as shown in the next section.
