Could you explain more about what you're doing? What are you using this potential to model? When you say, "the number of I's and J's will change," do you mean during a single simulation, or between different simulations? Also, am I correct in understanding that the constraints are not based on the distances between individual atoms, but rather on the average of a bunch of different pairwise distances? What is the typical size of an "I" or "J" group, and how many of those groups would typically exist at once?
Hi Peter,
Sorry for not being more clear. The potential is to model ambiguous NOE restraints from NMR. Typically you would specify two groups of atoms and the interaction would be based on an r^-6 average of all pairs of atoms between the two groups (the I's and J's). A particular restraint will always involve the same number of atoms during the course of the simulation.
We will sometimes have several thousand of these restraints in a single system (typically 1 per NOE peak from NMR). The values for #I and #J are varied. One NOE might only be explainable by a few pairs of atoms (#I and #J are both small), while others might be very ambiguous and the signal could be due to as many as several hundred atom pairs (#I and #J might be as large as 20).
For any particular value of #I and #J, I can easily write a CustomCompoundBond to compute the interaction. Just to be concrete, say (#I = 2, #J = 1) or (#I = 2, #J = 2). The code would look like:
Code: Select all
ambig_rest_2_1 = CustomCompoundBondForce(3, '0.5*kbond*( (distance(p1, p3)^-6.0 + distance(p2, p3)^-6.0) / 2.0) ) ^ (-1.0/6.0)^2.0' )
ambig_rest_2_2 = CustomCompoundBondForce(4, '0.5*kbond*( (distance(p1, p3)^-6.0 + distance(p1, p4)^-6.0 + distance(p2,p3)^-6.0 + distance(p2,p4)^-6.0) / 4.0) ) ^ (-1.0/6.0)^2.0' )
I could write similar code any arbitrary value of #I and #J. Obviously, this would get very tedious. But, I can easily write a python class that can generate a CustomCompoundBondForce for any particular #I and #J (meta-meta programming?). What I am worried about is that I will potentially end up with many different kernels, one for each combination of #I and #J. Some of these kernels will have many interactions to compute (which is good) because there are many NOEs from experiment that have a particular degree of ambiguity #I, #J. Other kernels might have a small number of interactions to compute (I assume this is bad) because that degree of ambiguity is rare in the dataset. I'm worried that the combination of many different kernels with varying degrees of occupancy might cause performance issues.
Here's another way to put it. I have a general interaction that I can write as a double sum over two sets of atoms (I and J). What I am proposing is to write an unrolled version of this loop for each particular combination of the sizes of the two sets (#I and #J). I can implement this unrolled loop using a CustomCompoundBondForce. To someone with little CUDA experience, but tons of python experience, this seems very easy and straightforward to do. But, I am worried about the performance of this set of unrolled kernels versus writing (from scratch) a custom kernel that actually loops over all interactions.