Skip to content

Commit

Permalink
Merge pull request #1788 from ReactionMechanismGenerator/terachem
Browse files Browse the repository at this point in the history
Terachem
  • Loading branch information
mjohnson541 authored Nov 22, 2019
2 parents e0a78cf + 1132931 commit bfffc5b
Show file tree
Hide file tree
Showing 55 changed files with 50,290 additions and 103 deletions.
100 changes: 97 additions & 3 deletions arkane/common.py
Original file line number Diff line number Diff line change
Expand Up @@ -41,6 +41,7 @@

import rmgpy.constants as constants
from rmgpy import __version__
from rmgpy.exceptions import InputError
from rmgpy.molecule.element import get_element
from rmgpy.molecule.translator import to_inchi, to_inchi_key
from rmgpy.pdep.collision import SingleExponentialDown
Expand All @@ -57,6 +58,8 @@

from arkane.pdep import PressureDependenceJob

################################################################################


# Add a custom string representer to use block literals for multiline strings
def str_repr(dumper, data):
Expand Down Expand Up @@ -353,9 +356,9 @@ def check_conformer_energy(energies, path):
energies = np.array(energies, np.float64)
e_diff = (energies[0] - np.min(energies)) * constants.E_h * constants.Na / 1000
if e_diff >= 2: # we choose 2 kJ/mol to be the critical energy
logging.warning('the species corresponding to {path} is different in energy from the lowest energy conformer '
'by {diff} kJ/mol. This can cause significant errors in your computed rate constants.'
.format(path=os.path.basename(path), diff=e_diff))
logging.warning(f'The species corresponding to {os.path.basename(path)} is different in energy from the '
f'lowest energy conformer by {e_diff:.2f} kJ/mol. This can cause significant errors in '
f'your computed thermodynamic properties and rate coefficients.')


def get_element_mass(input_element, isotope=None):
Expand Down Expand Up @@ -587,3 +590,94 @@ def get_element_mass(input_element, isotope=None):
'Lv': [[293, 293.20449]],
'Ts': [[292, 292.20746]],
'Og': [[294, 294.21392]]}


def get_center_of_mass(coords, numbers=None, symbols=None):
"""
Calculate and return the 3D position of the center of mass of the current geometry.
Either ``numbers`` or ``symbols`` must be given.
Args:
coords (np.array): Entries are 3-length lists of xyz coordinates for an atom.
numbers (np.array, list): Entries are atomic numbers corresponding to coords.
symbols (list): Entries are atom symbols corresponding to coords.
Returns:
np.array: The center of mass coordinates.
"""
if symbols is None and numbers is None:
raise IndexError('Either symbols or numbers must be given.')
if numbers is not None:
symbols = [symbol_by_number[number] for number in numbers]
center, total_mass = np.zeros(3, np.float64), 0
for coord, symbol in zip(coords, symbols):
mass = get_element_mass(symbol)[0]
center += mass * coord
total_mass += mass
center /= total_mass
return center


def get_moment_of_inertia_tensor(coords, numbers=None, symbols=None):
"""
Calculate and return the moment of inertia tensor for the current
geometry in amu*angstrom^2. If the coordinates are not at the center of mass,
they are temporarily shifted there for the purposes of this calculation.
Adapted from J.W. Allen: https://github.com/jwallen/ChemPy/blob/master/chempy/geometry.py
Args:
coords (np.array): Entries are 3-length lists of xyz coordinates for an atom.
numbers (np.array, list): Entries are atomic numbers corresponding to coords.
symbols (list): Entries are atom symbols corresponding to coords.
Returns:
np.array: The 3x3 moment of inertia tensor.
Raises:
InputError: If neither ``symbols`` nor ``numbers`` are given, or if they have a different length than ``coords``
"""
if symbols is None and numbers is None:
raise InputError('Either symbols or numbers must be given.')
if numbers is not None:
symbols = [symbol_by_number[number] for number in numbers]
if len(coords) != len(symbols):
raise InputError(f'The number of atoms ({len(symbols)}) is not equal to the number of '
f'atomic coordinates ({len(list(coords))})')
tensor = np.zeros((3, 3), np.float64)
center_of_mass = get_center_of_mass(coords=coords, numbers=numbers, symbols=symbols)
for symbol, coord in zip(symbols, coords):
mass = get_element_mass(symbol)[0]
cm_coord = coord - center_of_mass
tensor[0, 0] += mass * (cm_coord[1] * cm_coord[1] + cm_coord[2] * cm_coord[2])
tensor[1, 1] += mass * (cm_coord[0] * cm_coord[0] + cm_coord[2] * cm_coord[2])
tensor[2, 2] += mass * (cm_coord[0] * cm_coord[0] + cm_coord[1] * cm_coord[1])
tensor[0, 1] -= mass * cm_coord[0] * cm_coord[1]
tensor[0, 2] -= mass * cm_coord[0] * cm_coord[2]
tensor[1, 2] -= mass * cm_coord[1] * cm_coord[2]
tensor[1, 0] = tensor[0, 1]
tensor[2, 0] = tensor[0, 2]
tensor[2, 1] = tensor[1, 2]
return tensor


