[gmx-users] PBC treatment: need an explanation

[Justin A. Lemkul]

On 4/24/12 6:51 AM, Anna Marabotti wrote:
> Dear gmx-users,
> I know that this is one of the most frequent subjects in the gmx-users list,
> however please let me ask you for a direct answer, since it seems to me that
> this particular question was not treated before.
> I'm performing MD simulations on a dimeric protein, using a rhombic dodecahedric
> box. I made 3 simulations in which my system was subjected to different
> isotropic pressures (first simulation: room pressure; second simulation: small
> increase of pressure; third simulation: big increase of pressure). I run 50 ns
> simulation for each system, and at the end of simulations I checked for the
> visualization of the system with VMD and for the RMSD against the starting
> configuration.
> Using g_rms command, I checked for the backbone RMSD against starting structure
> (fullMD.tpr file). The first system stabilized after a few ns of simulation, and
> then the RMSD remained constant. The second system stabilized after a few ns of
> simulation, but with a quantity of "spikes". The third system stabilized after a
> few ns of simulation and then, at about 30 ns of simulation, the RMSD value
> jumped on from approx. 0.4 nm to > 4 nm and stayed stable on that new value
> until the end of simulation.
> I had a look at this trajectory with VMD, and saw that the dimeric protein
> separates into two monomers. This phenomenon is consistent with some
> experimental data about the protein, and it seems to me consistent also with the
> RMSD trend found on the trajectory. However, due to visualization problems with
> my rhombic system, I decided to apply trjconv -pbc nojump:
> trjconv -s fullMD.tpr -f fullMD.xtc -o fullMD_noj.xtc -pbc nojump
> (choosing System=0 as option)
> After this action, I re-calculated the RMSD of the simulations using the same
> options as before...and found that in the third simultion the RMSD is no longer
> jumping on to > 4 nm. The visualization of the trajectory shows the protein in
> form of a dimer that fluctuates into the zone of "spreaded" solvent (no longer a
> box).
> My question is: was the separation into two monomers a simple artifact of the
> simulation, corrected by trjconv, or is trjconv able to affect the results of
> the system in such a way that when monomers truly separate, trjconv is able to
> "force" them together again? How can I check for these two possibilities?

The real answer is "neither."  The results you obtained are not an artifact; 
they are normal behavior for a periodic cell.  The visual representation created 
by trjconv was simply a re-wrapping of periodic boundaries so that it was more 
convenient.

If the monomers truly separated, there is no way trjconv (through the simple use 
of -pbc options) can force them back together.  There are many ways in which the 
user can manipulate coordinates; this isn't one of them.

You can check the stability of the dimer interactions in a number of ways - 
hydrogen bonds (g_hbond), distances (g_dist), contacts (g_mindist), etc.

-Justin

-- 
========================================

Justin A. Lemkul
Ph.D. Candidate
ICTAS Doctoral Scholar
MILES-IGERT Trainee
Department of Biochemistry
Virginia Tech
Blacksburg, VA
jalemkul[at]vt.edu | (540) 231-9080
http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin

========================================