Gromacs umbrella sampling error for pull force calculation

GROMACS version: 2022
GROMACS modification: Yes/No

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.

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.

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.

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?

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.