def get_principal_moments_of_inertia(coords, numbers=None, symbols=None):
"""
Calculate and return the principal moments of inertia in amu*angstrom^2 in decending order
and the corresponding principal axes for the current geometry.
The moments of inertia are in translated to the center of mass. The principal axes have unit lengths.
Adapted from J.W. Allen: https://github.com/jwallen/ChemPy/blob/master/chempy/geometry.py
Args:
coords (np.array): Entries are 3-length lists of xyz coordinates for an atom.
numbers (np.array, list): Entries are atomic numbers corresponding to coords.
symbols (list): Entries are atom symbols corresponding to coords.
Returns:
tuple: The principal moments of inertia.
tuple: The corresponding principal axes.
"""
tensor0 = get_moment_of_inertia_tensor(coords=coords, numbers=numbers, symbols=symbols)
# Since tensor0 is real and symmetric, diagonalization is always possible
principal_moments_of_inertia, axes = np.linalg.eig(tensor0)
principal_moments_of_inertia, axes = zip(*sorted(zip(np.ndarray.tolist(principal_moments_of_inertia),
np.ndarray.tolist(axes)), reverse=True))
return principal_moments_of_inertia, axes
101 changes: 98 additions & 3 deletions arkane/commonTest.py
Original file line number Diff line number Diff line change
Expand Up @@ -46,7 +46,8 @@
from rmgpy.thermo import NASA, ThermoData

from arkane import Arkane, input
from arkane.common import ArkaneSpecies, get_element_mass
from arkane.common import ArkaneSpecies, get_element_mass, get_center_of_mass, \
get_moment_of_inertia_tensor, get_principal_moments_of_inertia
from arkane.input import job_list
from arkane.statmech import InputError, StatMechJob

Expand Down Expand Up @@ -443,9 +444,9 @@ def tearDownClass(cls):
shutil.rmtree(item_path)


class TestGetMass(unittest.TestCase):
class TestMomentOfInertia(unittest.TestCase):
"""
Contains unit tests of common.py
Contains unit tests for attaining moments of inertia from the 3D coordinates.
"""

def test_get_mass(self):
Expand All @@ -455,8 +456,102 @@ def test_get_mass(self):
self.assertEquals(get_element_mass('C', 13), (13.00335483507, 6)) # test specific isotope
self.assertEquals(get_element_mass('Bk'), (247.0703073, 97)) # test a two-element array (no isotope data)

def test_get_center_of_mass(self):
"""Test attaining the center of mass"""
symbols = ['C', 'H', 'H', 'H', 'H']
coords = np.array([[0.0000000, 0.0000000, 0.0000000],
[0.6269510, 0.6269510, 0.6269510],
[-0.6269510, -0.6269510, 0.6269510],
[-0.6269510, 0.6269510, -0.6269510],
[0.6269510, -0.6269510, -0.6269510]], np.float64)
center_of_mass = get_center_of_mass(coords=coords, symbols=symbols)
for cm_coord in center_of_mass:
self.assertEqual(cm_coord, 0.0)

symbols = ['O', 'C', 'C', 'H', 'H', 'H', 'H', 'H', 'H']
coords = np.array([[1.28706525, 0.52121353, 0.04219198],
[0.39745682, -0.35265044, -0.63649234],
[0.36441173, -1.68197093, 0.08682400],
[-0.59818222, 0.10068325, -0.65235399],
[0.74799641, -0.48357798, -1.66461710],
[0.03647269, -1.54932006, 1.12314420],
[-0.31340646, -2.38081353, -0.41122551],
[1.36475837, -2.12581592, 0.12433596],
[2.16336803, 0.09985803, 0.03295192]], np.float64)
center_of_mass = get_center_of_mass(coords=coords, symbols=symbols)
self.assertAlmostEqual(center_of_mass[0], 0.7201, 3)
self.assertAlmostEqual(center_of_mass[1], -0.4880, 3)
self.assertAlmostEqual(center_of_mass[2], -0.1603, 3)

