Dear all,
I am looking for a way to restrain the center of mass (COM) of two different molecules bathed in water, using GROMACS 2018. The only thing close I found on gromacs that does this is the pull module embedded in the mdp files. I want to clarify that I want to restrain each molecules (COM) to a different and specific coordinate in my simulated box. Indeed the pull module has an option - pull-coord1-geometry = direction, where I can set the origin (pull-coord1-origin= 0 0 0, for example) and a vector directing the movement of the COM (pull-coord1-vec= 1 1 1) with time (pull-coord1-rate = 0.01nm/ps). However I don’t want my COM to move in time, I want it stationary in a specific coordinate throughout my entire simulation time. I see that I can set a parameter - pull-coord1-start = no - so that it doesn’t start moving towards my vector and stays in what I labeled “origin”.
With all that said, when running my simulations, it looks like my specification in the mdp file lead to no restraints, and it looks like my molecules defuse freely without an energetic bias towards my specified origin coordinate. Here are my mdp settings:
pull = yes pull-ngroups = 2 pull-group1-name = NAH pull-group2-name = RAT
pull-ncoords = 2
pull-coord1-type = umbrella pull-coord1-geometry = direction pull-coord1-groups = 0 1 pull-coord1-dim = Y Y Y pull-coord1-origin = 3.75 2.5 2.5 pull-coord1-vec = 1.0 1.0 1.0 pull-coord1-k = 10000.0 pull-coord1-rate = 0.0 pull-coord1-init = 0 pull-coord1-start = no
pull-coord2-type = umbrella pull-coord2-geometry = direction pull-coord2-groups = 0 2 pull-coord2-dim = Y Y Y pull-coord2-origin = 3.75 2.5 2.5 pull-coord2-vec = 1.0 1.0 1.0 pull-coord2-k = 10000.0 pull-coord2-rate = 0.0 pull-coord2-init = 0 pull-coord2-start = no
where pull-coord2-init = 0 is supposedly irrelevant to the “direction” type of pulling.
If anybody has any experience with this pulling direction technique and might point me to where my settings are wrong, Id be extremely thankful.
Best to all, Ilan Shumilin
This problem is easier to be solved by PLUMED plug-in than options in mdp file.
Hi Ilan,
Using plumed is one option, as Quinglu pointed out, albeit for simple tasks like this you'll be better off with plain Gromacs solutions as it gives you slightly better performance.
With regard to your input, I see several options to try out:
(1) if the desired distance from the origin is 0, why not use "pull-coord1-geometry = distance"
(2) right now you're pulling both molecules to the same point in space (3.75 2.5 2.5), looks more like a bug than a feature
(3) if you want to keep geometry=direction, why not try with pull-coord1-start=yes, combined with pull-coord1-init = 0 it will try to keep the molecule in its initial position
(4) in an actual run, are you getting .xvg files with positions and forces? this might give you more insight into what Gromacs is seeing internally, and potentially troubleshoot more easily.
Best,
Miłosz
Fixing an atom (or the COM of a group of atoms) to a specific position in space (x,y,z) requires three constraints, not one. Right now, I think, you constrain your motion along the vector (1,1,1), meaning that the COM can still move freely along its perpendicular plane.
The same is true for the second group.
Matthias is right here, and seems that switching to geometry=distance (which already works in 3D with "pull-coord1-dim = Y Y Y") will fix that instantly, without the need to add more restraints.
I'm not 100% on this, but I don't think so. With the current parameters, if you change to distance, then it fixes the distance to the reference (3.75, 2.5, 2.5) in the box, but it can still move around this point in a sphere. You can't choose the reference as the starting position, since the distance has to be bigger then 0.
I'm reasoning that with "pull-coord1-init=0" and "pull-coord1-start=no", the molecule's COM is harmonically restrained at distance 0 from the reference (3.75, 2.5, 2.5). I don't see why distance would have to be bigger than 0, but maybe I'm missing something.
This might be a little out of topic now, but I checked the gromacs command and it throws an error if you do that (Fatal error: Distance for pull coordinate 1 is zero with constraint pulling, which is not allowed.). The reasoning is most probably, that at some point it would devide by zero (at r=0, the Jacobian is not defined).
Ok, that seems to be a fair concern, at least for free energy applications down the line. So I guess the most physical approach would be a flat-bottom COM restraint with a small but non-zero distance.
- Badges
- Science topic