Hi everyone,
So I am trying to do free energy calculations for a mutant (M205A) with the human Prion protein (1QLX). However I keep running into a lincs warning at the non equilirium phase of the free energy simulation. The simulation is executed through a framework I created but incorporates the same steps as this tutorial (Protein mutation — pmx documentation). Below is the output for the folded state non-equilibrium part (The whole output was too long for this post but essentially the commands are repeated per frame that is extracted), the MDP files
:-) GROMACS - gmx grompp, 2023.2 (-:
Executable: /usr/local/gromacs/bin/gmx_mpi
Data prefix: /usr/local/gromacs
Working dir: /home/21616019/git_projects/v5/PSA_framework/work/d5/eca632fcf7c6986caf04f9e062bb0a/scratch/nxf.lXHjGvJx22/frame4
Command line:
gmx_mpi grompp -f …/r_nonequil.mdp -c frame.gro -p …/newtop.top -o nonequil.tpr -maxwarn 1
Ignoring obsolete mdp entry ‘ns-type’
NOTE 1 [file …/r_nonequil.mdp]:
With PME there is a minor soft core effect present at the cut-off,
proportional to (LJsigma/rcoulomb)^6. This could have a minor effect on
energy conservation, but usually other effects dominate. With a common
sigma value of 0.34 nm the fraction of the particle-particle potential at
the cut-off at lambda=0.5 is around 0.0e+00, while ewald-rtol is 1.0e-05.
Generating 1-4 interactions: fudge = 1
Number of degrees of freedom in T-Coupling group System is 50060.00
There was 1 NOTE
GROMACS reminds you: “I admired Bohr very much. We had long talks together, long talks in which Bohr did practically all the talking.” (Paul Dirac)
Setting the LD random seed to -2097415
Generated 118802 of the 118828 non-bonded parameter combinations
Generated 88948 of the 118828 1-4 parameter combinations
Excluding 3 bonded neighbours molecule type ‘Protein_chain_P’
turning H bonds into constraints…
Excluding 2 bonded neighbours molecule type ‘SOL’
turning H bonds into constraints…
Excluding 1 bonded neighbours molecule type ‘NA’
turning H bonds into constraints…
Excluding 1 bonded neighbours molecule type ‘CL’
turning H bonds into constraints…
Analysing residue names:
There are: 104 Protein residues
There are: 7608 Water residues
There are: 49 Ion residues
Analysing Protein…
The largest distance between non-perturbed excluded atoms is 0.432 nm between atom 463 and 470
Determining Verlet buffer for a tolerance of 0.00025 kJ/mol/ps at 298.15 K
Calculated rlist for 1x1 atom pair-list as 1.205 nm, buffer size 0.005 nm
Set rlist, assuming 4x4 atom pair-list, to 1.201 nm, buffer size 0.001 nm
Note that mdrun will redetermine rlist based on the actual pair-list setup
Calculating fourier grid dimensions for X Y Z
Using a fourier grid of 60x60x60, spacing 0.116 0.116 0.116
Estimate for the relative computational load of the PME mesh part: 0.29
This run will generate roughly 7 Mb of data
GROMACS - gmx mdrun, 2023.2 (-:
Executable: /usr/local/gromacs/bin/gmx_mpi
Data prefix: /usr/local/gromacs
Working dir: /home/21616019/git_projects/v5/PSA_framework/work/d5/eca632fcf7c6986caf04f9e062bb0a/scratch/nxf.lXHjGvJx22/frame4
Command line:
gmx_mpi mdrun -s nonequil.tpr -deffnm nonequil -dhdl dgdl_rev4.xvg -v
Compiled SIMD is AVX2_256, but CPU also supports AVX_512 (see log).
Reading file nonequil.tpr, VERSION 2023.2 (single precision)
Changing nstlist from 10 to 100, rlist from 1.201 to 1.249
Update groups can not be used for this system because atoms that are (in)directly constrained together are interdispersed with other atoms
Using 1 MPI process
Using 72 OpenMP threads
starting mdrun ‘PMX MODEL in water’
40000 steps, 8.0 ps.
step 0
Step 1, time 0.0002 (ps) LINCS WARNING
relative constraint deviation after LINCS:
rms 2.918139, max 82.985596 (between atoms 1295 and 1298)
bonds that rotated more than 30 degrees:
atom 1 atom 2 angle previous, current, constraint length
1295 1296 90.0 0.1110 0.1288 0.1111
1295 1297 90.0 0.1110 0.1288 0.1111
1295 1298 81.3 0.1110 9.3308 0.1111
1284 1285 90.0 0.0996 0.9370 0.0997
Wrote pdb files with previous and current coordinates
Step 1, time 0.0002 (ps) LINCS WARNING
relative constraint deviation after LINCS:
rms 2.827562, max 80.210556 (between atoms 1295 and 1298)
bonds that rotated more than 30 degrees:
atom 1 atom 2 angle previous, current, constraint length
1295 1296 90.0 0.1110 0.5256 0.1111
1295 1297 90.0 0.1110 0.5243 0.1111
1295 1298 129.3 0.1110 9.0225 0.1111
1284 1285 90.0 0.0996 0.9369 0.0997
Back Off! I just backed up step1b.pdb to ./#step1b.pdb.1#
Back Off! I just backed up step1c.pdb to ./#step1c.pdb.1#
Wrote pdb files with previous and current coordinates
WARNING: Listed nonbonded interaction between particles 1291 and 1298
at distance 9.392 which is larger than the table limit 2.249 nm.
This is likely either a 1,4 interaction, or a listed interaction inside
a smaller molecule you are decoupling during a free energy calculation.
Since interactions at distances beyond the table cannot be computed,
they are skipped until they are inside the table limit again. You will
only see this message once, even if it occurs for several interactions.
IMPORTANT: This should not happen in a stable simulation, so there is
probably something wrong with your system. Only change the table-extension
distance in the mdp file if you are really sure that is the reason.
[n06:02481] *** Process received signal ***
[n06:02481] Signal: Segmentation fault (11)
[n06:02481] Signal code: Invalid permissions (2)
[n06:02481] Failing at address: 0x7f3792421d88
[n06:02481] [ 0] /lib/x86_64-linux-gnu/libc.so.6(+0x42520)[0x7f4b5f5bb520]
[n06:02481] [ 1] /usr/local/gromacs/bin/…/lib/libgromacs_mpi.so.8(+0xe1d67b)[0x7f4b6094967b]
[n06:02481] [ 2] /lib/x86_64-linux-gnu/libgomp.so.1(+0x1dc0e)[0x7f4b5f465c0e]
[n06:02481] [ 3] /lib/x86_64-linux-gnu/libc.so.6(+0x94ac3)[0x7f4b5f60dac3]
[n06:02481] [ 4] /lib/x86_64-linux-gnu/libc.so.6(clone+0x44)[0x7f4b5f69ea04]
[n06:02481] *** End of error message ***
/home/21616019/git_projects/v5/PSA_framework/work/d5/eca632fcf7c6986caf04f9e062bb0a/.command.sh: line 9: 2481 Segmentation fault (core dumped) gmx_mpi mdrun -s nonequil.tpr -deffnm nonequil -dhdl dgdl_rev$i.xvg -v.
MDP Non-equilibrium reverse:
;====================================================
; Non-equilibrium transition
;====================================================
; RUN CONTROL
;----------------------------------------------------
integrator = sd ; stochastic leap-frog integrator
nsteps = 40000 ; 2 * 40000 fs = 80 ps
dt = 0.0002 ; 2 fs
comm-mode = Linear ; remove center of mass translation
nstcomm = 1 ; frequency for center of mass motion removal
; OUTPUT CONTROL
;----------------------------------------------------
nstxout = 10000 ; save coordinates to .trr every 20 ps
nstvout = 10000 ; save velocities to .trr every 20 ps
nstfout = 10000 ; save forces to .trr every 20 ps
nstxout-compressed = 10000 ; xtc compressed trajectory output every 20 ps
compressed-x-precision = 1000 ; precision with which to write to the compressed trajectory file
nstlog = 5000 ; update log file every 10 ps
nstenergy = 1000 ; save energies every 2 ps
nstcalcenergy = 1 ; calculate energies at every step
; BONDS
;----------------------------------------------------
constraint_algorithm = lincs ; holonomic constraints
constraints = h-bonds ; constrain H-bonds
lincs-iter = 1 ; accuracy of LINCS (1 is default)
lincs-order = 6 ; also related to accuracy (4 is default)
lincs-warnangle = 30 ; maximum angle that a bond can rotate before LINCS will complain (30 is default)
continuation = yes ; formerly known as ‘unconstrained-start’ - useful for exact continuations and reruns
; NEIGHBOR SEARCHING
;----------------------------------------------------
cutoff-scheme = Verlet ; group or Verlet
ns-type = grid ; search neighboring grid cells
nstlist = 10 ; 20 fs (default is 10)
rlist = 1.2 ; short-range neighborlist cutoff (in nm)
pbc = xyz ; 3D PBC
; ELECTROSTATICS & EWALD
;----------------------------------------------------
coulombtype = PME ; Particle Mesh Ewald for long-range electrostatics
coulomb-modifier = Potential-shift-Verlet
rcoulomb = 1.2 ; short-range electrostatic cutoff (in nm)
ewald-geometry = 3d ; Ewald sum is performed in all three dimensions
pme-order = 4 ; interpolation order for PME (default is 4)
fourierspacing = 0.12 ; grid spacing for FFT
ewald-rtol = 1e-5 ; relative strength of the Ewald-shifted direct potential at rcoulomb
; VAN DER WAALS
;----------------------------------------------------
vdw-type = Cut-off ; potential switched off at rvdw-switch to reach zero at rvdw
vdw-modifier = force-switch
verlet-buffer-tolerance = 0.00025
rvdw = 1.2 ; van der Waals cutoff (in nm)
rvdw-switch = 1.0
DispCorr = no ; apply analytical long range dispersion corrections for Energy and Pressure
; TEMPERATURE COUPLING (Langevin)
;----------------------------------------------------
tc-grps = System
tau-t = 1.0
ref-t = 298.15
gen-vel = no ; Velocity generation (if gen-vel is ‘yes’, continuation should be ‘no’)
gen-seed = -1
; PRESSURE COUPLING
;----------------------------------------------------
pcoupl = C-rescale
pcoupltype = isotropic ; uniform scaling of box vectors
tau_p = 2.0 ; time constant (ps)
ref_p = 1.0 ; reference pressure (bar)
compressibility = 4.5e-05 ; isothermal compressibility of water (bar^-1)
refcoord-scaling = no
; FREE ENERGY
;----------------------------------------------------
free-energy = yes
init-lambda = 1 ; start from state B
delta-lambda = -2.5e-05 ; complete transition in this number of steps
sc-coul = yes ; use soft-core also for coulombic interactions
sc-alpha = 0 ; soft-core
sc-sigma = 0.34 ;
sc-power = 1 ;
nstdhdl = 1 ; write to dhdl at each step
The same also occurs for my tripeptide M2A.