carbonate system calculations

[5ir-Lancelot]

I want to calculate the DIC in an open carbonate system, assuming just NaHCO3 or NaOH is present. The input should be the pH, T, and pCO2.

So far I just modelled it step by step.

#shift the pH with NaOH (for alkaline range) or HCl (for acidic range)
 sol.change_ph(pH,'NaOH')
 
#concentration CO2
p=415 # ppm

#partial pressure pCO2 in atm
pCO2=p*1e-6            
            
#phreeqc always uses log10 values (NaOH will be used to increase and HCl used to lower the pH)
input_pCO2=np.log10(pCO2)
            
            
#equalize solution with CO2
#equilibrium with the atmosphere
sol.equalize(['CO2(g)'], [input_pCO2])

print(sol.pH)

The problem is that I basically want to change pH in equilibrium with CO2(g) directly. With the current approach, the resulting pH will always be much lower which makes sense.

Is there a function that I can use to get directly the water speciation for a given pH, pCO2 and T ? Assuming just pure carbonate system alkalinity is present?