# Steered MD Simulation error **URL:** \ **Category:** User discussions\ **Tags:** mdp-parameters, mdrun, simulation-setup\ **Created:** [January 25, 2026, 4:12am UTC](https://gromacs.bioexcel.eu/t/steered-md-simulation-error/13144 "2026-01-25T04:12:30Z")\ **Posts on this page:** 20\ **Page:** 1
**Author:** ![Aneesha.A](https://avatars.discourse-cdn.com/v4/letter/a/bc79bd/32.png) [@Aneesha.A](https://gromacs.bioexcel.eu/u/Aneesha.A)\ **Post date:** [January 25, 2026, 4:12am UTC](https://gromacs.bioexcel.eu/t/steered-md-simulation-error/13144/1 "2026-01-25T04:12:30Z")
GROMACS version: 2021.4/gpuvolta GROMACS modification: No Hi everyone, I am following Tutorial 3: Umbrella Sampling, but altering the .mdp files to a method section from a paper (unfortunately the authors are not much help - I’m building my own .mdp files from scratch). After running a steered MD simulation, the pullf.xvg file looks incorrect (below). My system consists of two alpha-helix structures each 35 residues in length and the goal of my simulation is to perform Umbrella Sampling to determine free energy. I appreciate any help, thank you! ![Screenshot 2026-01-25 at 2.54.33 pm](https://europe1.discourse-cdn.com/flex017/uploads/bioexcel1/original/2X/6/6895cffbe1b840dfb70db10a33da02680770e6bb.png) Also, during my Steered MD simulation I get these messages: ![Screenshot 2026-01-25 at 2.55.25 pm](https://europe1.discourse-cdn.com/flex017/uploads/bioexcel1/optimized/2X/c/cc47e405b4951b53de902724f5bb223c9396d232_2_690x23.png) Here is my md\_pull.mdp: title = Umbrella pulling simulation define = -DPOSRES ; Run parameters integrator = md dt = 0.002 tinit = 0 nsteps = 250000 ; 500 ps nstcomm = 10 ; Output parameters nstxout = 5000 ; Every 10 ps nstvout = 5000 nstfout = 500 nstxtcout = 500 nstenergy = 500 constraint\_algorithm = lincs lincs\_iter = 2 lincs\_order = 4 constraints = h-bonds continuation = yes pbc = xyz DispCorr = EnerPres ; PME electrostatics parameters coulombtype = PME fourierspacing = 0.12 fourier\_nx = 0 fourier\_ny = 0 fourier\_nz = 0 pme\_order = 4 ewald\_rtol = 1e-5 optimize\_fft = yes ; Pull code pull = yes pull\_ncoords = 1 ; Only one reaction coordinate pull\_ngroups = 2 ; Two groups defining one reaction coordinate pull\_group1\_name = Chain\_A pull\_group2\_name = Chain\_B pull\_group1\_pbcatom = 283 pull\_group2\_pbcatom = 867 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.01 ; 0.1 nm per ps = 50 nm per ns pull\_coord1\_k = 1000 ; kJ mol^-1 nm^-2 pull\_pbc\_ref\_prev\_step\_com = yes When I change pull\_coord1\_dim = Y Y Y I get: ![Screenshot 2026-01-25 at 3.11.01 pm](https://europe1.discourse-cdn.com/flex017/uploads/bioexcel1/original/2X/9/9ed802854fccc2f1388091339cd8b2cd37af34ed.png) ---
**Author:** ![hess](https://dub1.discourse-cdn.com/flex017/user_avatar/gromacs.bioexcel.eu/hess/32/416_2.png) [@hess](https://gromacs.bioexcel.eu/u/hess)\ **Post date:** [January 26, 2026, 12:00pm UTC](https://gromacs.bioexcel.eu/t/steered-md-simulation-error/13144/2 "2026-01-26T12:00:58Z")
You are pulling at an insanely high rate. I don’t expect the system could follow even with a strong spring constant. But with your rather weak spring you are simply extending the spring, leading to the results you show. ---
**Author:** ![Aneesha.A](https://avatars.discourse-cdn.com/v4/letter/a/bc79bd/32.png) [@Aneesha.A](https://gromacs.bioexcel.eu/u/Aneesha.A)\ **Post date:** [January 27, 2026, 6:06am UTC](https://gromacs.bioexcel.eu/t/steered-md-simulation-error/13144/3 "2026-01-27T06:06:44Z")
Hi Hess! Thanks for the reply. Actually, the paper I’m following says “A 500 ps COM pull sampling was performed with a consistent force approach (k = 1000 kJ/mol/nm) and a rate of 0.1 nm/ps between the COMs”. Unfortunately, the paper did not include their pullf.xvg graph, so every alteration I test for, goes until creating the pmf.xvg graph. When I change pull\_coord1\_rate to 0.0001, here are the pullf.xvg and pmf.xvg graphs I get, respectively. ![Screenshot 2026-01-27 at 1.53.38 pm](https://europe1.discourse-cdn.com/flex017/uploads/bioexcel1/original/2X/d/d565a75268f996bd8db82a35d335d081e77494c0.png) ![Screenshot 2026-01-27 at 1.54.25 pm](https://europe1.discourse-cdn.com/flex017/uploads/bioexcel1/original/2X/b/b3d7af75c6018a6f1f1a095d16351a2edc83b41f.jpeg) Is trial and error the only way to find a ‘good’ pull\_coord1\_rate? Also, how come in every simulation I do, the x-axis of my pmf.xvg graph spans 0.72 - 0.77 nm? The paper I am referring to has the following graph: ![Screenshot 2026-01-27 at 5.06.10 pm](https://europe1.discourse-cdn.com/flex017/uploads/bioexcel1/original/2X/b/b38d45559518bbbb08270ea57d8501ea07fd6b07.png) ---
**Author:** ![hess](https://dub1.discourse-cdn.com/flex017/user_avatar/gromacs.bioexcel.eu/hess/32/416_2.png) [@hess](https://gromacs.bioexcel.eu/u/hess)\ **Post date:** [January 27, 2026, 7:44am UTC](https://gromacs.bioexcel.eu/t/steered-md-simulation-error/13144/4 "2026-01-27T07:44:11Z")
You need to have some intuition/experience to see what would work. But I’d say that it’s rather obvious that two helices would never move apart at a speed of 0.1 nm/ps in reality. That’s far too fast. Maybe it’s a typo in the paper and it should be 0.1 nm/ns, which sounds reasonable. Then the rather force constant of 1000 kJ/mol/nm^2 (you write kJ/mol/nm) could work, although I would choose a higher value. Obviously you need to run much longer with a pull rate of 0.1 nm/ns. ---
**Author:** ![hess](https://dub1.discourse-cdn.com/flex017/user_avatar/gromacs.bioexcel.eu/hess/32/416_2.png) [@hess](https://gromacs.bioexcel.eu/u/hess)\ **Post date:** [January 27, 2026, 7:46am UTC](https://gromacs.bioexcel.eu/t/steered-md-simulation-error/13144/5 "2026-01-27T07:46:06Z")
A general note is that you will often not get a reasonable free-energy estimate from a single steered MD simulation, as there will be a lot of irreversible work included in the work. You would need run many simulations and used Jarzynski’s equality to get the free energy (but this only works with low variance). ---
**Author:** ![Sasha](https://avatars.discourse-cdn.com/v4/letter/s/54ee81/32.png) [@Sasha](https://gromacs.bioexcel.eu/u/Sasha)\ **Post date:** [January 29, 2026, 6:50am UTC](https://gromacs.bioexcel.eu/t/steered-md-simulation-error/13144/6 "2026-01-29T06:50:46Z")
> [@Aneesha.A](#): > > pull\_coord1\_rate = 0.01 ; 0.1 nm per ps = 50 nm per ns This is 0.01 nm/ps and also “0.1 nm per ps = 50 nm per ns“ is a false statement. Regardless, the above comments are correct. Also, it may be useful to convert these units to get a common sense feel for these numbers: 0.1 nm/ps = 100 m/s = 360 km/h, which is what a Bugatti does on a sunny day. That aside, it is possible that the paper you are trying to follow is nonsense, because free energy estimates don’t really come from direct pulls. Considerably slower (1e-5 to 1e-4 nm/ps) pulls can be useful if you wish to monitor the force and/or to get a Mickey Mouse estimate of the force integral over time or some carefully selected distance between the groups, but none of those things yield a reliable estimate of the free energy. To get the latter, I’d recommend following the Justin Lemkul’s tutorial: [Umbrella Sampling](http://www.mdtutorials.com/gmx/umbrella/index.html) ---
**Author:** ![Aneesha.A](https://avatars.discourse-cdn.com/v4/letter/a/bc79bd/32.png) [@Aneesha.A](https://gromacs.bioexcel.eu/u/Aneesha.A)\ **Post date:** [January 29, 2026, 10:33pm UTC](https://gromacs.bioexcel.eu/t/steered-md-simulation-error/13144/7 "2026-01-29T22:33:43Z")
Hi Sasha and Hess! Thank you both for your reply. Yes I am new to GROMACS and Umbrella Sampling, so I appreciate all the help I get. I ran another steered MD simulation (graph below), this time with a pull\_coord1\_rate of 0.0001. (0.1 nm/ns) and pull\_coord1\_k of 2000 kJ/mol/nm^2. I also ran the simulation for much longer (below parameters): ; Run parameters integrator = md dt = 0.002 tinit = 0 nsteps = 25000000 ; 50 ns nstcomm = 10 ; Output parameters nstxout = 500000 ; Every 10 ps nstvout = 500000 nstfout = 50000 nstxtcout = 50000 nstenergy = 50000 I’ve been following Justin Lemkul’s tutorial on Umbrella Sampling, but modifying the .mdp parameters to align with the paper’s method. Is there a better method for free energy estimates? ![Screenshot 2026-01-30 at 9.26.32 am](https://europe1.discourse-cdn.com/flex017/uploads/bioexcel1/original/2X/1/103800e56037aa1366436d8732205dadde2152b3.png) ---
**Author:** ![Sasha](https://avatars.discourse-cdn.com/v4/letter/s/54ee81/32.png) [@Sasha](https://gromacs.bioexcel.eu/u/Sasha)\ **Post date:** [January 29, 2026, 11:01pm UTC](https://gromacs.bioexcel.eu/t/steered-md-simulation-error/13144/8 "2026-01-29T23:01:42Z")
I think the group you’re pulling is just stuck in place, that’s all. The pull-code’s reference point moves 5 nm in 50 ns and now you have Hooke’s law (F = kx), where k = 2000 kJ/mol/nm^2 and x = 5 nm. At the end of the simulation, you end up with F = 10000kJ/mol/nm, which means whatever you’re pulling is not moving. I think you should look at your trajectory first. There are several ways to estimate free energies and the tutorial above is a very good way – even if it feels confusing at the moment. Direct pulling does not estimate free energies. ---
**Author:** ![Aneesha.A](https://avatars.discourse-cdn.com/v4/letter/a/bc79bd/32.png) [@Aneesha.A](https://gromacs.bioexcel.eu/u/Aneesha.A)\ **Post date:** [February 6, 2026, 11:18pm UTC](https://gromacs.bioexcel.eu/t/steered-md-simulation-error/13144/9 "2026-02-06T23:18:46Z")
Hi Sasha, Sorry for the late reply, I was unable to run simulations for some time. Your explanation about a group being pulled but stuck in place makes sense. I extracted the trajectory frames with the python script provided in the tutorial and here is what the “summary\_distances.dat” file looks like. I’m not sure what to make of it or how to fix it: ![Screenshot 2026-02-07 at 10.15.46 am](https://europe1.discourse-cdn.com/flex017/uploads/bioexcel1/original/2X/e/e6475b7a31a76728e5b37c645a01749016f3b256.png) Here is a full copy of my “md\_pull.mdp”: title = Umbrella pulling simulation define = -DPOSRES ; Run parameters integrator = md dt = 0.002 tinit = 0 nsteps = 25000000 ; 50 ns nstcomm = 10 ; Output parameters nstxout = 500000 ; Every 10 ps nstvout = 500000 nstfout = 50000 nstxtcout = 50000 nstenergy = 50000 constraint\_algorithm = lincs lincs\_iter = 2 lincs\_order = 4 constraints = h-bonds continuation = yes pbc = xyz DispCorr = EnerPres ; PME electrostatics parameters coulombtype = PME fourierspacing = 0.12 fourier\_nx = 0 fourier\_ny = 0 fourier\_nz = 0 pme\_order = 4 ewald\_rtol = 1e-5 optimize\_fft = yes ; Pull code pull = yes pull\_ncoords = 1 ; Only one reaction coordinate pull\_ngroups = 2 ; Two groups defining one reaction coordinate pull\_group1\_name = Chain\_A pull\_group2\_name = Chain\_B pull\_group1\_pbcatom = 283 pull\_group2\_pbcatom = 867 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.0001 pull\_coord1\_k = 2000 ; kJ mol^-1 nm^-2 pull\_pbc\_ref\_prev\_step\_com = yes Again, I appreciate all the help! ---
**Author:** ![hess](https://dub1.discourse-cdn.com/flex017/user_avatar/gromacs.bioexcel.eu/hess/32/416_2.png) [@hess](https://gromacs.bioexcel.eu/u/hess)\ **Post date:** [February 10, 2026, 8:08am UTC](https://gromacs.bioexcel.eu/t/steered-md-simulation-error/13144/10 "2026-02-10T08:08:27Z")
You likely need a much larger force constant. But using distance along one dimension is strange. Then I would use geometry=direction. And why are use using pull\_pbc\_ref\_prev\_step\_com=yes? ---
**Author:** ![Aneesha.A](https://avatars.discourse-cdn.com/v4/letter/a/bc79bd/32.png) [@Aneesha.A](https://gromacs.bioexcel.eu/u/Aneesha.A)\ **Post date:** [February 11, 2026, 5:26am UTC](https://gromacs.bioexcel.eu/t/steered-md-simulation-error/13144/11 "2026-02-11T05:26:50Z")
Hi Hess, appreciate the reply! The first change to my “md\_pull.mdp” file was my force constant to 8000 kJ mol^-1 nm^-2. Second, I realised I had a positional restraint on the whole system, “define = -DPOSRES”. I changed this to “define = -DPOSRES\_A” and accordingly changed “topol\_Protein\_chain\_A.itp” to include: #ifdef POSRES\_A #include “posre\_Protein\_chain\_A.itp” #endif I made this change with the logic of restraining one chain and leaving the second chain to freely move in the defined box, thus solving the problem of my chains remaining stationary. However, I received a “Segmentation Fault”. In addition to these new changes, I used “pull\_coord1\_geometry = direction”, which prompted me to use “pull\_coord1\_vec”, which I set to “0.0 0.0 4.5”. I chose 4.5 because I thought it will pull for 4.5 nm in the z axis only. Also, reading online forums I came to understand “pull\_pbc\_ref\_prev\_stom\_com=yes” reduces artifacts in pulling simulations. Actually, if I remove this line GROMACS provides the following error message: ”ERROR 1 [file md\_pull.mdp] The maximum distance from the chosen PBC atom (283) of pull group 1 to other atoms in the group is larger than 0.5 times half the box size. Set the pull-pbc-ref-prev-step-com option to yes.” ”ERROR 2 [file md\_pull.mdp] The maximum distance from the chosen PBC atom (867) of pull group 1 to other atoms in the group is larger than 0.5 times half the box size. Set the pull-pbc-ref-prev-step-com option to yes.” And if I remove the “pull\_group1\_pbcatom = 283” and “pull\_group2\_pbcatom = 867” line from the “md\_pull.mdp” file, I then recive this error: ”ERROR 1 [file md\_pull.mdp]: When the maximum distance from a pull group reference atom to other atoms in the group is larger than 0.5 times half the box size a centrally placed atom should be chosen as pbcatom. Pull group 1 is larger than that and does not have a specific atom selected as reference atom.” ”ERROR 2 [file md\_pull.mdp]: When the maximum distance from a pull group reference atom to other atoms in the group is larger than 0.5 times half the box size a centrally placed atom should be chosen as pbcatom. Pull group 2 is larger than that and does not have a specific atom selected as reference atom.” I really appreciate any help with fixing my Umbrella Sampling simulation, my end goal is to determine the binding energy between these chains ---
**Author:** ![hess](https://dub1.discourse-cdn.com/flex017/user_avatar/gromacs.bioexcel.eu/hess/32/416_2.png) [@hess](https://gromacs.bioexcel.eu/u/hess)\ **Post date:** [February 11, 2026, 8:25am UTC](https://gromacs.bioexcel.eu/t/steered-md-simulation-error/13144/12 "2026-02-11T08:25:41Z")
How can your pull group be larger than half the box? That doesn’t sound right. ---
**Author:** ![Aneesha.A](https://avatars.discourse-cdn.com/v4/letter/a/bc79bd/32.png) [@Aneesha.A](https://gromacs.bioexcel.eu/u/Aneesha.A)\ **Post date:** [February 15, 2026, 11:44pm UTC](https://gromacs.bioexcel.eu/t/steered-md-simulation-error/13144/13 "2026-02-15T23:44:19Z")
Hi Hess, I’m not sure why this is the case either. I’ve always run MD simulations with a defined unit cell like the tutorial “gmx editconf -f complex.gro -o newbox.gro -center 3.280 2.181 2.4775 -box 6.560 4.362 12”. After receiving this error, I did try enlarging the box with “gmx editconf -f complex.gro -o newbox.gro -center 3.280 2.181 2.4775 -box 8 8 12”. However, I received an error suggesting the system was not equilibrated well. What do you suggest the next step should be? Thanks! ---
**Author:** ![hess](https://dub1.discourse-cdn.com/flex017/user_avatar/gromacs.bioexcel.eu/hess/32/416_2.png) [@hess](https://gromacs.bioexcel.eu/u/hess)\ **Post date:** [February 16, 2026, 8:28am UTC](https://gromacs.bioexcel.eu/t/steered-md-simulation-error/13144/14 "2026-02-16T08:28:39Z")
But what is your pull group? If this is a whole protein, this could be larger than half the box size. ---
**Author:** ![Aneesha.A](https://avatars.discourse-cdn.com/v4/letter/a/bc79bd/32.png) [@Aneesha.A](https://gromacs.bioexcel.eu/u/Aneesha.A)\ **Post date:** [February 17, 2026, 6:49am UTC](https://gromacs.bioexcel.eu/t/steered-md-simulation-error/13144/15 "2026-02-17T06:49:07Z")
I’ve defined my pull group as two chains, since my system consists of two alpha-helix structures each 35 residues in length. From: “gmx make\_ndx -f npt.gro > ri 1-35 > ri 36-70 > name 19 Chain\_A > name 20 Chain\_B” The index.ndx file there looks like: ![Screenshot 2026-02-17 at 5.48.03 pm](https://europe1.discourse-cdn.com/flex017/uploads/bioexcel1/original/2X/8/8e2d37eeaf9e8fdf6bf347245b5967fc978c32a8.png) ---
**Author:** ![hess](https://dub1.discourse-cdn.com/flex017/user_avatar/gromacs.bioexcel.eu/hess/32/416_2.png) [@hess](https://gromacs.bioexcel.eu/u/hess)\ **Post date:** [February 17, 2026, 9:01am UTC](https://gromacs.bioexcel.eu/t/steered-md-simulation-error/13144/16 "2026-02-17T09:01:08Z")
I guess a helix could be longer than half the box size. If so, you can maybe select a pbc atom in the middle such that the distance of all atoms to this pbc atom is less than half the box size. But I would think the default pbc atom setting of the middle atom in the group should already result in a suitable atom. ---
**Author:** ![Aneesha.A](https://avatars.discourse-cdn.com/v4/letter/a/bc79bd/32.png) [@Aneesha.A](https://gromacs.bioexcel.eu/u/Aneesha.A)\ **Post date:** [February 17, 2026, 10:32am UTC](https://gromacs.bioexcel.eu/t/steered-md-simulation-error/13144/17 "2026-02-17T10:32:36Z")
So does this mean the “pull\_pbc\_ref\_prev\_step\_com=yes” with “pull\_group1\_pbcatom = 283 pull\_group2\_pbcatom = 867” setting in the md\_pull.mdp file should remain? If so, this takes me back to earlier simulations running into the same errors like “Segmentation Fault”. ---
**Author:** ![hess](https://dub1.discourse-cdn.com/flex017/user_avatar/gromacs.bioexcel.eu/hess/32/416_2.png) [@hess](https://gromacs.bioexcel.eu/u/hess)\ **Post date:** [February 17, 2026, 1:02pm UTC](https://gromacs.bioexcel.eu/t/steered-md-simulation-error/13144/18 "2026-02-17T13:02:45Z")
You should never get a segmentation fault. I would remove pull\_pbc\_ref\_prev\_step\_com=yes. I cannot say which pbc atom indices you should choose. ---
**Author:** ![Aneesha.A](https://avatars.discourse-cdn.com/v4/letter/a/bc79bd/32.png) [@Aneesha.A](https://gromacs.bioexcel.eu/u/Aneesha.A)\ **Post date:** [February 18, 2026, 1:04am UTC](https://gromacs.bioexcel.eu/t/steered-md-simulation-error/13144/19 "2026-02-18T01:04:46Z")
If not the middle atom in the group, how do I determine a suitable pbc atom? I mean to say is there a way of finding the best pbc atom index, or is it trial and error? ---
**Author:** ![hess](https://dub1.discourse-cdn.com/flex017/user_avatar/gromacs.bioexcel.eu/hess/32/416_2.png) [@hess](https://gromacs.bioexcel.eu/u/hess)\ **Post date:** [February 19, 2026, 3:24pm UTC](https://gromacs.bioexcel.eu/t/steered-md-simulation-error/13144/20 "2026-02-19T15:24:06Z")
You need to look for an atom that is geometrically in the middle. Often the middle atom work, but sometimes not. A worst case would be a beta-sheet consisting of two strands. The middle atom is then at one edge of the sheet.