Since molecular visualizers are able to display this system as a rectangle, I’m assuming there must be some way to generate an xtc file where the box dimensions are of the form x y z 0 0 0 0 0 0.
I would like to use gromacs-chain-coordinate to analyze solvent infiltration into the membrane, but it requires the box to have angles 90 90 90. If I try to analyze my current trajectory for the CV defined in this fork with gmx mdrun -v -rerun md.xtc -s pull.tpr, I get
Fatal error:
Your box does not have a rectangular x-y area (box[XX][YY] = 0 -- box[YY][XX]
= 9.35919)
We are not sure if slice pulling (chain coordinate) works with such a box
correctly.
gmx editconf is able to convert single frames into a massive cubic box, which does seem to be compatible with the collective variable, but it can’t handle .xtc files. Is there any way to generate a rectangular or cubic box trajectory from my current trajectory file?
This might be a long shot and probably not the correct way to do it, but you can try to use trjconv with the flag -sep to split the trajectory in separate gro files, use editconf to change the shape and write is as an .xtc again, and then rejoin the frames with trjcat. Lengthy procedure, not optimized at all, but might work.
I decided to see if I could code up something with MDAnalysis to convert the box to orthorhombic since both the editconf and -box options create boxes which are larger than the system. I have something that seems very close, but the potential energy of the system before and after transformation are not the same (after transformation, the potential is high, but there aren’t any clashes or stretched bonds bad enough for Gromacs to complain). I can’t tell if I missed something, or if that’s just imprecision caused by applying a bunch of transformations to an already low-precision xtc file.
Here’s the code:
import numpy as np
import MDAnalysis as mda
from MDAnalysis import transformations
import argparse
def triclinic_to_orthogonal_matrix(dims):
# Transform MDA's box representation to box matrix
box_matrix = mda.lib.mdamath.triclinic_vectors(dims)
v1, v2, v3 = box_matrix
# Original box volume
original_volume = mda.lib.mdamath.box_volume(dims)
# X & Z dimension: Correct from original file (ONLY FOR 90 90 X BOXES)
lx = np.linalg.norm(v1)
lz = np.linalg.norm(v3)
# Y dimension: volume needs to match the original volume
ly = original_volume / (lx * lz)
ortho_box = np.diag([lx, ly, lz])
# Get the transformation to go from the original box to the orthogonal box
box_inv = np.linalg.inv(box_matrix.T)
transform = ortho_box @ box_inv
return transform, ortho_box
def transform_trajectory(input_tpr, input_traj, output_traj):
"""
Transform a trajectory from triclinic to orthogonal box representation.
Parameters:
-----------
input_tpr : str
Path to topology file (.tpr, .pdb, .gro, etc.)
input_traj : str
Path to input trajectory file (.xtc, .trr, .gro, etc.)
output_traj : str
Path to output trajectory file
"""
# Load trajectory
u = mda.Universe(input_tpr, input_traj)
print(f"Processing {len(u.trajectory)} frames...")
# Process trajectory
with mda.Writer(output_traj, u.atoms.n_atoms) as W:
for i, ts in enumerate(u.trajectory):
if i == len(u.trajectory) - 1:
print("(Last frame stats)")
original_volume = mda.lib.mdamath.box_volume(ts.dimensions)
print(f"Original box: {ts.dimensions}")
print(f"Original box volume: {original_volume:.2f}")
# Calculate transformation
transform, ortho_matrix = triclinic_to_orthogonal_matrix(ts.dimensions)
# Apply transformations
u.atoms.positions = u.atoms.positions.copy() @ transform.T
ts.dimensions = mda.lib.mdamath.triclinic_box(*ortho_matrix)
transformations.wrap(u.atoms)
# Write frame
W.write(u.atoms)
if i == len(u.trajectory) - 1:
print(f"New orthogonal box: {ts.dimensions}")
print(f"Box volume preserved: {original_volume:.2f} -> {mda.lib.mdamath.box_volume(ts.dimensions):.2f}")
print(f"Transformation complete! Output written to {output_traj}")
def argument_parser():
parser = argparse.ArgumentParser(description="Transform trajectory to orthonormal box.")
parser.add_argument("topology", type=str, help="Path to the topology file (e.g., .tpr).")
parser.add_argument("trajectory", type=str, help="Path to the trajectory file (e.g., .xtc).")
parser.add_argument("output", type=str, help="Path to save the transformed trajectory.")
return parser
if __name__ == "__main__":
parser = argument_parser()
args = parser.parse_args()
transform_trajectory(args.topology, args.trajectory, args.output)
But you can’t change the unitcell of a simulation without affecting the physics! I thought you only wanted this to avoid an issue during visualization.
Correct me if I’m wrong, but isn’t there an isomorphism between a triclinic box with only one non-orthogonal angle and an orthorhombic box? I’m trying to run an analysis which calculates the positions of polar atoms relative to a membrane normal, but the way it was implemented, it only runs on orthorhombic boxes. So I’m trying to find the best way to convert my existing simulations to something which works with the calculation.
No, you can not convert these unit cells. What would work is duplicating the system along y. Then you can subtract the a-vector from the b-vector and get an orthorhombic box. But I don’t know if there is a tool that allows you to extend a trajectory of a system along a box vector.