[gmx-users] Getting non-bonded forces from GROMACS

[Justin Lemkul]

On 9/24/19 4:34 PM, James wrote:
> Hi Mark,
>
> Thanks a lot for the input. I understand what you mean about additive force
> fields, or at least I think I do, so let me summarize what I think is the
> point:
>
> In an additive force field, you have a lot of parameters that can be tuned,
> and the objective is to tune them in such a way that they reproduce
> experimental (or higher level of theory) results. Since adjusting one
> parameter can affect the tuning of others, you can't discretely map the
> force field parameters to real-world physical phenomena. So, attempting to
> decompose an additive force field to get, for example, VDW forces, is, at
> least mathematically speaking, incorrect. If that is the idea, let me make
> a couple observations, and then ask a question.
>
> 1) The force field does have parameters whose mathematical forms are chosen
> to be representative of real-world phenomenon, such as LJ potentials.
>
> 2) In some sense, some force field parameters may actually be used
> "separately" in certain simulation settings -- which, practically speaking,
> seems like an example of decomposition. For example, beyond the cutoff
> distance, VDW plays no role. But, if you bring two structures close to each
> other they experience VDW forces and react appropriately (or so we hope).
> The bonded forces didn't change (did they?), so it seems like if you can
> model something without VDW, and something with VDW, then it is legitimate
> to consider VDW a separate force that can be added or removed as
> appropriate.

You can't assume two species will simply interact as rigid bodies. As 
two molecules approach one another, all atoms within the cutoff will 
feel the LJ (and other) forces between the molecules. This can, and 
often does, change the internal structure of each of the molecules, thus 
affecting bonded forces.

The other instance in which decomposing forces this way doesn't really 
make sense are the 1-4 interactions. Many force fields apply a scaling 
factor to 1-4 LJ and electrostatic interactions, due to the fact that 
the molecular mechanics approximations describing such interactions are 
not very accurate for close distances. Dihedral terms are added 
effectively as a mathematical Band-Aid to balance out the forces (this 
is true for force fields that do not scale their 1-4 interactions). So 
in your proposal, you would be considering 1-4 interactions (because 
they are, in fact, nonbonded interactions) as part of your force 
evaluation, but those forces would not be offset by the corresponding 
dihedral term(s) that actually counterbalance these forces (and 
therefore conformational sampling) of your molecule. If you say you'll 
allow dihedral forces to be included to account for this phenomenon, 
then what about the angles that bend and the bonds that stretch as a 
function of that dihedral rotation and steric clashes among the relevant 
atoms? Everything is fundamentally linked.

> So, my question is this: While I understand that it may not be
> mathematically rigorous, does decomposing the force field really give
> meaningless answers (i.e., does not "correlate with something physical")? I
> could see it not being as accurate as desired, but for it to be useless
> seems strange to me.

There are no target data, theoretical or empirical, that can give you 
anything with which to compare. Sure, you can technically compute the 
quantity, but you have no context for it, and as the above examples 
demonstrate, you can't assume that the number is valid even in the 
context of the force field.

-Justin

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

Justin A. Lemkul, Ph.D.
Assistant Professor
Office: 301 Fralin Hall
Lab: 303 Engel Hall

Virginia Tech Department of Biochemistry
340 West Campus Dr.
Blacksburg, VA 24061

jalemkul at vt.edu | (540) 231-3129
http://www.thelemkullab.com

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