Commit db56b4e4 authored by eladnoor's avatar eladnoor

adding necessary function in PhasedReaction for the website

parent 88d9ecfc
Pipeline #48165563 (#130) failed with stage
in 17 minutes and 52 seconds
......@@ -30,7 +30,7 @@ del get_versions
import numpy as np
from equilibrator_cache import Compound, ureg, Q_, R
from equilibrator_cache import Compound
from equilibrator_cache import CompoundCache
from equilibrator_cache.thermodynamic_constants import (
FARADAY, POSSIBLE_REACTION_ARROWS, Q_, R, default_I, default_pH,
......
......@@ -27,6 +27,7 @@
from typing import Dict, List, Tuple
from uncertainties import ufloat
import numpy as np
from component_contribution import GibbsEnergyPredictor
......@@ -35,12 +36,12 @@ from . import (
FARADAY, Compound, R, default_I, default_pH, default_pMg, default_T, ureg)
from .phased_reaction import PhasedReaction
class ComponentContribution(object):
"""A wrapper class for GibbsEnergyPredictor.
Also holds default conditions for compounds in the different phases.
"""
predictor = GibbsEnergyPredictor()
@ureg.check(None, None, None, '[concentration]', '[temperature]')
def __init__(
......@@ -57,7 +58,6 @@ class ComponentContribution(object):
:param ionic_strength:
:param temperature:
"""
self.predictor = GibbsEnergyPredictor()
self.p_h = p_h
self.ionic_strength = ionic_strength
self.p_mg = p_mg
......@@ -82,19 +82,20 @@ class ComponentContribution(object):
U as the covariance of a Gaussian used for sampling
(where dG0_cc is the mean of that Gaussian).
"""
rxn_dg_pairs = map(lambda r: r.split_aqueous_phase(
rxn_dg_pairs = map(lambda r: r.separate_stored_dgf_values(
p_h=self.p_h, ionic_strength=self.ionic_strength,
temperature=self.temperature), reactions)
rxn_aquoues, dg_prime_nonaq = zip(*rxn_dg_pairs)
dg_prime_nonaq = np.array(dg_prime_nonaq)
standard_dg_prime, dg_sigma = self.predictor.standard_dg_prime_multi(
rxn_aquoues,
p_h=self.p_h,
ionic_strength=self.ionic_strength,
temperature=self.temperature
)
return standard_dg_prime + dg_prime_nonaq * self.RT, dg_sigma
residual_reactions, stored_dg_primes = zip(*rxn_dg_pairs)
stored_dg_primes = np.array(stored_dg_primes)
standard_dg_prime, dg_sigma = \
ComponentContribution.predictor.standard_dg_prime_multi(
residual_reactions,
p_h=self.p_h,
ionic_strength=self.ionic_strength,
temperature=self.temperature
)
return standard_dg_prime + stored_dg_primes, dg_sigma
def standard_dg_prime(
self,
......@@ -107,38 +108,36 @@ class ComponentContribution(object):
calculate the confidence interval, use the range -1.96 to 1.96 times
this value
"""
rxn_aquoues, dg_prime_nonaq = reaction.split_aqueous_phase(
residual_reaction, stored_dg_prime = reaction.separate_stored_dgf_values(
p_h=self.p_h,
ionic_strength=self.ionic_strength,
temperature=self.temperature
)
standard_dg_prime, dg_sigma = self.predictor.standard_dg_prime(
rxn_aquoues,
p_h=self.p_h,
ionic_strength=self.ionic_strength,
temperature=self.temperature
)
standard_dg_prime, dg_sigma = \
ComponentContribution.predictor.standard_dg_prime(
residual_reaction,
p_h=self.p_h,
ionic_strength=self.ionic_strength,
temperature=self.temperature
)
return standard_dg_prime + dg_prime_nonaq * self.RT, dg_sigma
return standard_dg_prime + stored_dg_prime, dg_sigma
def dg_prime(
self,
reaction: PhasedReaction,
compound_to_conc: Dict[Compound, float]
reaction: PhasedReaction
) -> Tuple[float, float]:
"""Calculate the dG'0 of a single reaction.
:param reaction: an object of type Reaction
:param compound_to_conc: a dictionary mapping Compound objects to
concentration in M (default is 1M)
:return: a tuple of (dG_r_prime, dG_uncertainty) where dG_r_prime is
the estimated Gibbs free energy of reaction and dG_uncertainty is the
standard deviation of estimation. Multiply it by 1.96 to get a 95%
confidence interval (which is the value shown on eQuilibrator).
"""
dg_prime, dg_uncertainty = self.standard_dg_prime(reaction)
dg_prime += self.RT * reaction.dg_correction(compound_to_conc)
dg_prime += self.RT * reaction.dg_correction()
return dg_prime, dg_uncertainty
def physiological_dg_prime(
......@@ -195,7 +194,61 @@ class ComponentContribution(object):
dG0_prime, dG0_uncertainty = self.standard_dg_prime(reaction)
E0_prime = -dG0_prime / (n_e*FARADAY) # in Volts
E0_uncertainty = dG0_uncertainty / (n_e*FARADAY) # in Volts
E0_uncertainty = dG0_uncertainty / abs(n_e*FARADAY) # in Volts
return E0_prime, E0_uncertainty
def physiological_e_prime(
self,
reaction: PhasedReaction
) -> Tuple[float, float]:
"""Calculate the E'0 of a single half-reaction.
:param reaction: a PhasedReaction object
:return: a tuple (E0_prime, E0_uncertainty) where E0_prime is
the estimated standard electrostatic potential of reaction and
E0_uncertainty is the standard deviation of estimation. Multiply it
by 1.96 to get a 95% confidence interval (which is the value shown on
eQuilibrator).
"""
n_e = reaction.check_half_reaction_balancing()
if n_e is None:
raise ValueError('reaction is not chemically balanced')
if n_e == 0:
raise ValueError('this is not a half-reaction, '
'electrons are balanced')
dG0_prime, dG0_uncertainty = self.physiological_dg_prime(reaction)
E0_prime = -dG0_prime / (n_e*FARADAY) # in Volts
E0_uncertainty = dG0_uncertainty / abs(n_e*FARADAY) # in Volts
return E0_prime, E0_uncertainty
def e_prime(
self,
reaction: PhasedReaction
) -> Tuple[float, float]:
"""Calculate the E'0 of a single half-reaction.
:param reaction: a PhasedReaction object
:return: a tuple (E0_prime, E0_uncertainty) where E0_prime is
the estimated standard electrostatic potential of reaction and
E0_uncertainty is the standard deviation of estimation. Multiply it
by 1.96 to get a 95% confidence interval (which is the value shown on
eQuilibrator).
"""
n_e = reaction.check_half_reaction_balancing()
if n_e is None:
raise ValueError('reaction is not chemically balanced')
if n_e == 0:
raise ValueError('this is not a half-reaction, '
'electrons are balanced')
dG0_prime, dG0_uncertainty = self.dg_prime(reaction)
E0_prime = -dG0_prime / (n_e*FARADAY) # in Volts
E0_uncertainty = dG0_uncertainty / abs(n_e*FARADAY) # in Volts
return E0_prime, E0_uncertainty
......@@ -208,12 +261,12 @@ class ComponentContribution(object):
:param reaction: a PhasedReaction object.
:return: the analysis results as a list of dictionaries
"""
return self.predictor.get_dg_analysis(reaction)
return ComponentContribution.predictor.get_dg_analysis(reaction)
def IsUsingGroupContributions(self, reaction: PhasedReaction) -> bool:
def is_using_group_contribution(self, reaction: PhasedReaction) -> bool:
"""Check whether group contribution is needed to get this reactions' dG.
:param reaction: a PhasedReaction object.
:return: true iff group contribution is needed
"""
return self.predictor.is_using_group_contribution(reaction)
return ComponentContribution.predictor.is_using_group_contribution(reaction)
This diff is collapsed.
This diff is collapsed.
......@@ -32,7 +32,7 @@ import numpy
import pytest
from equilibrator_api import Q_, ComponentContribution, Reaction, ccache
from equilibrator_api.phased_reaction import GAS_PHASE
from equilibrator_api.phased_compound import GAS_PHASE_NAME
from . import approx_unit
......@@ -49,16 +49,32 @@ def atp_hydrolysis() -> Reaction:
formula = "KEGG:C00002 + KEGG:C00001 = KEGG:C00008 + KEGG:C00009"
return Reaction.parse_formula(formula)
@pytest.fixture(scope="module")
def missing_h2o() -> Reaction:
"""Create a ATP hydrolysis reaction."""
formula = "KEGG:C00002 = KEGG:C00008 + KEGG:C00009"
return Reaction.parse_formula(formula)
@pytest.fixture(scope="module")
def fermentation_gas() -> Reaction:
"""Create a reaction object of glucose => 2 ethanol + 2 CO2(g)."""
formula = "KEGG:C00031 = 2 KEGG:C00469 + 2 KEGG:C00011"
rxn = Reaction.parse_formula(formula)
rxn.change_phase(ccache.get_compound("KEGG:C00011"), GAS_PHASE)
rxn.set_phase(ccache.get_compound("KEGG:C00011"), GAS_PHASE_NAME)
return rxn
def test_water_balancing(comp_contribution, missing_h2o):
assert not missing_h2o.is_balanced()
assert missing_h2o.is_balanced(ignore_atoms=("H", "O", "e-"))
balanced_rxn = missing_h2o.balance_with_compound(
ccache.get_compound('KEGG:C00001'),
ignore_atoms=("H")
)
assert balanced_rxn is not None
def test_atp_hydrolysis_physiological_dg(comp_contribution, atp_hydrolysis):
"""Test the dG adjustments for physiological conditions (with H2O)."""
warnings.simplefilter('ignore', ResourceWarning)
......@@ -66,12 +82,11 @@ def test_atp_hydrolysis_physiological_dg(comp_contribution, atp_hydrolysis):
assert float(atp_hydrolysis.physiological_dg_correction()) == \
pytest.approx(numpy.log(1e-3), rel=1e-3)
compound_to_conc = {ccache.get_compound('KEGG:C00001'): None,
ccache.get_compound('KEGG:C00002'): Q_("1uM"),
ccache.get_compound('KEGG:C00008'): Q_("1uM"),
ccache.get_compound('KEGG:C00009'): Q_("1uM")}
atp_hydrolysis.set_abundance(ccache.get_compound('KEGG:C00002'), Q_("1uM"))
atp_hydrolysis.set_abundance(ccache.get_compound('KEGG:C00008'), Q_("1uM"))
atp_hydrolysis.set_abundance(ccache.get_compound('KEGG:C00009'), Q_("1uM"))
assert float(atp_hydrolysis.dg_correction(compound_to_conc)) == \
assert float(atp_hydrolysis.dg_correction()) == \
pytest.approx(numpy.log(1e-6), rel=1e-3)
......@@ -86,10 +101,9 @@ def test_atp_hydrolysis_dg(comp_contribution, atp_hydrolysis):
"""Test the CC predictions for ATP hydrolysis."""
warnings.simplefilter('ignore', ResourceWarning)
compound_to_conc = {ccache.get_compound('KEGG:C00001'): None,
ccache.get_compound('KEGG:C00002'): Q_("1mM"),
ccache.get_compound('KEGG:C00008'): Q_("10mM"),
ccache.get_compound('KEGG:C00009'): Q_("10mM")}
atp_hydrolysis.set_abundance(ccache.get_compound('KEGG:C00002'), Q_("1mM"))
atp_hydrolysis.set_abundance(ccache.get_compound('KEGG:C00008'), Q_("10mM"))
atp_hydrolysis.set_abundance(ccache.get_compound('KEGG:C00009'), Q_("10mM"))
standard_dg_prime, dg_uncertainty = comp_contribution.standard_dg_prime(
atp_hydrolysis)
......@@ -97,7 +111,7 @@ def test_atp_hydrolysis_dg(comp_contribution, atp_hydrolysis):
approx_unit(standard_dg_prime, Q_("-25.8 kJ/mol"), abs=0.1)
approx_unit(dg_uncertainty, Q_("0.3 kJ/mol"), abs=0.1)
dg_prime, _ = comp_contribution.dg_prime(atp_hydrolysis, compound_to_conc)
dg_prime, _ = comp_contribution.dg_prime(atp_hydrolysis)
approx_unit(dg_prime, Q_("-31.5 kJ/mol"), abs=0.1)
physiological_dg_prime, _ = comp_contribution.physiological_dg_prime(
......@@ -148,3 +162,14 @@ def test_unresolved_reactions(comp_contribution):
_, dg_uncertainty = comp_contribution.standard_dg_prime(reaction)
assert float(dg_uncertainty / Q_("kJ/mol")) > 1e4
def test_reduction_potential(comp_contribution):
"""Test the CC predictions for a redox half-reaction."""
warnings.simplefilter('ignore', ResourceWarning)
# oxaloacetate = malate
formula = "KEGG:C00036 = KEGG:C00149"
reaction = Reaction.parse_formula(formula)
assert reaction.is_balanced(ignore_atoms=("H", "e-"))
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment