Skip to content

Commit 6aa0398

Browse files
KEHANGconnie
KEHANG
authored and
connie
committed
Reduced the unnecessary calls to calculateSymmetryNumber when calculating radical thermo.
calculateSymmetryNumber for saturated structure was called too frequently when calculating entropy of radical.
1 parent ae84b9b commit 6aa0398

File tree

1 file changed

+79
-59
lines changed

1 file changed

+79
-59
lines changed

‎rmgpy/data/thermo.py‎

Lines changed: 79 additions & 59 deletions
Original file line numberDiff line numberDiff line change
@@ -821,9 +821,6 @@ def estimateRadicalThermoViaHBI(self, molecule, stableThermoEstimator ):
821821
return None
822822
assert thermoData is not None, "Thermo data of saturated {0} of molecule {1} is None!".format(saturatedStruct, molecule)
823823

824-
# Undo symmetry number correction for saturated structure
825-
saturatedStruct.calculateSymmetryNumber()
826-
thermoData.S298.value_si += constants.R * math.log(saturatedStruct.symmetryNumber)
827824
# Correct entropy for symmetry number of radical structure
828825
molecule.calculateSymmetryNumber()
829826
thermoData.S298.value_si -= constants.R * math.log(molecule.symmetryNumber)
@@ -878,71 +875,94 @@ def estimateThermoViaGroupAdditivity(self, molecule):
878875
)
879876

880877
if molecule.getRadicalCount() > 0: # radical species
881-
return self.estimateRadicalThermoViaHBI(molecule, self.estimateThermoViaGroupAdditivity )
878+
return self.estimateRadicalThermoViaHBI(molecule, self.estimateThermoViaGroupAdditivityForSaturatedStruct)
882879

883880
else: # non-radical species
884-
cyclic = molecule.isCyclic()
885-
# Generate estimate of thermodynamics
886-
for atom in molecule.atoms:
887-
# Iterate over heavy (non-hydrogen) atoms
888-
if atom.isNonHydrogen():
889-
# Get initial thermo estimate from main group database
890-
try:
891-
self.__addGroupThermoData(thermoData, self.groups['group'], molecule, {'*':atom})
892-
except KeyError:
893-
logging.error("Couldn't find in main thermo database:")
894-
logging.error(molecule)
895-
logging.error(molecule.toAdjacencyList())
896-
raise
897-
# Correct for gauche and 1,5- interactions
898-
if not cyclic:
899-
try:
900-
self.__addGroupThermoData(thermoData, self.groups['gauche'], molecule, {'*':atom})
901-
except KeyError: pass
902-
try:
903-
self.__addGroupThermoData(thermoData, self.groups['int15'], molecule, {'*':atom})
904-
except KeyError: pass
905-
try:
906-
self.__addGroupThermoData(thermoData, self.groups['other'], molecule, {'*':atom})
907-
except KeyError: pass
908-
909-
# Do ring corrections separately because we only want to match
910-
# each ring one time
911-
912-
if cyclic:
913-
if molecule.getAllPolycyclicVertices():
914-
# If the molecule has fused ring atoms, this implies that we are dealing
915-
# with a polycyclic ring system, for which separate ring strain corrections may not
916-
# be adequate. Therefore, we search the polycyclic thermo group corrections
917-
# instead of adding single ring strain corrections within the molecule.
918-
# For now, assume only one polycyclic RSC can be found per molecule
919-
try:
920-
self.__addGroupThermoData(thermoData, self.groups['polycyclic'], molecule, {})
921-
except:
922-
logging.error("Couldn't find in polycyclic ring database:")
923-
logging.error(molecule)
924-
logging.error(molecule.toAdjacencyList())
925-
raise
926-
else:
927-
rings = molecule.getSmallestSetOfSmallestRings()
928-
for ring in rings:
929-
# Make a temporary structure containing only the atoms in the ring
930-
# NB. if any of the ring corrections depend on ligands not in the ring, they will not be found!
931-
try:
932-
self.__addGroupThermoData(thermoData, self.groups['ring'], molecule, {})
933-
except KeyError:
934-
logging.error("Couldn't find in ring database:")
935-
logging.error(ring)
936-
logging.error(ring.toAdjacencyList())
937-
raise
938-
881+
thermoData = self.estimateThermoViaGroupAdditivityForSaturatedStruct(molecule)
939882

940883
# Correct entropy for symmetry number
941884
molecule.calculateSymmetryNumber()
942885
thermoData.S298.value_si -= constants.R * math.log(molecule.symmetryNumber)
943886

944887
return thermoData
945888

889+
def estimateThermoViaGroupAdditivityForSaturatedStruct(self, molecule):
890+
"""
891+
Return the set of thermodynamic parameters corresponding to a given
892+
:class:`Molecule` object `molecule` by estimation using the group
893+
additivity values. If no group additivity values are loaded, a
894+
:class:`DatabaseError` is raised.
895+
"""
896+
# For thermo estimation we need the atoms to already be sorted because we
897+
# iterate over them; if the order changes during the iteration then we
898+
# will probably not visit the right atoms, and so will get the thermo wrong
899+
molecule.sortVertices()
900+
901+
# Create the ThermoData object
902+
thermoData = ThermoData(
903+
Tdata = ([300,400,500,600,800,1000,1500],"K"),
904+
Cpdata = ([0.0,0.0,0.0,0.0,0.0,0.0,0.0],"J/(mol*K)"),
905+
H298 = (0.0,"kJ/mol"),
906+
S298 = (0.0,"J/(mol*K)"),
907+
)
908+
909+
cyclic = molecule.isCyclic()
910+
# Generate estimate of thermodynamics
911+
for atom in molecule.atoms:
912+
# Iterate over heavy (non-hydrogen) atoms
913+
if atom.isNonHydrogen():
914+
# Get initial thermo estimate from main group database
915+
try:
916+
self.__addGroupThermoData(thermoData, self.groups['group'], molecule, {'*':atom})
917+
except KeyError:
918+
logging.error("Couldn't find in main thermo database:")
919+
logging.error(molecule)
920+
logging.error(molecule.toAdjacencyList())
921+
raise
922+
# Correct for gauche and 1,5- interactions
923+
if not cyclic:
924+
try:
925+
self.__addGroupThermoData(thermoData, self.groups['gauche'], molecule, {'*':atom})
926+
except KeyError: pass
927+
try:
928+
self.__addGroupThermoData(thermoData, self.groups['int15'], molecule, {'*':atom})
929+
except KeyError: pass
930+
try:
931+
self.__addGroupThermoData(thermoData, self.groups['other'], molecule, {'*':atom})
932+
except KeyError: pass
933+
934+
# Do ring corrections separately because we only want to match
935+
# each ring one time
936+
937+
if cyclic:
938+
if molecule.getAllPolycyclicVertices():
939+
# If the molecule has fused ring atoms, this implies that we are dealing
940+
# with a polycyclic ring system, for which separate ring strain corrections may not
941+
# be adequate. Therefore, we search the polycyclic thermo group corrections
942+
# instead of adding single ring strain corrections within the molecule.
943+
# For now, assume only one polycyclic RSC can be found per molecule
944+
try:
945+
self.__addGroupThermoData(thermoData, self.groups['polycyclic'], molecule, {})
946+
except:
947+
logging.error("Couldn't find in polycyclic ring database:")
948+
logging.error(molecule)
949+
logging.error(molecule.toAdjacencyList())
950+
raise
951+
else:
952+
rings = molecule.getSmallestSetOfSmallestRings()
953+
for ring in rings:
954+
# Make a temporary structure containing only the atoms in the ring
955+
# NB. if any of the ring corrections depend on ligands not in the ring, they will not be found!
956+
try:
957+
self.__addGroupThermoData(thermoData, self.groups['ring'], molecule, {})
958+
except KeyError:
959+
logging.error("Couldn't find in ring database:")
960+
logging.error(ring)
961+
logging.error(ring.toAdjacencyList())
962+
raise
963+
964+
return thermoData
965+
946966
def __addThermoData(self, thermoData1, thermoData2):
947967
"""
948968
Add the thermodynamic data `thermoData2` to the data `thermoData1`,

0 commit comments

Comments
 (0)