Reproducing simulations with truncated octahedron box & "Periodic box vectors must be in reduced form"

[austinweigle]

My issue is related to the truncated octahedron box conversion to triclinic, which has been discussed before. However, there are some issues that I think might be worth addressing for OpenMM. For background information, I have prepared and solvated my system in tleap, heated and equilibrated with AMBER22, and am now outputting the production traj to DCD format in OpenMM 8.1.1. This is my first time performing an OpenMM simulation with a non-cubic box.

When I first run a classical simulation in OpenMM 8.1.1, my topology information is read properly, and the output trajectory box type is "Rhombohedral". As it was discussed before (#3507), I am able to use the ucell option for trajin in cpptraj to output my trajectory with the proper visualization (i.e., back to looking like a truncated octahedron).

If I want to start a new simulation from my rhombohedral output, I can only do so in one of two ways:
(1) Simply take a frame from the Rhombohedral output as is and seed the new simulation. If I do this, I cannot recover the truncated octahedron box appearance for subsequent visualization.
(2) Use ucell to read in the trajectory, output the new seed state as an rst7, and then manually manipulate the alpha/beta/gamma angles from 109.471219 (AMBER tleap octahedron angles) to 109.471206. If I approach it this way, I can recover the truncated octahedron box appearance for subsequent visualization.

Otherwise, if I try to reseed a simulation directly from cpptraj'd openmm-ucell output, I get the error:
openmm.OpenMMException: Periodic box vectors must be in reduced form.
Note that this does not happen when the truncated octahedron output results from an AMBER simulation.

This issue has already been solved by Lester Hedges (@lohedges) (michellab/Sire#345), where he first proposed this makeshift conversion of angles to figure out a rounding issue in DCD format trajectory writing. However, it looks like this issue has not been resolved in the OpenMM base code?

Currently, without a fix, my options look like either redoing my system prep with a traditional cubic box, or going rogue with this rst7 tweaking strategy. I understand that this may seem like just a visualization issue (#3507). I can work around the visualization with stripped trajectories, but for anyone looking to visualize water or ion interactions with a solute, they may not be able to produce appropriate snapshots from OpenMM simulations with a non-cubic box. This can be problematic for people trying to work with longer trajectories on short HPC wallclock times, or people using OpenMM+truncated octahedrons for adaptive sampling regimes.

[lohedges]

Hi there. As you've found, the issue here is applying the reduction check on the basis of box vectors that might have been converted from angles (losing precision in the process) or read from a file with fixed-precision formatting. Depending on whether you are using the Python or C++ API, the appropriate checks are here or here.

Without reducing the tolerance of the assertions in the functions above it is essentially impossible to completely resolve this issue, since there will always be some combination of rounding that would fall foul of the checks. Our workaround is to pre-reduce the box vectors using the internal OpenMM reduction algorithm, so by the time they are checked via the functions above they are already guaranteed to be in a compliant form. (Obviously this only makes sense for box vectors that you know are reduced within numerical precision, i.e. not arbitrary box vectors.) This can be done as follows using the function here.

from openmm.app import AmberPrmtopFile
from openmm.app.internal.unitcell import reducePeriodicBoxVectors

# Load a prmtop file containing triclinic box vectors in reduced form
# (to within numerical precision).
prmtop = AmberPrmtopFile("system.prmtop")

# Now generate the reduced box vectors.
box_vectors = reducePeriodicBoxVectors(prmtop.box_vectors)

# Create the simulation object...
...

# Set the box vectors in the context.
simulation.context.setPeriodicBoxVectors(*box_vectors)

[austinweigle]

Thank you very much @lohedges for the concise fix posted here! I can confirm that I was able to (1) perform a simulation with an unmodified rst7 (angles as 109.471219 in the last line) file as my seed; (2) visualize the output trajectory with no problems across different trajectory formats; and (3) repeat this process again from trajectory output in Step 1 here. I am able to consistently get a triclinic periodicity without any visualization artifacts.

Just an FYI since I saw you discuss this on your Sire threads - you can use parm as well instead of just the AmberPrmtop/Inpcrd commands:

# Load system
parm=LoadParm('system.prmtop', 'system.rst7')

# Now generate the reduced box vectors
box_vectors = reducePeriodicBoxVectors(parm.box_vectors)
.
.
# Set the box vectors in the context.
simulation.context.setPeriodicBoxVectors(*box_vectors)

This could be useful for other people that have this issue when handling OpenMM-derived files that do not necessarily originate from AMBER formatted inputs (topology, coordinates). Thanks again! Really appreciate the help.

[lohedges]

Great, glad to hear that this worked for you. Triclinic space representations and conversions are the bane of my existence (at least from an interoperability perspective) so I'm pleased this helped. Yes, there is nothing AMBER specific about this (other than it being a common input format to OpenMM) so it can certainly be used as a general approach for any format (even using ParmEd to read too).