[gmx-users] FEP and equilibration problems

[João Martins]

Hi,

I'm having problems with running a FEP decoupling procedure for my POPE,
protein and ligand system. I'm following the Alchemistry.org Gromacs 4.6
tutorial (
http://www.alchemistry.org/wiki/GROMACS_4.6_example:_n-phenylglycinonitrile_binding_to_T4_lysozyme)
while changing the necessary parameters for my system, but while the
restraint-lambdas and coul-lambda steps are working correctly, the
vdw-lambdas steps are presenting some problems.

My procedure is as follows:

Starting with a pre-equilibrated membrane + protein + ligand system, I run
a water-restrained minimization step for the specific init-lambda-state.

free-energy              = yes
> sc-alpha                   = 0.5
> sc-power                  = 1
> sc-sigma                 = 0.3
> init-lambda-state        = 23
> restraint-lambdas        = 0.0 0.01 0.025 0.05 0.075 0.1 0.2 0.35 0.5 0.75
> 1.0 1.00 1.0 1.00 1.0 1.00 1.0 1.0 1.0 1.0 1.0 1.0 1.00 1.0 1.00 1.0 1.00
> 1.0 1.00 1.0
> coul-lambdas             = 0.0 0.00 0.000 0.00 0.000 0.0 0.0 0.00 0.0 0.00
> 0.0 0.25 0.5 0.75 1.0 1.00 1.0 1.0 1.0 1.0 1.0 1.0 1.00 1.0 1.00 1.0 1.00
> 1.0 1.00 1.0
> vdw-lambdas              = 0.0 0.00 0.000 0.00 0.000 0.0 0.0 0.00 0.0 0.00
> 0.0 0.00 0.0 0.00 0.0 0.05 0.1 0.2 0.3 0.4 0.5 0.6 0.65 0.7 0.75 0.8 0.85
> 0.9 0.95 1.0
> nstdhdl                  = 10
> calc-lambda-neighbors    = -1


The following step is the actual production run with the same FEP
parameters and for the same init-lambda-state. While this works correctly
for any lambda relating with restraints and coulombic forces, I'm getting
errors in steps 22 to 27. The type of error depends on the run, these two
being the ones I'm getting:

Fatal error:
> 1 particles communicated to PME node 11 are more than 2/3 times the
> cut-off out of the domain decomposition cell of their charge group in
> dimension x.
> This usually means that your system is not well equilibrated.
> For more information and tips for troubleshooting, please check the GROMACS
> website at http://www.gromacs.org/Documentation/Errors


and

Fatal error:
> 2 of the 518891 bonded interactions could not be calculated because some
> atoms involved moved further apart than the multi-body cut-off distance (1
> nm) or the two-body cut-off distance (1 nm), see option -rdd, for pairs and
> tabulated bonds also see option -ddcheck
> For more information and tips for troubleshooting, please check the GROMACS
> website at http://www.gromacs.org/Documentation/Errors


The first one seems to stop appearing when I decrease the time-step to 1 fs
from 2 fs. However, I can't seem to overcome this second error, which has
left me stumped so far. For the sake of completeness, this is the
production MDP I'm using:

define                   =
> integrator               = md
> dt                       = 0.001
> nsteps                   = 1000000
> nstxout                  = 5000
> nstvout                  = 5000
> nstlog                   = 5000
> nstenergy                = 5000
> nstxtcout                = 5000
> nstcalcenergy            = 5
> nstlist                  = 5
> nstcomm                  = 5
> comm_mode                = Linear
> comm-grps                = Water_and_ions Protein_&_!AMA_Lipid Ligand
> ns_type                  = grid
> rlist                    = 1.0
> rlistlong                = 1.4
> rvdw_switch              = 0.8
> vdwtype                  = Switch
> coulombtype              = pme
> rcoulomb                 = 1.0
> rcoulomb_switch          = 0.0
> rvdw                     = 1.0
> fourierspacing           = 0.15
> pme_order                = 6
> tcoupl                   = V-rescale
> nhchainlength            = 1
> tc-grps                  = POPE Water_and_ions_Ligand Protein_&_!Ligand
> tau_t                    = 0.1  0.1  0.1
> ref_t                    = 310.15  310.15  310.15
> Pcoupl                   = berendsen
> Pcoupltype               = semiisotropic
> tau_p                    = 5.0
> compressibility          = 4.5e-5       4.5e-5
> ref_p                    = 1.0          1.0
> refcoord-scaling         = all
> pbc                      = xyz
> gen_vel                  = yes
> gen_temp                 = 310.15
> optimize_fft             = yes
> constraints              = all-bonds
> continuation             = no
> constraint_algorithm     = Lincs
> lincs-order              = 4            ; 8 is needed for BD with large
> time-steps.
> lincs-iter               = 1            ; 1 is fine for normal
> simulations, but use 2 to conserve energy in NVE runs.
> ; Free energy control stuff =
> free-energy              = yes
> sc-alpha                 = 0.5
> sc-power                 = 1
> sc-sigma                 = 0.3
> init-lambda-state        = 23
> restraint-lambdas        = 0.0 0.01 0.025 0.05 0.075 0.1 0.2 0.35 0.5 0.75
> 1.0 1.00 1.0 1.00 1.0 1.00 1.0 1.0 1.0 1.0 1.0 1.0 1.00 1.0 1.00 1.0 1.00
> 1.0 1.00 1.0
> coul-lambdas             = 0.0 0.00 0.000 0.00 0.000 0.0 0.0 0.00 0.0 0.00
> 0.0 0.25 0.5 0.75 1.0 1.00 1.0 1.0 1.0 1.0 1.0 1.0 1.00 1.0 1.00 1.0 1.00
> 1.0 1.00 1.0
> vdw-lambdas              = 0.0 0.00 0.000 0.00 0.000 0.0 0.0 0.00 0.0 0.00
> 0.0 0.00 0.0 0.00 0.0 0.05 0.1 0.2 0.3 0.4 0.5 0.6 0.65 0.7 0.75 0.8 0.85
> 0.9 0.95 1.0
> nstdhdl                  = 10
> calc-lambda-neighbors    = -1
> pull           = umbrella
> pull_geometry  = distance
> pull_dim       = Y Y Y
> pull_start     = no
> pull_init1     = 0.610
> pull_ngroups   = 1
> pull_group0    = C10_Ligand
> pull_group1    = Ca_SER
> pull_k1        = 0.0   ; kJ*mol^(-1)*nm^(-2)
> pull_kB1       = 4184  ; kJ*mol^(-1)*nm^(-2)



Does anyone have any recommendations on making this calculation work?

Best,
João Martins

Best Regards,

*Joao Martins*

PhD Fellow


Structural Biology and NMR laboratory (SBiNLab)



*University of Copenhagen*,

*Faculty of Science*

Department of Biology

Ole Maaløes Vej 5

2200 Copenhagen N

Denmark


TEL +45 35 32 37 10

DIR +45 35 32 22 59


*joao.martins at bio.ku.dk <serlendsson at bio.ku.dk>*