numbers = [6, 6, 8, 1, 1, 1, 1, 1, 1]
coords = np.array([[1.1714680, -0.4048940, 0.0000000],
[0.0000000, 0.5602500, 0.0000000],
[-1.1945070, -0.2236470, 0.0000000],
[-1.9428910, 0.3834580, 0.0000000],
[2.1179810, 0.1394450, 0.0000000],
[1.1311780, -1.0413680, 0.8846660],
[1.1311780, -1.0413680, -0.8846660],
[0.0448990, 1.2084390, 0.8852880],
[0.0448990, 1.2084390, -0.8852880]], np.float64)
center_of_mass = get_center_of_mass(coords=coords, numbers=numbers)
self.assertAlmostEqual(center_of_mass[0], -0.0540, 3)
self.assertAlmostEqual(center_of_mass[1], -0.0184, 3)
self.assertAlmostEqual(center_of_mass[2], -0.0000, 3)

def test_get_moment_of_inertia_tensor(self):
"""Test calculating the moment of inertia tensor"""
symbols = ['O', 'C', 'C', 'H', 'H', 'H', 'H', 'H', 'H']
coords = np.array([[1.28706525, 0.52121353, 0.04219198],
[0.39745682, -0.35265044, -0.63649234],
[0.36441173, -1.68197093, 0.08682400],
[-0.59818222, 0.10068325, -0.65235399],
[0.74799641, -0.48357798, -1.66461710],
[0.03647269, -1.54932006, 1.12314420],
[-0.31340646, -2.38081353, -0.41122551],
[1.36475837, -2.12581592, 0.12433596],
[2.16336803, 0.09985803, 0.03295192]], np.float64)
tensor = get_moment_of_inertia_tensor(coords=coords, symbols=symbols)
expected_tensor = [[50.24197604, -15.43600683, -3.07977736],
[-15.43600683, 22.20416597, 2.5935549],
[-3.07977736, 2.5935549, 55.49144794]]
np.testing.assert_almost_equal(tensor, expected_tensor)

def test_get_principal_moments_of_inertia(self):
"""Test calculating the principal moments of inertia"""
numbers = [6, 6, 8, 1, 1, 1, 1, 1, 1]
coords = np.array([[1.235366, -0.257231, -0.106315],
[0.083698, 0.554942, 0.046628],
[-1.210594, -0.239505, -0.021674],
[0.132571, 1.119728, 0.987719],
[0.127795, 1.278999, -0.769346],
[-1.272620, -0.962700, 0.798216],
[-2.074974, 0.426198, 0.055846],
[-1.275744, -0.785745, -0.965493],
[1.241416, -0.911257, 0.593856]], np.float64)
principal_moments_of_inertia = get_principal_moments_of_inertia(coords=coords, numbers=numbers)[0]
expected_principal_moments_of_inertia = [60.98026894, 53.83156297, 14.48858465]
for moment, expected_moment in zip(principal_moments_of_inertia, expected_principal_moments_of_inertia):
self.assertAlmostEqual(moment, expected_moment)

symbols = ['N', 'O', 'O'] # test a linear molecule
coords = np.array([[0.000000, 0.000000, 1.106190],
[0.000000, 0.000000, -0.072434],
[0.000000, 0.000000, -1.191782]], np.float64)
with self.assertRaises(InputError):
get_principal_moments_of_inertia(coords=coords, numbers=numbers)
principal_moments_of_inertia, axes = get_principal_moments_of_inertia(coords=coords, symbols=symbols)
expected_principal_moments_of_inertia = [39.4505153, 39.4505153, 0.0]
expected_axes = [[1, 0, 0], [0, 1, 0], [0, 0, 1]]
for moment, expected_moment in zip(principal_moments_of_inertia, expected_principal_moments_of_inertia):
self.assertAlmostEqual(moment, expected_moment)
for axis, expected_axis in zip(axes, expected_axes):
for entry, expected_entry in zip(axis, expected_axis):
self.assertAlmostEqual(entry, expected_entry)
self.assertIsInstance(principal_moments_of_inertia, tuple)
self.assertIsInstance(axes, tuple)

################################################################################


if __name__ == '__main__':
unittest.main(testRunner=unittest.TextTestRunner(verbosity=2))
File renamed without changes.
File renamed without changes.
File renamed without changes.
File renamed without changes.
File renamed without changes.
File renamed without changes.
File renamed without changes.
File renamed without changes.
File renamed without changes.
File renamed without changes.
File renamed without changes.
File renamed without changes.
File renamed without changes.
File renamed without changes.
8 changes: 8 additions & 0 deletions arkane/data/terachem/ethane_coords.xyz
Original file line number Diff line number Diff line change
@@ -0,0 +1,8 @@
6
coordinates
C 0.66409651 0.00395265 0.07100793
C -0.66409647 -0.00395253 -0.07100790
H 1.24675866 0.88983869 -0.16137840
H 1.19483972 -0.87530680 0.42244414
H -1.19483975 0.87530673 -0.42244421
H -1.24675868 -0.88983873 0.16137844
Loading

0 comments on commit bfffc5b

Please sign in to comment.