How to pull a molecule across a bilayer with fixed boundaries?

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

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

Best regards
Alessandra

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.

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 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):"

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

the above link contain the gro file

Hi,

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.

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 cylinder.

Best regards
Alessandra

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

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

\Alessandra

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?

Hi,
did you try different values for pull_coord1_rate and pull_coord1_k ?
\Alessandra

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.

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

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.

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