Page 1 of 2
Density-dependent potential
Posted: Thu Jul 20, 2017 9:34 am
by npschafer
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
Re: Density-dependent potential
Posted: Fri Jul 21, 2017 9:43 am
by peastman
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?
Re: Density-dependent potential
Posted: Fri Jul 21, 2017 10:55 am
by npschafer
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?
Re: Density-dependent potential
Posted: Fri Jul 21, 2017 11:24 am
by peastman
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.
Re: Density-dependent potential
Posted: Fri Jul 21, 2017 11:43 am
by npschafer
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?
Re: Density-dependent potential
Posted: Fri Jul 21, 2017 11:46 am
by peastman
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.
Re: Density-dependent potential
Posted: Fri Jul 21, 2017 11:49 am
by npschafer
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?
Re: Density-dependent potential
Posted: Fri Jul 21, 2017 11:50 am
by peastman
Makes sense.
Re: Density-dependent potential
Posted: Sat Jul 22, 2017 7:47 am
by npschafer
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
Re: Density-dependent potential
Posted: Mon Jul 24, 2017 8:55 am
by peastman
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.