diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml
index d1d3080..b8ac249 100644
--- a/.github/workflows/CI.yml
+++ b/.github/workflows/CI.yml
@@ -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 }}
diff --git a/README.md b/README.md
index 02e00c0..dadf3e2 100644
--- a/README.md
+++ b/README.md
@@ -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 |
diff --git a/devtools/conda-envs/test_env.yaml b/devtools/conda-envs/test_env.yaml
index b383cf7..ae1a89a 100644
--- a/devtools/conda-envs/test_env.yaml
+++ b/devtools/conda-envs/test_env.yaml
@@ -7,3 +7,5 @@ dependencies:
- openff-toolkit
- pytest
- pytest-cov
+ - pandas
+ - click
diff --git a/docs/water-models.md b/docs/water-models.md
index 99fd3cc..b2b31c9 100644
--- a/docs/water-models.md
+++ b/docs/water-models.md
@@ -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.
diff --git a/openforcefields/data/README.md b/openforcefields/data/README.md
index d8f0155..59f899b 100644
--- a/openforcefields/data/README.md
+++ b/openforcefields/data/README.md
@@ -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
diff --git a/openforcefields/data/ionslm_126_opc.csv b/openforcefields/data/ionslm_126_opc.csv
new file mode 100644
index 0000000..5f2a33a
--- /dev/null
+++ b/openforcefields/data/ionslm_126_opc.csv
@@ -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
diff --git a/openforcefields/offxml/opc-1.0.0.offxml b/openforcefields/offxml/opc-1.0.0.offxml
new file mode 100644
index 0000000..ebd2fdf
--- /dev/null
+++ b/openforcefields/offxml/opc-1.0.0.offxml
@@ -0,0 +1,139 @@
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
\ No newline at end of file
diff --git a/openforcefields/offxml/opc.offxml b/openforcefields/offxml/opc.offxml
new file mode 100644
index 0000000..ebd2fdf
--- /dev/null
+++ b/openforcefields/offxml/opc.offxml
@@ -0,0 +1,139 @@
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
\ No newline at end of file
diff --git a/openforcefields/tests/test_water_models.py b/openforcefields/tests/test_water_models.py
index 26df48e..d9db2ad 100644
--- a/openforcefields/tests/test_water_models.py
+++ b/openforcefields/tests/test_water_models.py
@@ -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):
@@ -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
@@ -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:
@@ -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')
\ No newline at end of file
+ ForceField('openff-2.0.0.offxml', 'tip3p_fb.offxml')
+ ForceField('openff-2.0.0.offxml', 'opc.offxml')
\ No newline at end of file
diff --git a/scripts/get_ion_nb_params_from_ambertools_frcmod.py b/scripts/get_ion_nb_params_from_ambertools_frcmod.py
new file mode 100644
index 0000000..6743f7e
--- /dev/null
+++ b/scripts/get_ion_nb_params_from_ambertools_frcmod.py
@@ -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()
+
+
diff --git a/scripts/write_opc.py b/scripts/write_opc.py
new file mode 100644
index 0000000..07ce71c
--- /dev/null
+++ b/scripts/write_opc.py
@@ -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"))
+