GROMACS version: 2021.3
GROMACS modification: No
Here post your question:
I ran an MD with a peptide in presence of POPG and POPE, the system was built using CHARMM-GUI with box-type HEXAGONAL, ran MD with Gromacs (2021.3). I am following instructions from the manual to design an optimal workflow for PBC connection. Few lipid molecules are still breaking and moving when I visualise in VMD.
Workflow I am using:
echo 0 | gmx trjconv -s step7_1.tpr -f step7_1.trr -o analysis2/step7_1_pbc_0_whole_trial7.xtc -pbc whole -tu ns
gmx make_ndx -f step5_input.pdb -o index2/index_POPG_POPE.ndx
“POPG” | “POPE” # and q to save and quit
echo 18 0 | gmx trjconv -s step7_1.tpr -f analysis2/step7_1_pbc_0_whole_trial7.xtc -o analysis2/step7_1_pbc_4_whole_cluster_trial7.xtc -n index2/index_POPG_POPE.ndx -pbc cluster -tu ns
echo 0 | gmx trjconv -s step7_1.tpr -f analysis2/step7_1_pbc_4_whole_cluster_trial7.xtc -o analysis2/step7_1_pbc_4_whole_nojump_cluster_trial7.xtc -pbc nojump -tu ns
4. centering using -pbc mol
echo 1 0 | gmx trjconv -s step7_1.tpr -f analysis2/step7_1_pbc_4_whole_nojump_cluster_trial7.xtc -o analysis2/step7_1_pbc_4_whole_nojump_cluster_center_compact_mol_trial7.xtc -pbc mol -center -ur compact -tu ns
vmd step7_1.gro analysis2/step7_1_pbc_4_whole_nojump_cluster_center_compact_mol_trial7.xtc
I am not able to figure out a way to keep those lipids intact.
I had the same issue with pure bilayers in small rectangular boxes. Something which worked for me:
- Center an atom in the middle of the bilayer with -pbc atom
- Center the bilayer with -pbc mol
For your case the workflow would be:
Make an index group with an atom in the middle of the bilayer (Last atom in an acyl chain for example) → e.g. CenterAtom.ndx
gmx trjconv -s step7_1.tpr -f step7_1.trr -o analysis2/step7_1_pbc_0_atom_trial7.xtc -pbc atom -tu ns -center -n CenterAtom.ndx → Choose first the group with the atom and then the group “System”
Make an index group with the whole bilayer → e.g. Bilayer.ndx
gmx trjconv -s step7_1.tpr -f analysis2/step7_1_pbc_0_atom_trial7.xtc -o analysis2/step7_1_pbc_0_mol_trial7.xtc -pbc mol -tu ns -center -n Bilayer.ndx → Choose first the group with the bilayer and then the group “System”
Hope that helps! :)
Thank you for your quick response.
I tried your suggestions for the workflow, it took me some time to understand making an index file for the last atom in the acyl chain.
Just wanted to check with you if I am doing it right. I did the following:
1. Making index group with an atom in the middle of the bilayer
gmx make_ndx -f …/step5_input.pdb -o CenterAtom.ndx
a C218 | a C316 # (selected atom C218 | atom C316 (which are at the end of POPE and POPG))
2. Centering for the selected atoms
gmx trjconv -s …/step7_1.tpr -f …/step7_1.trr -o step7_1_pbc_0_atom_trial10.xtc -pbc atom -tu ns -center -n CenterAtom.ndx # selected new group created by index above (C218_C316) for centering and ‘system’ for output
3. Making an index group with the whole bilayer e.g. Bilayer.ndx
gmx make_ndx -f …/step5_input.pdb -o Bilayer.ndx
“POPG” | “POPE” # Selecting both POPG and POPE
4. Centering for the selected mol (i.e. bilayer)
gmx trjconv -s …/step7_1.tpr -f step7_1_pbc_0_atom_trial10.xtc -o step7_1_pbc_0_mol_trial10.xtc -pbc mol -tu ns -center -n Bilayer.ndx # Selected POPG_POPE group created above
vmd …/step7_1.gro step7_1_pbc_0_mol_trial10.xtc
And this worked, this time bilayer (POPG and POPE) are together during the whole run.
The only problem left now is that protein is jumping around, I guess need to use nojump using the step7_1_pbc_0_mol_trial10.xtc, without any index file?
Again, thanks a lot for your help. This was driving me crazy.
nice to hear that it worked! :)
To center the protein in the box I would suggest to use again center to “fix” your protein at the center of your box:
gmx trjconv -s …/step7_1.tpr -f step7_1_pbc_0_mol_trial10.xtc -o step7_1_pbc_0_mol_protein_trial10.xtc -tu ns -center → Choose first the Protein Group and then the system.
Hope that works, too. It could fail maybe if your box in z-direction is too small…
I would not suggest using no jump since your lipids and the protein will then “diffuse” out of the box. If you are familiar with Python you could take a look at the trajectory transformation tools of the MDAnalysis library.
You were right again. nojump did not work when applied directly.
Following your suggestions, centred whole thing around the ‘protein’ and selected ‘system’ for output. It gave better results, but then the membrane was jumping. So, created an index file for protein and bilayer combined and centred for the combined index, but there was some jumps in the protein. It looks much better though. I guess need to use ‘rotate and translate’ for it to centre around the protein and stabilise, I will try.
I have limited knowledge of Python, but learning it and exploring MDAnalysis library as well, looks very promising.
I was wondering, how much of trajectory transformation is required for downstream analysis, I understand it may depend on the analysis. I intend to see the influence of peptide on bilayer, so probably LiPyphilic, which is based on MDAnalysis would be helpful, but was unsure about how much trajectory transformation (for removing PBC) before beginning those analysis.
Thanks a lot for the help again.