[gmx-users] Re: Water molecules cannot be settled, why?
John Ladasky
blind.watchmaker at yahoo.com
Wed Jul 11 10:12:07 CEST 2012
Hello, everyone,
I am reviving a discussion I started about a month back.
http://gromacs.5086.n6.nabble.com/Water-molecules-cannot-be-settled-why-tt4998059.html
I thought that I would try to work around the problem by only proceeding
with those simulations that didn't crash, but I've reached a point of
diminishing returns. I even got one SETTLE error during my
position-restrained preprocessing step. I have to fix this problem.
Jun 04, 2012; 1:36pm (Justin A. Lemkul):
> The fact that the water molecule cannot satisfy the SETTLE algorithm
> does not necessarily indicate that water is a problem - other clashes and
> particles moving rapidly across the system can slam into an
unsuspecting water
> molecule and cause issues.
OK, I'll keep this in mind. I am making chimeric molecules by
manipulating PDB files. It's possible that I'm making bad joints when I
build my chimeras. I will be surprised to learn that the minor
alterations I've made to bond angles could impart so much energy to the
system, but I'll take your word for it that I should be concerned.
If there is a standard way to make chimeras with GROMACS, I am unaware
of it. I hardly see how there could be, given that there are many
arbitrary choices to make. The user would want a lot of control over
bond angles and the positions of the fragments that are joined together.
I have tried two other software tools to build my chimeras. PyMol's
"fuse" command does something that I do not want. It makes nice bond
angles at the joint, but it refuses to adjust bond angles anywhere
inside the fragments. And so, PyMol has no compunction about putting
one fragment straight through the middle of another (a condition which,
I think, precludes the use of GROMACS pull code to correct it). If I
want to work with what PyMol gives me, I have to untangle the two halves
of the chimera by selecting and rotating various bonds on the peptide
backbone. My choices are arbitrary and chosen by eye. I may be
introducing strain into the system when I do this.
Because I was not satisfied with PyMol, I also wrote a custom chimera
program (based on Biopython) which assumed that the peptide fragments
were in approximately the correct orientation before joining. I still
need to bend a bond angle or two 20 or 30 degrees here and there, in
order to expose the terminal amino acids before I start. I superimpose
a template of a trans-peptide bond at equilibrium on the C-terminus of
the N-terminal peptide, then use the coordinates of the idealized
peptide bond to translate (but not rotate) the C-terminal peptide into
position. The results, by eye at least, look believable.
> Have you watched the
> trajectory to see if anything weird is going on?
I haven't seen anything egregious, but I'm not watching the solvent
molecules, only the peptide. The final two frames that GROMACS dumps
when the simulation crashes don't tell me much, though my own lack of
understanding may be at fault for that.
> We would need to see an .mdp file to be able to provide any help. MD is
> chaotic, and there is no guarantee that surviving the first few (hundred,
> thousand...) steps that the simulation will be stable.
> What is your minimization and equilibration protocol?
Here is my entire protocol, please pardon the length.
Shell commands for setting up the initial energy minimization:
pdb2gmx -f <pdb file> -ff gromos45a3 -water spc -ignh
editconf -f <pdb2gmx output> -bt cubic -d 1.2
genbox -cp <editconf output> -cs spc216.gro
grompp -f em.mdp -c <genbox output>
echo SOL | genion -s <genbox output> -neutral -conc 0.1
grompp -f em.mdp -c <genion output>
mpirun -np 5 mdrun_mpi -s <genion output> -o em
Here's my em.mdp file:
; em.mdp
; molecular dynamics parameters for the initial energy minimization,
; just after adding solvent
integrator = steep ; Algorithm (steepest descent minimization)
emtol = 1000.0 ; Stop when max force < 1000 kJ/mol/nm
emstep = 0.01 ; Energy step size
nsteps = 10000 ; Maximum number of minimization steps to perform.
; In practice, no more than 1000 cycles are used
; before the emtol force minimum is satisfied.
nstlist = 1 ; Frequency to update neighbor list, long-range
forces
ns_type = grid ; Method to determine neighbor list (simple, grid)
coulombtype = PME ; Treatment of long range electrostatic interactions
rlist = 1.4 ; Cut-off for neighbor list (short-range forces)
rcoulomb = 1.4 ; Short-range electrostatic cut-off (nm)
rvdw = 1.4 ; Short-range Van der Waals cut-off (nm)
pbc = xyz ; Periodic Boundary Conditions (in all axes)
My pr.mdp file (used in the next ground of grompp and mpirun, shell
commands not shown):
; pr.mdp
; molecular dynamics parameters for the second energy minimization,
; settles the solvent further while holding the protein stationary.
define = -DPOSRES ; position-restrain the protein
; Run parameters
integrator = md ; leap-frog integrator
dt = 0.002 ; 2 fs
nsteps = 50000 ; 2 fs * 25000 = 100 ps
; Output control
; Notice no nstfout, we aren't interested in the details of this run.
; In truth, we might not need the .trr file from this run at all.
; Only the last frame is of interest, unless perhaps for debugging.
nstxout = 2500 ; save coordinates every 5 ps
nstvout = 2500 ; save velocities every 5 ps
nstenergy = 2500 ; save energies every 5 ps
nstlog = 2500 ; update log file every 5 ps
; Bond parameters
continuation = no ; Initial simulation
constraint_algorithm = lincs ; holonomic constraints
constraints = all-bonds ; (even heavy atom-H bonds)
lincs_iter = 1 ; accuracy of LINCS
lincs_order = 4 ; also related to accuracy
; Neighbor searching
ns_type = grid ; search neighboring grid cels
nstlist = 5 ; update neighbor list, long range forces
; every 5 cycles = 10 fs
rlist = 1.4 ; short-range neighbor list cutoff (in nm)
rcoulomb = 1.4 ; short-range electrostatic cutoff (in nm)
rvdw = 1.4 ; short-range van der Waals cutoff (in nm)
; Electrostatics
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 ; Weak coupling for initial equilibration
tc-grps = Protein Non-Protein ; two coupling groups, more accurate
tau_t = 0.1 0.1 ; time constant, in ps
ref_t = 310 310 ; reference temperature, one for each
; group, in Kelvins
; Pressure coupling is on
pcoupl = Berendsen ; Pressure coupling on in NPT, also weak
coupling
pcoupltype = isotropic ; uniform scaling of x-y-z box vectors
tau_p = 2.0 ; time constant, in ps
ref_p = 1.0 1.0 ; reference pressure (in bar)
compressibility = 4.5e-5 ; isothermal compressibility, bar^-1
; Periodic boundary conditions
pbc = xyz ; 3-D PBC
; Dispersion correction
DispCorr = EnerPres ; account for cut-off vdW scheme
; Velocity generation
gen_vel = yes ; Velocity generation is on
gen_temp = 310 ; temperature for velocity generation (Kelvin)
gen_seed = -1 ; random seed
; COM motion removal
; These options remove COM motion of the system
comm-mode = Linear
comm-grps = System
And finally, my md.mdp file (used in the third round with grompp and
mpirun, shell commands not shown):
; md.mdp
; This is the MDP file which will be used for the main
; run, after the completion of preprocessing.
; Run parameters
integrator = md
dt = 0.002 ; ps, which is 2.0 fs
nsteps = 500000 ; 1.0 ns. MAKE THIS A VARIABLE!
nstcomm = 5 ; frequency for center of mass motion removal,
; default is 10, in the original file it was 1.
; tinit = 0 ; NOT NEEDED I THINK, this is the default.
; HOWEVER this might be a useful variable if
; runs are extended?
; Trajectory (.trr) file output parameters
nstxout = 5000 ; positions, every 10 ps
nstvout = 5000 ; velocities
nstfout = 1000 ; forces, every 2 ps
nstenergy = 1000 ; energies
; Frequency to write coordinates to .xtc compressed trajectory file
nstxtcout = 1000 ; every 2 ps
; Periodic boundary conditions are on in all directions
pbc = xyz
; Single-range cutoff scheme
nstlist = 5 ; update the neighbor list every 5 steps
ns_type = grid
rlist = 1.4
rcoulomb = 1.4
rvdw = 1.4
; PME electrostatics parameters
; I don't understand these yet, I just copied them!
coulombtype = PME
fourierspacing = 0.12
fourier_nx = 0
fourier_ny = 0
fourier_nz = 0
pme_order = 4
ewald_rtol = 1e-5
optimize_fft = yes
; Nose-Hoover temperature coupling is on in two groups
Tcoupl = Nose-Hoover
tc_grps = Protein Non-Protein
ref_t = 310 310 ; 310 Kelvin = 37 C
tau_t = 0.2 0.2 ; picoseconds. Temperature coupling time
; constant. I changed this from 0.1, due to
; this warning: "For proper integration of
; the Nose-Hoover thermostat, tau_t (0.1)
; should be at least 20 times larger than
; nsttcouple*dt (0.01). dt = 2.0 fs, and
; I'm not changing that. nsttcouple is
; set to its default value, equal to
; nstlist when it is not specified.
; Pressure coupling is on
Pcoupl = Parrinello-Rahman
pcoupltype = isotropic
tau_p = 1.0 ; ps. Pressure coupling time constant.
compressibility = 4.5e-5 ; bar-1. should be a function of temperature?
; this is the compressibility at 300K, and
; I'm at 310K. Does that matter?
ref_p = 1.0 ; bar.
; Generate random velocities in the initial frame is off, useful for
; reproducing runs?
gen_vel = no
; Bond parameters
constraint_algorithm = lincs
constraints = all-bonds ; apparently you should always use the
; all-bonds option with LINCS.
continuation = yes ; do not apply constraints to the start
; configuration and do not reset shells.
; Useful for exact continuations, reruns.
; Long-range dispersion correction
DispCorr = EnerPres
I'm using GROMACS 4.5.4 on Linux. I have a multi-core CPU, so I use
mpirun. Simulation temperature is 37C, solvent is 0.1 M NaCl. The
proteins in question are soluble, so I'm using isotropic conditions for
all parameters.
I have changed as little as possible from the .mdp files which have
served me well in the past. I don't understand many of the parameters
and of course, it's best not to fiddle with what you don't understand.
As you can see from the notes I wrote in md.mdp, I got a novel error
message when I switched over from GROMACS 4.0, so I made what I hope was
a minor adjustment.
Convergence during the energy minimization (steep) preprocessing doesn't
appear to be a problem. I thought that was a good sign that I could
proceed safely, I guess not?
Thanks for any advice you can share.
More information about the gromacs.org_gmx-users
mailing list