[gmx-users] Re: Water molecules cannot be settled, why?
Tsjerk Wassenaar
tsjerkw at gmail.com
Wed Jul 11 10:32:35 CEST 2012
Hi John,
Check where the unsettling water molecule is placed. If it's in the
protein. that may be the cause of the problem. Otherwise, it's some of
the other stuff you're doing, but rule out the simple things first.
Cheers,
Tsjerk
On Wed, Jul 11, 2012 at 10:12 AM, John Ladasky
<blind.watchmaker at yahoo.com> wrote:
> 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.
>
> --
> gmx-users mailing list gmx-users at gromacs.org
> http://lists.gromacs.org/mailman/listinfo/gmx-users
> * Only plain text messages are allowed!
> * Please search the archive at
> http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
> * Please don't post (un)subscribe requests to the list. Use the www
> interface or send it to gmx-users-request at gromacs.org.
> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
--
Tsjerk A. Wassenaar, Ph.D.
post-doctoral researcher
Molecular Dynamics Group
* Groningen Institute for Biomolecular Research and Biotechnology
* Zernike Institute for Advanced Materials
University of Groningen
The Netherlands
More information about the gromacs.org_gmx-users
mailing list