Add Lipid21 and update documentation (#5066)

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

@@ -487,6 +487,7 @@ File                                 Parameters
 :file:`amber19/protein.ff19ipq.xml`  Protein (alternative)
 :file:`amber19/DNA.OL21.xml`         DNA\ :cite:`Zgarbova2021`
 :file:`amber14/RNA.OL3.xml`          RNA
+:file:`amber19/lipid21.xml`          Lipid
 :file:`amber14/GLYCAM_06j-1.xml`     Carbohydrates and glycosylated proteins\ :cite:`Kirschner2007`
 :file:`amber19/tip3p.xml`            TIP3P water model\ :cite:`Jorgensen1983` and ions
 :file:`amber19/tip3pfb.xml`          TIP3P-FB water model\ :cite:`Wang2014` and ions
@@ -498,8 +499,8 @@ File                                 Parameters
 ===================================  ===========================================
 
 As a convenience, the file :file:`amber19-all.xml` can be used as a shortcut to
-include :file:`amber19/protein.ff19SB.xml`, :file:`amber19/DNA.OL21.xml`, and
-:file:`amber14/RNA.OL3.xml`.  In most cases, you can simply include that file,
+include :file:`amber19/protein.ff19SB.xml`, :file:`amber19/DNA.OL21.xml`,
+:file:`amber14/RNA.OL3.xml`, and :file:`amber19/lipid21.xml`.  In most cases, you can simply include that file,
 plus one of the water models, such as :file:`amber19/tip3pfb.xml` for the
 TIP3P-FB water model and ions\ :cite:`Wang2014`:
 ::
@@ -530,15 +531,6 @@ to them.
          you run the risk of having :class:`ForceField` throw an exception since
          :file:`tip3p.xml` will be missing parameters for ions in your system.
 
-.. warning::
-   The updated Lipid21 lipid force field is not yet supported in this port of
-   Amber19, as it makes use of Amber features not yet supported in
-   `ParmEd <https://github.com/parmed/parmed>`_.  Amber19 should be preferred
-   over Amber14 for simulations not requiring a lipid force field, but Amber14
-   should be used if the Lipid17 force field is desired.  Alternatively, to use
-   Amber19 with Lipid21, you can prepare your system with AmberTools_ before
-   loading it into OpenMM, as described in Section :numref:`using_amber_files`.
-
 The converted parameter sets come from the `AmberTools 24 release <http://ambermd.org/AmberTools.php>`_
 and were converted using the openmmforcefields_ package and `ParmEd <https://github.com/parmed/parmed>`_.
 

wrappers/python/openmm/app/data/amber19-all.xml

@@ -2,4 +2,5 @@
  <Include file="amber19/protein.ff19SB.xml"/>
  <Include file="amber19/DNA.OL21.xml"/>
  <Include file="amber14/RNA.OL3.xml"/>
+ <Include file="amber19/lipid21.xml"/>
 </ForceField>

wrappers/python/tests/TestModeller.py

@@ -1234,47 +1234,48 @@ class TestModeller(unittest.TestCase):
         """Test adding a membrane to a realistic system."""
 
         mol = PDBxFile('systems/gpcr.cif')
-        modeller = Modeller(mol.topology, mol.positions)
-        ff = ForceField('amber14-all.xml', 'amber14/tip3p.xml')
+        for ff_files in [['amber14-all.xml', 'amber14/tip3p.xml'], ['amber19-all.xml', 'amber19/opc3.xml']]:
+            modeller = Modeller(mol.topology, mol.positions)
+            ff = ForceField(*ff_files)
 
-        # Add a membrane around the GPCR
-        modeller.addMembrane(ff, minimumPadding=1.1*nanometers, ionicStrength=1*molar)
+            # Add a membrane around the GPCR
+            modeller.addMembrane(ff, minimumPadding=1.1*nanometers, ionicStrength=1*molar)
 
-        # Make sure we added everything correctly
-        resCount = defaultdict(int)
-        for res in modeller.topology.residues():
-            resCount[res.name] += 1
+            # Make sure we added everything correctly
+            resCount = defaultdict(int)
+            for res in modeller.topology.residues():
+                resCount[res.name] += 1
 
-        self.assertEqual(16, resCount['ALA'])
-        self.assertEqual(226, resCount['POP'])  # 2x128 - overlapping
-        self.assertTrue(resCount['HOH'] > 1)
+            self.assertEqual(16, resCount['ALA'])
+            self.assertEqual(226, resCount['POP'])  # 2x128 - overlapping
+            self.assertTrue(resCount['HOH'] > 1)
 
-        deltaQ = resCount['CL'] - resCount['NA']
-        self.assertEqual(deltaQ, 10)  # protein net q: +10
+            deltaQ = resCount['CL'] - resCount['NA']
+            self.assertEqual(deltaQ, 10)  # protein net q: +10
 
-        # Check _addIons did the right thing.
-        expected_ion_fraction = 1.0*molar/(55.4*molar)
+            # Check _addIons did the right thing.
+            expected_ion_fraction = 1.0*molar/(55.4*molar)
 
-        total_water = resCount['HOH']
-        total_water_ions = resCount['HOH'] + resCount['CL'] + resCount['NA']
+            total_water = resCount['HOH']
+            total_water_ions = resCount['HOH'] + resCount['CL'] + resCount['NA']
 
-        # total_water_ions - protein charge
-        expected_sodium = math.floor((total_water_ions-10)*expected_ion_fraction+0.5)
-        expected_chlorine = expected_sodium + 10
+            # total_water_ions - protein charge
+            expected_sodium = math.floor((total_water_ions-10)*expected_ion_fraction+0.5)
+            expected_chlorine = expected_sodium + 10
 
-        self.assertEqual(resCount['CL'], expected_chlorine)
-        self.assertEqual(resCount['NA'], expected_sodium)
+            self.assertEqual(resCount['CL'], expected_chlorine)
+            self.assertEqual(resCount['NA'], expected_sodium)
 
-        # Check lipid numbering for repetitions
-        lipidIdList = [(r.chain.id, r.id) for r in modeller.topology.residues()
-                       if r.name == 'POP']
-        self.assertEqual(len(lipidIdList), len(set(lipidIdList)))
+            # Check lipid numbering for repetitions
+            lipidIdList = [(r.chain.id, r.id) for r in modeller.topology.residues()
+                        if r.name == 'POP']
+            self.assertEqual(len(lipidIdList), len(set(lipidIdList)))
 
-        # Check dimensions to see if padding was respected
-        originalSize = max(mol.positions) - min(mol.positions)
-        newSize = modeller.topology.getUnitCellDimensions()
-        for i in range(3):
-            self.assertTrue(newSize[i] >= originalSize[i]+1.1*nanometers)
+            # Check dimensions to see if padding was respected
+            originalSize = max(mol.positions) - min(mol.positions)
+            newSize = modeller.topology.getUnitCellDimensions()
+            for i in range(3):
+                self.assertTrue(newSize[i] >= originalSize[i]+1.1*nanometers)
 
     def test_bondTypeAndOrderPreserved(self):
         """ Check that bond type and order are preserved across multiple operations.