AmberPrmtopFile accepts box vectors as a constructor argument (#4094)

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

@@ -190,7 +190,7 @@ Using AMBER Files
 OpenMM can build a system in several different ways.  One option, as shown
 above, is to start with a PDB file and then select a force field with which to
 model it.  Alternatively, you can use AmberTools_ to model your system.  In that
-case, you provide a :class:`prmtop` file and an :class:`inpcrd` file.  OpenMM loads the files and
+case, you provide a :file:`prmtop` file and an :file:`inpcrd` file.  OpenMM loads the files and
 creates a :class:`System` from them.  This is illustrated in the following script.  It can be
 found in OpenMM’s :file:`examples` folder with the name :file:`simulateAmber.py`.
 
@@ -202,15 +202,13 @@ found in OpenMM’s :file:`examples` folder with the name :file:`simulateAmber.p
         from openmm.unit import *
         from sys import stdout
 
-        prmtop = AmberPrmtopFile('input.prmtop')
         inpcrd = AmberInpcrdFile('input.inpcrd')
+        prmtop = AmberPrmtopFile('input.prmtop', periodicBoxVectors=inpcrd.boxVectors)
         system = prmtop.createSystem(nonbondedMethod=PME, nonbondedCutoff=1*nanometer,
                 constraints=HBonds)
         integrator = LangevinMiddleIntegrator(300*kelvin, 1/picosecond, 0.004*picoseconds)
         simulation = Simulation(prmtop.topology, system, integrator)
         simulation.context.setPositions(inpcrd.positions)
-        if inpcrd.boxVectors is not None:
-            simulation.context.setPeriodicBoxVectors(*inpcrd.boxVectors)
         simulation.minimizeEnergy()
         simulation.reporters.append(PDBReporter('output.pdb', 1000))
         simulation.reporters.append(StateDataReporter(stdout, 1000, step=True,
@@ -225,12 +223,12 @@ This script is very similar to the previous one.  There are just a few
 significant differences:
 ::
 
-    prmtop = AmberPrmtopFile('input.prmtop')
     inpcrd = AmberInpcrdFile('input.inpcrd')
+    prmtop = AmberPrmtopFile('input.prmtop', periodicBoxVectors=inpcrd.boxVectors)
 
-In these lines, we load the prmtop file and inpcrd file.  More precisely, we
-create :class:`AmberPrmtopFile` and :class:`AmberInpcrdFile` objects and assign them to the
-variables :code:`prmtop` and :code:`inpcrd`\ , respectively.  As before,
+In these lines, we load the :file:`inpcrd` file and :file:`prmtop` file.  More precisely, we
+create :class:`AmberInpcrdFile` and :class:`AmberPrmtopFile` objects and assign them to the
+variables :code:`inpcrd` and :code:`prmtop`\ , respectively.  As before,
 you can change these lines to specify any files you want.  Be sure to include
 the single quotes around the file names.
 
@@ -241,6 +239,19 @@ the single quotes around the file names.
     example files that are in the "old-style" :file:`prmtop` format. These "old-style" files will
     not run in OpenMM.
 
+Notice that when we load the :file:`prmtop` file we include the argument :code:`periodicBoxVectors=inpcrd.boxVectors`\ .
+AMBER stores information about the periodic box in the :file:`inpcrd` file.  To let
+:class:`AmberPrmtopFile` create a :class:`Topology` object, we therefore need to
+tell it the periodic box vectors that were loaded from the :file:`inpcrd` file.  You
+only need to do this if you are simulating a periodic system.  For implicit
+solvent simulations, it usually can be omitted.
+
+.. note::
+
+    For historical reasons, :file:`prmtop` files also have the ability to store
+    periodic box information, but it is deprecated.  It is always better to get
+    the box vectors from the :file:`inpcrd` file instead.
+
 Next, the :class:`System` object is created in a different way:
 ::
 
@@ -259,22 +270,7 @@ directly.
 
 Notice that we now get the topology from the :file:`prmtop` file and the atom positions
 from the :file:`inpcrd` file.  In the previous section, both of these came from a PDB
-file, but AMBER puts the topology and positions in separate files.  We also add the
-following lines:
-::
-
-    if inpcrd.boxVectors is not None:
-        simulation.context.setPeriodicBoxVectors(*inpcrd.boxVectors)
-
-For periodic systems, the :file:`prmtop` file specifies the periodic box vectors, just
-as a PDB file does.  When we call :meth:`createSystem`, it sets those as the default
-periodic box vectors, to be used automatically for all simulations.  However, the
-:file:`inpcrd` may *also* specify periodic box vectors,
-and if so we want to use those ones instead.  For example, if the system has been
-equilibrated with a barostat, the box vectors may have changed during equilibration.
-We therefore check to see if the :file:`inpcrd` file contained box vectors.  If so,
-we call :meth:`setPeriodicBoxVectors` to tell it to use those ones, overriding the
-default ones provided by the :class:`System`.
+file, but AMBER puts the topology and positions in separate files.
 
 .. _using_gromacs_files:
 
@@ -315,16 +311,10 @@ with the name :file:`simulateGromacs.py`.
 
 This script is nearly identical to the previous one, just replacing
 :class:`AmberInpcrdFile` and :class:`AmberPrmtopFile` with :class:`GromacsGroFile` and :class:`GromacsTopFile`.
-Note that when we create the :class:`GromacsTopFile`, we specify values for two extra
-options.  First, we specify
-:code:`periodicBoxVectors=gro.getPeriodicBoxVectors()`\ .  Unlike OpenMM and
-AMBER, which can store periodic unit cell information with the topology, Gromacs
-only stores it with the coordinates.  To let :class:`GromacsTopFile` create a :class:`Topology`
-object, we therefore need to tell it the periodic box vectors that were loaded
-from the :file:`gro` file.  You only need to do this if you are simulating a periodic
-system.  For implicit solvent simulations, it usually can be omitted.
-
-Second, we specify :code:`includeDir='/usr/local/gromacs/share/gromacs/top'`\ .  Unlike AMBER,
+As with AMBER files, when we create the :class:`GromacsTopFile` we specify
+:code:`periodicBoxVectors=gro.getPeriodicBoxVectors()` to tell it the periodic
+box vectors that were loaded from the :file:`gro` file.  In addition, we specify
+:code:`includeDir='/usr/local/gromacs/share/gromacs/top'`\ .  Unlike AMBER,
 which stores all the force field parameters directly in a :file:`prmtop` file, Gromacs just stores
 references to force field definition files that are installed with the Gromacs
 application.  OpenMM needs to know where to find these files, so the

examples/simulateAmber.py

@@ -3,14 +3,12 @@ from openmm import *
 from openmm.unit import *
 from sys import stdout
 
-prmtop = AmberPrmtopFile('input.prmtop')
 inpcrd = AmberInpcrdFile('input.inpcrd')
+prmtop = AmberPrmtopFile('input.prmtop', periodicBoxVectors=inpcrd.boxVectors)
 system = prmtop.createSystem(nonbondedMethod=PME, nonbondedCutoff=1*nanometer, constraints=HBonds)
 integrator = LangevinMiddleIntegrator(300*kelvin, 1/picosecond, 0.004*picoseconds)
 simulation = Simulation(prmtop.topology, system, integrator)
 simulation.context.setPositions(inpcrd.positions)
-if inpcrd.boxVectors is not None:
-    simulation.context.setPeriodicBoxVectors(*inpcrd.boxVectors)
 simulation.minimizeEnergy()
 simulation.reporters.append(PDBReporter('output.pdb', 1000))
 simulation.reporters.append(StateDataReporter(stdout, 1000, step=True, potentialEnergy=True, temperature=True))

wrappers/python/openmm/app/amberprmtopfile.py

@@ -82,11 +82,26 @@ def _strip_optunit(thing, unit):
 class AmberPrmtopFile(object):
     """AmberPrmtopFile parses an AMBER prmtop file and constructs a Topology and (optionally) an OpenMM System from it."""
 
-    def __init__(self, file):
-        """Load a prmtop file."""
+    def __init__(self, file, periodicBoxVectors=None, unitCellDimensions=None):
+        """Load a prmtop file.
+
+        file : str
+            the name of the file to load
+        periodicBoxVectors : tuple of Vec3=None
+            the vectors defining the periodic box
+        unitCellDimensions : Vec3=None
+            the dimensions of the crystallographic unit cell.  For
+            non-rectangular unit cells, specify periodicBoxVectors instead.
+        """
         ## The Topology read from the prmtop file
         self.topology = top = Topology()
         self.elements = []
+        if periodicBoxVectors is not None:
+            if unitCellDimensions is not None:
+                raise ValueError("Specify either periodicBoxVectors or unitCellDimensions, but not both")
+            top.setPeriodicBoxVectors(periodicBoxVectors)
+        else:
+            top.setUnitCellDimensions(unitCellDimensions)
 
         # Load the prmtop file
 
@@ -151,7 +166,7 @@ class AmberPrmtopFile(object):
 
         # Set the periodic box size.
 
-        if prmtop.getIfBox():
+        if top.getPeriodicBoxVectors() is None and prmtop.getIfBox():
             box = prmtop.getBoxBetaAndDimensions()
             top.setPeriodicBoxVectors(computePeriodicBoxVectors(*(box[1:4] + box[0:1]*3)))
 

wrappers/python/openmm/app/internal/amber_file_parser.py

@@ -12,7 +12,7 @@ Simbios, the NIH National Center for Physics-Based Simulation of
 Biological Structures at Stanford, funded under the NIH Roadmap for
 Medical Research, grant U54 GM072970. See https://simtk.org.
 
-Portions copyright (c) 2012-2022 Stanford University and the Authors.
+Portions copyright (c) 2012-2023 Stanford University and the Authors.
 Authors: Randall J. Radmer, John D. Chodera, Peter Eastman
 Contributors: Christoph Klein, Michael R. Shirts, Jason Swails, Kye Won Wang
 
@@ -890,7 +890,7 @@ def readAmberSystem(topology, prmtop_filename=None, prmtop_loader=None, shake=No
     # Add nonbonded interactions.
     if verbose: print("Adding nonbonded interactions...")
     force = mm.NonbondedForce()
-    if (prmtop.getIfBox() == 0):
+    if topology.getPeriodicBoxVectors() is None and prmtop.getIfBox() == 0:
         # System is non-periodic.
         if nonbondedMethod == 'NoCutoff':
             force.setNonbondedMethod(mm.NonbondedForce.NoCutoff)
@@ -904,9 +904,12 @@ def readAmberSystem(topology, prmtop_filename=None, prmtop_loader=None, shake=No
     else:
         # System is periodic.
         # Set periodic box vectors for periodic system
-        (boxBeta, boxX, boxY, boxZ) = prmtop.getBoxBetaAndDimensions()
-        xVec, yVec, zVec = computePeriodicBoxVectors(boxX, boxY, boxZ, boxBeta, boxBeta, boxBeta)
-        system.setDefaultPeriodicBoxVectors(xVec, yVec, zVec)
+        if topology.getPeriodicBoxVectors() is None:
+            (boxBeta, boxX, boxY, boxZ) = prmtop.getBoxBetaAndDimensions()
+            xVec, yVec, zVec = computePeriodicBoxVectors(boxX, boxY, boxZ, boxBeta, boxBeta, boxBeta)
+            system.setDefaultPeriodicBoxVectors(xVec, yVec, zVec)
+        else:
+            system.setDefaultPeriodicBoxVectors(*topology.getPeriodicBoxVectors())
 
         # Set cutoff.
         if nonbondedCutoff is None:

wrappers/python/tests/TestAmberPrmtopFile.py

@@ -7,16 +7,16 @@ from openmm import *
 from openmm.unit import *
 import openmm.app.element as elem
 
+inpcrd3 = AmberInpcrdFile('systems/ff14ipq.rst7')
+inpcrd4 = AmberInpcrdFile('systems/Mg_water.inpcrd')
+inpcrd7 = AmberInpcrdFile('systems/18protein.rst7')
 prmtop1 = AmberPrmtopFile('systems/alanine-dipeptide-explicit.prmtop')
 prmtop2 = AmberPrmtopFile('systems/alanine-dipeptide-implicit.prmtop')
-prmtop3 = AmberPrmtopFile('systems/ff14ipq.parm7')
-prmtop4 = AmberPrmtopFile('systems/Mg_water.prmtop')
+prmtop3 = AmberPrmtopFile('systems/ff14ipq.parm7', periodicBoxVectors=inpcrd3.boxVectors)
+prmtop4 = AmberPrmtopFile('systems/Mg_water.prmtop', periodicBoxVectors=inpcrd4.boxVectors)
 prmtop5 = AmberPrmtopFile('systems/tz2.truncoct.parm7')
 prmtop6 = AmberPrmtopFile('systems/gaffwat.parm7')
-prmtop7 = AmberPrmtopFile('systems/18protein.parm7')
-inpcrd3 = AmberInpcrdFile('systems/ff14ipq.rst7')
-inpcrd4 = AmberInpcrdFile('systems/Mg_water.inpcrd')
-inpcrd7 = AmberInpcrdFile('systems/18protein.rst7')
+prmtop7 = AmberPrmtopFile('systems/18protein.parm7', periodicBoxVectors=inpcrd7.boxVectors)
 
 class TestAmberPrmtopFile(unittest.TestCase):
 
@@ -177,7 +177,6 @@ class TestAmberPrmtopFile(unittest.TestCase):
         # 'working', and the system is plenty small so this won't be too slow
         sim = Simulation(prmtop3.topology, system, integrator, Platform.getPlatformByName('Reference'))
         # Check that the energy is about what we expect it to be
-        sim.context.setPeriodicBoxVectors(*inpcrd3.boxVectors)
         sim.context.setPositions(inpcrd3.positions)
         ene = sim.context.getState(getEnergy=True, enforcePeriodicBox=True).getPotentialEnergy()
         ene = ene.value_in_unit(kilocalories_per_mole)
@@ -208,7 +207,6 @@ class TestAmberPrmtopFile(unittest.TestCase):
         # 'working', and the system is plenty small so this won't be too slow
         sim = Simulation(prmtop3.topology, system, integrator, Platform.getPlatformByName('Reference'))
         # Check that the energy is about what we expect it to be
-        sim.context.setPeriodicBoxVectors(*inpcrd3.getBoxVectors())
         sim.context.setPositions(inpcrd3.getPositions())
         ene = sim.context.getState(getEnergy=True, enforcePeriodicBox=True).getPotentialEnergy()
         ene = ene.value_in_unit(kilocalories_per_mole)
@@ -254,12 +252,19 @@ class TestAmberPrmtopFile(unittest.TestCase):
         self.assertTrue(has_nonbond_force)
         self.assertTrue(has_custom_nonbond_force)
         self.assertEqual(nonbond_exceptions, custom_nonbond_exceptions)
+        # Make sure the periodic box vectors match the ones in the inpcrd file, which are different from
+        # the ones in the prmtop file.
+        systemBox = system.getDefaultPeriodicBoxVectors()
+        topologyBox = prmtop4.topology.getPeriodicBoxVectors()
+        for i in range(3):
+            for j in range(3):
+                self.assertEqual(inpcrd4.boxVectors[i][j], systemBox[i][j])
+                self.assertEqual(inpcrd4.boxVectors[i][j], topologyBox[i][j])
         integrator = VerletIntegrator(1.0*femtoseconds)
         # Use reference platform, since it should always be present and
         # 'working', and the system is plenty small so this won't be too slow
         sim = Simulation(prmtop4.topology, system, integrator, Platform.getPlatformByName('Reference'))
         # Check that the energy is about what we expect it to be
-        sim.context.setPeriodicBoxVectors(*inpcrd4.boxVectors)
         sim.context.setPositions(inpcrd4.positions)
         ene = sim.context.getState(getEnergy=True, enforcePeriodicBox=True).getPotentialEnergy()
         ene = ene.value_in_unit(kilocalories_per_mole)
@@ -358,7 +363,6 @@ class TestAmberPrmtopFile(unittest.TestCase):
 
         simulation = Simulation(prmtop.topology, system, integrator)
         simulation.context.setPositions(inpcrd.positions)
-        simulation.context.setPeriodicBoxVectors(*inpcrd.boxVectors)
 
         fname = tempfile.mktemp(suffix='.dcd')
         simulation.reporters.append(DCDReporter(fname, 1))  # This is an explicit test for the bugs in issue #850
@@ -435,8 +439,7 @@ class TestAmberPrmtopFile(unittest.TestCase):
 
         simulation = Simulation(prmtop.topology, system, integrator)
         simulation.context.setPositions(inpcrd.positions)
-        simulation.context.setPeriodicBoxVectors(*inpcrd.boxVectors)
-        
+
         for i, force in enumerate(system.getForces()):
             force.setForceGroup(i)