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.