Limited range for reference group detection in cylinder pulling
The cylinder pulling is mostly used with membranes and it works great until the membrane is more or less planar. However, if one works with complex curved membranes the problems begin.
Imagine that one is simulating the liposome and want to perform cylinder pulling of some drug through the membrane on top of it along Z axis with "pull-coord1-vec = 0 0 1". Current implementation will search for atoms of the reference group (membrane) across the whole box and will include two patches of the membrane - on top and on bottom of the liposome. However, only top one is needed.
It is possible to workaround this issue by assigning different group names for top and bottom halves of the liposome, but this won't work in long simulations where complete mixing of lipids occurs.
Currently we are solving this problem by naive modification of the pullutil.cpp. Instead of the line
if (dr2_rel < 1)
if (dr2_rel < 1 && abs(axialLocation)< dist)
where dist is hardcoded distance, which is smaller than the diameter of liposome.
Is it possible to introduce additional parameter for cylindrical pulling i.e. "pull-cylinder-detection-max-dist" which will limit the range for detecting reference group atoms along the axis?
Have hard cutoffs is dangerous. That leads to delta forces when taking the gradient of the potential, which usually disappear in MD. We could use a continuous switching function.
But how should we determine the reference z coordinate for determining this distance?
#2 Updated by Semen Yesylevskyy 9 months ago
This could be done with a third group serving as a separate reference for distance computations. For example we have the groups "vesicle" and "ligand" and want to use only top part of vesicle. The user knows that the redius of vesicle is ~10 nm and defines something like:
pull-coord1-vec = 0 0 1 ; detect distance along this vetor
; dist.ref.group cyl.group pull.group
pull-coord1-grps = vesicle vesicle ligand
; only take atoms from cyl.group in limits (min,max) from COM of dist.ref.group along vec
pull-coord1-cylinder-detection-dist-min = 5
pull-coord1-cylinder-detection-dist-max = 15
It works like this. The com of the whole vesicle (which is dist.ref.group) is computed. The axis of cylinder is defined along vec. Only atoms of vesicle (which is also a cyl.group) within (min,max) along vec are included using switch function. After that the distance between com of selected atoms and com of ligand is computed and the forces applied as usually.
dist.ref.group could be distinct from cyl.group providing better flexibility.