Lab 05: File I/O and Plotting#

Examples (by instructor)#

Example 1: Ideal Gas Law Calculation with File I/O#

Define a function to compute the pressure of an ideal gas using PV = nRT. The function takes n (mol), V (L), T (K), and R as arguments and returns pressure in atm. Loop over several temperatures from a list. Print the results on the screen and save them into the file “ideal_gas_results.dat”. Use R = 0.082057 L·atm/(K·mol), n = 1.0 mol, and V = 10.0 L.

# Constants
R = 0.082057   # L·atm/(K·mol)
n = 1.0        # mol
V = 10.0       # L

# Temperatures to test (K)
temperatures = [200, 250, 300, 350, 400, 450, 500]

Example 2: Plotting and Saving a Figure#

Using the data from Example 1, create a plot of Pressure vs. Temperature. Add axis labels with units, a title, and a grid. Save the plot as a PNG file with dpi=300.

import matplotlib.pyplot as plt

# Recalculate pressures for plotting
R = 0.082057
n = 1.0
V = 10.0
temperatures = [200, 250, 300, 350, 400, 450, 500]

Practice Problems (by students)#

Problem 1: Van der Waals Equation with File Input#

This part was done in the previous lab. Using the following van der Waals equation of states, define the Python function returning the pressure of a gas (p), given its volume (V), amount (n) and temperature (T). Use 1.0 mol of CO₂ at different temperatures from the file (“temperatures.dat”) confined to a volume of 0.5 L. Use R = 0.082057 [L atm K⁻¹ mol⁻¹]. Take the van der Waals parameters a = 3.640 L² bar/mol² and b = 0.04267 L/mol.

Van der Waals equation:

\[ p = \frac{nRT}{V - nb} - \frac{n^2a}{V^2} \]

First, create the file “temperatures.dat” containing the following:

# Temperature, K
50
70
120
150
200
250
# Create the temperatures.dat file
with open('temperatures.dat', 'w') as f:
    f.write("# Temperature, K\n")
    f.write("50\n")
    f.write("70\n")
    f.write("120\n")
    f.write("150\n")
    f.write("200\n")
    f.write("250\n")

print("File 'temperatures.dat' created successfully")
File 'temperatures.dat' created successfully
# Constants
R = 0.082057     # L·atm/(K·mol)
a = 3.640        # L²·bar/mol² → convert to L²·atm/mol²
a = a / 1.01325  # 1 bar = 0.986923 atm, so a (L²·atm/mol²)
b = 0.04267      # L/mol
n = 1.0          # mol
V = 0.5          # L

Problem 2: P-V Diagrams for Van der Waals Gas#

Your goal is to use the above van der Waals equation of state Python function to graph p-V diagrams for the temperatures defined in the file “temperatures.dat”. You will need to create one plot with P-V diagrams for all temperatures, as well as one output datafile containing p-V data for each temperature. Proceed as follows.

a) In addition to the math library (import math), import the numerical Python library called NumPy using the shortcut name np:

import numpy as np

Also, import Matplotlib visualization library using the shortcut name plt:

import matplotlib.pyplot as plt

b) Read the temperature data from the “temperatures.dat” file into a list as follows:

T = np.loadtxt('temperatures.dat', comments='#', dtype='float')

Here, the “comments” keyword indicates to omit lines starting with “#” sign (i.e., treat them as comments), and assign the datatype of temperature values to be floats. Check that you have correctly read the values: print the whole list; print any element from the list.

c) To provide 1000 equidistant volume points in the range from 0.05 to 0.3, use the following command:

volume = np.linspace(0.05, 0.3, 1000)

You will need to loop over this space of volumes using for-loop and save the results into one plot. In this plot each line will correspond to a different temperature. Thus, you will have 6 pressure lists (P1, P2, etc.) for each temperature. Start with empty lists such as P1 = [] and then populate them using P1.append method. Create the p-V plot for all six temperatures from the file “temperatures.dat”. Add x and y labels. Add the title “The van der Waals equation of state for CO2”. Save the plot into a file using the following command:

savefig('vdW_PV_diagram.png')

Use two different settings for figure quality: dpi=300 and dpi=50. Save two files.

d) Save the p-V data into six files each one for a given temperature. E.g, the file “vdW_PV_50K.dat” would correspond to 50 K. The file should begin with the following two lines:

# T = 50 K
# P, V

Use the format keyword to create such filenames.

Hint 1: It is recommended to implement two for-loops here, one for the temperature range and one for volume.

Hint 2: It is better to create a list of lists for pressures (i.e., creating a two-dimensional matrix). This way you can access elements with two indices as P[i][j].

# Problem 2: P-V Diagrams for Van der Waals Gas
# Step a) Import libraries
import numpy as np
import matplotlib.pyplot as plt

# Constants (same as Problem 1)
R = 0.082057
a = 3.640 / 1.01325  # L²·bar/mol² → L²·atm/mol²
b = 0.04267
n = 1.0


# Step b) Read temperature data from file



# Step c) Create volume array
volume = np.linspace(0.05, 0.3, 1000)


# plt.savefig('vdW_PV_diagram_low.png', dpi=50)

Problem 3: Reading and Plotting P-V Data#

Now, read the P-V data you just created from a file corresponding to 150 K temperature and plot P vs V. Save the plot into a png file. You can read the data in one go as follows:

np.loadtxt('vdW_PV_150.0K.dat', comments='#', dtype='float', unpack=True)
# Problem 3: Reading and Plotting P-V Data for 150 K

Problem 4: Lennard-Jones Potential for Argon#

Implement the function to compute the Lennard-Jones interatomic potential for argon:

\[ V(r) = 4\varepsilon\left[\left(\frac{\sigma}{r}\right)^{12} - \left(\frac{\sigma}{r}\right)^6\right] \]

where r is the interatomic distance, ε = 3.4 kJ/mol is the depth of the potential well and σ = 3.4 Angstroms is the separation distance at which the potential is zero.

Plot V(r) at a grid of r values, e.g., using np.linspace(3.5, 10, 1000) similar to Problem 2c). Add the titles for the figure and the axes. Set appropriate plot limits. Save V-r values into the external file like Problem 2d).

# Problem 4 (Bonus): Lennard-Jones Potential for Argon
import numpy as np
import matplotlib.pyplot as plt

# Lennard-Jones parameters for Argon
epsilon = 3.4    # kJ/mol (depth of potential well)
sigma = 3.4      # Angstroms (distance where V = 0)

# Grid of interatomic distances
r = np.linspace(3.0, 10, 1000)