I’m trying to make an analysis package of mine compatible with triclinic cells, particularly with Gromacs dodecahedal and octahedral cells.

Yet, I’m not a regular Gromacs user and I struggling to obtain something that is probably very simple (so simple that tutorials do not follow that route).

What I want to obtain is a box of ~10k Neon atoms (available at the OPLS force field), in octahedral and dodecahedral boxes, and obtain a frame of an equilibrated trajectory of such system.

My goal is to use that as a test file for some routines to compute distances, energies, etc, in my analysis package, for when reading Gromacs trajectories.

I’m failing miserably in obtaining such a simple Gromacs setup. I have a PDB file with the Neon atoms. Is there some tutorial anyone knows of in the simplest setup possible is explained, starting from a given PDB file with coordinates?

After running the above tutorial, I’m almost there in the implementation of my tests. I have perhaps a more specific question.

I’m trying to compute the Lennard-Jones potential energy from the trajectory (a simulation of 100 Argon atoms), and I’m almost there.

Specifically, I get:

u = -0.8340485862647321

While Gromacs reports, for the same frame:

LJ (SR) Coulomb (SR) Potential Kinetic En. Total Energy
-8.28121e-01 0.00000e+00 -8.28121e-01 1.27260e+02 1.26432e+02

Thus, I’m off by 1.5%.

That is close, but I do not think it is justifiable by numerical inaccuracies (with the same code I can reproduce much precisely energies obtained with NAMD). I’m using the following parameters in the mdp file, to simplify as much as possible the LJ computation:

; long-range cut-off for switched potentials
cutoff-scheme = Verlet
; NEIGHBORSEARCHING PARAMETERS
; nblist update frequency
nstlist = 1
; Periodic boundary conditions: xyz (default), no (vacuum)
; or full (infinite systems only)
pbc = xyz
; nblist cut-off
rlist = 1.3
; Method for doing Van der Waals
vdw-type = Cut-off
rvdw = 1.2
vdw-modifier = None
; Apply long range dispersion corrections for Energy and Pressure
DispCorr = no

Any hint on what I may be missing to be off by those 1.5% in the computation of the LJ energy? Is there any additional modification of the LJ parameters and interactions I may be missing yet, such that the computation is a bare-bones applications of the energy function for each pair?

The parameters for Argon are simply:

[ defaults ]
; nbfunc comb-rule gen-pairs fudgeLJ fudgeQQ
1 1 no 1.0 1.0
[ atomtypes ]
AR 39.948 0.0 A 0.00622127 9.69576e-06

and the function I use to compute the energy from the simulation (for each pair - this is computed for every pair, considering PBC and the appropriate cutoff) is just:

function lj_Argon_Gromacs(d2, u)
c6 = 0.00622127e6
c12 = 9.69576
u += c12 / d2^6 - c6 / d2^3
end

(this is a small Julia code). The only conversion here is that of the units, because the trajectory (xtc file) is given in Angstroms instead of nm.