Unable to correct PBC effect of MD trajectory of a tetrameric protein + ligand

[GuhLelouch]

GROMACS version: 2023.2
GROMACS modification: No

Dear GMX users,

I need to remove the PBC effect of my MD trajectory to perform MM/PB(GB)SA calculations. My system is made of a tetrameric protein (dimer of dimers) with one of the monomers bound to a ligand (for the binding energy calculation). Here you can see the complex:

The issue is that I have tried several approaches suggested in this forum, but none of them has worked. For instance, here is how the unprocessed trajectory looks like (I open the .gro file in PyMOL and then drop the .xtc file into it):

When I try to correct the PBC effect with

gmx_mpi trjconv -s md_0_1.tpr -f md_0_1.xtc -o md_0_1_noPBC.xtc -pbc mol -center -ur compact
1 (protein) 0 (system)

I get:

With

gmx_mpi trjconv -s md_0_1.tpr -f md_0_1.xtc -o md_0_1_noPBC2.xtc -pbc mol -center
1 (protein) 0 (system)

I get:

I created an index file with a group containing only the atoms of one of the chains, and then I centered using this group:

gmx_mpi trjconv -s md_0_1.tpr -f md_0_1.xtc -o md_0_1_noPBC.xtc -n index_pbc.ndx -pbc mol -center -ur compact
1 (Chain_A) 0 (system)

This improved the situation, but the complex was not correctly centralized in the box:

I tried mdvwhole using the command

mdvwhole -f md_0_1.gro -x md_0_1.xtc -o whole.xtc

but it didn’t help:

And with

mdvwhole -fmd_0_1.tpr -x md_0_1.xtc -mol True -o whole_sm.xtc

the issue persisted:

I would appreciate it if anyone could suggest a possible solution for this issue.

Best regards,

Gustavo

Edit:

I also tried to re-center (not sure which of the options below is appropriate) the trajectory previously “centered” using a single chain:

gmx_mpi trjconv -s md_0_1.tpr -f md_0_1_noPBC3.xtc -o md_0_1_noPBC7.xtc -center

gmx_mpi trjconv -s md_0_1.tpr -f md_0_1_noPBC3.xtc -o md_0_1_noPBC8.xtc -pbc mol -center

…but no success.

I also tried the “-pbc cluster” in two different ways:

gmx_mpi trjconv -s md_0_1.tpr -f md_0_1.xtc -o md_0_1_noPBC9.xtc -pbc cluster -center -ur compact
Cluster: 1 (Protein)
Center: 1 (Protein)
Output: 0 (System)

gmx_mpi trjconv -s md_0_1.tpr -f md_0_1.xtc -o md_0_1_noPBC10.xtc -n index.ndx -pbc cluster -center -ur compact
Cluster: 21 (Protein+Ligand)
Center: 1 (Protein)
Output: 0 (System)

…however, this also failed.

[milosz.wieczor]

What often helps in such cases is to (1) first -center using a single chosen atom that is as close to the middle point of the molecule as possible so that the complex gets translated to fit inside the box, and (2) run with -pbc mol that will put each molecule in its nearest periodic image.

It might take a while to find the right reference atom, but as long the complex doesn’t change its shape much, it should work eventually.

[GuhLelouch]

Dear milosz.wieczor,

Thanks for the input.

Just to clarify: in (1), do you mean a single atom that is near the center of the complex and that is part of a single protein chain?

If that is the case, then I would need to do (starting from the unprocessed trajectory):

1. Create an index file with a group containing only that atom.

gmx_mpi trjconv -s md_0_1.tpr -f md_0_1.xtc -o md_0_1_centered.xtc -n index.ndx -center
21 (Atom_X) 0 (system)

gmx_mpi trjconv -s md_0_1.tpr -f md_0_1_centered.xtc -o md_0_1_centered_noPBC.xtc -pbc mol
0 (system)

Does that make sense?

Should “-ur compact” be present in any of these commands to recover the dodecahedron? Or is it something that I have to try out if the previous commands are not sufficient?

Best regards,

Gustavo

Edit:

The commands I proposed above didn’t work, so I tried

gmx_mpi trjconv -s md_0_1.tpr -f md_0_1.xtc -o md_0_1_noPBC11.xtc -n index_atom.ndx -pbc mol -center -ur compact
21 (Atom_X) 0 (system)

This appears to have solved the issue (see the system below).

Here, “Atom_X” was a single atom near the center of the complex, in a single protein chain.

Now, I have the following questions: How do I know whether or not the complex is properly centralized? Does the centralization have to be perfect for correct analysis?

[milosz.wieczor]

Glad to see it worked!

If you plan to run an MM/PBSA calculation, then you’d usually want remove the solvent, and so the layout of the water box is irrelevant. For the solute, I’d say the most basic prerequisite for any energy calculation software is that it should be translation- and rotation- (or SE(3)-) invariant, so no, the details of centering shouldn’t affect your result at all.