marrink-lab / polyply_1.0

Generate input parameters and coordinates for atomistic and coarse-grained simulations of polymers, ssDNA, and carbohydrates
Apache License 2.0
119 stars 21 forks source link

Implementation of MC bending sampling #344

Closed fgrunewald closed 8 months ago

fgrunewald commented 10 months ago

This code implements a local bending stiffness to the polymers. The idea is that one can set a bending constant for any triplet of residues. Subsequently, during the RW protocol the angle if a trail is accepted or rejected according to a MC-like procedure. First, the normalized probability of an angle is computed using:

$p(\theta) = \huge \frac{e^{\frac{lp \times \theta}{180}}}{\frac{180 \times e^{lp}-1}{lp}}$

This probability distribution is about equal for low values of lp and has an increased probability at higher angles with increasing lp. Note lp is a dimensionless bending constant not per se the persistence length. Once the probability is computed we check if the value is higher than the previous bending probability. If yes the move is accepted and if not we compare it to the probability of a uniform sampling of the distribution.

So far bending constant can be given as tuples of three residue names. This means local bending stiffness can be accounted for. We could therefore use it to really fine-tune the DNA parameters. However, it also means the bending term is non-symmetric which means for combinations of many different monomers this might become cumbersome. Ideally, we could set a default value for all monomers.

Points to discuss:

@jan-stevens let me know your thoughts.

I did some testing of this implementation and in general works quite well to reproduce the expected end-to-end distance by setting an appropriately large bending constant. Any relation to the experimental value of lp is, however, to me not obvious. Probably the estimated volume of the residue itself contributes to the stiffness to some extent so one potentially has to determine a qualitative bending constant by trail and error. For DNA, we can simply figure it out and report. The main limitation of this implementation is that the distributions are not as peaked for cases where contour-length >~ lp. In those cases the WCM distributions are probably not too good either.

fgrunewald commented 10 months ago

some observation:

stiff polymers where contour-length >~ lp need a low number of rewind steps (e.g. setting -rw 1). For single polymers this helps to get relatively straight chain conformations.

fgrunewald commented 9 months ago

summoning @pckroon for a technical code review