[gmx-users] Re: FEP problem

[David Mobley]

Dear Qin,

I believe the GROMACS manual recommends doing constant volume
equilibration prior to doing constant pressure equilibration. This is
what I always do. My understanding is that the Berendsen barostat
tends to be rather "picky" and crash (as you're seeing) when the
system is changing rapidly, as typically happens in the early stages
of equilibration. So typically one needs to run a very small amount
(i.e. 10 ps) of constant volume equilibration prior to going to
constant pressure or one gets exactly the sorts of problems you're
seeing.

I am CC-ing the GROMACS users list (which is a better place for
questions like this) for the record.

Best wishes,
David


On 10/15/07, Qin Wang <gwangq at gmail.com> wrote:
> Dear Sir,
>        I am a grad student from the Goran Krilov's group in Boston College,
> and when we asked questions to Prof. Michael Shirts, he recommended you, who
> is doing great work in GROMACS, to us and here I ask you for help. Thank you
> very much.
>        Here is the problem.
>       I am trying to do a FEP. In the equilibrium step, it fails after
> about 500 steps with a warnning:  *pressure scaling more than 1%*, mu:
> 1.06718 1.06718 1.06718. I searched the old discussion about it and then
> try to increase the tau_p to 2 till 10, or use Berendsen temperature
> coupling (with tau_t=0.3). And it just goes futher but still meets the
> problem at the end.
>
>     The system is a protein binding with one ligand and the ligand
> mutates to another one. The problem is found when lamda=1. First of all
> I run a energy minimization using an initial structure from the product
> of lamda=0.9, or even lamda=0.95. Then the problem happens when I run a
> equilibrium.
> Here is the mdp option:
> cpp                      =/usr/bin/cpp
> integrator               = md
>
>
> tinit                    = 0
> dt                       = 0.002
> nsteps                   = 50000
> nstcomm                  = 1
> nstxout                  = 0
> nstvout                  = 0
> nstfout                  = 0
> nstlog                   = 500
> nstenergy                =100
>
> nstxtcout                = 5000
> xtc-precision            = 1000
> nstlist                  = 1
> ns_type                  = grid
>
>
> Pcoupl                   = berendsen
> tau_p                    = 1.0
> compressibility          = 4.5e-05
> ref_p                    = 1.0
>
>
> tcoupl                   = nose-hoover
> tc_grps                  = system
> tau_t                    = 0.1
> ref_t                    = 300
>
> constraints              = none
> constraint-algorithm     =lincs
>
> shake-tol                = 0.0001
> lincs-order              = 12
> lincs-warnangle          =30
>
> coulombtype              =pme
>
> rcoulomb                 =1.0
>
> epsilon-r                =1
> vdw-type                 =switch
> rvdw-switch              = 1.0
> rvdw                     = 1.3
> DispCorr                  =EnerPres
>
> fourierspacing           =0.08
>
> fourier_nx               = 0
> fourier_ny               = 0
> fourier_nz               = 0
>
>
> pme_order                = 6
> ewald_rtol               = 1e-05
> epsilon_surface          = 0
> optimize_fft             = yes
>
>
> free_energy              = yes
> init_lambda              = 1.0
> delta_lambda             = 0
> sc_alpha                 = 0.5
> sc-power                 =1.0
> sc-sigma                 = 0.3
>
>         In purpose to solve the problem, I was trying to calculate the
> free energy difference from PTEN-LIG4 to PTEN-LIG6 scaling at lamda=0.0
> instead of Lig6 to Lig 4 at lamda=1.0. But I met the same problem.
>          Thank you very much if you have any ideas on how to solve the
> problem.
>
>          Best wishes,
>
>  Sincerely,
>  Qin Wang
>
>
>
>