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

[Abhishek Acharya]

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.
-------------------------------------------------------
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.

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

Sincerely,
Abhishek Acharya