How should the box size be increased ?
Since I am very new to Gromacs, kindly help…
I find we have to use editconf command with -d option…
Should I change the box vectors ?
If this is the actual command,
gmx editconf -f complex.gro -o complex_newbox.gro -c -box 6.41840 6.44350 6.59650
changing it to
gmx editconf -f complex.gro -o complex_newbox.gro -c -d 1.5 -box 6.41840 6.44350 6.59650
is correct ?
How should the box size be increased ?
So you have this error when youuse gmx grompp???
What is your rvdw in mdp file??? In most force field it is 1.0 or 1.2 etc
It’s very strange that you have this error when you have 6x6x6 box.
For example if you have in you mdp file rvdw=1. Your box should be min 2x2x2 (each dimension two times bigger than rvdw)
If you want have a bigger box just define a bigger box.
gmx editconf -f complex.gro -o complex_newbox.gro -c -box 6.41840 6.44350 6.59650
gmx editconf -f complex.gro -o complex_newbox.gro -c -box 8.41840 8.44350 8.59650
and you will have bigger box (2 nm longer in x,y,z dimension).
To be honest I never used -d option. You can read that
Distance between the solute and the box
-d and a
triclinic box the size of the system in the x -, y -, and z -directions is used. With
octahedron boxes, the dimensions are set to the diameter of the system (largest distance between atoms) plus twice the specified distance
I did as you said…
For me rvdw = 1.2, so I added 2.4 to the existing box vectors and applied in editconf command but the same error appeared during shrinking iterations…
I tried with -d 1.5, again the same error after 4 th iteration …so am really confused as I couldnt find the root of the error
Please copy and paste the actual output from the
gmx editconf command. There is some other information that is printed out by the command that should assist with working out what is going on here.
This is the command and output …Its a prot-lig system and the lig topology is gen using PRODRG server.
gmx editconf -f complex.gro -o complex1_newbox.gro -c -box 8.81840 8.84350 8.99650
Read 7589 atoms
Volume: 681.651 nm^3, corresponds to roughly 306700 electrons
No velocities found
system size : 6.618 8.981 11.469 (nm)
center : 5.097 5.597 6.025 (nm)
box vectors : 6.618 8.981 11.468 (nm)
box angles : 90.00 90.00 90.00 (degrees)
box volume : 681.65 (nm^3)
shift : -0.688 -1.175 -1.526 (nm)
new center : 4.409 4.422 4.498 (nm)
new box vectors : 8.818 8.844 8.997 (nm)
new box angles : 90.00 90.00 90.00 (degrees)
new box volume : 701.60 (nm^3)
.mdp file is given below
minim.mdp - used as input into grompp to generate em.tpr
; Parameters describing what to do, when to stop and what to save
define = -DSTRONG_POSRES -DSTRONG_POSRES_LIG ; Prevent protein AND ligand from moving
integrator = steep ; Algorithm (steep = steepest descent minimization)
emtol = 1000.0 ; Stop minimization when the maximum force < 1000.0 kJ/mol/nm
emstep = 0.01 ; Energy step size
nsteps = 50000 ; Maximum number of (minimization) steps to perform
; Parameters describing how to find the neighbors of each atom and how to calculate the interactions
nstlist = 1 ; Frequency to update the neighbor list and long range forces
cutoff-scheme = Verlet ; Buffered neighbor searching
ns_type = grid ; Method to determine neighbor list (simple, grid)
rlist = 1.2 ; Cut-off for making neighbor list (short range forces)
coulombtype = cutoff ; Treatment of long range electrostatic interactions
rcoulomb = 1.2 ; Short-range electrostatic cut-off
rvdw = 1.2 ; Short-range Van der Waals cut-off
pbc = xyz ; Periodic Boundary Conditions
Also, while generating the system.gro, I pasted the new box vectors at the bottom…that is 888 dimension and removed the 666 vectors…
Thank you in advance…
A few things here.
Your assigned box vectors are inconsistent with the actual size of the system. Clearly your system is much larger along the z-axis than you’re trying to assign it, and that will just create problems.
You should not need to do any manual assignment or changing of box vectors when using InflateGRO, which I’m assuming is what you’re doing. Start over. Post the original error message you get and don’t manipulate the box. We should be able to help you sort out the original error without going down another unproductive rabbit hole.
PRODRG is an unreliable source for ligand parameters. Your simulation will be deeply flawed at best, and entirely wrong at worst. If you want GROMOS parameters, use ATB.
Thank you so much for the timely response…
I will try as suggested and post the message.
Regarding usage of ATB server
I tried using ATB server (already read various responses on how to use ATB files) where I have couple of queries like,
The format to be used …GROMOS or GROMACS ?
Bcz in GROMACS format, we have united atom ITP and united atom PDB(original geo) files and in
GROMOS96 format, we have united atom MTB and united atom PDB (original geo) files.
I read that MTB file should not be used on account of some charge calculation disparity…
Which should be chosen ?
Also, I get “atom type not found” while using the ATB gen files ?
Pardon me if am wrong in any of the above statements and hope to seek a clarity on these …
The following is the error that I get when I run the iterations (third run) without changing any of the box vectors…Here I didnt use the bash script for automation…
In a previous run of the same file with automation, I got the output but when I take it for NVT, I got the error “segmentation fault, core dumped”, with too many LINCS warning…I was trying to sort out that error when I ended up on reasons like 1) energy not minimised properly or 2) fault in ligand topology or 3) error in nvt.mdp file…
gmx grompp -f minim_inflategro.mdp -c system_inflated3.gro -p topol.top -r system_inflated3.gro -o system_inflated3_em.tpr
NOTE 1 [file minim_inflategro.mdp]:
With Verlet lists the optimal nstlist is >= 10, with GPUs >= 20. Note
that with the Verlet scheme, nstlist has no effect on the accuracy of
Setting the LD random seed to -1227739524
Generated 813 of the 2346 non-bonded parameter combinations
Excluding 3 bonded neighbours molecule type ‘Protein’
Excluding 3 bonded neighbours molecule type ‘PDB’
Excluding 3 bonded neighbours molecule type ‘DPPC’
NOTE 2 [file topol.top, line 47383]:
System has non-zero total charge: -19.000000
Total charge should normally be an integer. See
for discussion on how close it should be to an integer.
ERROR 1 [file topol.top, line 47383]:
ERROR: The cut-off length is longer than half the shortest box vector or
longer than the smallest box diagonal element. Increase the box size or
Removing all charge groups because cutoff-scheme=Verlet
There were 2 notes
Program: gmx grompp, version 2018.8
Source file: src/gromacs/gmxpreprocess/grompp.cpp (line 1991)
There was 1 error in input file(s)
For more information and tips for troubleshooting, please check the GROMACS
website at http://www.gromacs.org/Documentation/Errors
Thank you for being patient …
Also, I didn’t face any error for my apoprotein simulation, I mean I used bash script automation there…
You need GROMACS-formatted files to use in GROMACS. If you have missing atom types, consult the ATB FAQ, I bet they address this issue…
Regarding the error from
grompp, please share your
system_inflated3.gro files. You may need to post to Google Drive, Dropbox, etc. and provide us the link.
Here are the links …
The links are restricted. Please make them accessible.
I hope the files can now be accessed.
system_inflated3.gro in a text editor. It is completely misformatted. Not just the box vectors but the coordinates themselves. I have no idea how that would have happened but you should back up and start over without any manual manipulations or changes. If you have further issues, please post a complete accounting of all your commands, exactly as entered, especially the InflateGRO commands including the scaling factors and relevant screen output from running InflateGRO.
I startef from first and did E minimisation and once I put the system for nvt equilibration, I got the following errors…
Can anyone pls. tell where exactly this error points to ?
Forcefield used : Gromos96
Membrane prot-lig system in dppc
I used bash script automation for shrinking till the total area is 70 Ang sq.units and I got energy minimisation as
Steepest Descents converged to Fmax < 1000 in 5118 steps
Potential Energy = -1.9396520e+05
Maximum force = 9.2956146e+02 on atom 3723
Norm of force = 3.6939106e+01
gmx mdrun -v -deffnm nvt
Compiled SIMD: None, but for this host/run AVX2_256 might be better (see log).
The current CPU can measure timings more accurately than the code in
gmx mdrun was configured to use. This might affect your simulation
speed as accurate timings are needed for load-balancing.
Please consider rebuilding gmx mdrun with the GMX_USE_RDTSCP=ON CMake option.
Reading file nvt.tpr, VERSION 2018.8 (single precision)
Changing nstlist from 12 to 100, rlist from 1.2 to 1.252
Using 1 MPI thread
Using 8 OpenMP threads
WARNING: Using the slow plain C kernels. This should
not happen during routine usage on supported platforms.
Setting the maximum number of constraint warnings to 2147483647
starting mdrun ‘protcut - hbond-opt in water’
50000 steps, 100.0 ps.
step 100, will finish Fri Dec 4 12:59:44 2020
WARNING: Listed nonbonded interaction between particles 2174 and 2179
at distance 2.994 which is larger than the table limit 2.252 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.
My nvt.mdp file is as below (Pls correct me wherever wrong)
title = NVT equilibration for DPPC_Protein_LIG
define = -DPOSRES -DPOSRES_LIG ; position restrain the protein
; Run parameters
integrator = md ; leap-frog integrator
nsteps = 50000 ; 2 * 50000 = 100 ps
dt = 0.002 ; 2 fs
; Output control
nstxout = 100 ; save coordinates every 0.2 ps
nstvout = 100 ; save velocities every 0.2 ps
nstenergy = 100 ; save energies every 0.2 ps
nstlog = 100 ; update log file every 0.2 ps
; Bond parameters
continuation = no ; first dynamics run
constraint_algorithm = lincs ; holonomic constraints
constraints = h-bonds ; H bonds constrained
lincs_iter = 1 ; accuracy of LINCS
lincs_order = 4 ; also related to accuracy
cutoff-scheme = Verlet
ns_type = grid ; search neighboring grid cels
nstlist = 10 ; 10 fs
rlist = 1.2 ; short-range neighborlist cutoff (in nm)
rcoulomb = 1.4 ; short-range electrostatic cutoff (in nm)
rvdw = 1.4 ; short-range van der Waals cutoff (in nm)
coulombtype = PME ; Particle Mesh Ewald for long-range electrostatics
pme_order = 4 ; cubic interpolation
fourierspacing = 0.16 ; grid spacing for FFT
; Temperature coupling is on
tcoupl = V-rescale ; modified Berendsen thermostat
tc-grps = DPPC_Protein_LIG Water_and_ions ; two coupling groups - more accurate
tau_t = 0.1 0.1 ; time constant, in ps
ref_t = 323 323 ; reference temperature, one for each group, in K
; Pressure coupling is off
pcoupl = no ; no pressure coupling in NVT
; Periodic boundary conditions
pbc = xyz ; 3-D PBC
; Dispersion correction
DispCorr = EnerPres ; account for cut-off vdW scheme
; Velocity generation
gen_vel = yes ; assign velocities from Maxwell distribution
gen_temp = 323 ; temperature for Maxwell distribution
gen_seed = -1 ; generate a random seed
; COM motion removal
; These options remove motion of the protein/bilayer relative to the solvent/ions
nstcomm = 100
comm-mode = Linear
comm-grps = DPPC_Protein_LIG Water_and_ions
Hope to get some help,
Thank you in advance
After the above warnings, the following appeared…
step 200: One or more water molecules can not be settled.
Check for bad contacts and/or reduce the timestep if appropriate.
Wrote pdb files with previous and current coordinates
step 200, will finish Fri Dec 4 12:59:36 2020Segmentation fault (core dumped)
What is the source of your ligand topology? How have you validated it?
Follow the steps here to diagnose: