# How to pull a molecule across a bilayer with fixed boundaries? **URL:** \ **Category:** User discussions\ **Tags:** pull-code\ **Created:** [April 21, 2021, 10:26am UTC](https://gromacs.bioexcel.eu/t/how-to-pull-a-molecule-across-a-bilayer-with-fixed-boundaries/1999 "2021-04-21T10:26:58Z")\ **Posts on this page:** 15\ **Page:** 1
**Author:** ![drrahulsuresh](https://avatars.discourse-cdn.com/v4/letter/d/e68b1a/32.png) [@drrahulsuresh](https://gromacs.bioexcel.eu/u/drrahulsuresh)\ **Post date:** [April 21, 2021, 10:26am UTC](https://gromacs.bioexcel.eu/t/how-to-pull-a-molecule-across-a-bilayer-with-fixed-boundaries/1999/1 "2021-04-21T10:26:58Z")
GROMACS version: 2016.4 GROMACS modification: No Dear Users I am trying to pull a small molecule across a lipid bilayer. The bilayer is constructed using charmm-gui and the box is fixed the size of the lipid bilayer. During pull process, I end up with error “Distance between pull groups 1 and 2 (5.763902 nm) is larger than 0.49 times the box size (4.446700).” So in case I have to increase the box size, can I do it before minimization and continue with the steps? Or is there anything which I can add to the script to pull the molecule with the same dimensions as now? Pull code ; Pull code pull = yes pull\_ncoords = 1 ; only one reaction coordinate pull\_ngroups = 2 ; two groups defining one reaction coordinate pull\_group1\_name = POPC pull\_group2\_name = POL pull\_coord1\_type = umbrella ; harmonic potential pull\_coord1\_geometry = distance ; simple distance increase pull\_coord1\_dim = N N Y pull\_coord1\_groups = 1 2 pull\_coord1\_start = yes ; define initial COM distance \> 0 pull\_coord1\_rate = 0.005 ; 0.01 nm per ps = 10 nm per ns pull\_coord1\_k = 1000 ; kJ mol^-1 nm^-2 The Bilayer also consists of sterols. Any suggestions would be really helpful ---
**Author:** ![alevilla](https://dub1.discourse-cdn.com/flex017/user_avatar/gromacs.bioexcel.eu/alevilla/32/439_2.png) [@alevilla](https://gromacs.bioexcel.eu/u/alevilla)\ **Post date:** [April 22, 2021, 10:24am UTC](https://gromacs.bioexcel.eu/t/how-to-pull-a-molecule-across-a-bilayer-with-fixed-boundaries/1999/2 "2021-04-22T10:24:46Z")
Hi, you can use the mdp option `pull-pbc-ref-prev-step-com` when there are large pull groups (like bilayer) see more here (under Definition of the center of mass) [https://manual.gromacs.org/2021/reference-manual/special/pulling.html?highlight=pull%20code%20umbrella](https://manual.gromacs.org/2021/reference-manual/special/pulling.html?highlight=pull%20code%20umbrella) Best regards Alessandra ---
**Author:** ![drrahulsuresh](https://avatars.discourse-cdn.com/v4/letter/d/e68b1a/32.png) [@drrahulsuresh](https://gromacs.bioexcel.eu/u/drrahulsuresh)\ **Post date:** [April 26, 2021, 9:50am UTC](https://gromacs.bioexcel.eu/t/how-to-pull-a-molecule-across-a-bilayer-with-fixed-boundaries/1999/3 "2021-04-26T09:50:24Z")
Hi @alevilla The command doesn’t help. The molecule is on the surface of lipid leaflet and the aim of the simulation is to calculate PMF by pulling the molecule from one side to another. Any suggestions will be really helpful. ---
**Author:** ![alevilla](https://dub1.discourse-cdn.com/flex017/user_avatar/gromacs.bioexcel.eu/alevilla/32/439_2.png) [@alevilla](https://gromacs.bioexcel.eu/u/alevilla)\ **Post date:** [April 27, 2021, 8:13am UTC](https://gromacs.bioexcel.eu/t/how-to-pull-a-molecule-across-a-bilayer-with-fixed-boundaries/1999/4 "2021-04-27T08:13:30Z")
Hi, Do you mean that you have tried all the options suggested in the link above? That are pull-pbc-ref-prev-step-com, `pull-group?-pbcatom, and pull-coord?-geometry = cylinder \Alessandra Extract from the link above "When there are large pull groups, such as a lipid bilayer, `pull-pbc-ref-prev-step-com` can be used to avoid potential large movements of the center of mass in case that atoms in the pull group move so much that the reference atom is too far from the intended center of mass. With this option enabled the center of mass from the previous step is used, instead of the position of the reference atom, to determine the reference position. The position of the reference atom is still used for the first step. For large pull groups it is important to select a reference atom that is close to the intended center of mass, i.e. do not use `pull-group?-pbcatom = 0` . For a layered system, for instance a lipid bilayer, it may be of interest to calculate the PMF of a lipid as function of its distance from the whole bilayer. The whole bilayer can be taken as reference group in that case, but it might also be of interest to define the reaction coordinate for the PMF more locally. The [mdp](https://manual.gromacs.org/2021/reference-manual/file-formats.html#mdp) option `pull-coord?-geometry = cylinder` does not use all the atoms of the reference group, but instead dynamically only those within a cylinder with radius `pull-cylinder-r` around the pull vector going through the pull group. This only works for distances defined in one dimension, and the cylinder is oriented with its long axis along this one dimension. To avoid jumps in the pull force, contributions of atoms are weighted as a function of distance (in addition to the mass weighting):" ---
**Author:** ![drrahulsuresh](https://avatars.discourse-cdn.com/v4/letter/d/e68b1a/32.png) [@drrahulsuresh](https://gromacs.bioexcel.eu/u/drrahulsuresh)\ **Post date:** [April 28, 2021, 3:29pm UTC](https://gromacs.bioexcel.eu/t/how-to-pull-a-molecule-across-a-bilayer-with-fixed-boundaries/1999/5 "2021-04-28T15:29:31Z")
Yes I have tried all the options. ; Pull code pull = yes pull\_print\_com = no pull\_print\_ref\_value = yes pull\_print\_components = Yes pull\_nstxout = 500 pull\_nstfout = 500 pull\_ncoords = 1 ; only one reaction coordinate pull\_ngroups = 2 ; two groups defining one reaction coordinate pull\_group1\_name = POL pull\_group2\_name = POPC pull-group1-pbcatom = 1 pull\_coord1\_type = umbrella ; harmonic potential pull\_coord1\_geometry = cylinder pull\_coord1\_vec = 0.0 0.0 1.0 pull\_coord1\_dim = N N Y pull\_coord1\_groups = 1 2 pull\_coord1\_start = yes ; define initial COM distance \> 0 pull\_coord1\_rate = 0.005 ; 0.01 nm per ps = 10 nm per ns pull\_coord1\_k = 1000 ; kJ mol^-1 nm^-2 Fatal error: The total potential energy is -nan, which is not finite. The LJ and electrostatic contributions to the energy are 1768.28 and -742206, respectively. A non-finite potential energy can be caused by overlapping interactions in bonded interactions or very large or NaN coordinate values. Usually this is caused by a badly or non-equilibrated initial configuration or incorrect interactions or parameters in the topology. For the option pull-pbc-ref-prev-step-com Warning Unknown left-hand ‘pull-pbc-ref-prev-step-com’ in parameter file > **[step6.6\_equilibration.gro](https://drive.google.com/file/d/1qg_M1hf0y2jxlZVNGBiwSM-r1jBbiwSx/view?usp=sharing)** > > Google Drive file. the above link contain the gro file ---
**Author:** ![alevilla](https://dub1.discourse-cdn.com/flex017/user_avatar/gromacs.bioexcel.eu/alevilla/32/439_2.png) [@alevilla](https://gromacs.bioexcel.eu/u/alevilla)\ **Post date:** [April 29, 2021, 1:28pm UTC](https://gromacs.bioexcel.eu/t/how-to-pull-a-molecule-across-a-bilayer-with-fixed-boundaries/1999/6 "2021-04-29T13:28:04Z")
Hi, > [@drrahulsuresh](#): > > Warning > Unknown left-hand ‘pull-pbc-ref-prev-step-com’ in parameter file Correct, you use a old version, where this option is not implemented. If you have the possibility, it is always better to use the most recent version of the program. > [@drrahulsuresh](#): > > Fatal error: > The total potential energy is -nan, which is not finite. The LJ and > electrostatic contributions to the energy are 1768.28 and -742206, > respectively. A non-finite potential energy can be caused by overlapping > interactions in bonded interactions or very large or NaN coordinate values. > Usually this is caused by a badly or non-equilibrated initial configuration or > incorrect interactions or parameters in the topology. when in time occurs such a error? It seems that you are pulling too fast. Since you pull a molecule from one side to the other side of the membrane, the lipids need the time to rearrange and make space for the molecule to pass through. Note that pull\_group1-pbcatom = 1 is not used with [`pull-coord1-geometry`](https://manual.gromacs.org/current/user-guide/mdp-options.html#mdp-pull-coord1-geometry) cylinder. Best regards Alessandra ---
**Author:** ![drrahulsuresh](https://avatars.discourse-cdn.com/v4/letter/d/e68b1a/32.png) [@drrahulsuresh](https://gromacs.bioexcel.eu/u/drrahulsuresh)\ **Post date:** [April 29, 2021, 4:01pm UTC](https://gromacs.bioexcel.eu/t/how-to-pull-a-molecule-across-a-bilayer-with-fixed-boundaries/1999/7 "2021-04-29T16:01:15Z")
> [@drrahulsuresh](#): > > ; Pull code > pull = yes > pull\_print\_com = no > pull\_print\_ref\_value = yes > pull\_print\_components = Yes > pull\_nstxout = 500 > pull\_nstfout = 500 > pull\_ncoords = 1 ; only one reaction coordinate > pull\_ngroups = 2 ; two groups defining one reaction coordinate > pull\_group1\_name = POL > pull\_group2\_name = POPC > pull\_group1-pbcatom = 1 > pull\_coord1\_type = umbrella ; harmonic potential > pull\_coord1\_geometry = cylinder > pull\_coord1\_vec = 0.0 0.0 1.0 > pull\_coord1\_dim = N N Y > pull\_coord1\_groups = 1 2 > pull\_coord1\_start = yes ; define initial COM distance \> 0 > pull\_coord1\_rate = 0.005 ; 0.01 nm per ps = 10 nm per ns > pull\_coord1\_k = 1000 ; kJ mol^-1 nm^-2 [https://drive.google.com/file/d/1n8y7ZG-s6nRZqykoPHPZE7EXdiG9K5Io/view?usp=sharing](https://drive.google.com/file/d/1n8y7ZG-s6nRZqykoPHPZE7EXdiG9K5Io/view?usp=sharing) Is the link to the mdp file. I have 2016.4 version gromacs. Yes I will try to simulate with 2018 or higher version but I am curious to understand the error here. Thank you ---
**Author:** ![alevilla](https://dub1.discourse-cdn.com/flex017/user_avatar/gromacs.bioexcel.eu/alevilla/32/439_2.png) [@alevilla](https://gromacs.bioexcel.eu/u/alevilla)\ **Post date:** [April 30, 2021, 7:23am UTC](https://gromacs.bioexcel.eu/t/how-to-pull-a-molecule-across-a-bilayer-with-fixed-boundaries/1999/8 "2021-04-30T07:23:56Z")
Hi, In my previous answer, I wrote what I thought was the problem. Maybe you have missed the second part of the post. Here I repost > [@alevilla](#): > > when in time occurs such a error? It seems that you are pulling too fast. Since you pull a molecule from one side to the other side of the membrane, the lipids need the time to rearrange and make space for the molecule to pass through. > > Note that pull\_group1-pbcatom = 1 is not used with [`pull-coord1-geometry`](https://manual.gromacs.org/current/user-guide/mdp-options.html#mdp-pull-coord1-geometry) cylinder. \Alessandra ---
**Author:** ![drrahulsuresh](https://avatars.discourse-cdn.com/v4/letter/d/e68b1a/32.png) [@drrahulsuresh](https://gromacs.bioexcel.eu/u/drrahulsuresh)\ **Post date:** [June 1, 2021, 6:35am UTC](https://gromacs.bioexcel.eu/t/how-to-pull-a-molecule-across-a-bilayer-with-fixed-boundaries/1999/9 "2021-06-01T06:35:50Z")
Sorry for the late reply @alevilla Lockdown is quite severe here in town Regarding the query you asked, I am not sure when it came because I couldn’t reproduce it. however, I overcame the issue. ; Pull code pull = yes pull\_ncoords = 1 ; only one reaction coordinate pull\_ngroups = 2 ; two groups defining one reaction coordinate pull\_group1\_name = lipids pull\_group2\_name = POL pull\_coord1\_type = umbrella ; harmonic potential pull\_coord1\_dim = N N Y ; pull along z pull\_coord1\_groups = 1 2 ; pull\_coord1\_start = yes ; define initial COM distance \> 0 pull-pbc-ref-prev-step-com = yes pull-coord1-geometry = cylinder pull\_coord1\_vec = 0.0 0.0 1.0 pull\_coord1\_rate = 0.0001 ; pull\_coord1\_k = 1000 ; kJ mol^-1 nm^-2 The above is the pull code I used for simulation. The ligand has to be pulled across the bilayer, i.e. from one end to another but the molecule is stopped at the center of the bilayer. How do I pull further? ---
**Author:** ![alevilla](https://dub1.discourse-cdn.com/flex017/user_avatar/gromacs.bioexcel.eu/alevilla/32/439_2.png) [@alevilla](https://gromacs.bioexcel.eu/u/alevilla)\ **Post date:** [June 15, 2021, 8:43am UTC](https://gromacs.bioexcel.eu/t/how-to-pull-a-molecule-across-a-bilayer-with-fixed-boundaries/1999/10 "2021-06-15T08:43:49Z")
Hi, did you try different values for pull\_coord1\_rate and pull\_coord1\_k ? \Alessandra ---
**Author:** ![drrahulsuresh](https://avatars.discourse-cdn.com/v4/letter/d/e68b1a/32.png) [@drrahulsuresh](https://gromacs.bioexcel.eu/u/drrahulsuresh)\ **Post date:** [June 15, 2021, 5:53pm UTC](https://gromacs.bioexcel.eu/t/how-to-pull-a-molecule-across-a-bilayer-with-fixed-boundaries/1999/11 "2021-06-15T17:53:30Z")
Hi @alevilla Yes I have tried different values and 1000 was my highest and 0.0001 was the least. The motion of the molecule was clear with this rate. ---
**Author:** ![alevilla](https://dub1.discourse-cdn.com/flex017/user_avatar/gromacs.bioexcel.eu/alevilla/32/439_2.png) [@alevilla](https://gromacs.bioexcel.eu/u/alevilla)\ **Post date:** [June 16, 2021, 9:25am UTC](https://gromacs.bioexcel.eu/t/how-to-pull-a-molecule-across-a-bilayer-with-fixed-boundaries/1999/12 "2021-06-16T09:25:08Z")
Hi, I guess you have simulate long enough to get the molecule from one side to the other, knowing that the rate is in nm/ps. Alessandra ---
**Author:** ![drrahulsuresh](https://avatars.discourse-cdn.com/v4/letter/d/e68b1a/32.png) [@drrahulsuresh](https://gromacs.bioexcel.eu/u/drrahulsuresh)\ **Post date:** [June 16, 2021, 4:16pm UTC](https://gromacs.bioexcel.eu/t/how-to-pull-a-molecule-across-a-bilayer-with-fixed-boundaries/1999/13 "2021-06-16T16:16:09Z")
Hi Yes, I have simulated it for 10ns. the molecule reach the center of the bilayer before 3ns and after that the molecule vibrate in the center of the bilayer. ---
**Author:** ![alevilla](https://dub1.discourse-cdn.com/flex017/user_avatar/gromacs.bioexcel.eu/alevilla/32/439_2.png) [@alevilla](https://gromacs.bioexcel.eu/u/alevilla)\ **Post date:** [June 17, 2021, 9:27am UTC](https://gromacs.bioexcel.eu/t/how-to-pull-a-molecule-across-a-bilayer-with-fixed-boundaries/1999/14 "2021-06-17T09:27:32Z")
Hi in 10 ns (10^4 ps) with a rate of 10^(-4) nm/ps the molecule will `move` of 1 nm from the begin. Or? Is this what you observe? \Alessandra ---
**Author:** ![drrahulsuresh](https://avatars.discourse-cdn.com/v4/letter/d/e68b1a/32.png) [@drrahulsuresh](https://gromacs.bioexcel.eu/u/drrahulsuresh)\ **Post date:** [June 25, 2021, 7:48am UTC](https://gromacs.bioexcel.eu/t/how-to-pull-a-molecule-across-a-bilayer-with-fixed-boundaries/1999/15 "2021-06-25T07:48:43Z")
Yes. The molecule does not cross the center of the Bilayer and remain in the middle of the bilayer for nearly 6ns of the simulation. But the module do move from its position. Only thing is I couldn’t pull it further.