OPC and OPC3 water (#3654)

* Add benchmarks from Amber20 benchmark suite to standard benchmark script

* Add ensemble option; don't change hydrogen mass in amber input files

* Download and extract .tar.gz using pure python code, no wget/tar dependencies

* Rename amber tests

* add opc and opc3 models

* update to match https://bioinformatics.cs.vt.edu/~izadi/OPC_Gromacs/opc.top



* opc box, converted from ambertools-22.0-py38h6177452_1/dat/leap/lib/opcbox.off

* change values to make serialized system match one created from prmtop as close as possible

* unit test for opc water

* opc - final values, match frcmod.opc; derivation in comments

* opc3 water - final values, shows derivation

* opc3box made from ambertools 22 dat/leap/lib/opc3box.off

* add opc3 water test

* add opc and opc3 to docs

* move tests to TestForceField.py

* move opc tests out of amoeba tests, oops

* move opcbox and opc3box pdb files
Co-authored-by: Alex Izvorski <alex@genesistherapeutics.ai>

docs-source/usersguide/application/02_running_sims.rst

@@ -756,6 +756,8 @@ File                 Water Model
 :code:`tip5p.xml`    TIP5P water model\ :cite:`Mahoney2000`
 :code:`spce.xml`     SPC/E water model\ :cite:`Berendsen1987`
 :code:`swm4ndp.xml`  SWM4-NDP water model\ :cite:`Lamoureux2006`
+:code:`opc.xml`      OPC water model\ :cite:`Izadi2014`
+:code:`opc3.xml`     OPC3 water model\ :cite:`Izadi2016`
 ===================  ============================================
 
 Small molecule parameters

docs-source/usersguide/references.bib

@@ -671,3 +671,27 @@
     doi = {10.1021/acs.jpca.9b02771},
     type = {Journal Article}
 }
+
+@article{Izadi2014,
+   title={Building water models: a different approach},
+   author={Izadi, Saeed and Anandakrishnan, Ramu and Onufriev, Alexey V},
+   journal={The journal of Physical Chemistry Letters},
+   volume={5},
+   number={21},
+   pages={3863-3871},
+   year={2014},
+   doi = {10.1021/jz501780a},
+   type = {Journal Article}
+}
+
+@article{Izadi2016,
+   title={Accuracy limit of rigid 3-point water models},
+   author={Izadi, Saeed and Onufriev, Alexey V},
+   journal={The Journal of Chemical Physics},
+   volume={145},
+   number={7},
+   pages={074501},
+   year={2016},
+   doi = {10.1063/1.4960175},
+   type = {Journal Article}
+}

wrappers/python/openmm/app/data/opc.xml

+<ForceField>
+ <Info>
+  <DateGenerated>2022-06-21</DateGenerated>
+  <Reference>Izadi S, Anandakrishnan R, Onufriev AV. Building Water Models: A Different Approach. The Journal of Physical Chemistry Letters. 2014 Nov 6;5(21):3863-71. DOI:10.1021/jz501780a</Reference>
+ </Info>
+ <AtomTypes>
+  <Type name="opc-O" class="OW" element="O" mass="16.00000"/>
+  <Type name="opc-H" class="HW" element="H" mass="1.00800"/>
+  <Type name="opc-M" class="MW" mass="0"/>
+ </AtomTypes>
+ <Residues>
+  <Residue name="HOH">
+   <Atom name="O" type="opc-O"/>
+   <Atom name="H1" type="opc-H"/>
+   <Atom name="H2" type="opc-H"/>
+   <Atom name="M" type="opc-M"/>
+   <!-- weight2 = 0.15939833/(2*(cos(angle/2) * 0.87243313)) -->
+   <VirtualSite type="average3" index="3" atom1="0" atom2="1" atom3="2" weight1="0.7045587810585758" weight2="0.1477206094707121" weight3="0.1477206094707121"/>
+   <Bond from="0" to="1"/>
+   <Bond from="0" to="2"/>
+  </Residue>
+ </Residues>
+ <HarmonicBondForce>
+  <Bond class1="OW" class2="HW" length="0.087243313" k="502416.0"/>
+ </HarmonicBondForce>
+ <HarmonicAngleForce>
+  <!-- angle = arccos(1 - (1.37120510**2 / (2*0.87243313**2))) -->
+  <Angle class1="HW" class2="OW" class3="HW" angle="1.8081424254418306" k="628.02"/>
+ </HarmonicAngleForce>
+ <NonbondedForce coulomb14scale="0.833333" lj14scale="0.5">
+  <!-- sigma = (1.777167268 * 2 / (2**(1/6)))*unit.angstrom/unit.nanometer -->
+  <!-- epsilon = (0.2128008130*unit.kilocalorie)/(unit.kilojoule) -->
+  <Atom type="opc-O" charge="0" sigma="0.3166552081964338" epsilon="0.8903586015920001"/>
+  <Atom type="opc-H" charge="0.679142" sigma="0.08908987181403394" epsilon="0"/>
+  <Atom type="opc-M" charge="-1.358284" sigma="0.08908987181403394" epsilon="0"/>
+ </NonbondedForce>
+</ForceField>

