[gmx-users] R: Re: g_mindist on rhombic dodecahedron system

[Anna Marabotti]

Dear Justin, dear all,
following your suggestion I used the command:
trjconv -f prot_boxdodfull.xtc -s prot_boxdodfull.tpr -pbc mol -ur compact
-o prot_boxdodfull_mol.xtc
to convert my simulations, and as you anticipated all went OK: now my
trajectories are without spikes, the protein is entirely in the rhombic
dodecahedric box and the minimum distance is never lower than 2.5 nm, so I
think that the simulations went properly.
 
Now I have another doubt. In order to analyze my data (e.g. with g_rmsf,
g_rms etc), which is the more correct reference file? Usually I use as
reference file the .tpr input file of the full MD trajectory. I used this
file also this time with the command:
g_rms -f prot_boxdodfull_mol.xtc -s prot_boxdodfull.tpr -o
prot_boxdodfull_rms.xvg

Apparently, there are no main irregular behaviours (the RMS value oscillates
between 4.2 and 4.25 nm, so I think I can assume it is quite stable). The
absolute value of RMS is however very high with respect to the ones I used
to see (that are generally lower than 1 nm). I assume that the important
thing in this kind of analysis is the variation of the value, not the value
itself; however, I would like to know if I use the correct reference file or
if I have to create another reference file in which I "trjconv'ed" (how?)
also the reference structure. Could you please give me some suggestion about
my question?

Many thanks and best regards
Anna

-----Messaggio originale-----

Message: 5
Date: Tue, 05 Jul 2011 12:24:11 -0400
From: "Justin A. Lemkul" <jalemkul at vt.edu>
Subject: Re: [gmx-users] Re: g_mindist on rhombic dodecahedron system
To: Discussion list for GROMACS users <gmx-users at gromacs.org>
Message-ID: <4E133AAB.1030105 at vt.edu>
Content-Type: text/plain; charset=ISO-8859-1; format=flowed



Anna Marabotti wrote:
> Dear Tsjerk,
> thank you very much for your answer. I completely re-analyzed my
simulations
> and I'm telling you all the problems I found.
> First of all, please remember that this is a dimeric protein, formed by
two
> identical subunits, with no covalent interactions (no disulphide bridges
and
> so on).
> Here the steps I used for the simulation.
> 
> pdb2gmx -f prot.pdb -o prot.gro -p prot.top
> 
> I used editconf to create the rhombic dodecahedric box and to center the
> protein in the box. Here the exact string of commands:
> editconf -f prot.gro -o prot_boxdod.gro -bt dodecahedron -d 1.5 -c
> genbox -cp prot_boxdod.gro -cs -o prot_boxdodwat.gro -p prot.top
> If I look at the prot_boxdodwat.gro file in VMD, I see a cubic box in
which
> the dimeric protein is in a corner (and with a little part of the protein
> out of the box).
> What's wrong with this command? Why I don't see a rhombic box? Why the
> protein is not centered in the box? Why a part of the protein is outside
the
> box?
> 

The default representation is a triclinic cell.  You won't see a
dodecahedral 
representation unless you use trjconv -pbc mol -ur compact with a suitable
.tpr 
file.

> Using
> grompp -f em.mdp -c prot_boxdodwat.gro -o prot_boxdodwat.tpr -p prot.top
> trjconv -f prot_boxdodwat.gro -s prot_boxdodwat.tpr -pbc whole -ur compact
> -o prot_boxdodwat_whole.gro
> and looking at prot_boxdodwat_whole.gro with VMD I see exactly the same
> thing as before. Why I cannot see a rhombic box?
> 

Use -pbc mol, not whole.

> I neutralized the system:
> genion -s prot_boxdodwat.tpr -o prot_boxdodneu.gro -np 6 -pname NA -p
> prot.top
> and used
> grompp -f em.mdp -c prot_boxdodneu.gro -o prot_boxdodneu.tpr -p prot.top
> trjconv -f prot_boxdodneu.gro -s prot_boxdodneu.tpr -pbc whole -ur compact
> -o prot_boxdodneu_whole.gro
> and looking at prot_boxdodneu_whole.gro with VMD I see exactly the same
> thing: cubic box with dimeric protein in a corner (and a little part of
the
> protein out of box)
> 

See above.

> I minimized the system
> mdrun -s prot_boxdodneu.tpr -deffnm prot_boxdodmin
> (converged to Emtol < 500)
> and obtained the system after minimization. 
> Looking at prot_boxdodmin.gro file with VMD, now I see a rhombic box, but
my
> protein now is splitted into two monomers, one into the rhombic box (with
a
> little part outside) and the other monomer completely out of the box.
> Using 
> trjconv -f prot_boxdodmin.gro -s prot_boxdodneu.tpr -pbc whole -ur compact
> -o prot_boxdodmin_whole.gro
> and looking at the prot_boxdodmin_whole.gro file with VMD, nothing is
> changed.
> 

See above.

> Continuing with NVT and NPT, nothing changes: always I see the protein
> splitted into two monomers. Please remember that I NEVER see any error
> message before/during/after calculations.
> 

Apparent separation of monomers would not trigger an error, because it's not
an 
error.  It's simply a normal consequence of PBC and is especially common
with 
multimeric assemblies.

> After the full MD, instead, I look at the prot_dodboxfull.gro file and I
can
> see the dimeric protein into the rhombic box (with a little part outside
the
> box), but the trajectory behaves as I told before: spikes for the first 10
> ns and then nothing.
> The reason why I obtained a <1 distance with g_mindist is probably due to
> the reference I used, i.e. the .tpr file at the start of the simulation
> (obtained from the npt trajectory, with the monomers splitted and one of
> them outside the box).
> 
> I really didn't face a similar behaviour before, and I really don't know
how
> to cope with this protein. Any help will be appreciated.
> Sorry for the length of the message, I tried to be as precise as possible.

This has been a very useful and complete message and hopefully you'll now
solve 
your issues.

The apparent problem likely comes from the dimer splitting across periodic 
boundaries.  When both subunits are within the central unit cell (and not on

opposite sides), then your distances are what you would expect, on the order
of 
3 nm.  When the monomers appear to split, then the minimum periodic distance
is 
thus affected, because protein-protein distances appear to become smaller. 
Likely there is absolutely nothing wrong.  A simple test would be to
generate a 
reference configuration in which both protein subunits are centered in the
box 
and then follow suit with the trajectory (trjconv -center); a creative index

group may be necessary such that the dimers don't still split.
Geometrically, 
if monomer A is on one side and monomer B is on the other, the protein is
still 
centered but not useful.  Run g_mindist on the centered trajectory and you 
should have a consistent periodic distance and avoid the artifacts of the 
monomers jumping back and forth.

-Justin

> Anna
> 
> 
>