Maintainer triage (June 29, 2026)
- Status: Verified bug against the current
master source in this fork.
- Severity: High.
- Area labels: area:kinetics, area:crystallizers.
- Tackle batch: P4 - high-impact numerical/model correctness.
The source review confirmed that the cited code path and failure mode are still present. The original audit details are preserved below for reproduction notes and suggested fixes.
Original audit details
Part of the model-correctness audit — see #17.
Severity: high | Category: logic | Location: PharmaPy/Kinetics.py lines 959-967
What's wrong
In deriv_cryst.dmech_dparam, the non-reformulated branch computes the derivative of a mechanism rate w.r.t. the activation energy E (params[1]) as -mech * np.exp(e_act/gas_ct/temp) (line 966). The rate model in cryst_mechanism is mech = kexp(-E/(RT))S^exp (line 35, 37). The correct partial derivative is d(mech)/dE = mech * (-1/(RT)) = -mech/(gas_cttemp). The code instead multiplies by exp(+E/(RT)), which is neither dimensionally correct (it has units of 1/[mech] inverse) nor numerically equal to -mech/(R*T). For typical E ~ 4e4 J/mol, T=300 K, exp(E/RT) ~ exp(16) ~ 9e6, so the reported derivative is off by many orders of magnitude and has the wrong functional form. The reformulated branch (lines 952-957) is correct, confirming the non-reformulated branch is the erroneous one.
Why it's incorrect
mech = k * exp(-E/(RT)) * S^exp. Partial derivative w.r.t. E: dmech/dE = k * S^exp * exp(-E/(RT)) * (-1/(RT)) = -mech/(RT). The code's -mech*exp(+E/(R*T)) replaces the factor 1/(RT) with exp(+E/(RT)), which is dimensionally and numerically wrong. This corrupts any gradient-based parameter estimation / sensitivity that consumes deriv_cryst.
Suggested regression test
Compare deriv_cryst output element [1] against a finite-difference estimate d(mech)/dE = (mech(E+dE)-mech(E))/dE for a single mechanism (e.g. growth) with reformulate_kin=False. The finite difference will match -mech/(gas_cttemp), not -mechexp(E/(gas_ct*temp)).
Verification: hand-confirmed by direct code reading.
Maintainer triage (June 29, 2026)
mastersource in this fork.The source review confirmed that the cited code path and failure mode are still present. The original audit details are preserved below for reproduction notes and suggested fixes.
Original audit details
Part of the model-correctness audit — see #17.
Severity: high | Category: logic | Location:
PharmaPy/Kinetics.pylines 959-967What's wrong
In deriv_cryst.dmech_dparam, the non-reformulated branch computes the derivative of a mechanism rate w.r.t. the activation energy E (params[1]) as
-mech * np.exp(e_act/gas_ct/temp)(line 966). The rate model in cryst_mechanism is mech = kexp(-E/(RT))S^exp (line 35, 37). The correct partial derivative is d(mech)/dE = mech * (-1/(RT)) = -mech/(gas_cttemp). The code instead multiplies by exp(+E/(RT)), which is neither dimensionally correct (it has units of 1/[mech] inverse) nor numerically equal to -mech/(R*T). For typical E ~ 4e4 J/mol, T=300 K, exp(E/RT) ~ exp(16) ~ 9e6, so the reported derivative is off by many orders of magnitude and has the wrong functional form. The reformulated branch (lines 952-957) is correct, confirming the non-reformulated branch is the erroneous one.Why it's incorrect
mech = k * exp(-E/(RT)) * S^exp. Partial derivative w.r.t. E: dmech/dE = k * S^exp * exp(-E/(RT)) * (-1/(RT)) = -mech/(RT). The code's
-mech*exp(+E/(R*T))replaces the factor 1/(RT) with exp(+E/(RT)), which is dimensionally and numerically wrong. This corrupts any gradient-based parameter estimation / sensitivity that consumes deriv_cryst.Suggested regression test
Compare deriv_cryst output element [1] against a finite-difference estimate d(mech)/dE = (mech(E+dE)-mech(E))/dE for a single mechanism (e.g. growth) with reformulate_kin=False. The finite difference will match -mech/(gas_cttemp), not -mechexp(E/(gas_ct*temp)).
Verification: hand-confirmed by direct code reading.