Worked Examples: Comparisons and Flow Control#
These worked solutions correspond to the exercises on the Comparisons and Flow Control page.
How to use this notebook:
Try each exercise yourself first before looking at the solution
The code cells show both the code and its output
Download this notebook if you want to run and experiment with the code yourself
Your solution might look different - that’s fine as long as it gives the correct answer!
Exercise 1: pH with Activity Coefficient#
Problem: The correct definition of pH is
where \(a_\mathrm{H}^+\) is the activity of hydrogen ions, which is related to the concentration \([\mathrm{H}^+]\) by the activity coefficient, \(\gamma_\mathrm{H}^+\):
Write an updated version of the calculate_pH() function from the earlier example so that it calculates pH correctly. Your function should take two arguments: the concentration of hydrogen ions and their activity coefficient.
Check that for \(\gamma_\mathrm{H}^+ = 1\) you get the same answer as before.
Now rewrite your function so that the activity coefficient is an optional argument, with a default value of \(\gamma_\mathrm{H}^+ = 1\). Check that if you call this function in your code as calculate_pH(h_plus) you still get the same result.
Solution
Part 1: Function with Two Required Arguments
First, let’s write the function with both the hydrogen ion concentration and activity coefficient as required arguments:
import math
def calculate_pH(h_plus, gamma):
"""
Calculates pH.
Args:
h_plus (float): Hydrogen ion concentration (M).
gamma (float): Activity coefficient (dimensionless).
Returns:
(float): pH value.
"""
activity = gamma * h_plus
return -math.log10(activity)
Let’s test this function with an example.
If \([\mathrm{H}^+] = 1.0 \times 10^{-5}\) M and \(\gamma_\mathrm{H}^+ = 1.0\):
h_plus = 1.0e-5
gamma = 1.0
pH = calculate_pH(h_plus, gamma)
print(f"pH = {pH:.2f}")
pH = 5.00
Verification: When \(\gamma_\mathrm{H}^+ = 1\), we get the same result as the simplified formula \(\mathrm{pH} = -\log_{10}[\mathrm{H}^+]\).
Part 2: Function with Optional Argument
Now let’s rewrite the function so that the activity coefficient is an optional argument with a default value of 1:
def calculate_pH(h_plus, gamma=1.0):
"""
Calculates pH.
Args:
h_plus (float): Hydrogen ion concentration (M).
gamma (float, optional): Activity coefficient (dimensionless).
Defaults to 1.0 (ideal behaviour).
Returns:
(float): pH value.
"""
activity = gamma * h_plus
return -math.log10(activity)
Now we can call this function in three different ways:
h_plus = 1.0e-5
# Method 1: Using the default gamma = 1.0
pH_default = calculate_pH(h_plus)
print(f"pH (default γ=1.0) = {pH_default:.2f}")
# Method 2: Explicitly specifying gamma = 1.0
pH_explicit = calculate_pH(h_plus, 1.0)
print(f"pH (explicit γ=1.0) = {pH_explicit:.2f}")
# Method 3: Using a different activity coefficient
pH_nonideal = calculate_pH(h_plus, 0.85)
print(f"pH (γ=0.85) = {pH_nonideal:.2f}")
pH (default γ=1.0) = 5.00
pH (explicit γ=1.0) = 5.00
pH (γ=0.85) = 5.07
Exercise 2: Using the Generic quadratic_roots() Function#
Problem: Rewrite your calculate_h_plus() function to call the generic quadratic_roots() function. You will need to pass in appropriate values for the arguments a, b, and c. Remember that the quadratic_roots() function returns both roots: you want the positive root.
Solution
Background: The Original Function
From the earlier example in the Functions section, the original calculate_h_plus() function was:
def calculate_h_plus(concentration, Ka):
"""Calculate the hydrogen ion concentration of a weak acid solution."""
return (-Ka + math.sqrt(Ka**2 + 4*Ka*concentration)) / 2
This function directly implements the positive root of the quadratic equation:
where \(x = [\mathrm{H}^+]\), derived from the equilibrium expression for weak acid dissociation.
The Generic quadratic_roots() Function
The generic function solves \(ax^2 + bx + c = 0\):
def quadratic_roots(a, b, c):
"""
Calculate both roots of a quadratic equation ax^2 + bx + c = 0.
Args:
a (float): Coefficient of x^2.
b (float): Coefficient of x.
c (float): Constant term.
Returns:
(tuple): A tuple containing both roots (r1, r2).
"""
d = b**2 - 4*a*c
r1 = (-b + math.sqrt(d)) / 2 / a
r2 = (-b - math.sqrt(d)) / 2 / a
return r1, r2
Rewritten calculate_h_plus() Function
Comparing our quadratic equation $\(x^2 + K_a x - K_a C = 0\)\( with the standard form \)\(ax^2 + bx + c = 0\)\( we see we want the following coefficients: \)\(a = 1\)\( \)\(b = K_a\)\( \)\(c = -K_a \times \text{concentration}\)$
Here’s the rewritten function:
def calculate_h_plus(concentration, Ka):
"""
Calculate the hydrogen ion concentration of a weak acid solution.
Args:
concentration (float): Initial concentration of weak acid (M).
Ka (float): Acid dissociation constant.
Returns:
(float): Equilibrium hydrogen ion concentration (M).
"""
# Set up coefficients for ax^2 + bx + c = 0
a = 1
b = Ka
c = -Ka * concentration
# Get both roots
root1, root2 = quadratic_roots(a, b, c)
# Return the positive root
# (root1 uses +sqrt, root2 uses -sqrt)
return root1
Testing the Function
Let’s test this with an example. For acetic acid with \(K_a = 1.8 \times 10^{-5}\) and initial concentration \(C = 0.1\) M:
# Example: acetic acid
Ka = 1.8e-5
concentration = 0.1
# Calculate H+ concentration
h_plus = calculate_h_plus(concentration, Ka)
print(f"[H+] = {h_plus:.6f} M")
# Calculate pH for verification
pH = -math.log10(h_plus)
print(f"pH = {pH:.2f}")
[H+] = 0.001333 M
pH = 2.88
Verification: Comparing with Original Function
Let’s verify that both versions give the same result:
# Original version (direct implementation)
def calculate_h_plus_original(concentration, Ka):
return (-Ka + math.sqrt(Ka**2 + 4*Ka*concentration)) / 2
# Test both functions
h_plus_new = calculate_h_plus(concentration, Ka)
h_plus_original = calculate_h_plus_original(concentration, Ka)
print(f"New function: [H+] = {h_plus_new:.10e} M")
print(f"Original function: [H+] = {h_plus_original:.10e} M")
New function: [H+] = 1.3326709731e-03 M
Original function: [H+] = 1.3326709731e-03 M
Benefits of this approach:
Reusability: We can now use quadratic_roots() for any quadratic equation
Clarity: The code explicitly shows the quadratic equation being solved
Maintainability: If we improve quadratic_roots() (e.g., error handling), all functions using it benefit
Alternative Solution: More Robust Approach
The solution above assumes we know that root1 (using the \(+\sqrt{\cdots}\) term) is the positive root. A more robust approach doesn’t rely on knowing the order the two roots are returned in quadratic_roots(). But then we need an approach that will give us the positive root whether this is root1 or root2. We can do this by using max(roots) to get the maximum root (which is guaranteed to be the positive root in our case).
def calculate_h_plus(concentration, Ka):
"""
Calculate the hydrogen ion concentration of a weak acid solution.
Args:
concentration (float): Initial concentration of weak acid (M).
Ka (float): Acid dissociation constant.
Returns:
(float): Equilibrium hydrogen ion concentration (M).
"""
# Set up coefficients for ax^2 + bx + c = 0
a = 1
b = Ka
c = -Ka * concentration
# Get both roots
roots = quadratic_roots(a, b, c)
# Return the positive root (the physically meaningful one)
return max(roots)
Exercise 3: Distance Between Atoms#
Problem: Write a function to calculate the distance between two atoms. Then rewrite the exercise in the loops section to utilise this function. Make sure to include a docstring describing the purpose of the function and outlining the arguments.
Solution
The distance between two atoms in 3D space is calculated using the Euclidean distance formula: $\(d = \sqrt{(x_2-x_1)^2 + (y_2-y_1)^2 + (z_2-z_1)^2}\)$ Here’s the function implementation:
def distance(atom1, atom2):
"""
Calculate the distance between two atoms in 3D space.
Args:
atom1 (list): [x, y, z] coordinates of atom 1 (Å).
atom2 (list): [x, y, z] coordinates of atom 2 (Å).
Returns:
(float): Distance between the two atoms (Å).
"""
dx = atom1[0] - atom2[0]
dy = atom1[1] - atom2[1]
dz = atom1[2] - atom2[2]
return math.sqrt(dx**2 + dy**2 + dz**2)
Testing the function:
# Example: distance between two atoms
atom_a = [0.0, 0.0, 0.0]
atom_b = [0.0, 0.94, 0.38]
d = distance(atom_a, atom_b)
print(f"Distance: {d:.4f} Å")
Distance: 1.0139 Å
Using the Function in Nested Loops:
Now let’s use this function to calculate all pairwise distances between three atoms. From the loops section, we want to avoid calculating the same distance twice (e.g., distance 1\(\to\)2 is the same as 2\(\to\)1) and avoid calculating the distance from an atom to itself.
# Ammonia (NH3) atomic coordinates in Ångströms
atom_N = [0.0, 0.0, 0.0] # Nitrogen at origin
atom_H1 = [0.0, 0.94, 0.38] # Hydrogen 1
atom_H2 = [0.81, -0.47, 0.38] # Hydrogen 2
atom_H3 = [-0.81, -0.47, 0.38] # Hydrogen 3
atoms = [atom_N, atom_H1, atom_H2, atom_H3]
atom_labels = ['N', 'H1', 'H2', 'H3']
# Calculate all unique pairwise distances using nested loops
print("Unique pairwise distances in NH3")
print("=" * 40)
for i in range(len(atoms)):
for j in range(i + 1, len(atoms)):
d = distance(atoms[i], atoms[j])
print(f"Distance from {atom_labels[i]:3s} to {atom_labels[j]:3s}: {d:.4f} Å")
Unique pairwise distances in NH3
========================================
Distance from N to H1 : 1.0139 Å
Distance from N to H2 : 1.0106 Å
Distance from N to H3 : 1.0106 Å
Distance from H1 to H2 : 1.6261 Å
Distance from H1 to H3 : 1.6261 Å
Distance from H2 to H3 : 1.6200 Å
Benefits of using a function:
Reusability: We can calculate any pairwise distance with a simple function call:
distance(atoms[i], atoms[j])Readability: The code in the loop is much clearer than writing out the full distance formula
Maintainability: If we need to change the distance calculation (e.g., add periodic boundary conditions), we only update one function
Testing: We can test the distance function independently before using it in the loop
Alternative Solution: Using enumerate() and zip()
Here’s a more Pythonic approach that works directly with atom coordinates and labels rather than indices:
# Calculate all unique pairwise distances using nested loops
print("Unique pairwise distances in NH3")
print("=" * 40)
for i, (atom_i, label_i) in enumerate(zip(atoms, atom_labels)):
for atom_j, label_j in zip(atoms[i+1:], atom_labels[i+1:]):
d = distance(atom_i, atom_j)
print(f"Distance from {label_i:3s} to {label_j:3s}: {d:.4f} Å")
Unique pairwise distances in NH3
========================================
Distance from N to H1 : 1.0139 Å
Distance from N to H2 : 1.0106 Å
Distance from N to H3 : 1.0106 Å
Distance from H1 to H2 : 1.6261 Å
Distance from H1 to H3 : 1.6261 Å
Distance from H2 to H3 : 1.6200 Å
Benefits of this approach:
Works directly with values rather than indices
More readable inside the loops:
distance(atom_i, atom_j)instead ofdistance(atoms[i], atoms[j])Less error-prone: no risk of mixing up indices inside the loops
Exercise 4: Einstein Model for Heat Capacity#
Problem: The Einstein model for the heat capacity of a crystal is:
where the Einstein temperature, \(\Theta_\mathrm{E}\) is a constant for a given material.
Write a function to calculate \(C_{V,m}\) at 300 K for:
(a) sodium (\(\Theta_\mathrm{E}\) = 192 K)
(b) diamond (\(\Theta_\mathrm{E}\) = 1450 K)
Solution:
def einstein_heat_capacity(T, theta_E):
"""
Calculate the molar heat capacity using the Einstein model.
Args:
T (float): Temperature (K).
theta_E (float): Einstein temperature for the material (K).
Returns:
(float): Molar heat capacity at constant volume, C_V,m (J/(mol·K)).
"""
R = 8.314 # Gas constant in J/(mol·K)
# Calculate the ratio theta_E / T
x = theta_E / T
# Calculate exp(theta_E / T)
exp_x = math.exp(x)
# Apply the Einstein formula
C_V = 3 * R * (x**2) * exp_x / (exp_x - 1)**2
return C_V
# Part (a): Sodium at 300 K
T = 300 # K
theta_E_sodium = 192 # K
C_V_sodium = einstein_heat_capacity(T, theta_E_sodium)
print(f"Sodium at {T} K:")
print(f"θ_E = {theta_E_sodium} K")
print(f"C_V,m = {C_V_sodium:.2f} J/(mol·K)")
print()
# Part (b): Diamond at 300 K
theta_E_diamond = 1450 # K
C_V_diamond = einstein_heat_capacity(T, theta_E_diamond)
print(f"Diamond at {T} K:")
print(f"θ_E = {theta_E_diamond} K")
print(f"C_V,m = {C_V_diamond:.2f} J/(mol·K)")
Sodium at 300 K:
θ_E = 192 K
C_V,m = 24.11 J/(mol·K)
Diamond at 300 K:
θ_E = 1450 K
C_V,m = 4.71 J/(mol·K)
Breaking down the formula:
To make the calculation clearer and avoid errors, we break it into parts:
Define \(x = \frac{\Theta_E}{T}\)
Calculate \(e^x = \exp(x)\)
Apply the formula: \(C_{V,m} = 3R \cdot x^2 \cdot \frac{e^x}{(e^x - 1)^2}\)
Why we use intermediate variables:
x = theta_E / T- makes the formula more readableexp_x = math.exp(x)- we use this value twice, so calculate it once
This approach reduces errors and makes debugging easier