[gmx-users] Protein in vacum

[Justin Lemkul]

On 1/28/13 4:47 AM, Shima Arasteh wrote:
> Dear users,
>
>
> I have done simulation with my protein in water earlier successfully using CHARMM36 FF. Now, I am trying to simulate protein in vacum similar to the last simulation of protein-water. When I go to the NPT equilibration, I get a fatal error.
> The commands I ran:
>
> 1.pdb2gmx -f monomer.pdb -o monomer.gro  -ter
> 2.editconf -f monomer.gro -o monomer_newbox.gro -c -d 1.0 -bt cubic
> 3.grompp -f EM.mdp -c monomer_solv_ions.gro -p topol.top -o EM.tpr
> 4.mdrun -v -deffnm EM
> Steepest Descents converged to Fmax < 100 in 1632 steps
> Potential Energy  = -2.3169937e+02
> Maximum force     =  9.9520912e+01 on atom 360
> Norm of force     =  1.9172571e+01
>
>
> 5.grompp -f NVT.mdp -c EM.gro -p topol.top -o NVT.tpr
> 6.mdrun -deffnm NVT -v
>
> 7.grompp -f NPT.mdp -c NVT.gro -t NVT.cpt -p topol.top -o NPT.tpr
> 8.mdrun -v -deffnm NPT
> Fatal error:
> The X-size
> of the box (3.202905) times the triclinic skew factor (1.000000) is smaller than the number of DD cells (4) times the smallest allowed cell size (0.800718)
>
>
> The npt.mdp file settings are:
> ; Parameters describing neighbors searching and details about interaction calculations
> ns_type        = grid        ; Neighbor list search method (simple, grid)
> nstlist        = 5        ; Neighbor list update frequency (after every given number of steps)
> rlist        = 1.2        ; Neighbor list search cut-off distance (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     = 0.8
> ; Parameters for treating bonded interactions
> continuation    = no        ; Whether a fresh start or a continuation from a previous run (yes/no)
> constraint_algorithm = LINCS    ; Constraint algorithm (LINCS / SHAKE)
> constraints    = all-bonds    ; Which bonds/angles to constrain (all-bonds / hbonds / none / all-angles / h-angles)
> lincs_iter    = 1        ; Number of iterations to correct for rotational lengthening in LINCS (related to accuracy)
> lincs_order    = 4
>      ; Highest order in the expansion of the constraint coupling matrix (related to accuracy)
>
> ; Parameters for treating electrostatic interactions
> coulombtype    = PME        ; Long range electrostatic interactions treatment (cut-off, Ewald, PME)
> pme_order    = 4        ; Interpolation order for PME (cubic interpolation is represented by 4)
> fourierspacing    = 0.16        ; Maximum grid spacing for FFT grid using PME (nm)
>
> ; Temperature coupling parameters
> tcoupl        = Nose-Hoover            ; Modified Berendsen thermostat using velocity rescaling
> tc-grps        = Protein     ; Define groups to be coupled separately to temperature
> bath
> tau_t        = 0.5            ; Group-wise coupling time constant (ps)
> ref_t        = 310             ; Group-wise reference temperature (K)
>
> ; Pressure coupling parameters
> pcoupl        = Parrinello-Rahman        ; Pressure coupler used under NPT conditions
> pcoupltype    = isotropic            ; Isotropic scaling in the x-y direction, independent of the z direction
> tau_p        = 5.0                ; Coupling time constant (ps)
> ref_p        = 1.0                ; Reference pressure
> for coupling, x-y, z directions (bar)
> compressibility = 4.5e-5    4.5e-5        ; Isothermal compressibility (bar^-1)
> refcoord_scaling = com
> ; Miscellaneous control parameters
> ; Dispersion correction
> DispCorr    = EnerPres        ; Dispersion corrections for Energy and Pressure for vdW cut-off
> ; Initial Velocity Generation
> gen_vel        = no            ; Velocity is read from the previous run
> ; Centre of mass (COM) motion removal relative to the specified groups
> nstcomm        = 1            ; COM removal frequency (steps)
> comm_mode    = Linear        ; Remove COM translation (linear / angular / no)
> comm_grps    = Protein
>      ; COM removal relative to the specified groups
>
>
> Would you please help me? How would it be possible to solve my problem?
>

NPT is incompatible with a vacuum simulation.  The error message gives you a 
clue that the unit cell is rapidly shrinking as the system collapses in around 
the protein, leaving you with a pseudo-crystalline state that is full of 
periodicity artifacts.  Please consult the list archive for practical 
considerations for running in vacuo; the topic has been discussed several times, 
even within the last few days/weeks.

-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

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