PBCs in a slightly complex system

Hi everyone,

I have a system made up of:

  1. Two protein chains (not interacting, that was the purpose of the simulation to see if they do).
  2. ~ 100 molecules of hexanediol.
  3. water and ions.

after the simulatioin is complete, I was trying to remove the artifacts caused by the PBCs. I tried different tricks and combinations, but nothing worked as intended (we need everything to be in its original shape/form, no diffusion, no any sort of PBCs). There is always a problem left with something. Here is what I did

; I extracted a sample to try things on

gmx trjconv -s production_run.tpr -f production_run.xtc -o sample.xtc -skip 10
0 (system)

First, I tried
gmx trjconv -s production_run.tpr -f sample.xtc -o mol.xtc -pbc mol -ur compact -center -trans 10 10 10
0
(protein chains still jumping)
gmx trjconv -s production_run.tpr -f sample.xtc -o mol.xtc -pbc res -ur compact -center -trans 10 10 10
0
(protein chains still stretching)

; Then I tried a different approach based on the suggested workflow (here), I made it whole

gmx trjconv -s production_run.tpr -f sample.xtc -o made_whole.xtc -pbc whole
0

up to this stage, water molecules, hexanediol, and protein are fine without stretching or any funny appearance. The only thing is that the protein chains are still jumping. Here is what I tried

; extracted the first frame
gmx trjconv -s production_run.tpr -f sample.xtc -o first_frame.pdb -dump 0
0

; I invoked the nojump option:

gmx trjconv -s first_frame.pdb -f made_whole.xtc -o whole_no_jump.xtc -pbc nojump
0

It destroys the box (it’s now a sphere instead of the 10x10x10 cube), the proteins, the hexanediol molecules, and the water molecules. Everything, except ions now has artefacts.

I tried to wrap it back to a box
gmx trjconv -s production_run.tpr -f whole_no_jump.xtc -o whole_no_jump_centered.xtc -center -n index.ndx
X (I tried centered on the system, the hexanediol, the protein, a single chain even though they are not interacting), nothing worked!

The option -pbc mol -ur compact introduces the jumping again. Even the other option for -pbc (res).

I even tried

Any recommendation is highly appreciated.

Kind regards,
Jameel

You cannot both have everything in the box and not jumping.

If you are lucky, you might be able to get a non-jumping pair of proteins. But that is not guaranteed, as it will depend on how the proteins move. For that you need to first run -pbc nojump, without the whole step before (assuming molecules are whole in your tpr file). Then center the result on the COM of the two proteins. And then run a third step with -pbc mol -ur compact.

Thank you @hess so much for your reply!

I tried that before and now after receiving your reply. To double check that the molecules are not broken in the production_run.tpr file, I invoked

gmx editconf -f production_run.tpr -o prod.gro

and visualized the structure in PyMOL, nothing is broken at all.

I executed the commands you recommend in the correct order, but I still get PBC artifacts on water molecules from all sides of the box that gradually become severe as the trajectory continue. I tried the -pbc mol -ur compact with and without -trans 10 10 10 but the outcome is the same as in the image below.

Thanks again and kind regards,

Jameel

Your last step cannot have been -pbc mol, as then molecules would not be broken over PBC.

Dear @hess ,

Thanks for your quick responses!

The command indeed included -pbc mol! These were the invoked command:
gmx trjconv -s production_run.tpr -f sample.xtc -o nojump.xtc -pbc nojump
0 (system)
gmx trjconv -s production_run.tpr -f nojump.xtc -o nojump_centered.xtc -center
1 (protein)
gmx trjconv -s production_run.tpr -f nojump_centered.xtc -o nojump_centered_mol.xtc -pbc mol -ur compact
0 (system)

If you would like to test it on your side, here is a link to a zip file that contains a sample of the raw .xtc, .tpr, and .gro file.

Thank you in advance and kind regards,

Jameel

It took me some time to figure out the issue. The problem is that your trajectory has a time step of 1 ns. And in 1 ns atoms will now and then move more than half a both length. So you cannot use the nojump option. You should still be able to first do -pbc mol, then center and then -pbc mol again. Maybe you get the result that you want.

Dear @hess ,

Thank you so much for your time, help, and support.

I think there is a some sort of misunderstanding or my message was not clear enough. The sample.xtc was created by skipping 10 frames from the original trajectory (as seen in the first command in my first message in this thread.

The time step for the simulation is actually 0.002 ps. This value is from the prod_run.mdp prod_run.mdp and from the .tpr file via

gmx dump -s production_run.tpr | grep -i "dt"

Since the nstxout-compressed = 50000 in the mdp file, I believe I am saving a frame every 100 ps (0.002 * 50000) not every 1000 ps. So when I extracted the sample.xtc, the spacing got magnified by 10 which made it to looks like as if the time step is 1 ns.

Anyway, I ran the following commands on the original trajectory (not the sample) and I still see the jumping.

gmx trjconv -s production_run.tpr -f production_run.xtc -o mol.xtc -pbc mol
0 (system)
gmx trjconv -s production_run.tpr -f mol.xtc -o mol_centered.xtc -center
1 (protein)- I also tried system and hexanediol
gmx trjconv -s production_run.tpr -f mol_centered.xtc -o mol_centered_again.xtc -pbc mol -ur compact -trans 10 10 10
0 (system)

Although I know this has been an annoyingly long thread, and sorry about that, I have created another continuous sample from the original xtc and zipped it here if you still interested in looking into the issue.

This continuous sample was extracted via
gmx trjconv -s production_run.tpr -f production_run.xtc -o continuous_sample.xtc -b 50000 -e 60000

Thanks again and sorry for any inconvenience caused.

Kind regards,
Jameel

As I wrote, there is no guarantee that you can get what you want. The proteins can move apart by more than half a box length and then interact with a different periodic copy. Such jumps can not be removed. I think this is the case in your trajectory.