[gmx-users] restraints on water oxygen atoms

[Shima Arasteh]

Thanks for your reply.
 
The system I am trying to equilibrate is composed of popc/ peptide/ions/water. I built the system by InflateGRO methodology as you wrote in kalp-dppc tutorial. The last gro file which gave me the acceptable APL, was used as the starting configuration in next steps. But before going on, I found some overlaps between lipids acyl chains and peptide structure. So tried to move the  lipid residues caused the overlaps. Then there are some gaps between peptide and near lipid chains. As a result, I see the water cleavage in these regions.
Would you please give me any suggestions? How would I make a better packed system without disturbing overlaps or crashes?
Would you please help me?

Thanks for your suggestions and help. 


 
Sincerely,
Shima


----- Original Message -----
From: Justin Lemkul <jalemkul at vt.edu>
To: Shima Arasteh <shima_arasteh2001 at yahoo.com>; Discussion list for GROMACS users <gmx-users at gromacs.org>
Cc: 
Sent: Wednesday, June 5, 2013 9:25 PM
Subject: Re: [gmx-users] restraints on water oxygen atoms



On 6/4/13 11:47 PM, Shima Arasteh wrote:
> Dear gmx users,
>
> I have a POPC/peptide/water/ions system. I ran NVT and then NPT on my system. I'd prefer to run the equilibrium steps with position restraints on water oxygen atoms, because the water molecules penetrate the lipid bilayer when running the equilibrium and I don't want it to happen.
> I tried the NVT with position restraints on water by adding -DPOSRES_WATER to nvt.mdp file and editing the top file by changing the fc to 100000.
>
> #ifdef POSRES_WATER
> ; Position restraint for each water oxygen
> [ position_restraints ]
> ;  i funct       fcx        fcy        fcz
>     1    1       100000       100000       100000
> #endif
>
> This edition turned into a better result.
>
> Now I tried to put such a restraint on npt but the gromacs does not allow it by turning a fatal error:
> A charge group moved too far between two domain decomposition steps.
>
>
> npt.mdp file:
> ;NPT equlibration Dimer-POPC-Water - CHARMM36
> define        = -DPOSRES_LIPID -DPOSRES -DPOSRES_WATER    ; P headgroups of POPC and Protein is position restrained (uses the posres.itp file information)
> integrator      = md            ; leap-frog integrator
> nsteps          =1000000         ; 1 * 1000000 = 1000 ps
> dt              = 0.001         ; 1 fs
> ; Output control
> nstxout         = 2000           ; save coordinates every 0.2 ps
> nstvout         = 1000           ; save velocities every 0.2 ps
> nstenergy       = 1000           ; save energies every 0.2 ps
> nstlog          = 1000           ; update log file every 0.2 ps
>
> continuation    = yes            ; first dynamics run
> constraint_algorithm = lincs    ; holonomic constraints
> constraints     = h-bonds     ; all bonds (even heavy atom-H bonds) constrained
> lincs_iter      = 1             ; accuracy of LINCS
> lincs_order     = 4             ; also related to accuracy
> ; Neighborsearching
> ns_type         = grid          ; search neighboring grid cells
> nstlist         = 5             ; 10 fs
> rlist           = 1.2           ; short-range neighborlist cutoff (in nm)
> rlistlong       = 1.4
> rcoulomb        = 1.2           ; short-range electrostatic cutoff (in nm)
> rvdw            = 1.2           ; short-range van der Waals cutoff (in nm)
> vdwtype         = switch
> rvdw_switch     = 1.0
> ; 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          = Nose-Hoover     ; modified Berendsen thermostat
> tc-grps         = Protein_POPC Water_Ces_CL        ; two coupling groups - more accurate
> tau_t           = 0.5   0.5       ; time constant, in ps
> ref_t           = 310   310    ; reference temperature, one for each group, in K
> pcoupl          = Berendsen            ; no pressure coupling in NVT
> pcoupltype      = semiisotropic
> tau_p           = 4
> ref_p           = 1.01325 1.01325
> compressibility = 4.5e-5 4.5e-5
>
> ; Periodic boundary conditions
> pbc             = xyz           ; 3-D PBC
> ; Velocity generation
> gen_vel         = no           ; assign velocities from Maxwell distribution
> ;gen_temp        = 310           ; temperature for Maxwell distribution
> ;gen_seed        = -1            ; generate a random seed
> nstcomm         = 1
> comm_mode       = Linear
> comm_grps       = Protein_POPC Water_Ces_CL
> refcoord_scaling    = com
>
>
>
> I am wondering how it is possible to prevent penetrating the water molecules through equilibrium? And how I can put the restraint in npt as well as nvt? Would you please help me in this issue please?
>

Restraints in all directions, especially with pressure coupling, will 
undoubtedly lead to nasty atomic clashes and the failure that you're seeing. 
Normally, one does not need to apply restraints in any other dimension aside 
from z, to prevent "vertical" diffusion into the membrane.  The structure of the 
membrane should be reasonable after NVT; if it's not, it suggests a poor 
starting configuration.  Occasional water leakage is one thing, having to 
completely restrain water molecules to prevent it from happening suggests 
something is very bad.

-Justin

-- 
========================================

Justin A. Lemkul, Ph.D.
Research Scientist
Department of Biochemistry
Virginia Tech
Blacksburg, VA
jalemkul[at]vt.edu | (540) 231-9080
http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin

========================================