[gmx-users] 1-4 interactions free energy calculations

[Justin Lemkul]

On 6/24/13 1:07 PM, Sonia Aguilera wrote:
> Hi all!
>
> I have been running a series of simulations for coupling a molecule (free
> energy calculations). All my settings and parameters in the mdp file worked
> great for one system, but then I change the molecule to a bigger one and I
> started to have warnings about 1-4 interactions:
>
> Warning: 1-4 interaction between 1 and 143 at distance 2.579 which is larger
> than the 1-4 table size 2.500 nm
> These are ignored for the rest of the simulation
> This usually means your system is exploding,
> if not, you should increase table-extension in your mdp file
> or with user tables increase the table size
>
> I made severals tests to find the error: minimization, nvt, mdp and md
> without the Free Energy settings and It ran without the warning. So I guess
> the warning must be related to the FE implementation. However, I'm clueless
> about what can I change to stop having that warning. I already read other
> posts but I'm still the same. Should I just ignore it? Also, the initial
> structure is the result of a md simulation, so I don't know what can be
> wrong in the topology or mdp file.  Should I change the continuation
> parameter (continuation=yes) for the minimization and the nvt?
>
> Energy minimization converged but not to machine precision:
>
> Converged to machine precision,
> but not to the requested precision Fmax < 100
>
> Double precision normally gives you higher accuracy.
> You might need to increase your constraint accuracy, or turn
> off constraints alltogether (set constraints = none in mdp file)
>
> writing lowest energy coordinates.
>
> Steepest Descents converged to machine precision in 481 steps,
> but did not reach the requested Fmax < 100.
> Potential Energy  = -7.8070088e+05
> Maximum force     =  7.6036969e+02 on atom 68
> Norm of force     =  2.3164215e+01
>
>
> This is the copy of the mdp file (based os Justin Lemluk's tutorial), I'm
> using GROMACS 4.5.4.
>
> title                    = steepest descent energy minimization
> ; Run control
> integrator               = steep
> nsteps                   = 30000
> ; EM criteria and other stuff
> emtol                    = 100
> emstep                   = 0.01
> niter                    = 20
> nbfgscorr                = 10
> ; Output control
> nstlog                   = 1
> nstenergy                = 1
> ; Neighborsearching and short-range nonbonded interactions
> nstlist                  = 1
> ns_type                  = grid
> pbc                      = xyz
> periodic-molecules       = yes
> rlist                    = 1.5
> ; Electrostatics
> coulombtype              = PME
> rcoulomb                 = 1.5
> ; van der Waals
> vdw-type                 = switch
> rvdw-switch              = 0.8
> rvdw                     = 0.9
> ; Apply long range dispersion corrections for Energy and Pressure
> DispCorr                  = EnerPres
> ; Spacing for the PMEPPPM FFT grid
> fourierspacing           = 0.12
> ; EWALDPMEPPPM parameters
> pme_order                = 6
> ewald_rtol               = 1e-06
> epsilon_surface          = 0
> optimize_fft             = no
> ; Temperature and pressure coupling are off during EM
> tcoupl                   = no
> pcoupl                   = no
> ; Free energy control stuff
> free_energy              = yes
> init_lambda              = 1.0
> delta_lambda             = 0
> foreign_lambda           = 0.95
> sc-alpha                 = 0
> sc-power                 = 0
> sc-sigma                 = 0
> couple-moltype           = Protein_chain_A   ; name of moleculetype to
> decouple
> couple-lambda0           = vdw              ;  only vdw interactions
> couple-lambda1           = vdw-q            ; all interactions
> couple-intramol          = no
> nstdhdl                  = 10
> ; Generate velocities to start
> gen_vel                  = no
> ; options for bonds
> constraints              = none  ; we only have C-H bonds here
> ; Type of constraint algorithm
> constraint-algorithm     = lincs
> ; Do not constrain the starting configuration
> continuation             = no
> ; Highest order in the expansion of the constraint coupling matrix
> lincs-order              = 12
>

The problem almost certainly comes from the use of "couple-intramol = no," which 
for large molecules can cause issues (which have been reported before).  Larger 
molecules are more flexible, and thus the interaction distances can vary a lot. 
  The manual suggests using "couple-intramol = yes" for large molecules. 
Whether or not this is an appropriate setting for whatever you're trying to do 
is up to you and perhaps table-extension could be increased in this case, as it 
is not likely a crash but also something you cannot simply ignore.

-Justin

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

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

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