Simulate octahedral/dodecaedral boxes of simple Neon atoms

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?

Thank you all.

Edit: I found now this tutorial simulating Argon, which probably will solve my issue. Just need to convert the boxes: Introduction to molecular dynamics simulations

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

; 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

(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.

Just for the records, it was a simple unit mistake. I was missing a e6 in the c12 parameter, fixed above.