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

[Justin A. Lemkul]

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
> 
> 
> 
> ------------------------------
> 
> Message: 4
> Date: Tue, 5 Jul 2011 15:06:11 +0200
> From: Tsjerk Wassenaar <tsjerkw at gmail.com>
> Subject: Re: [gmx-users] R: g_mindist on rhombic dodecahedron system
> To: Discussion list for GROMACS users <gmx-users at gromacs.org>
> Message-ID:
> 	<CABzE1Sg7n3EdyHLuTbJ=hWNdx0+1kr003E0gbuiC6RS78Nd6Pg at mail.gmail.com>
> Content-Type: text/plain; charset=ISO-8859-1
> 
> Hi Anna,
> 
> The spikes going down to 0.1 nm can only be explained by a fragmented
> representation of your protein in the frame. Mind that Gromacs does
> not write molecules in one piece as it used to. You therefore have to
> make all molecules whole prior to calling g_mindist. You mention doing
> so, which raises the question what your actual command lines were, in
> particular what you used as a reference structure. Or did you
> accidentally use the original trajectory again, rather than the one
> made whole? Please keep a record of the commands you used and paste
> them in your mail so that we can better understand what may be
> happening, rather than making (educated ;)) guesses.
> 
> Cheers,
> 
> Tsjerk
> 
> On Tue, Jul 5, 2011 at 2:49 PM, Anna Marabotti
> <anna.marabotti at isa.cnr.it> wrote:
>> Dear users,
>> can anybody give me suggestions about my questions below?
>> Thank you very much
>> Anna
>> ________________________________
>> Da: Anna Marabotti [mailto:anna.marabotti at isa.cnr.it]
>> Inviato: lunedl 4 luglio 2011 17.43
>> A: 'gmx-users at gromacs.org'
>> Oggetto: g_mindist on rhombic dodecahedron system
>>
>> Dear users,
>> following Tsjerk' suggestions, I simulate the protein which had some
>> problems in periodic distance violations in a triclinic box, using instead
> a
>> rhombic dodecahedron box created with the options:
>>
>> editconf -f prot.gro -o prot_boxdod.gro -bt dodecahedron -d 1.5 -c
>>
>> Minimization, PR-NVT and PR-NPT went OK (no error messages, potential
> energy
>> after minimization of about -1e+6). I run 30 ns full MD (no PR).
>> At the end of the production MD, I would like to repeat the g_mindist
>> analysis, so I can be sure my system does not suffer again of a periodic
>> distance violation.
>> I performed g_mindist on 4 different systems, and these are my results:
>> 1) on the resulting trajectory: on average, the distance between two
>> periodic distances calculated on the entire protein (option = 1) is > 3
> nm.
>> However, until 10 ns I can see several "spikes" going down to < 0.1 nm
>> 2) on the resulting trajectory, applying trjconv -pbc nojump before
>> g_mindist: on average, the distance is < 0.1 nm
>> 3) on the resulting trajectory, applying trjconv -pbc whole -ur compact
>> before g_mindist: the same as 1) (perfectly superimposable)
>> 4) on the resulting trajectory, applying trjconv -pbc whole -ur compact,
>> then trjconv -pbc nojump, before g_mindist calculation: the same as 2)
>> (perfectly superimposable)
>> My questions are:
>> 1) do I have to transform in some (other) way the trajectory before
>> calculating g_mindist? Do I have to calculate g_mindist in another way
> (eg.
>> using a different group for calculations)?
>> 2) given these results, can I consider this trajectory suitable for
>> analysis? (maybe excluding the first 10 ns)
>> 3) If not, what can I do more to have a suitable trajectory?
>>
>> Many thanks for any suggestions, I'm quite frustrating for not
> understanding
>> where the problem comes from...
>> Anna
>>
>> __________________________________________________________________
>> Anna Marabotti, Ph.D.
>> Laboratory of Bioinformatics and Computational Biology
>> Institute of Food Science - CNR
>> Via Roma, 64
>> 83100 Avellino
>> Phone: +39 0825 299651
>> Fax: +39 0825 781585
>> E-mail: amarabotti at isa.cnr.it
>> Skype account: annam1972
>> Web site: http://bioinformatica.isa.cnr.it/anna/anna.htm
>>
>> "When a man with a gun meets a man with a pen, the man with the gun is a
>> dead man"
>> (Roberto Benigni, about Roberto Saviano)
>>
>> --
>> gmx-users mailing list    gmx-users at gromacs.org
>> http://lists.gromacs.org/mailman/listinfo/gmx-users
>> Please search the archive at
>> http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
>> Please don't post (un)subscribe requests to the list. Use the
>> www interface or send it to gmx-users-request at gromacs.org.
>> Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>>
> 
> 
> 

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

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

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