Skip to content

Commit

Permalink
WIP
Browse files Browse the repository at this point in the history
  • Loading branch information
@hwpang
hwpang committed Jan 30, 2024
1 parent a34e5f6 commit d0ebc42
Show file tree
Hide file tree
Showing 2 changed files with 78 additions and 62 deletions.
24 changes: 21 additions & 3 deletions arkane/encorr/isodesmic.py
View file
Expand Up @@ -291,6 +291,8 @@ def __eq__(self, other):
f"GenericConstraint object has no __eq__ defined for other object of "
f"type {type(other)}"
)
def __repr__(self):
return self.constraint_str


def bond_centric_constraints(
Expand All @@ -307,8 +309,14 @@ def bond_centric_constraints(


def _buerger_rc2(bond: Bond) -> BondConstraint:
atom1 = AtomConstraint(label=bond.atom1.symbol)
atom2 = AtomConstraint(label=bond.atom2.symbol)
atom1 = bond.atom1
atom2 = bond.atom2

if atom1.symbol > atom2.symbol:
atom1, atom2 = atom2, atom1

atom1 = AtomConstraint(label=atom1.symbol)
atom2 = AtomConstraint(label=atom2.symbol)

return BondConstraint(atom1=atom1, atom2=atom2, bond_order=int(bond.order))

Expand All @@ -317,16 +325,25 @@ def _buerger_rc3(bond: Bond) -> BondConstraint:
atom1 = bond.atom1
atom2 = bond.atom2

if atom1.symbol > atom2.symbol:
atom1, atom2 = atom2, atom1

atom1 = AtomConstraint(label=f"{atom1.symbol}{atom1.connectivity1}")
atom2 = AtomConstraint(label=f"{atom2.symbol}{atom2.connectivity1}")

return BondConstraint(atom1=atom1, atom2=atom2, bond_order=int(bond.order))


def _buerger_rc4(bond: Bond) -> BondConstraint:
atom1 = bond.atom1
atom2 = bond.atom2

if atom1.symbol > atom2.symbol:
atom1, atom2 = atom2, atom1

atoms = []

for atom in [bond.atom1, bond.atom2]:
for atom in [atom1, atom2]:
connections = []
for a, b in atom.bonds.items():
ac = AtomConstraint(label=f"{a.symbol}{a.connectivity1}")
Expand Down Expand Up @@ -412,6 +429,7 @@ def _get_all_constraints(
if self.conserve_ring_size:
features += self._get_ring_constraints(species)

features.sort(key=lambda x: x.__repr__())
return features

def _enumerate_constraints(self, full_constraints_list):
Expand Down
116 changes: 57 additions & 59 deletions test/arkane/encorr/isodesmicTest.py
View file
Expand Up @@ -31,7 +31,6 @@
This script contains unit tests of the :mod:`arkane.isodesmic` module.
"""


import numpy as np

from rmgpy.molecule import Molecule
Expand Down Expand Up @@ -137,12 +136,11 @@ def test_initializing_constraint_map(self):
"""
Test that the constraint map is properly initialized when a SpeciesConstraints object is initialized
"""
caffeine_consts = SpeciesConstraints(self.caffeine, [self.butane, self.benzene])
assert set(caffeine_consts.constraint_map.keys()) == {
"H",
"C",
"O",
"N",
consts = SpeciesConstraints(self.caffeine, [self.butane, self.benzene])
caffeine_features = consts._get_all_constraints(self.caffeine)
caffeine_constraint_list = [feat.__repr__() for feat in caffeine_features]

assert set(caffeine_constraint_list) == {
"C=O",
"C-N",
"C-H",
Expand All @@ -154,55 +152,69 @@ def test_initializing_constraint_map(self):
}

no_rings = SpeciesConstraints(self.caffeine, [self.butane, self.benzene], conserve_ring_size=False)
assert set(no_rings.constraint_map.keys()) == {"H", "C", "O", "N", "C=O", "C-N", "C-H", "C=C", "C=N", "C-C"}

atoms_only = SpeciesConstraints(self.caffeine, [self.butane], conserve_ring_size=False, conserve_bonds=False)
assert set(atoms_only.constraint_map.keys()) == {"H", "C", "O", "N"}
caffeine_features = no_rings._get_all_constraints(self.caffeine)
caffeine_constraint_list = [feat.__repr__() for feat in caffeine_features]
assert set(caffeine_constraint_list) == {"C=O", "C-N", "C-H", "C=C", "C=N", "C-C"}

def test_enumerating_constraints(self):
"""
Test that a SpeciesConstraints object can properly enumerate the constraints of a given ErrorCancelingSpecies
"""
spcs_consts = SpeciesConstraints(self.benzene, [])
assert set(spcs_consts.constraint_map.keys()) == {"C", "H", "C=C", "C-C", "C-H", "6_ring"}

# Now that we have confirmed that the correct keys are present, overwrite the constraint map to set the order
spcs_consts.constraint_map = {
"H": 0,
"C": 1,
"C=C": 2,
"C-C": 3,
"C-H": 4,
"6_ring": 5,
}
benzene_features = spcs_consts._get_all_constraints(self.benzene)
benzene_constraint_list = [feat.__repr__() for feat in benzene_features]
assert set(benzene_constraint_list) == {"C=C", "C-C", "C-H", "6_ring"}

target_constraints, _ = spcs_consts._enumerate_constraints([benzene_features])
benzene_constraints = target_constraints

assert np.array_equal(
spcs_consts._enumerate_constraints(self.propene),
np.array([6, 3, 1, 1, 6, 0]),
benzene_constraints,
np.array([1, 3, 6, 3]),
)

spcs_consts.all_reference_species = [self.propene]
propene_features = spcs_consts._get_all_constraints(self.propene)
_, reference_constraints = spcs_consts._enumerate_constraints([benzene_features, propene_features])
propene_constraints = reference_constraints[0]
assert np.array_equal(
spcs_consts._enumerate_constraints(self.butane),
np.array([10, 4, 0, 3, 10, 0]),
propene_constraints,
np.array([0, 1, 6, 1]),
)

spcs_consts.all_reference_species = [self.butane]
butane_features = spcs_consts._get_all_constraints(self.butane)
_, reference_constraints = spcs_consts._enumerate_constraints([benzene_features, butane_features])
butane_constraints = reference_constraints[0]
assert np.array_equal(
spcs_consts._enumerate_constraints(self.benzene),
np.array([6, 6, 3, 3, 6, 1]),
butane_constraints,
np.array([0, 3, 10, 0]),
)

# Caffeine and ethyne should return None since they have features not found in benzene
assert spcs_consts._enumerate_constraints(self.caffeine) is None
assert spcs_consts._enumerate_constraints(self.ethyne) is None
# Caffeine and ethyne should return empty list since they have features not found in benzene
spcs_consts.all_reference_species = [self.caffeine]
caffeine_features = spcs_consts._get_all_constraints(self.caffeine)
_, reference_constraints = spcs_consts._enumerate_constraints([benzene_features, caffeine_features])
assert len(reference_constraints) == 0

spcs_consts.all_reference_species = [self.ethyne]
ethyne_features = spcs_consts._get_all_constraints(self.ethyne)
_, reference_constraints = spcs_consts._enumerate_constraints([benzene_features, ethyne_features])
assert len(reference_constraints) == 0

def test_calculating_constraints(self):
"""
Test that a SpeciesConstraints object can properly return the target constraint vector and the constraint matrix
"""
spcs_consts = SpeciesConstraints(self.caffeine, [self.propene, self.butane, self.benzene, self.ethyne])
assert set(spcs_consts.constraint_map.keys()) == {
"H",
"C",
"O",
"N",
caffeine_features = spcs_consts._get_all_constraints(self.caffeine)
propene_features = spcs_consts._get_all_constraints(self.propene)
butane_features = spcs_consts._get_all_constraints(self.butane)
benzene_features = spcs_consts._get_all_constraints(self.benzene)
ethyne_features = spcs_consts._get_all_constraints(self.ethyne)

caffeine_feature_list = [feat.__repr__() for feat in caffeine_features]
assert set(caffeine_feature_list) == {
"C=O",
"C-N",
"C-H",
Expand All @@ -213,36 +225,20 @@ def test_calculating_constraints(self):
"6_ring",
}

# Now that we have confirmed that the correct keys are present, overwrite the constraint map to set the order
spcs_consts.constraint_map = {
"H": 0,
"C": 1,
"O": 2,
"N": 3,
"C=O": 4,
"C-N": 5,
"C-H": 6,
"C=C": 7,
"C=N": 8,
"C-C": 9,
"5_ring": 10,
"6_ring": 11,
}

target_consts, consts_matrix = spcs_consts.calculate_constraints()

# First, test that ethyne is not included in the reference set
assert spcs_consts.reference_species == [self.propene, self.butane, self.benzene]

# Then, test the output of the calculation
assert np.array_equal(target_consts, np.array([10, 8, 2, 4, 2, 10, 10, 1, 1, 1, 1, 1]))
assert np.array_equal(target_consts, np.array([ 1, 1, 1, 10, 10, 1, 1, 2, 0, 8, 10, 4, 2]))
assert np.array_equal(
consts_matrix,
np.array(
[
[6, 3, 0, 0, 0, 0, 6, 1, 0, 1, 0, 0],
[10, 4, 0, 0, 0, 0, 10, 0, 0, 3, 0, 0],
[6, 6, 0, 0, 0, 0, 6, 3, 0, 3, 0, 1],
[ 0, 0, 1, 6, 0, 1, 0, 0, 0, 3, 6, 0, 0],
[ 0, 0, 3, 10, 0, 0, 0, 0, 0, 4, 10, 0, 0],
[ 0, 1, 3, 6, 0, 3, 0, 0, 0, 6, 6, 0, 0],
]
),
)
Expand Down Expand Up @@ -278,6 +274,8 @@ def test_creating_error_canceling_schemes(self):
scheme = ErrorCancelingScheme(
self.propene,
[self.butane, self.benzene, self.caffeine, self.ethyne],
"rc2",
True,
True,
True,
)
Expand Down Expand Up @@ -328,7 +326,7 @@ def test_multiple_error_canceling_reactions(self):
)

reaction_list = scheme.multiple_error_canceling_reaction_search(n_reactions_max=20)
assert len(reaction_list) == 20
assert len(reaction_list) == 6
reaction_string = reaction_list.__repr__()
# Consider both permutations of the products in the reaction string
rxn_str1 = "<ErrorCancelingReaction 1*C=CC + 1*CCCC <=> 1*CCC + 1*C=CCC >"
Expand Down Expand Up @@ -361,9 +359,9 @@ def test_calculate_target_enthalpy(self):
)

target_thermo, rxn_list = scheme.calculate_target_enthalpy(n_reactions_max=3, milp_software=["lpsolve"])
assert target_thermo.value_si == 115000.0
assert target_thermo.value_si == 110000.0
assert isinstance(rxn_list[0], ErrorCancelingReaction)

if self.pyo is not None:
target_thermo, _ = scheme.calculate_target_enthalpy(n_reactions_max=3, milp_software=["pyomo"])
assert target_thermo.value_si == 115000.0
assert target_thermo.value_si == 110000.0

0 comments on commit d0ebc42

Please sign in to comment.