[gmx-users] simulating protein in gas-phase

[Justin Lemkul]

On 1/6/13 3:09 PM, Sanku M wrote:
> Dear Gromacs users,
>    I am trying to perform a MD simulation in gromacs 4.5.4 using a protein and ligand in gas phase.
> I have previously run the simulation in water without any problem. Now, I want to compare the result to see the effect of solvent.
>
> But for running the same simulation after stripping the waters from the simulation box I am looking for some advise.
> I am having following issues:
> a) If I use 'md' i.e usual leap frog integrator with nstcomm=10 and comm-mode=Linear( this was the setup I had when I ran the simulation in water), I find after running the system of protein and ligand in gas phase for some time ( about 1 ns ), the protein starts rotating i.e it generates a very high angular momentum.
> b) So, I thought of using comm-mode=angular,,, But then grompp gives me a warning that removing rotation is not a problem only if I have only 1 molecule in the system ...which is not the case here as I have ligand and some ions . Any suggestion is highly appreciated. I can change nstcomm=1 using comm-mode=Linear but I am not sure it will stop the rotation of the protein about its own axis.

Use comm-mode = angular.

> I am running the system using periodic boundary condition and using PME for electrostatics.
>

You're certain to get periodicity artifacts.  You should run without PBC and 
with infinite cutoffs for nonbonded interactions to simulate the molecules in vacuo.

-Justin

> c) another option I thought of was using 'sd' integrator but considering my simulation in water having been performed using md integrator, I was looking to keep my .mdp file option as similar as possible in both cases.
>
> Any help on how to perform the simulations in gas-phase will be highly appreciated.
>
> Here is my current .mdp options:
> title       = Umbrella pulling simulation
> ; Run parameters
> integrator  = md
> dt          = 0.002
> tinit       = 0
> nsteps      = 1500000     ; 800 ps
> nstcomm     = 10
> ; Output parameters
> nstxout     = 0000      ; every 10 ps
> nstvout     = 0000
> nstfout     = 0000
> nstxtcout   = 250
> nstenergy   = 250
> ; Bond parameters
> constraint_algorithm    = lincs
> constraints             = hbonds
> continuation            = yes
> ; Single-range cutoff scheme
> nstlist     = 5
> ns_type     = grid
> rlist       = 1.4
> rcoulomb    = 1.4
> rvdw        = 1.4
> ; PME electrostatics parameters
> coulombtype     = PME
> fourierspacing  = 0.12
> fourier_nx      = 0
> fourier_ny      = 0
> fourier_nz      = 0
> pme_order       = 4
> ewald_rtol      = 1e-5
> optimize_fft    = yes
> ; Berendsen temperature coupling is on in two groups
> Tcoupl      = Nose-Hoover
> tc_grps     =  System
> tau_t       = 0.5
> ref_t       = 300
> ; Pressure coupling is on
> Pcoupl          = no
> pcoupltype      = isotropic
> tau_p           = 1.0
> compressibility = 4.5e-5
> ref_p           = 1.0
> ; Generate velocities is on
> gen_vel     =  yes
> ; Periodic boundary conditions are on in all directions
> pbc     = xyz
> ; Long-range dispersion correction
> DispCorr    = EnerPres
>
> Thanks
> Sanku
>

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

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

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