Used Python to solve for acid base titration. Can any chemists see if this code is accurate with HCl volume output?
I ran an experiment in a lab and the output of my code double the amount of HCl I actually put in.
I checked the concentrations of the three carbon species and they seemed accurate according to the Bjerrum plot based off my initial molarity.
I did not take into account the CO2 that can escape into the air since the test I conducted was not in a closed system.
Any constraints that are not explicitly mentioned are as follows: this titration is performed within a bioreactor
The code will be posted below the parenthesis, if anyone has tips on optimizing the recursion or anything else in general I'd love to hear! Thank you for your time.
# This code is to simulate pH decrease as acid is added into a solution
# Based on the solution of sodium bicarbonate and hydrochloric acid to carbonic acid and sodium chloride
# NOTE: Calculations derived from the formula pH = pKa - log(base/acid)
print(f"NOTE: Do not enter a pH of greater than 13 and/or less than 1\n "
f"Do not enter a pH titration of less than 0.1\n "
f"Do not enter a pH titration of greater than the pH difference\n")
hcl_list = []
# List keeps track of volumes of titration per pH change
h2co3_mol_list = []
hco3_mol_list = []
co3_mol_list = []
total_volume = 0
# Variable for summation of pH change for function
total_hcl = 0
# Placeholder for summed HCl moles
pka_1 = round(6.37, 3)
# Physical constant for bicarbonate to carbonic acid
pka_2 = round(10.33, 3)
# Physical constant for bicarbonate to carbonate
full_bicarbonate_ph = round((pka_1 + pka_2) / 2, 3)
# pH of maximum bicarbonate ~ 8.35
mol_mass = round(84.007, 3)
# Physical constant for sodium bicarbonate molar mass
delta_h_1 = round(-950.8, 3)
delta_h_2 = round(-1130.8, 3)
# Units: kJ/mol
r = round(8.314, 3)
standard_temperature = round(273.15 + 25, 3)
initial_ph_input = round(float(input("Enter the initial pH:\n")), 3)
final_ph_input = round(float(input("\nEnter the final pH:\n")), 3)
increments_input = input("\nEnter the increment of pH change: ")
if increments_input == "":
increments_input = round(float(abs(initial_ph_input - final_ph_input)), 3)
print(f"Singular Titration of {round(increments_input, 3)} pH change")
# User inputs for initial pH, final pH, and change in pH increments
elif increments_input != "":
increments_input = float(increments_input)
increments_input = round(increments_input, 3)
# User can choose to leave titration step empty for single titration from initial pH to final pH
volume_input = round(float(input("\nEnter the container volume in Liters:\n")), 3)
gram_input = round(float(input("\nEnter the number of grams per Liter of sodium bicarbonate:\n")), 3)
initial_hcl_moles_input = round(float(input("\nEnter the acid molarity:\n")), 3)
initial_hcl_volume_input = round(float(input("\nEnter the acid volume in Liters:\n")), 3)
# User inputs for container volume, solution concentration, and acid molarity
initial_hco3_molarity = round(gram_input / mol_mass, 3)
initial_co3_molarity = round(0, 3)
initial_h2co3_molarity = round(0, 3)
initial_hydronium_molarity = 10 ** -initial_ph_input
# Initial molarity
ph_range = []
# pH range for graph depiction
remaining_volume = round(2 - volume_input, 3)
if remaining_volume == 0:
print(f"\nNo room for HCl acid addition, try again")
sys.exit(0)
# Exits program if initial volume is max capacity for container
if initial_ph_input > full_bicarbonate_ph:
adjustment_a = round(10 ** (initial_ph_input - pka_2), 3)
initial_co3_molarity = round((adjustment_a * initial_hco3_molarity) / (adjustment_a + 1), 3)
initial_hco3_molarity -= initial_co3_molarity
initial_hco3_molarity = round(initial_hco3_molarity, 3)
pka_2 = round(math.log10(initial_hco3_molarity / (initial_hydronium_molarity * initial_co3_molarity)), 3)
elif initial_ph_input <= full_bicarbonate_ph:
adjustment_b = round(10 ** (initial_ph_input - pka_1), 3)
initial_h2co3_molarity = round(initial_hco3_molarity / (adjustment_b + 1), 3)
initial_hco3_molarity -= initial_h2co3_molarity
initial_hco3_molarity = round(initial_hco3_molarity, 3)
pka_1 = round(math.log10(initial_h2co3_molarity / (initial_hydronium_molarity * initial_hco3_molarity)), 3)
# Accounts for proportionate concentrations for initial conditions and adjusted pka
def temperature_correction(pka, delta_h):
new_pka = abs(round(math.log10((10 ** -pka) * math.exp(-1 * (delta_h / r) * ((1 / temperature_input)
- (1 / standard_temperature)))), 3))
print(f"At {temperature_input} Kelvin, the new pKa is {new_pka}\n")
return new_pka
# Corrects pKa for temperature
# Equation: pka_t2 = -log10( (10 ^ -pka_t1) * e ^ ( -(H / R) * ( (1 / T2) - (1 / T1)) )
# NOTE: Only accounts for CONSTANT Temperature, it is not a dynamic temperature reading
temperature_check = input(f"\nIs there a new temperature to account for? (Yes/No):\n")
if temperature_check == "Yes" or temperature_check == "yes" or temperature_check == "Y" or temperature_check == "y":
temperature_input = round(float(input("\nEnter the container temperature in Celsius:\n")) + 273.15, 3)
pka_1 = temperature_correction(pka_1, delta_h_1)
pka_2 = temperature_correction(pka_2, delta_h_2)
# Corrects pKa for temperature
def find_volume(start, increments, goal, hco3_molarity, co3_molarity,
h2co3_molarity, container_volume, total_acid_volume):
ph = round(start, 3)
# pH variable starts off at the initial pH
h2co3_moles = round(h2co3_molarity * volume_input, 3)
hco3_moles = round(hco3_molarity * volume_input, 3)
co3_moles = round(co3_molarity * volume_input, 3)
new_h2co3_molarity = round(0, 3)
new_hco3_molarity = round(0, 3)
new_co3_molarity = round(0, 3)
total_hcl_moles = round(0, 3)
if ph == goal:
return
# Ends function after pH variable reaches desired value
else:
if ph > full_bicarbonate_ph:
# co3_moles = hco3_moles * (10 ** (ph - pka_2))
# Formula: pH = pKa + log(base / acid) -> base = acid * 10 ^ (pH - pKa)
adjustment_2 = 10 ** (ph - (pka_2 + increments))
# Adjustment for initial run = 10 ^ (10 - (6.37 + 0.1)) = 10 ^ (10 - 6.47) = 10 ^ 3.53
hcl_moles_2 = (co3_moles - (adjustment_2 * hco3_moles)) / (adjustment_2 + 1)
# - total_hcl_moles
hcl_molarity_2 = round(hcl_moles_2 / container_volume, 3)
# Used to calculate the new ratio of acid and base for next titration
new_hco3_molarity = round(hco3_molarity + hcl_molarity_2, 3)
new_co3_molarity = round(co3_molarity - hcl_molarity_2, 3)
# print(hcl_moles_2, new_co3_molarity, co3_molarity, new_hco3_molarity, hco3_molarity)
total_hcl_moles = hcl_moles_2
# Add HCl for next titration down
dic = round(sum([hco3_moles, co3_moles, h2co3_moles]), 3)
new_dic = round(sum([new_hco3_molarity, new_co3_molarity, new_h2co3_molarity])
* volume_input / container_volume, 3)
print(f"Initial Dissolved Carbon: {dic}, Final Dissolved Carbon: {new_dic}")
elif ph <= full_bicarbonate_ph:
# h2co3_moles = hco3_moles / (10 ** (ph - pka_1))
# Formula: pH = pKa + log(base / acid) -> acid = base / 10 ^ (pH - pKa)
adjustment_1 = 10 ** (ph - (pka_1 + increments))
# Adjustment for initial run = 10 ^ (10 - (6.37 + 0.1)) = 10 ^ (10 - 6.47) = 10 ^ 3.53
hcl_moles_1 = (hco3_moles - (adjustment_1 * h2co3_moles)) / (adjustment_1 + 1)
# - total_hcl_moles
hcl_molarity_1 = round(hcl_moles_1 / container_volume, 3)
# Used to calculate the new ratio of acid and base for next titration
new_h2co3_molarity = round(h2co3_molarity + hcl_molarity_1, 3)
new_hco3_molarity = round(hco3_molarity - hcl_molarity_1, 3)
# print(hcl_molarity_1, new_h2co3_molarity, h2co3_molarity, new_hco3_molarity, hco3_molarity)
total_hcl_moles = hcl_moles_1
# Add HCl for next titration down
dic = round(sum([hco3_moles, co3_moles, h2co3_moles]) / container_volume, 3)
new_dic = round(sum([new_hco3_molarity, new_co3_molarity, new_h2co3_molarity])
* volume_input / container_volume, 3)
print(f"Initial Dissolved Carbon: {dic}, Final Dissolved Carbon: {new_dic}\n\n")
# HCl is calculated from the following formula:
# 1) pH - increment = pKa + log(base - (x + t) / acid + (x + t)
# NOTE: x is the HCl moles needed for titration
# NOTE: t is the total HCl moles already present in the system so far
# NOTE: t is reset during each recursion/iteration -> high volume of HCl required
# 2) 10^(pH - (pKa + increment)) = base - (x + t) / acid + (x + t)
# NOTE: Adjustment = 10^(pH - (pKa + increment))
# Adjustment = (base - x - t) / (acid + x + t)
# 3) (Adjustment * x) + (Adjustment * (acid + t)) = base - x - t
# 4) ((Adjustment + 1) * x) = base - (Adjustment * acid) - ((Adjustment + 1) * t)
# 5) x = (base - (Adjustment * acid) - ((Adjustment + 1) * t)) / (Adjustment + 1)
hcl_volume = total_hcl_moles * initial_hcl_volume_input / initial_hcl_moles_input
# Acid volume calculated from the following formula: C1V1 = C2V2
# Final_Acid_Molarity / Final_Acid_Volume = Initial_Acid_Molarity / Initial_Acid_Volume
# Final_Acid_Volume = Final_Acid_Molarity * Initial_Acid_Volume / Initial_Acid_Molarity
new_volume = hcl_volume + container_volume
# Takes total volume of acid + initial container volume
total_acid_volume += hcl_volume
if total_acid_volume >= remaining_volume:
print(f"\nAny pH value below exceeds container capacity\n")
return
# Placeholder variable for logging total acid volume per pH change
# Exits program once total acid volume reaches container capacity
hcl_list.append(total_acid_volume)
# Adds value of acid into a list for graph visual
h2co3_mol_list.append(h2co3_moles)
hco3_mol_list.append(hco3_moles)
co3_mol_list.append(co3_moles)
ph_range.append(int(start * 10))
new_ph = round(ph - increments, 2)
print(f"\nFrom {round(ph, 1)} to {new_ph},\n"
f"Moles of HCl required: {round(total_hcl_moles, 3)}, "
f"Volume of HCl required [mL]: {round(hcl_volume * 1000, 3)}, "
f"Total HCl required [mL]: {round(total_acid_volume * 1000, 3)}\n"
f"Carbonic Acid Concentration [Moles]: {round(h2co3_moles, 3)}, "
f"Bicarbonate Concentration [Moles]: {round(hco3_moles, 3)}, "
f"Carbonate Concentration [Moles]: {round(co3_moles, 3)}\n")
return find_volume(ph - increments, increments, goal, new_hco3_molarity,
new_co3_molarity, new_h2co3_molarity, new_volume, total_acid_volume)
find_volume(initial_ph_input, increments_input, final_ph_input, initial_hco3_molarity,
initial_co3_molarity, initial_h2co3_molarity, volume_input, total_volume)
dec_range = [initial_ph_input - (int(x) / 10) for x in ph_range]
plt.plot(dec_range, hcl_list)
plt.xlabel(f'Change in pH starting at {initial_ph_input}')
plt.ylabel('HCl [L]')
plt.title('Total HCl Volume to pH change')
plt.show()
carbon_ph_range = [(int(x) / 10) for x in ph_range]
# print(co3_mol_list, hco3_mol_list, h2co3_mol_list)
plt.plot(carbon_ph_range, h2co3_mol_list, color='r', label='H2CO3')
plt.plot(carbon_ph_range, hco3_mol_list, color='g', label='HCO3')
plt.plot(carbon_ph_range, co3_mol_list, color='b', label='CO3')
plt.xlabel(f'pH')
plt.ylabel(f'Carbon Species [Moles]')
plt.title('Carbon Species to pH ratio')
plt.legend()
plt.show()