Skip to content

Commit

Permalink
Merge pull request #65 from openforcefield/opc_water
Browse files Browse the repository at this point in the history
Add OPC water model
  • Loading branch information
j-wags authored May 3, 2023
2 parents adde054 + f9e8e71 commit 18cfdc3
Show file tree
Hide file tree
Showing 11 changed files with 580 additions and 6 deletions.
5 changes: 3 additions & 2 deletions .github/workflows/CI.yml
Original file line number Diff line number Diff line change
Expand Up @@ -77,8 +77,9 @@ jobs:
- name: Regenerate water models
if: ${{ matrix.toolkit == true }}
run:
ls scripts/*py | xargs -n 1 -P 3 python
run: |
python scripts/get_ion_nb_params_from_ambertools_frcmod.py $AMBERHOME/dat/leap/parm/frcmod.ionslm_126_opc
ls scripts/write_*.py | xargs -n 1 -P 3 python
- name: Error out if any OFFXML files changed
if: ${{ matrix.toolkit == true }}
Expand Down
1 change: 1 addition & 0 deletions README.md
Original file line number Diff line number Diff line change
Expand Up @@ -26,6 +26,7 @@ Details for each force field in this repository can be found in the following ta

| Filename | DOI | FF line | Release Date | Major format changes? |
|---------------------------------------------|-----------------------------------------------------------------------------------------------------------|---------|---------------|-----------------------|
| `opc-1.0.0.offxml` | see `docs/water-models.md` | Ports | May 1, 2023 | No |
| `openff-2.1.0-rc.1.offxml` | | Sage | Apr 10, 2023 | No |
| `openff_unconstrained-2.1.0-rc.1.offxml` | | Sage | Apr 10, 2023 | No |
| `tip3p_fb-1.0.0.offxml` | see `docs/water-models.md` | Ports | Feb 27, 2023 | No |
Expand Down
2 changes: 2 additions & 0 deletions devtools/conda-envs/test_env.yaml
Original file line number Diff line number Diff line change
Expand Up @@ -7,3 +7,5 @@ dependencies:
- openff-toolkit
- pytest
- pytest-cov
- pandas
- click
3 changes: 3 additions & 0 deletions docs/water-models.md
Original file line number Diff line number Diff line change
Expand Up @@ -4,9 +4,12 @@ The following water models are distributed with this package:

* [TIP3P](https://doi.org/10.1063/1.445869)
* [TIP3P-FB](https://doi.org/10.1021/jz500737m)
* [OPC](https://doi.org/10.1021/jz501780a)

The TIP3P file includes monovalent ion parameters from [Juong and Cheatham](https://dx.doi.org/10.1021/jp8001614).

The OPC file includes the 12-6 hydration free energy (HFE) monovalent ion parameters from [Sengupta et al. 2021](https://doi.org/10.1021/acs.jcim.0c01390), the 12-6 compromise (CM) divalent ion parameters from [Li et al. 2020](https://doi.org/10.1021/acs.jctc.0c00194), and the 12-6 ion-oxygen distance (IOD) trivalent and tetravalent ion parameters from [Li et al. 2021](https://doi.org/10.1021/acs.jctc.0c01320).

These files are not original force fields made by Open Force Field.
Instead they are ports of existing force fields into the [SMIRNOFF](https://openforcefield.github.io/standards/standards/smirnoff/) force field format.

Expand Down
1 change: 1 addition & 0 deletions openforcefields/data/README.md
Original file line number Diff line number Diff line change
@@ -1 +1,2 @@
* `jc.csv`: Data taken from Table 5 in https://pubs.acs.org/doi/10.1021/jp8001614
* `ionslm_126_opc.csv` taken from $AMBERHOME/dat/leap/parm/frcmod.ionslm_126_opc in ambertools-22.0-py310h206695f
61 changes: 61 additions & 0 deletions openforcefields/data/ionslm_126_opc.csv
Original file line number Diff line number Diff line change
@@ -0,0 +1,61 @@
element,opc_rmin2 (A),opc_eps (kcal/mol)
Li+,1.242,0.00216058
Na+,1.467,0.02960343
K+,1.702,0.13953816
Rb+,1.818,0.2279546
Cs+,1.96,0.35308749
Tl+,1.66,0.11250137
Cu+,1.18,0.00078213
Ag+,1.316,0.00602547
F-,1.84,0.24650465
Cl-,2.36,0.6787887
Br-,2.499,0.75971103
I-,2.9,0.90300541
Be2+,0.921,1.28e-06
Cu2+,1.178,0.0007549
Ni2+,1.197,0.00104974
Pt2+,1.219,0.00150903
Zn2+,1.219,0.00150903
Co2+,1.263,0.00294683
Pd2+,1.269,0.00321068
Ag2+,1.305,0.00523385
Cr2+,1.316,0.00602547
Fe2+,1.323,0.00657749
Mg2+,1.33,0.0071693
V2+,1.334,0.00752608
Mn2+,1.381,0.0128646
Hg2+,1.381,0.0128646
Cd2+,1.389,0.01400886
Yb2+,1.602,0.08034231
Ca2+,1.608,0.08337961
Sn2+,1.625,0.09235154
Pb2+,1.707,0.14295367
Eu2+,1.753,0.17618319
Sr2+,1.766,0.18612361
Sm2+,1.766,0.18612361
Ba2+,1.96,0.35308749
Ra2+,1.96,0.35308749
Al3+,1.25,0.00243637
Fe3+,1.4,0.01570749
Cr3+,1.336,0.00770969
In3+,1.44,0.02322071
Tl3+,1.483,0.03393126
Y3+,1.575,0.06751391
La3+,1.7,0.13818331
Ce3+,1.725,0.15557763
Pr3+,1.717,0.14990448
Nd3+,1.662,0.11371963
Sm3+,1.638,0.09957472
Eu3+,1.646,0.10417397
Gd3+,1.6,0.07934493
Tb3+,1.608,0.08337961
Dy3+,1.583,0.07117158
Er3+,1.575,0.06751391
Tm3+,1.575,0.06751391
Lu3+,1.558,0.06014121
Hf4+,1.461,0.02808726
Zr4+,1.488,0.03537062
Ce4+,1.683,0.12693448
U4+,1.683,0.12693448
Pu4+,1.658,0.11129023
Th4+,1.704,0.14089951
139 changes: 139 additions & 0 deletions openforcefields/offxml/opc-1.0.0.offxml

Large diffs are not rendered by default.

139 changes: 139 additions & 0 deletions openforcefields/offxml/opc.offxml

Large diffs are not rendered by default.

36 changes: 32 additions & 4 deletions openforcefields/tests/test_water_models.py
Original file line number Diff line number Diff line change
Expand Up @@ -91,6 +91,34 @@ def test_tip3p_fb(water_molecule):
},
)

# TODO: The interchange that results from applying the OPC water model yields a sigma that is exactly
# twice what it should be on the hydrogens. This problem is only cosmetic as the epsilon is 0.
def test_opc(water_molecule):
from openmm.app import Modeller
omm_water = water_molecule.to_openmm()
omm_ff = OpenMMForceField("opc.xml")
mod = Modeller(omm_water, water_molecule.get_positions().to_openmm())
mod.addExtraParticles(omm_ff)
reference = omm_ff.createSystem(
mod.getTopology(),
constraints=HAngles,
rigidWater=True,
)

interchange = ForceField("opc-1.0.0.offxml").create_interchange(
water_molecule,
)
system = interchange.to_openmm()
compare_water_systems(
reference,
system,
{
"charge": 1e-10 * openmm.unit.elementary_charge,
#"sigma": 1e-10 * openmm.unit.nanometer,
"sigma": .1 * openmm.unit.nanometer,
"epsilon": 1e-10 * openmm.unit.kilojoule_per_mole,
},
)

@pytest.mark.skip(reason="Skipping in first pass")
def test_tip4p_ew(water_molecule):
Expand Down Expand Up @@ -132,7 +160,7 @@ def test_tip5p(water_molecule):

@pytest.mark.parametrize(
"water_model,pattern",
[("tip3p", "^tip3p(?!.*fb)"), ("tip3p_fb", "^tip3p_fb")],
[("tip3p", "^tip3p(?!.*fb)"), ("tip3p_fb", "^tip3p_fb"), ("opc", "^opc")],
)
def test_most_recent_version_match(water_model, pattern):
import re
Expand All @@ -151,8 +179,7 @@ def test_most_recent_version_match(water_model, pattern):
assert len(matched_files) > 0, f"Failed to match any files for pattern {pattern}!"

for file in matched_files:
split = re.split(pattern + "-", file.strip(".offxml"))

split = re.split(pattern + "-", file.replace(".offxml", ""))
if len(split) == 2:
found_verison = version.Version(split[1])
if found_verison > maximum_version:
Expand Down Expand Up @@ -247,4 +274,5 @@ def test_water_model_is_compatible_with_mainline():
# Since we don't have a way to get the most recent mainline FF, be sure
# to occasionally update the first FF listed here
ForceField('openff-2.0.0.offxml', 'tip3p.offxml')
ForceField('openff-2.0.0.offxml', 'tip3p_fb.offxml')
ForceField('openff-2.0.0.offxml', 'tip3p_fb.offxml')
ForceField('openff-2.0.0.offxml', 'opc.offxml')
44 changes: 44 additions & 0 deletions scripts/get_ion_nb_params_from_ambertools_frcmod.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,44 @@
import click
import pandas
from pathlib import Path

@click.command()
@click.argument(
"amber_frcmod_path",
type = click.Path(),
)
def main(amber_frcmod_path):

ion_nb_params = list()

nonbonded = False

with open(Path(amber_frcmod_path), "r") as amber_frcmod_file:
for line in amber_frcmod_file:

if line == '\n':
nonbonded = False

if nonbonded:

fields = line.split()

ion_nb_params.append({
"element": fields[0],
"opc_rmin2 (A)": float(fields[1]),
"opc_eps (kcal/mol)": float(fields[2]),
})

if line.startswith('NONBON'):
nonbonded = True

output_path = Path(
"openforcefields", "data", f"{Path(amber_frcmod_path).suffix[1:]}.csv"
)
pandas.DataFrame(ion_nb_params).set_index("element").to_csv(output_path)


if __name__ == "__main__":
main()


155 changes: 155 additions & 0 deletions scripts/write_opc.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,155 @@
from openff.toolkit.typing.engines.smirnoff.forcefield import ForceField
from openff.toolkit.typing.engines.smirnoff.parameters import (
ConstraintHandler,
ElectrostaticsHandler,
LibraryChargeHandler,
VirtualSiteHandler,
vdWHandler,
)
from openff.units import unit
from openff.units.elements import SYMBOLS
from packaging import version
import pandas
from pathlib import Path


# Data taken from ambertools-22.0-py310h206695f
# Water Lennard-Jones parameters and geometry from $AMBERHOME/dat/leap/parm/frcmod.opc
# Water hydrogen charge from $AMBERHOME/dat/leap/lib/solvents.lib
# Ion Lennard-Jones parameters from $AMBERHOME/dat/leap/parm/frcmod.ionslm_126_opc

VERSION = version.Version("1.0.0")
OFFXML_PATH = Path("openforcefields", "offxml")

ion_nb_params_df = pandas.read_csv("openforcefields/data/ionslm_126_opc.csv")

opc = ForceField()

opc_electrostatics = ElectrostaticsHandler(version=0.4, scale14=0.8333333333)
opc_library = LibraryChargeHandler(version=0.3)
opc_vdw = vdWHandler(version=0.3)
opc_constraints = ConstraintHandler(version=0.3)
opc_virtual_sites = VirtualSiteHandler(version=0.3)

opc_vdw.add_parameter(
{
"smirks": "[#1]-[#8X2H2+0:1]-[#1]",
"epsilon": unit.Quantity(0.2128008130, unit.kilocalorie_per_mole),
"rmin_half": unit.Quantity(1.777167268, unit.angstrom),
"id": "n-opc-O",
}
)
opc_vdw.add_parameter(
{
"smirks": "[#1:1]-[#8X2H2+0]-[#1]",
"epsilon": unit.Quantity(0.0, unit.kilocalorie_per_mole),
"rmin_half": unit.Quantity(1.0, unit.angstrom),
"id": "n-opc-H",
}
)

opc_library.add_parameter(
{
"smirks": "[#1]-[#8X2H2+0:1]-[#1]",
"charge1": unit.Quantity(0.0, unit.elementary_charge),
"id": "q-opc-O",
}
)
opc_library.add_parameter(
{
"smirks": "[#1:1]-[#8X2H2+0]-[#1]",
"charge1": unit.Quantity(0.0, unit.elementary_charge),
"id": "q-opc-H",
}
)

opc_virtual_sites.add_parameter(
{
"type": "DivalentLonePair",
"smirks": "[#1:2]-[#8X2H2+0:1]-[#1:3]",
"match": "once",
"name": "EP",
"distance": unit.Quantity(-0.15939833, unit.angstrom),
"rmin_half": unit.Quantity(1.0, unit.angstrom),
"epsilon": unit.Quantity(0.0, unit.kilocalories_per_mole),
"outOfPlaneAngle": unit.Quantity(0.0, unit.degree),
"charge_increment1": unit.Quantity(0.0, unit.elementary_charge),
"charge_increment2": unit.Quantity(0.679142, unit.elementary_charge),
"charge_increment3": unit.Quantity(0.679142, unit.elementary_charge),
}
)

opc_constraints.add_parameter(
{
"smirks": "[#1:1]-[#8X2H2+0:2]-[#1]",
"id": "c-opc-H-O",
"distance": unit.Quantity(0.87243313 , unit.angstrom),
}
)
opc_constraints.add_parameter(
{
"smirks": "[#1:1]-[#8X2H2+0]-[#1:2]",
"id": "c-opc-H-O-H",
"distance": unit.Quantity(1.37120510, unit.angstrom),
}
)

# Construct dict of element symbol to atomic number
_SYMBOLS_TO_ATOMIC_NUMBER = {
symbol: atomic_number for atomic_number, symbol in SYMBOLS.items()
}

for _, row in ion_nb_params_df.iterrows():

ion_name = row["element"]

# Handle variable length element names
if any(charge_str in ion_name for charge_str in ["2+", "3+", "4+"]):

atomic_number = _SYMBOLS_TO_ATOMIC_NUMBER[ion_name[:-2]]
charge_sign = "+"
charge_magnitude = ion_name[-2]

else:

atomic_number = _SYMBOLS_TO_ATOMIC_NUMBER[ion_name[:-1]]
charge_sign = ion_name[-1]
charge_magnitude = "1"

smirks = f"[#{atomic_number}X0{charge_sign}{charge_magnitude}:1]"

opc_vdw.add_parameter(
{
"smirks": smirks,
"rmin_half": unit.Quantity(row["opc_rmin2 (A)"], unit.angstrom),
"epsilon": unit.Quantity(
row["opc_eps (kcal/mol)"],
unit.kilocalorie_per_mole,
),
"id": f"n-ionslm-126-opc-{ion_name}"
}
)

opc_library.add_parameter(
{
"smirks": smirks,
"charge1": unit.Quantity(
int(f"{charge_sign}{charge_magnitude}"), unit.elementary_charge
),
"id": f"q-ionslm-126-opc-{ion_name}"
}
)

for handler in [
opc_vdw,
opc_library,
opc_electrostatics,
opc_constraints,
opc_virtual_sites,
]:

opc.register_parameter_handler(handler)

opc.to_file(Path(OFFXML_PATH, "opc.offxml"))
opc.to_file(Path(OFFXML_PATH, f"opc-{VERSION}.offxml"))

0 comments on commit 18cfdc3

Please sign in to comment.