[gmx-users] Free Energy Calculations: Error during minimization step.

[Wes Barnett]

On Sat, Sep 9, 2017 at 2:09 AM, Abhishek Acharya <abhi117acharya at gmail.com>
wrote:

> Dear GROMACS users,
>
> I am trying to estimate the free energy of solvation for a ion. But, I am
> facing a problem while running simulations for the deltaG_LJ calculation.
> The simulation at vdw_lambda=1.0 crashes during the steepest-descent
> minimization step, with the following error.
>
> -------------------------------------------------------
> Program gmx mdrun, VERSION 5.1.2
> Source code file:
> /home/bp-lab/Downloads/gromacs-5.1.2/src/gromacs/mdlib/constr.cpp, line:
> 555
>
> Fatal error:
>
> step 10: Water molecule starting at atom 120 can not be settled.
> Check for bad contacts and/or reduce the timestep if appropriate.
>
> For more information and tips for troubleshooting, please check the GROMACS
> website at http://www.gromacs.org/Documentation/Errors
> -------------------------------------------------------
>
> The free energy parameters I used for this run are as follows:
> ************************************************************
> *************************************************************
> free_energy              = yes
> couple-moltype           = Protein
> couple-lambda0           = vdw
> couple-lambda1           = none
> couple-intramol          = no
> init-lambda              = -1
> init_lambda_state        = 20
> delta_lambda             = 0
> nstdhdl                  = 10
> fep-lambdas              =
> mass_lambdas             = 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
> 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
> coul_lambdas             = 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
> 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
> vdw_lambdas              = 0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35 0.40
> 0.45 0.50 0.55 0.60 0.65 0.70 0.75 0.80 0.85 0.90 0.95 1.00
> bonded_lambdas           = 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
> 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
> restraint_lambdas        = 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
> 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
> temperature_lambdas      = 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
> 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
> calc_lambda_neighbors    = 1
> init-lambda-weights      =
> dhdl-print-energy        = no
> sc-alpha                 = 0.5
> sc-power                 = 1
> sc-r-power               = 6
> sc-sigma                 = 0.3
> sc-coul                  = no
> separate-dhdl-file       = yes
> dhdl-derivatives         = yes
> dh_hist_size             = 0
> dh_hist_spacing          = 0.1
> ************************************************************
> *************************************************************
>
> Out of curiosity, I also ran simulations for a reverse transformation with
> the following parameters (and I hope they are correct):
>
> ************************************************************
> *************************************************************
> free_energy              = yes
> couple-moltype           = Protein
> couple-lambda0           = none
> couple-lambda1           = vdw
> couple-intramol          = no
> init-lambda              = -1
> init_lambda_state        = 0
> delta_lambda             = 0
> nstdhdl                  = 10
> fep-lambdas              =
> mass_lambdas             = 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
> 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
> coul_lambdas             = 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
> 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
> vdw_lambdas              = 0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35 0.40
> 0.45 0.50 0.55 0.60 0.65 0.70 0.75 0.80 0.85 0.90 0.95 1.00
> bonded_lambdas           = 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
> 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
> restraint_lambdas        = 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
> 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
> temperature_lambdas      = 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
> 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
> calc_lambda_neighbors    = 1
> init-lambda-weights      =
> dhdl-print-energy        = no
> sc-alpha                 = 0.5
> sc-power                 = 1
> sc-r-power               = 6
> sc-sigma                 = 0.3
> sc-coul                  = no
> separate-dhdl-file       = yes
> dhdl-derivatives         = yes
> dh_hist_size             = 0
> dh_hist_spacing          = 0.1
> ************************************************************
> *************************************************************
> Intriguingly, I this case I got the exact error as above, but this time at
> lambda=0.
>

That make sense, since in both cases the vdw interactions are off in both
situations (couple-lambda0 or couple-lambda1 = none).



> -------------------------------------------------------
> Program gmx mdrun, VERSION 5.1.2
> Source code file:
> /home/bp-lab/Downloads/gromacs-5.1.2/src/gromacs/mdlib/constr.cpp, line:
> 555
>
> Fatal error:
>
> step 10: Water molecule starting at atom 120 can not be settled.
> Check for bad contacts and/or reduce the timestep if appropriate.
>
> For more information and tips for troubleshooting, please check the GROMACS
> website at http://www.gromacs.org/Documentation/Errors
> -------------------------------------------------------
>
> LINCS has been used to constrain all-bonds with lincs_order = 12.
>
> Analysis of step10b.pdb and step10c.pdb did not reveal any obvious clashes
> between water molecules.
>
> I have also run simulations for charge annihilation step and I did not face
> such a problem.
>

Has the charge been turned off for this series of simulations? You should
generally turn electrostatics off through a series of free energy
simulations and then turn off VDW while the charge is off.


>
> Looking for the some advice on how to proceed from here.
>

I do something like this:

vdw-lambdas              = 0.00 0.00 0.00 0.00 0.00 0.10 0.20 0.30 0.40
0.50 0.60 0.70 0.80 0.90 1.00
coul-lambdas             = 0.00 0.25 0.50 0.75 1.00 1.00 1.00 1.00 1.00
1.00 1.00 1.00 1.00 1.00 1.00
couple-lambda0           = vdw-q
couple-lambda1           = none

See:

http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin/gmx-tutorials/free_energy/index.html
http://wbarnett.us/tutorials/4_methane_fe/

-- 
James "Wes" Barnett
Postdoctoral Research Scientist
Department of Chemical Engineering
Kumar Research Group <http://www.columbia.edu/cu/kumargroup/>
Columbia University
w.barnett at columbia.edu
http://wbarnett.us