Skip to content

Commit

Permalink
Feature: Add non_element_types support for atomtyping (#648)
Browse files Browse the repository at this point in the history 
* Feature: Add non_element_types support for atomtyping

* Minor fixes to imports/ comments

* Address review comments, couple of more asserts

* WIP- few more asserts

* WIP- Empty branch trigger for GHA

* Fix test case
  • Loading branch information
@umesh-timalsina
umesh-timalsina authored Apr 28, 2022
1 parent 822247d commit 1818e46
Show file tree
Hide file tree
Showing 3 changed files with 156 additions and 0 deletions.
16 changes: 16 additions & 0 deletions gmso/core/forcefield.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -5,6 +5,7 @@

from lxml import etree

from gmso.core.element import element_by_symbol
from gmso.exceptions import MissingPotentialError
from gmso.utils._constants import FF_TOKENS_SEPARATOR
from gmso.utils.ff_utils import (
Expand Down Expand Up @@ -106,6 +107,21 @@ def __init__(self, xml_loc=None, strict=True, greedy=True):
self.combining_rule = "geometric"
self.units = {}

@property
def non_element_types(self):
"""Get the non-element types in the ForceField."""
non_element_types = set()

for name, atom_type in self.atom_types.items():
element_symbol = atom_type.get_tag(
"element"
) # FixMe: Should we make this a first class citizen?
if element_symbol:
element = element_by_symbol(element_symbol)
non_element_types.add(element_symbol) if not element else None

return non_element_types

@property
def atom_class_groups(self):
"""Return a dictionary of atomClasses in the Forcefield."""
Expand Down
104 changes: 104 additions & 0 deletions gmso/tests/files/non-element-type-ff.xml
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,104 @@
<?xml version='1.0' encoding='UTF-8'?>
<ForceField name="Trappe_CH3_CH2s" version="0.0.1">
<FFMetaData electrostatics14Scale="0" nonBonded14Scale="0">
<Units energy="kJ/mol" mass="amu" charge="elementary_charge" distance="nm"/>
</FFMetaData>
<AtomTypes expression="4 * epsilon * ((sigma/r)**12 - (sigma/r)**6)">
<ParametersUnitDef parameter="epsilon" unit="kJ/mol"/>
<ParametersUnitDef parameter="sigma" unit="nm"/>
<AtomType name="CH3_sp3" atomclass="CH3" element="_CH3" charge="0.0" mass="15.03500" definition="[_CH3;X1][_CH3,_CH2]" doi="10.1021/jp972543+">
<Parameters>
<Parameter name="epsilon" value="0.814817"/>
<Parameter name="sigma" value="0.375"/>
</Parameters>
</AtomType>
<AtomType name="CH2_sp3" atomclass="CH2" element="_CH2" charge="0.0" mass="14.02700" definition="[_CH2;X2]([_CH3,_CH2])[_CH3,_CH2]" doi="10.1021/jp972543+">
<Parameters>
<Parameter name="epsilon" value="0.382465"/>
<Parameter name="sigma" value="0.395"/>
</Parameters>
</AtomType>
</AtomTypes>
<BondTypes expression="k / 2 * (r-r_eq)**2">
<ParametersUnitDef parameter="k" unit="kJ/mol/nm**2"/>
<ParametersUnitDef parameter="r_eq" unit="nm"/>
<BondType name="BondType_1" class1="CH3" class2="CH3">
<Parameters>
<Parameter name="k" value="502416.0"/>
<Parameter name="r_eq" value="0.1540"/>
</Parameters>
</BondType>
<BondType name="BondType_2" class1="CH3" class2="CH2">
<Parameters>
<Parameter name="k" value="502416.0"/>
<Parameter name="r_eq" value="0.1540"/>
</Parameters>
</BondType>
<BondType name="BondType_3" class1="CH2" class2="CH2">
<Parameters>
<Parameter name="k" value="502416.0"/>
<Parameter name="r_eq" value="0.1540"/>
</Parameters>
</BondType>
</BondTypes>
<AngleTypes expression="k / 2 * (theta - theta_eq)**2">
<ParametersUnitDef parameter="k" unit="kJ/mol/rad**2"/>
<ParametersUnitDef parameter="theta_eq" unit="rad"/>
<AngleType name="AngleType_1" class1="CH3" class2="CH2" class3="CH3">
<Parameters>
<Parameter name="k" value="519.65389"/>
<Parameter name="theta_eq" value="1.98967"/>
</Parameters>
</AngleType>
<AngleType name="AngleType_2" class1="CH3" class2="CH2" class3="CH2">
<Parameters>
<Parameter name="k" value="519.65389"/>
<Parameter name="theta_eq" value="1.98967"/>
</Parameters>
</AngleType>
<AngleType name="AngleType_3" class1="CH2" class2="CH2" class3="CH2">
<Parameters>
<Parameter name="k" value="519.65389"/>
<Parameter name="theta_eq" value="1.98967"/>
</Parameters>
</AngleType>
</AngleTypes>
<DihedralTypes expression="c0 + c1 * cos(phi)**1 + c2 * cos(phi)**2 + c3 * cos(phi)**3 + c4 * cos(phi)**4 + c5 * cos(phi)**5">
<ParametersUnitDef parameter="c5" unit="kJ"/>
<ParametersUnitDef parameter="c4" unit="kJ"/>
<ParametersUnitDef parameter="c3" unit="kJ"/>
<ParametersUnitDef parameter="c2" unit="kJ"/>
<ParametersUnitDef parameter="c1" unit="kJ"/>
<ParametersUnitDef parameter="c0" unit="kJ"/>
<DihedralType name="Dihedral_1" class1="CH3" class2="CH2" class3="CH2" class4="CH3">
<Parameters>
<Parameter name="c0" value="8.39736"/>
<Parameter name="c1" value="16.78632"/>
<Parameter name="c2" value="1.13393"/>
<Parameter name="c3" value="-26.31760"/>
<Parameter name="c4" value="0.0"/>
<Parameter name="c5" value="0.0"/>
</Parameters>
</DihedralType>
<DihedralType name="Dihedral_2" class1="CH3" class2="CH2" class3="CH2" class4="CH2">
<Parameters>
<Parameter name="c0" value="8.39736"/>
<Parameter name="c1" value="16.78632"/>
<Parameter name="c2" value="1.13393"/>
<Parameter name="c3" value="-26.31760"/>
<Parameter name="c4" value="0.0"/>
<Parameter name="c5" value="0.0"/>
</Parameters>
</DihedralType>
<DihedralType name="Dihedral_3" class1="CH2" class2="CH2" class3="CH2" class4="CH2">
<Parameters>
<Parameter name="c0" value="8.39736"/>
<Parameter name="c1" value="16.78632"/>
<Parameter name="c2" value="1.13393"/>
<Parameter name="c3" value="-26.31760"/>
<Parameter name="c4" value="0.0"/>
<Parameter name="c5" value="0.0"/>
</Parameters>
</DihedralType>
</DihedralTypes>
</ForceField>
36 changes: 36 additions & 0 deletions gmso/tests/test_forcefield.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -29,6 +29,10 @@ def opls_ethane_foyer(self):
get_path(filename=get_path("oplsaa-ethane_foyer.xml"))
)

@pytest.fixture(scope="session")
def non_element_ff(self):
return ForceField(get_path(filename="non-element-type-ff.xml"))

def test_ff_name_version_from_xml(self, ff):
assert ff.name == "ForceFieldOne"
assert ff.version == "0.4.1"
Expand Down Expand Up @@ -560,3 +564,35 @@ def test_get_improper_type_missing(self, opls_ethane_foyer):
opls_ethane_foyer._get_improper_type(
["opls_359", "opls_600", "opls_700", "opls_800"], warn=True
)

def test_non_element_types(self, non_element_ff, opls_ethane_foyer):
assert "_CH3" in non_element_ff.non_element_types
assert "_CH2" in non_element_ff.non_element_types
assert opls_ethane_foyer.non_element_types == set()
assert len(opls_ethane_foyer.atom_types) > 0

assert (
non_element_ff.get_potential(
group="atom_type", key="CH2_sp3"
).charge
== 0
)
assert (
non_element_ff.get_potential(
group="atom_type", key="CH3_sp3"
).charge
== 0
)

assert (
non_element_ff.get_potential(
group="atom_type", key="CH3_sp3"
).definition
== "[_CH3;X1][_CH3,_CH2]"
)
assert (
non_element_ff.get_potential(
group="atom_type", key="CH2_sp3"
).definition
== "[_CH2;X2]([_CH3,_CH2])[_CH3,_CH2]"
)

0 comments on commit 1818e46

Please sign in to comment.