simtk.org

Hi,

Is there a functionality, which would allow me to extract ALL force contributions (electrostatic, vdw, bonded, angle, etc.) acting on atom A, resulting from interactions with atom B?

If not, is it possible to change to code just a little bit (and not restructure the whole package), so that I could get the information I want?

Cheers,
MD

Hi Maciej,

Could you give more information about what you want to do? The force on atom A from only interactions with atom B isn't well defined, since some of the force field terms involve more than two atoms. For example, there could be several different angles that all involve atoms A and B, but each also involves a different third atom. Do you want all interactions that involve A and B, regardless of what other atoms they also involve? Or only the ones that involve no other atoms? Or reading between the lines, perhaps what you really want is the second derivative of the energy with respect to the positions of atoms A and B?

Also, how will this be used? For example, will you want to calculate this for all possible pairs of atoms? Or only for a single pair of atoms, but doing it repeatedly for many different configurations? Knowing a bit more about your application would help me to suggest an appropriate solution.

Peter

Sorry for the delay.

You asked the right questions, and to be honest, I am not quite sure yet what I want (at least not in terms of pairwise forces between A and B), but I do know how I might use this information later on.

I'll try to clarify:
I'm interested in all pairs of atoms, and how they act on each other (in terms of pairwise forces), for a given configuration.
Let's say we have atoms $A, B$ and $C_1,...,C_N$, and that the part of the potential involving $A$ and $B$ looks like this:
$$V_2(A,B) + \sum_i V_3(A,B,C_i) + \sum_{i>j} V_4(A,B,C_i,C_j)$$.
The force acting, say, on $A$ can be decomposed into:
$$F_2(A,B) + \sum_i F_3(A,B,C_i) + \sum_{i>j} F_4(A,B,C_i,C_j)$$.
These forces are partial derivatives of V with respect to the position of A, and these are the forces I am interested in.

Now, $F_2$ is easy -- it's the force acting on $A$, 'coming from' $B$. But I'm not quite sure what is the proper way of decomposing the $F_3(A,B,C_i)$ into contributions $A$--$B$ and $A$--$C_i$. The simplest way would be to say that 1/2 of this force comes from $B$, and the other 1/2 from $C_i$. This ensures additivity, but is it `proper'? I'm just in the middle of figuring this out.

I want to create a connected graph, in which nodes would correspond to atoms, and weights would somehow depend on these pairwise forces (and how they `influence' a given reaction coordinate).Then, I want to carry out a clustering of this graph, that would say "atoms from the group ONE tend to cooperate, and move the transition forward, whereas atoms in the group TWO do the opposite". Next, I want to operate on many such clustering, acquired for configurations produced through a simulation, spanning the whole range of the reaction coordinate. This signal will be noisy and will require some statistical magic, but that's the easy part. I think.

I'd suggest using the Jacobian. I think that tells you what you really care about: how does the position of each atom influence the force on each other atom? It's quite easy to calculate the whole matrix through finite differencing. At one point I even had some code to do that. I could try to find it, if you like.

Peter

Hi Peter,

The Jacobian that you are suggesting will give the matrix of second derivatives in the energy, i.e. the force constants. I think this is likely to produce a graph where 2, 3, and 4-tuples of bonded atoms might be connected, but also nonbonded contacts would be connected. Incidentally this is the first step towards doing a normal mode analysis to figure out which degrees of freedom are coupled to one another.

I imagine the original request might have been to extract the pairwise-additive parts of the nonbonded forces, at least that is how I interpreted it. The 3-body and 4-body terms in the force field cannot be decomposed into pairwise contributions in the sense that they cannot be expressed as the gradient of several pairwise-additive potentials.

If you have a constant force vector acting on several atoms (e.g. the vector function evaluated at a particular geometry), you might be able to decompose it into several pairwise terms, but I don't think this is unique (the number of pairwise terms is n(n+1)/2, where the number of atoms is n) and it won't work for any other geometry.

Thanks,

- Lee-Ping

But he doesn't actually need an additive decomposition of the force. He wants to do a clustering, and for that he just needs a measure of how atoms affect the forces on other atoms.

Peter

Yeah, grouping the atoms based on the Jacobian elements would be interesting. I'm not familiar with the literature here but it sounds like a reasonable thing to do.

The Jacobian - that's a good idea.
It would be interesting to cluster it, yes, but also for my purposes, it might suggest a reasonable decomposition of the multibody force into pairwise components.

I imagine the original request might have been to extract the pairwise-additive parts of the nonbonded forces, at least that is how I interpreted it. The 3-body and 4-body terms in the force field cannot be decomposed into pairwise contributions in the sense that they cannot be expressed as the gradient of several pairwise-additive potentials.
If you have a constant force vector acting on several atoms (e.g. the vector function evaluated at a particular geometry), you might be able to decompose it into several pairwise terms, but I don't think this is unique (the number of pairwise terms is n(n+1)/2, where the number of atoms is n) and it won't work for any other geometry.

Actually, I would like to get not only non-bonded forces, but also (among other) the force acting on A, as a result of interactions with B and C (for a 3-body term, as an example). In other words, dV(A,B,C)/dA, where the "d"s are partial derivatives. I can't decompose this into a pairwise force acting on A from C, but I don't want to. I will consider a dot product of this force with the gradient of the reaction coordinate. And from that I want to get the weights in my graph. My worry is that these weights are not definable in a unique way (which is related to the fact, that the pairwise components for multibody terms cannot be specified).

I've checked and it's impossible to get these multibody forces from NAMD. I've managed to get them from openbabel (which required a minor change in the C++ part of the program), but it's too slow to actually run a simulation.

In order to get a pairwise force decomposition from higher-body interactions, another method which might be helpful for you is the multiscale coarse graining (aka force matching) approach.

Basically, you fit an effective pair potential to the forces in all of your trajectory snapshots, and then your pairwise decomposition would come from simply evaluating the individual terms in the potential. The form of the pair potential can be highly flexible like a cubic spline, and you can have a different cubic spline for each pair of atom types. Greg Voth has applied this method extensively to build coarse grained models, and your problem could be an interesting application for this method. Since the forces are linear in the spline parameters, you can get the optimal fit without doing any nonlinear iterations; it is essentially a matrix inversion. The drawback is that since you have a pair potential that is trying to describe the higher-body interactions, the fit will not be perfect.

Incidentally, fitting a cubic spline pair potential using force matching is something my code (ForceBalance) should be able to do, but I haven't tried it yet.

Thanks,

- Lee-Ping

That might be an interesting solution.
I've just downloaded forcebalance-1.1, I'll give it a spin.