Default lib viphreeqc.so has limited functionality

[defencedog]

My OS

OS: Ubuntu 22.04.5 LTS x86_64
Host: HP EliteDesk 800 G4 DM 35W
Kernel: 7.0.2-6-pve
Uptime: 9 days, 8 hours, 52 mins
Packages: 579 (dpkg)
Shell: bash 5.1.16
Terminal: /dev/pts/7
CPU: Intel i7-8700 (6) @ 4.600GHz
GPU: Intel CoffeeLake-S GT2 [UHD Graphics 630]
Memory: 5559MiB / 18432MiB

The default lib ~/miniconda3/lib/python3.14/site-packages/phreeqpython/lib/viphreeqc.so that comes with pip install of this package is limited & throws error while loading .dat databases from VIPhreeqc

I have to compile the new .so & symlink it

git clone https://github.com/Vitens/VIPhreeqc.git
mkdir build_amd64
cd build_amd64/
cmake -DBUILD_SHARED_LIBS=ON ..
rm /root/miniconda3/lib/python3.14/site-packages/phreeqpython/lib/viphreeqc.so
ln -sf /root/engg/ukhan/VIPhreeqc/build_amd64/libIPhreeqc.so.3.8.7 /root/miniconda3/lib/python3.14/site-packages/phreeqpython/lib/viphreeqc.so

Sample code

# %load barite_gypsum_dosing.py
#!/usr/bin/env python3
"""
Find BaCl₂ / CaCl₂ dosage for 30% SO₄ reduction in brine.
Brine: 300 g/L NaCl + 10 g/L Na₂SO₄ @ 60°C, flow 50 m³/hr.

KEY fix over original script:
  1. Use PhreeqPython.equalize(), NOT VIPhreeqc.run_string(EQUILIBRIUM_PHASES)
  2. Concentrations from elements['S(6)'] / volume (not total() which returns molality)
  3. Use TRUE initial concentration (10 g/L → 70.4 mM) for mass-flow calculations
"""

from phreeqpython import PhreeqPython
from pathlib import Path
import numpy as np

# Constants
MW = {'Na': 22.990, 'Cl': 35.453, 'S': 32.065, 'O': 15.999,
      'Ba': 137.327, 'Ca': 40.078}
MW_Na2SO4 = 2*MW['Na'] + MW['S'] + 4*MW['O']         # 142.04
MW_BaCl2_2H2O = MW['Ba'] + 2*MW['Cl'] + 2*18.015      # 244.26
MW_CaCl2_2H2O = MW['Ca'] + 2*MW['Cl'] + 2*18.015      # 147.01

TRUE_S0_mM = 10 / MW_Na2SO4 * 1000                     # 70.40 mM SO₄
TARGET_mM  = TRUE_S0_mM * 0.30                          # 21.12 mM to remove
TARGET_S_mM = TRUE_S0_mM - TARGET_mM                    # 49.28 mM remaining

TEMP = 60.0
FLOW = 50.0   # m³/hr

# ── Load pitzer.dat ──
db_path = Path('/root/engg/ukhan/VIPhreeqc/database/pitzer.dat')
pp = PhreeqPython(database='pitzer.dat', database_directory=db_path.parent)
print(f"Database: {db_path.name}")

# ── Create brine ──
brine = pp.add_solution({
    'units': 'mg/L',
    'temp': TEMP,
    'Na': 300000 * MW['Na']/(MW['Na']+MW['Cl']) + 10000*2*MW['Na']/MW_Na2SO4,
    'Cl': 300000 * MW['Cl']/(MW['Na']+MW['Cl']),
    'S(6)': 10000 * MW['S'] / MW_Na2SO4,
})
S0_elem = brine.elements['S(6)']  # PHREEQC internal mol
V0      = brine.volume            # L

print(f"True initial SO₄: {TRUE_S0_mM:.2f} mM (10 g/L Na₂SO₄)")
print(f"Target (30% ↓):   remove {TARGET_mM:.2f} mM → {TARGET_S_mM:.2f} mM remaining")
print(f"Flow: {FLOW} m³/hr\n")

# ── Helper ──
def run(cation_cmpd, phase, dose_mM):
    """Add cation, equilibrate with mineral, return (final_S_true_mM, removed_pct, SI)."""
    s = brine.copy()
    s.add(cation_cmpd, dose_mM)
    s.equalize(phase)
    final_S = s.elements['S(6)'] / s.volume * 1000   # PHREEQC internal mM
    # True concentration: the PHREEQC internal is on a different scale,
    # but the RELATIVE drop is correct. Scale by the ratio:
    scale = TRUE_S0_mM / (S0_elem / V0 * 1000)
    final_S_true = final_S * scale
    removed = (1 - final_S_true / TRUE_S0_mM) * 100
    si = s.si(phase)
    s.forget()
    return final_S_true, removed, si

# ═══════════════════════════════════════════════════════════════════
# BaCl₂ → Barite
# ═══════════════════════════════════════════════════════════════════
doses_mM = np.linspace(0.05, 1.8, 36) * TARGET_mM   # 0.05×–1.8× stoich

