SMD simulation - protein pulling and resulting rotation

User discussions
, ,
1

GROMACS version: 2025.2
GROMACS modification: Yes

Hello everyone,

I am having some questions about my pull code.

The idea is to hold a protein via two springs on the two terminals and to pull one terminal along the x axis (the reference point should be pulled actively and the connected terminal will be dragged along).

I tried to simulate this by creating two individual reference points for each terminal and to simulate a spring via pull-coord1-type = umbrella. In this example, I want to pull in the -x axis, which is why I am setting pull-coord1-vec = -1.0 0.0 0.0.

Do I understand correctly, that pull-coord1-groups = 0 1 creates the setup the way, that the pulling force is directly applied to the chosen reference?

When looking at my protein (it is placed into an ionized and equilibrated water box) it is being realistically pulled apart, but my protein is changing its position a lot in y and z axes (especially y axis). About +4 nm in y axis and +2-3 nm in z axis. It seems to “fly away” and rotate a lot.

I do understand, that I do not have a fixed reference point in the protein (the pull-groupx-name referred as “fixed” or “pulled” in the protein are atoms, defined in index.ndx, that should be connected to the springs). However, why would it be “flying away” in y axis and what could be done about it ? ( I know about the freeze option and also about creating a posre_fixed.itp file, but I would not like to put additional forces onto one specific atom, as the springs should be enough (theoretically)).

Thanks a lot for your help and any suggestions for improvement!

Best,

AL

;---------------code--------------
pull = yes
pull-ngroups = 2
pull-ncoords = 2

pull-group1-name = pulled
pull-group2-name = fixed

pull-coord1-type = umbrella
;pull-pbc-ref-prev-step-com = yes

pull-coord1-geometry = direction
pull-coord1-origin = 0.8 4 4
pull-coord1-vec = -1.0 0.0 0.0
pull-coord1-groups = 0 1
pull-coord1-start = yes
pull-coord1-dim = Y N N
pull-coord1-init = 0
pull-coord1-rate = 0.0001
pull-coord1-k = 660

pull-coord2-type = umbrella
pull-coord2-geometry = direction
pull-coord2-vec = 1.0 0.0 0.0
pull-coord2-groups = 0 2
pull-coord2-start = yes
pull-coord2-dim = Y Y Y
pull-coord2-origin = 17 4 4
;pull-coord2-init = 5.0
pull-coord2-rate = 0.0
pull-coord2-k = 660

2

Since both springs are only applying a force along the x-axis, it seems natural that motion would still be occurring along the other 2 axis. Considering how the initial positions of the pulling groups have different initial coordinates along the y and z axis, it also seems natural that pulling them only along the x-axis could drive additional rotation.

Is the goal here just to pull the termini apart from each other along the x-direction while keeping one terminus fixed? It seems to me that it might be easier to set just one pull group including the 2 terminal ends and specify the geometry as distance while keeping one end fixed with constraints along all axis if you want to keep things fixed in place. Also, setting the groups to include the entire terminal residue rather than individual atoms might give more physically realistic results.

3

By far the simplest is using two pull groups, one for each terminus, and pull those apart. But you will run into issues when the distance becomes larger than half the box size. One could use an elongated box and pull along one dimension. Then one can use geometry=distance-periodic to allow for distances larger than half the box size. One would then have to orient the protein to avoid large initial rotations.