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

Justin A. Lemkul jalemkul at vt.edu
Mon Jul 11 16:29:51 CEST 2011

Anna Marabotti wrote:
> 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?

If you've manipulated the trajectory with trjconv, then you need a corresponding 
reference frame for position-dependent quantities.  For RMSD, RMSF, etc I 
usually do something like:

editconf -f start.tpr -o 0ns.gro
trjconv -s start.tpr -f 0ns.gro -pbc mol -ur compact -o 0ns_fix.gro
g_rms -s 0ns_fix.gro -f traj_fix.xtc

Other manipulations may be necessary if the protein is a dimer, etc.


> 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


Justin A. Lemkul
Ph.D. Candidate
ICTAS Doctoral Scholar
Department of Biochemistry
Virginia Tech
Blacksburg, VA
jalemkul[at]vt.edu | (540) 231-9080


More information about the gromacs.org_gmx-users mailing list