wrappers/python/openmm/app/data/opc3.xml

+<ForceField>
+ <Info>
+  <DateGenerated>2022-06-21</DateGenerated>
+  <Reference>Izadi S, Onufriev AV. Accuracy limit of rigid 3-point water models. The Journal of Chemical Physics 145, 074501 (2016); DOI:10.1063/1.4960175</Reference>
+ </Info>
+ <AtomTypes>
+  <Type name="opc3-O" class="OW" element="O" mass="16.00000"/>
+  <Type name="opc3-H" class="HW" element="H" mass="1.00800"/>
+ </AtomTypes>
+ <Residues>
+  <Residue name="HOH">
+   <Atom name="O" type="opc3-O"/>
+   <Atom name="H1" type="opc3-H"/>
+   <Atom name="H2" type="opc3-H"/>
+   <Bond atomName1="O" atomName2="H1"/>
+   <Bond atomName1="O" atomName2="H2"/>
+  </Residue>
+ </Residues>
+ <HarmonicBondForce>
+  <Bond class1="OW" class2="HW" length="0.978882" k="502416.0"/>
+ </HarmonicBondForce>
+ <HarmonicAngleForce>
+  <!-- angle = arccos(1 - (1.598507**2 / (2*0.978882**2))) -->
+  <Angle class1="HW" class2="OW" class3="HW" angle="1.9106321528999624" k="628.02"/>
+ </HarmonicAngleForce>
+ <NonbondedForce coulomb14scale="0.833333" lj14scale="0.5">
+  <!-- sigma = (1.7814990 * 2 / (2**(1/6)))*unit.angstrom/unit.nanometer -->
+  <!-- epsilon = (0.163406*unit.kilocalorie)/(unit.kilojoule) -->
+  <Atom type="opc3-O" charge="-0.895170" sigma="0.31742703509365927" epsilon="0.683690704"/>
+  <Atom type="opc3-H" charge="0.447585" sigma="1" epsilon="0"/>
+ </NonbondedForce>
+</ForceField>

wrappers/python/tests/TestForceField.py

@@ -1264,6 +1264,67 @@ self.scriptExecuted = True
         energy2 = context.getState(getEnergy=True).getPotentialEnergy().value_in_unit(kilojoules_per_mole)
         self.assertAlmostEqual(energy1, energy2)
 
+    def test_OpcEnergy(self):
+        pdb = PDBFile('systems/opcbox.pdb')
+        topology, positions = pdb.topology, pdb.positions
+        self.assertEqual(len(positions), 864)
+        forcefield = ForceField('opc.xml')
+        system = forcefield.createSystem(
+            topology,
+            nonbondedMethod=PME,
+            nonbondedCutoff=0.7*nanometer,
+            constraints=HBonds,
+            rigidWater=True,
+        )
+
+        integrator = LangevinIntegrator(300*kelvin, 2.0/picoseconds, 2.0*femtoseconds)
+        simulation = Simulation(topology, system, integrator)
+        context = simulation.context
+        context.setPositions(positions)
+
+        # Compare to values computed with Amber (sander).
+        energy_amber = -2647.6233 # kcal/mol
+        energy_tolerance = 1.0
+
+        state = context.getState(getEnergy=True)
+        energy1 = state.getPotentialEnergy().value_in_unit(kilocalorie_per_mole)
+        # -2647.2222697324237
+        self.assertTrue(abs(energy1 - energy_amber) < energy_tolerance)
+
+        context.applyConstraints(1e-12)
+        state = context.getState(getEnergy=True)
+        energy2 = state.getPotentialEnergy().value_in_unit(kilocalorie_per_mole)
+        # -2647.441600693312
+        self.assertTrue(abs(energy1 - energy_amber) < energy_tolerance)
+        self.assertTrue(abs(energy1 - energy2) < energy_tolerance)
+
+    def test_Opc3Energy(self):
+        pdb = PDBFile('systems/opc3box.pdb')
+        topology, positions = pdb.topology, pdb.positions
+        self.assertEqual(len(positions), 648)
+        forcefield = ForceField('opc3.xml')
+        system = forcefield.createSystem(
+            topology,
+            nonbondedMethod=PME,
+            nonbondedCutoff=0.7*nanometer,
+            constraints=HBonds,
+            rigidWater=True,
+        )
+
+        integrator = LangevinIntegrator(300*kelvin, 2.0/picoseconds, 2.0*femtoseconds)
+        simulation = Simulation(topology, system, integrator)
+        context = simulation.context
+        context.setPositions(positions)
+
+        # Compare to values computed with Amber (sander).
+        energy_amber = -2532.1414 # kcal/mol
+        energy_tolerance = 1.0
+
+        state = context.getState(getEnergy=True)
+        energy1 = state.getPotentialEnergy().value_in_unit(kilocalorie_per_mole)
+        # -2532.4862082354407
+        self.assertTrue(abs(energy1 - energy_amber) < energy_tolerance)
+
 class AmoebaTestForceField(unittest.TestCase):
     """Test the ForceField.createSystem() method with the AMOEBA forcefield."""