[gmx-users] unable to equilibrate protein in membrane with NPT

[Justin A. Lemkul]

Peter C. Lai wrote:
> On 2011-04-08 08:16:19AM -0500, Justin A. Lemkul wrote:
>>
>> Peter C. Lai wrote:
>>> On 2011-04-08 08:02:14AM -0500, Justin A. Lemkul wrote:
>>>> Peter C. Lai wrote:
>>>>> Hello again
>>>>>
>>>>> In my protein-membrane-water-ion system (inserted via g_membed) I have run 
>>>>> a 1ns NVT equilibration with the protein restrained and now I am trying to 
>>>>> equilibrate with NPT and LINCS/mdrun is crashing after about 5000 iterations.
>>>>>
>>>>> Here is my NPT mdp file:
>>>>>
>>>>> define                  =-DPOSRES
>>>>> integrator               = md
>>>>> ; Start time and timestep in ps
>>>>> tinit                    = 0
>>>>> dt                       = 0.002
>>>>> nsteps                   = 500000
>>>>> init_step                = 0
>>>>> comm-mode                = Linear
>>>>> nstcomm                  = 1
>>>>> comm-grps                = Protein_POPC SOL_CL
>>>>> nstxout         = 100           ; save coordinates every 0.2 ps
>>>>> nstvout         = 100           ; save velocities every 0.2 ps
>>>>> nstenergy       = 100           ; save energies every 0.2 ps
>>>>> nstlog          = 100           ; update log file every 0.2 ps
>>>>> continuation    = yes
>>>>> constraint_algorithm = lincs    ; holonomic constraints
>>>>> constraints     = all-bonds     ; all bonds (even heavy atom-H bonds)
>>>>> lincs_iter      = 1             ; accuracy of LINCS
>>>>> lincs_order     = 4             ; also related to accuracy
>>>>> 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     = 0.8
>>>>> coulombtype     = PME
>>>>> pme_order       = 4             ; cubic interpolation
>>>>> fourierspacing  = 0.16          ; grid spacing for FFT
>>>>> tcoupl          = Nose-Hoover
>>>>> tc-grps         = Protein POPC SOL_CL
>>>>> tau-t                    = 0.5 0.5 0.5
>>>>> ref-t                    = 300 300 300
>>>>> gen-vel                  = no
>>>>> pcoupl          = Parrinello-Rahman         ; Pressure coupling on in NPT
>>>>> pcoupltype      = semiisotropic             ;
>>>>> tau_p           = 5.0                           ; time constant, in ps
>>>>> ref_p           = 1.01325 1.01325
>>>>> compressibility = 4.5e-5        4.5e-5
>>>>>
>>>>> Any suggestions?
>>>> Try using the Berendsen barostat.  P-R allows for wider oscillations that can 
>>>> lead to instability in incompletely equilibrated systems.
>>>>
>>>> -Justin
>>>>
>>> Switching to V-rescale and tau_t to 0.1 0.1 0.1 also crashes.
>> My recommendation was for the barostat, not the thermostat.
>>
>>> LINCS is complaining about C-H interactions in both the mobile lipid chains
>>> as well as atoms restrained by posre.itp
>>>
>>> Maybe I will run without LINCS all-bond?
>> No, that will make your life worse.  Just because the constraints fail doesn't 
>> mean that you need to disable LINCS.
>>
>> http://www.gromacs.org/Documentation/Errors#LINCS.2fSETTLE.2fSHAKE_warnings
>> http://www.gromacs.org/Documentation/Terminology/Blowing_Up
> 
> So basically, just keep running NVT until I can get NPT stabilized?
> Thanks.

At some point you need to equilibrate under NPT.  Indefinite use of NVT serves 
no purpose.  If your NVT conditions stabilized, then proceeding to NPT is 
logical.  If NVT didn't stabilize, then yes, backtrack and finish that 
equilibration properly.

The main point is if your system still has some unfavorable geometry somewhere, 
the application of a barostat that can cause big pressure oscillations is going 
to treat your system very nastily.  Some systems require a bit more finesse, so 
when it's time to introduce pressure coupling to the system, use a weak coupling 
(i.e. Berendsen) method so the system relaxes.  Then, once that's stable, switch 
to P-R to get a proper NPT ensemble.

-Justin

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

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

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