# Gromacs umbrella sampling error for pull force calculation **URL:** \ **Category:** User discussions\ **Tags:** mdp-parameters, mdrun, umbrella-sampling, pull-code\ **Created:** [October 30, 2024, 12:18pm UTC](https://gromacs.bioexcel.eu/t/gromacs-umbrella-sampling-error-for-pull-force-calculation/10453 "2024-10-30T12:18:15Z")\ **Posts on this page:** 5\ **Page:** 1
**Author:** ![pratyush](https://avatars.discourse-cdn.com/v4/letter/p/cc9497/32.png) [@pratyush](https://gromacs.bioexcel.eu/u/pratyush)\ **Post date:** [October 30, 2024, 12:18pm UTC](https://gromacs.bioexcel.eu/t/gromacs-umbrella-sampling-error-for-pull-force-calculation/10453/1 "2024-10-30T12:18:16Z")
GROMACS version: 2022 GROMACS modification: Yes/No ![Screenshot 2024-10-30 172314](https://europe1.discourse-cdn.com/flex017/uploads/bioexcel1/original/2X/a/a717cd602d2f99b4928c4df3a2602280f5032170.jpeg) So I guess this will be silly question but I am very much clueless, so please help me out. I have a system where the ligand is complexed with the dimer at their interface. While most of the interactions come from chain A, the chain B has a single interaction with the ligand. I want to calculate the energy required for dissociating the ligand from the dimer but I have some issue. The pullf.xvg seems fine upon plotting. but when I generate the snapshot of the very last time frame of the pull production, its parts of chain A that is being pulled instead of the ligand. The chain B as well as the ligand are exactly at their place and not moving though. I have done pull force calculation for protein dimers before for which I had to pull one of them from the other and had no such problem. So the problem here arises from me trying to pull the ligand from the interface of two chains. I am unable to figure out what is going wrong. I tried bunch of suggestions that I found from here and other forums but the problem remains exactly the same as before as the snapshot is being pulled in the very same way. Here is my pull code: pull = yes pull\_ncoords = 1 ; only one reaction coordinate pull\_ngroups = 2 ; two groups defining one reaction coordinate pull\_group1\_name = UNK pull\_group2\_name = chain\_A pull\_coord1\_type = umbrella ; harmonic potential pull\_coord1\_geometry = distance ; simple distance increase pull\_coord1\_dim = Y N N pull\_coord1\_groups = 1 2 pull\_coord1\_start = yes ; define initial COM distance \> 0 pull\_coord1\_rate = 0.01 ; 0.01 nm per ps = 10 nm per ns pull\_coord1\_k = 1000 ; kJ mol^-1 nm^-2 pull-group2-pbcatom = 16130 pull-pbc-ref-prev-step-com = yes The ligand (UNK) here is the mobile entity and the chain\_A is the stationary entity. The last two lines were added from suggestions were followed. The atom number 16130 is a central carbon atom of the ligand. And to put constraints on both the protein chains I used the following: title = Umbrella pulling simulation define = -DPOSRES\_A -DPOSRES\_B Please help me out in solving this issue. Thank you in advance for any heads up. ---
**Author:** ![MagnusL](https://dub1.discourse-cdn.com/flex017/user_avatar/gromacs.bioexcel.eu/magnusl/32/2380_2.png) [@MagnusL](https://gromacs.bioexcel.eu/u/MagnusL)\ **Post date:** [October 30, 2024, 12:38pm UTC](https://gromacs.bioexcel.eu/t/gromacs-umbrella-sampling-error-for-pull-force-calculation/10453/2 "2024-10-30T12:38:04Z")
I guess you mean that atom number 16130 is a central carbon of the protein, right? Otherwise you should change it to `pull-group1-pbcatom`, but I assume that the ligand does not have over 16000 atoms. Does it make any difference at all if you set `pull_coord1_groups = 2 1`? I.e., use the protein as reference group. ---
**Author:** ![pratyush](https://avatars.discourse-cdn.com/v4/letter/p/cc9497/32.png) [@pratyush](https://gromacs.bioexcel.eu/u/pratyush)\ **Post date:** [October 30, 2024, 1:14pm UTC](https://gromacs.bioexcel.eu/t/gromacs-umbrella-sampling-error-for-pull-force-calculation/10453/3 "2024-10-30T13:14:26Z")
Thank you for your quick reply. I realized the mistake and corrected it, upon which my code looks like this now: ; Pull code pull = yes pull\_ncoords = 1 ; only one reaction coordinate pull\_ngroups = 2 ; two groups defining one reaction coordinate pull\_group1\_name = UNK pull\_group2\_name = chain\_A pull\_coord1\_type = umbrella ; harmonic potential pull\_coord1\_geometry = distance ; simple distance increase pull\_coord1\_dim = Y N N pull\_coord1\_groups = 1 2 pull\_coord1\_start = yes ; define initial COM distance \> 0 pull\_coord1\_rate = 0.01 ; 0.01 nm per ps = 10 nm per ns pull\_coord1\_k = 1000 ; kJ mol^-1 nm^-2 pull-group1-pbcatom = 16130 pull-pbc-ref-prev-step-com = yes However, the problem remains. I am getting the same kind of xvg plot and the same parts of chain A being pulled. Also for your second question, I thought 1 2 should be okay (I honestly do not understand that particular option well). But now, I am going to try keeping it 2 1. I shall update. Thanks again. ---
**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:** [October 30, 2024, 1:19pm UTC](https://gromacs.bioexcel.eu/t/gromacs-umbrella-sampling-error-for-pull-force-calculation/10453/4 "2024-10-30T13:19:30Z")
With geometry=distance the order of the two pull groups does not matter. I would say that your results shows that is it easier to detach part of chain A then to detach the ligand, which could be the case. Or? ---
**Author:** ![MagnusL](https://dub1.discourse-cdn.com/flex017/user_avatar/gromacs.bioexcel.eu/magnusl/32/2380_2.png) [@MagnusL](https://gromacs.bioexcel.eu/u/MagnusL)\ **Post date:** [October 30, 2024, 1:31pm UTC](https://gromacs.bioexcel.eu/t/gromacs-umbrella-sampling-error-for-pull-force-calculation/10453/5 "2024-10-30T13:31:53Z")
I would recommend choosing a centrally located atom (remember that the indexing is within the molecule, not a global atom number) for group2 as well.