implement a cosine-based repulsive force

[golobor]

I found a nice cosine-based repulsive force which may be worth testing:

Not urgent, but I just wanted to document it somewhere.

[golobor]

btw, @mimakaev , I implemented this force - it's ~10% slower than polynomial_repulsive, but seems a bit more stable. (1) It led to quicker energy equilibration and (2) I could use it with lower friction /higher error tolerance. Do you think it's worth including into polychrom? If it would make the library a bit more messy, we could also have a separate file for "exotic_forces".

[mimakaev]

Yep, cosine is a difficult function to compute. But it is probably the best function as it has the least maximum gradient (at least up to ~4-5 kT). (1-x)^2 * (1+x)^2 has the same curve, and is easier to compute.

I agree about "Exotic" forces. Right now we have legacy forces that weren't fully ported. It we add exotic forces, we can slowly migrate legacy forces to exotic.

@anton: what was the truncation, and what was the polymer "wiggle distance" in your simulations?

Below there is a plot to look at. It suggests that if you're using polymer bonds with "wiggle distance" of 0.05, then it doesn't matter what nonbonded force you use: max derivative of polynomial force is larger in the "relevant" region of U less than 6-8kT. If you go to wiggle distance of 0.1, then polynomial force (labeled as "Anton") is better than Grosberg. But if you look at 6-8kT (which you would frequently reach), it still is on par with polynomial force with truncation=10 (which is as much as one should go, right?).

If you soften the polynomial force to 0.15, then you are probably not needing strong repulsion anyways, so 4kT polynomial repulsive force would be sufficient and should be used.

If you crank up truncation to 25, then old Grosberg force becomes on par with polynomial repulsive and with polymer bonds at wiggle=0.05. At truncation of 50 one should use grosberg force or cosine force, but I'm not sure why would anyone use such high repulsive energy.

[mimakaev]

Yep, cosine is a difficult function to compute. But it is probably the best function as it has the least maximum gradient (at least up to ~4-5 kT). (1-x)^2 * (1+x)^2 has the same curve, and is easier to compute.

I agree about "Exotic" forces. Right now we have legacy forces that weren't fully ported. It we add exotic forces, we can slowly migrate legacy forces to exotic.

@anton: what was the truncation, and what was the polymer "wiggle distance" in your simulations?

Below there is a plot to look at. It suggests that if you're using polymer bonds with "wiggle distance" of 0.05, then it doesn't matter what nonbonded force you use: max derivative of polynomial force is larger in the "relevant" region of U less than 6-8kT. If you go to wiggle distance of 0.1, then polynomial force (labeled as "Anton") is better than Grosberg. But if you look at 6-8kT (which you would frequently reach), it still is on par with polynomial force with truncation=10 (which is as much as one should go, right?).

If you soften the polynomial force to 0.15, then you are probably not needing strong repulsion anyways, so 4kT polynomial repulsive force would be sufficient and should be used.

If you crank up truncation to 25, then old Grosberg force becomes on par with polynomial repulsive and with polymer bonds at wiggle=0.05. At truncation of 50 one should use grosberg force or cosine force, but I'm not sure why would anyone use such high repulsive energy.

[golobor]

(a) (1+x)^2 * (1-x)^2 is a really pretty function, wow!
I just tested it in the example.py and it is very similar to cosine, both in terms of speed (10% fewer SPS), stability (similarly lower kinetic energy at the same time point comparing to the polynomial force); all three lead to ~same dt picked by the variableLangevin. In my understanding, cos must be a special function, whose calculation is optimized by GPUs?..

[golobor]

Your plots look super interesting, but I would need your help to digest it. Why are you plotting bond and repulsive potentials on the same plot? To compare their stiffness (i.e. for computational stability), or for some other reason?

[mimakaev]

wow, that could be! fewer SPS could be related to the size of the neighbors list! It is a very "soft" function, so you would have a lot of neighbors in the neighbors list.

Also, a fair comparison should be to compare cosine with r_cutoff = 1.5, polynomial with cutoff=1, and Lennard-Jones with the cutoff of 0.8-0.9 - that would put the "up" part of the potential in the same place, the forces would take roughly the same total volume, etc. I'm suspicious that LJ may actually win due to the steepness of potential that starts right away.

We should not forget about the density and the total volume occupied by our excluded volume interactions. Switch to polynomial force made density drop (and I often set r to 1.1 for it now to combat that). Cosine would leave even more space for the chains.

[mimakaev]

Our messages crossed. Replying to it now. What matters to the integrator is the largest realistically achievable gradient. I'm not sure of the true practical definition, but I would define it to be as a maximum gradient up to 6kT-8kT-10kT. ln(0.001) ~ 7, meaning one particle in ~1000 would have an energy of 7kT. Which is why I'm thinking it may be up to 10kT, but I"m really not sure.

So the maximum realistically achievable gradient could be realized at any force: bondforce, nonbonded and external. External is out of question - we rarely use more than 5kT/mon gradients. But the bondforce is evil - it is very steep! As you can see, gradient goes up to 100kt/mon and even larger. So an easy way to speed things up if you're doing a "relaxed" simulation with 1.5kT truncation energy and no fancy forces is to set wiggle_dist to 0.1 or even 0.15. There are my theoretical expectations... I can do some experiments to verify them.

[golobor]

(a) re: SPS,
(a1) I would expect that the size of the neighbour list would only depend on the cutOff distance and the density. In my test, I follow example.py, where density is 0.85 and keep the cutoff distance as 1.0 for all forces.
(a2) In a separate test, I also compared the polynomial force at trunc 4.0 with the other two forces at 10.0 and the results were the same. If anything, the difference in dt became even more noticeable.
(a3) also, the polynomial force is shitty for simulations with attraction-induced-condensation, because it's very flat at small r and particles can sometimes collapse onto each other.
So, yeah, in my simulations I will go with either of the two new ones.

(b) Yeah, agreed with the other point re: bond force. I've been also playing with duplicating bonds with "constraints", which are mathematical constraints on the bond lengths. The results are interesting - it sometimes breaks and it leads to massive drops in SPS, but the speed up in dt pretty much compensates it. I think these constraints could be useful for topological simulations.

[mimakaev]

re: SPS especially at density of 0.85, the size of the neighbors list would depend on how many monomers can sneak in "under" the cutoff radius. It is large for attractive forces, because many monomers reside in the potential well. If your force is repulsive and very steep, then the force will only be calculated for a thin sphere of monomers. For the cosine force alone, this shell would be pretty thick. But if you combine it with attraction, then the difference would probably vanish because the cutoff distance would be much larger than where the cosine (or (1-x^2) * (1+x^2)) ends