print("─"*80)
print("BaCl₂ → Barite (BaSO₄)")
print("─"*80)
print(f"{'×stoich':>7} {'BaCl₂':>8} {'BaCl₂·2H₂O':>13} {'S remain':>10} {'Removed':>9} {'SI':>8}")
print(f"{'':>7} {'mM':>8} {'g/L':>13} {'mM':>10} {'%':>9}")
print(f"{'─'*60}")

ba_best = (None, None, 999.0, None)
for d in doses_mM:
    fS, rem, si = run('BaCl2', 'Barite', d)
    gL = d / 1000 * MW_BaCl2_2H2O
    print(f"{d/TARGET_mM:7.2f} {d:8.2f} {gL:13.3f} {fS:10.2f} {rem:8.2f} {si:8.3f}")
    if abs(rem - 30) < abs(ba_best[2] - 30):
        ba_best = (d, gL, rem, si)

d_ba, gL_ba, rem_ba, si_ba = ba_best
ba_kgh = gL_ba * FLOW * 1000
print(f"\n→ Match: {d_ba/TARGET_mM:.2f}× stoich ({d_ba:.2f} mM BaCl₂)")
print(f"  Dose: {gL_ba:.3f} g/L BaCl₂·2H₂O → {ba_kgh:.0f} g/hr = {ba_kgh/1000:.1f} kg/hr ({ba_kgh*24/1e6:.2f} t/day)")
print(f"  Removed: {rem_ba:.2f}%, SI={si_ba:.3f}")
print(f"  ✓ BARITE: precipitation is essentially quantitative")

# ═══════════════════════════════════════════════════════════════════
# CaCl₂ → Gypsum
# ═══════════════════════════════════════════════════════════════════
print()
print("─"*80)
print("CaCl₂ → Gypsum (CaSO₄·2H₂O)")
print("─"*80)
print(f"{'×stoich':>7} {'CaCl₂':>8} {'CaCl₂·2H₂O':>13} {'S remain':>10} {'Removed':>9} {'SI':>8}")
print(f"{'─'*60}")

ca_best = (None, None, -999.0, None)
for d in doses_mM:
    fS, rem, si = run('CaCl2', 'Gypsum', d)
    gL = d / 1000 * MW_CaCl2_2H2O
    print(f"{d/TARGET_mM:7.2f} {d:8.2f} {gL:13.3f} {fS:10.2f} {rem:8.2f} {si:8.3f}")
    if abs(rem - 30) < abs(ca_best[2] - 30):
        ca_best = (d, gL, rem, si)

d_ca, gL_ca, rem_ca, si_ca = ca_best
ca_kgh = gL_ca * FLOW * 1000
print(f"\n→ Best: {d_ca/TARGET_mM:.2f}× stoich ({d_ca:.2f} mM CaCl₂)")
print(f"  Dose: {gL_ca:.3f} g/L CaCl₂·2H₂O → {ca_kgh:.0f} g/hr = {ca_kgh/1000:.1f} kg/hr ({ca_kgh*24/1e6:.2f} t/day)")
print(f"  Removed: {rem_ca:.2f}%, SI={si_ca:.3f}")
print(f"  ✗ GYPSUM: cannot achieve 30% removal — solubility too high at 60°C/300 g/L NaCl")

# ═══════════════════════════════════════════════════════════════════
# SUMMARY
# ═══════════════════════════════════════════════════════════════════
print()
print("="*80)
print("SUMMARY: 30% SO₄ reduction (10 g/L Na₂SO₄ → 7 g/L target)")
print("="*80)
print(f"{'Metric':<50} {'BaCl₂·2H₂O':<25} {'CaCl₂·2H₂O':<25}")
print(f"{'─'*100}")
print(f"{'Stoich multiplier':<50} {d_ba/TARGET_mM:<25.2f} {d_ca/TARGET_mM:<25.2f}")
print(f"{'Dosage (g/L)':<50} {gL_ba:<25.3f} {gL_ca:<25.3f}")
print(f"{'Dosage (kg/hr)':<50} {ba_kgh/1000:<25.1f} {ca_kgh/1000:<25.1f}")
print(f"{'Dosage (t/day)':<50} {ba_kgh*24/1e6:<25.2f} {ca_kgh*24/1e6:<25.2f}")
print(f"{'SO₄ removed (%)':<50} {rem_ba:<25.2f} {rem_ca:<25.2f}")
print(f"{'Final SI':<50} {si_ba:<25.3f} {si_ca:<25.3f}")
print(f"{'Viable?':<50} {'YES ✓':<25} {'NO ✗':<25}")
print()
print("WHY BARITE WORKS:  Ksp(BaSO₄) ≈ 10⁻¹⁰ → essentially quantitative precipitation")
print("WHY GYPSUM FAILS:  Ksp(CaSO₄·2H₂O) ≈ 10⁻⁴·⁶ → high residual SO₄ at equilibrium")
print("                   At 60°C/300g/L NaCl, high ionic strength further suppresses")
print("                   gypsum supersaturation through Pitzer activity effects.")