Density-dependent potential

[Nicholas Schafer]

Hi Peter -

Based on my reading of the CustomGBForce documentation, I don't think that it is suitable for a particular type of interaction that I am interested in, but I wanted to check with you to see if my understanding is correct.

I want to implement a pairwise interaction V_ij that depends on both the distance between the atoms and a pairwise quantity that, in turn, depends on the local densities of the interacting atoms, such that you have something like this:

Code: Select all

V_ij(sigma_ij(rho_i({r}), rho_j({r})), r_ij)

The way I guess this would be done would be to compute rho for each atom, which depends on the position of all of the other atoms in the system, then compute sigma for each pair of particles, then compute the energy. I think that CustomGBForce is not suitable in at least two ways, however. First, single particle quantities (like rho) are not allowed to depend on the positions of the other particles. Second, pair quantities like sigma must be computed before single particle quantities like rho.

Does that sound right to you, i.e., CustomGBForce will not be suitable for what I want to do?

Thanks,
Nick

[Peter Eastman]

I can't tell from your description. What are the various functions? Especially, how does rho depend on the positions of other particles? Can it be formulated as a sum of pairwise terms? Also, does sigma really need to be computed and saved separately, or can it just be made part of the definition of V?

[Nicholas Schafer]

Thanks for helping me try to figure this out. The details are as follows:

Code: Select all

theta(rmin, rmax, r_ij)=0.25(1+tanh[eta*(r_ij-rmin)])(1+tanh[eta(rmax-r_ij)]) # a smooth step function = 1 between rmin and rmax
rho_i=\sum_j theta(.45, .65, r_ij) # the local density of residue i
sigma_ij=0.25(1-tanh([eta*(rho_i-rho0)])(1-tanh[eta(rho_j)])*gamma_ij # a product of step functions involving the local densities multiplied by a residue-pair specific weighting parameter
V_ij=-k*theta(.65,.95,r_ij)*sigma_ij # The potential

rho_i is a sum of functions that depend on all r_ij. sigma_ij could almost be included as part of V except that it has to be multiplied by a residue pair-specific weighting parameter that depends on both i and j.

What do you think? Still possible?

[Peter Eastman]

I don't see any difficulties with implementing that. Use a particle pair computed value for rho, and a particle pair energy term for V. Is gamma the weighting parameter you mentioned? You can use a lookup table for that. Store it in a Discrete2DFunction, and have a per-particle parameter containing the index into the table for each particle.

[Nicholas Schafer]

That's great! I'm going to try it out.

Could a Discrete2DFunction also be used to exclude some pairs of residues from the density calculation based on their sequence separation?

[Peter Eastman]

Sure. Create a function that's 0 for ones that should be excluded, 1 for others, and multiply by it. Though exclusions might be a more efficient way of doing the same thing.

[Nicholas Schafer]

The density has a different sequence separation criterion than the pairwise interaction. I plan to use the exclusions to handle the fact that the interaction should only be applied to pairs of Cbeta atoms, but the density calculation includes pairs of residues that are not included in the pairwise interaction. In that case, does using exclusions for the pairwise interactions and a Discrete2DFunction for the density calculation make sense?

[Peter Eastman]

Makes sense.

[Nicholas Schafer]

This is great. I think that I got the density-dependent 1- and 2-body terms working. I'm also interested in implementing a density-dependent hydrogen bonding potential. The way this is supposed to work is that there are interactions between all pairs of residues (i, i+4) and the energy is the product of a density-dependent term and a term that depends on two pair distances from a triplet of atoms, one in residue i and two in residue i+4.

Code: Select all

V(i,i+4)=sigma(i,i+4)*f(r_{Oi->Ni+4}, r{Oi->Hi+4})

where, similar to before

Code: Select all

sigma_ij=0.25(1-tanh([eta*(rho_i-rho0)])(1-tanh[eta(rho_j)])
rho_i=\sum_j theta(.45, .65, r_ij) # the local density of residue i
theta(rmin, rmax, r_ij)=0.25(1+tanh[eta*(r_ij-rmin)])(1+tanh[eta(rmax-r_ij)]) # a smooth step function = 1 between rmin and rmax

Having the function "f" depend on two pair distances helps to ensure directionality as well as proximity, but I think that it also puts it just beyond the reach of CustomGBForce because it can't be expressed as a pairwise energy term anymore. Does that sound right to you?

I was thinking that I could maybe make this work, in principle, as a CustomCompoundBondForce and including all of the (required) particles in the system. Do you think that building up a huge string to compute rho_i and rho_i+4 would break the expression parsing?

Can you maybe think of a better way to do this?

Much appreciated,
Nick

[Peter Eastman]

You won't break it, but it could be very slow. CustomCompoundBondForce is designed for cases where each bond only involves a small number of particles. As the number gets large, it becomes very inefficient.

Unfortunately, none of the custom forces is really designed to handle this case. The "right" way to do it is to write a plugin implementing the force you want. That will be a lot more work though, especially if you want to support GPUs.