I am currently working with a trajectory obtained from a GROMACS simulation involving 200 hydrogen molecules solvated in TIP4P water. The simulation data includes both positions and velocities and is saved in the .xtc and .trr formats.
My objective is to calculate the pressure exerted by the hydrogen atoms within the system. Specifically, I am interested in utilizing LAMMPS’s compute stress/atom feature to perform these pressure calculations directly from the existing GROMACS trajectory data, thereby avoiding the need to rerun the simulation in LAMMPS.
Is it feasible to:
Import a GROMACS trajectory (.xtc or .trr) into LAMMPS for post-processing?
Apply the compute stress/atom command (or similar LAMMPS features) to the imported trajectory to accurately calculate the pressure profiles of the hydrogen atoms?
If this approach is possible, I would greatly appreciate any guidance on briefly mentioning the required commands for data conversion.
LAMMPS does not natively read the two file formats you mention, but can be configured to read them though an interface to VMD molfile plugins. It should be possible to do this with a “rerun”, but I have not tried this for a long time.
To use compute stress/atom, you need to have the entire system definition, forcefield parameters and partial changes, topology, i.e. “everything” translated to LAMMPS, since the compute doesn’t do any computation, but rather makes data available that is tallied and stored during force field calculations.
To convert the virial stress to pressure, you need the volume over which the stress is summed, which is ill defined for individual atoms.
I’ve previously been able to initialise a LAMMPS system (by defining the force field’s pair and other styles and initialising the topology with read_data from a LAMMPS data file) and subsequently post-process trajectory coordinates from a GROMACS XTC file with rerun (and MOLFILE plugins).
So that’s possible and not particularly difficult. However:
LAMMPS will not exactly replicate GROMACS forces:
GROMACS neighbourlisting is different;
Long-range electrostatics are slightly different in GROMACS (PME) vs LAMMPS (PPPM);
LAMMPS doesn’t implement GROMACS’s LINCS for rigid bonds
While reading coordinates from an XTC file, you can’t also then read velocities “in parallel” from a TRR file.
Your best bet may be to process the configurational stress components first with the XTC coordinates, then process the kinetic stress component separately from the TRR file and somehow add those together.
Thank you all for your responses. I saved the positions and velocities in .trr format and then converted that to lammpstrj format using MDAnalysis and proper unit conversion. A few frames are attached here for your reference.
I also made a data file which is a system of hydrogen gas (treated as beads instead of molecules, Buch model, ref [1]) solvated in TIP4P water (4 point geometry).
Now when I rerun the calculations, I get this warning:
Also the temperature is ~150K instead of the original 300K and pressure is very low : ~-229000 atm. Any suggestions or insights into what I might be missing or need to adjust would be immensely helpful.
Here’s my input script for rerun command:
# LAMMPS Input Script with Bonds and Angles
units real
atom_style full
boundary p p p
# ---------------------- Read Data File ----------------------
read_data simulation.data
# ---------------------- Force Fields ----------------------
pair_style lj/cut/coul/long 10.0
kspace_style pppm 1.0e-4
# Pair Coefficients
pair_coeff 1 1 0.1550 3.1536 # OW-O interactions
pair_coeff 2 2 0.0 1.0 # HW1-HW1 & HW2-HW2 interactions
pair_coeff 3 3 0.0 1.0 # MW-MW interactions
pair_coeff 4 4 0.06796 2.96 # H₂ beads interactions
pair_modify mix arithmetic
# ---------------------- Bond and Angle Styles ----------------------
bond_style zero
bond_coeff 1 0.9574
angle_style zero
angle_coeff 1 104.52
# ---------------------- Define Groups ----------------------
group water type 1 2 3
group hydrogen_beads type 4
# ---------------------- Computes and Variables ----------------------
# Compute center of mass of hydrogen beads
compute com_hyd hydrogen_beads com
# Define variables for COM coordinates
variable comx equal c_com_hyd[1]
variable comy equal c_com_hyd[2]
variable comz equal c_com_hyd[3]
# Initialize computes by performing a zero-step run
run 0
# Compute chunking of hydrogen atoms into spherical shells around COM
compute chunk_hyd all chunk/atom bin/sphere ${comx} ${comy} ${comz} 0.0 30.0 30 discard yes
# Compute per-atom stress tensor for hydrogen beads only
compute peratom_hyd all stress/atom NULL
# Define a per-atom variable for the trace of the stress tensor
variable stress_trace atom "c_peratom_hyd[1] + c_peratom_hyd[2] + c_peratom_hyd[3]"
# Fix to average and write per-chunk stress trace every 1000 timesteps for hydrogen beads
fix pressure_profile all ave/chunk 1 1 1 chunk_hyd v_stress_trace file pressure_profile.txt norm none
# ---------------------- Thermo Settings ----------------------
thermo_style custom step temp press etotal pe density
thermo 1
# ---------------------- Rerun Command ----------------------
timestep 2.0
rerun trajectory.lammpstrj dump x y z vx vy vz box no
[1] Tsimpanogiannis, I. N., Maity, S., Celebi, A. T., & Moultos, O. A. (2021). Engineering model for predicting the intradiffusion coefficients of hydrogen and oxygen in vapor, liquid, and supercritical water based on molecular dynamics simulations. Journal of Chemical & Engineering Data , 66 (8), 3226-3244.
@srtee is correct, since you are trying to use an explicit TIP4P rigid model, the temperature should only be computed using the relevant rotational degrees of freedom of the rigid bodies. It might be easier to retrieve them using unwrapped coordinates so that LAMMPS can reconstruct them using fix rigid/small molecule and compute the relevant degrees of freedom correctly. Using fix shake might require adding other bonds with the massless particle, though I am unsure about that.
This would also affect the computation of the pressure.
With TIP4P you would have 4\cdot3 = 12 DOFs for the 4 atoms, but making it rigid, you are left with 6 DOFs (3 for rotation, 3 translation). So you could just remove 6 DOFs per molecule and the temperature compute should return the expected value.
Since this is going to be a rerun job, there is no need to issue a fix (fix shake won’t work, BTW) to remove the DOFs.
Update: what you do want to do in this case, is to remove the intra-molecular force computation inside the TIP4P modules by either using neigh_modify or by adding bonds with bond style zero between all internal interactions.