[gmx-users] First frame already out of box, getting very large RMSD

yunshi11 . yunshi09 at gmail.com
Mon Mar 2 03:03:40 CET 2015 

On Sun, Mar 1, 2015 at 5:42 PM, Justin Lemkul <jalemkul at vt.edu> wrote:

>
>
> On 3/1/15 8:26 PM, yunshi11 . wrote:
>
>> On Sun, Mar 1, 2015 at 4:36 PM, Justin Lemkul <jalemkul at vt.edu> wrote:
>>
>>
>>>
>>> On 3/1/15 3:47 PM, yunshi11 . wrote:
>>>
>>>  On Sun, Mar 1, 2015 at 11:43 AM, Justin Lemkul <jalemkul at vt.edu> wrote:
>>>>
>>>>
>>>>
>>>>> On 3/1/15 1:21 PM, yunshi11 . wrote:
>>>>>
>>>>>   On Sat, Feb 28, 2015 at 4:40 PM, Justin Lemkul <jalemkul at vt.edu>
>>>>> wrote:
>>>>>
>>>>>>
>>>>>>
>>>>>>
>>>>>>  On 2/28/15 7:38 PM, yunshi11 . wrote:
>>>>>>>
>>>>>>>    On Sat, Feb 28, 2015 at 4:21 PM, Justin Lemkul <jalemkul at vt.edu>
>>>>>>> wrote:
>>>>>>>
>>>>>>>
>>>>>>>>
>>>>>>>>
>>>>>>>>   On 2/28/15 7:17 PM, yunshi11 . wrote:
>>>>>>>>
>>>>>>>>>
>>>>>>>>>     On Sat, Feb 28, 2015 at 3:03 PM, Justin Lemkul <
>>>>>>>>> jalemkul at vt.edu>
>>>>>>>>> wrote:
>>>>>>>>>
>>>>>>>>>
>>>>>>>>>
>>>>>>>>>>
>>>>>>>>>>    On 2/28/15 6:00 PM, yunshi11 . wrote:
>>>>>>>>>>
>>>>>>>>>>
>>>>>>>>>>>      Dear all,
>>>>>>>>>>>
>>>>>>>>>>>
>>>>>>>>>>>   I am running MD for a protein-ligand complex in a dodecahedron
>>>>>>>>>>> box
>>>>>>>>>>>
>>>>>>>>>>>> and
>>>>>>>>>>>> followed the "Suggested trjconv workflow" from
>>>>>>>>>>>> http://www.gromacs.org/Documentation/Terminology/
>>>>>>>>>>>> Periodic_Boundary_Conditions
>>>>>>>>>>>> .
>>>>>>>>>>>>
>>>>>>>>>>>> Now I wonder how to remove jumps (across periodic boxes) when
>>>>>>>>>>>> the
>>>>>>>>>>>> first
>>>>>>>>>>>> frame (actually the .gro file from equilibration run) is already
>>>>>>>>>>>> jumping
>>>>>>>>>>>> across two (or more) boxes (according to visualization in VMD).
>>>>>>>>>>>>
>>>>>>>>>>>> I understand that this is just an visualizing artifact, but it
>>>>>>>>>>>> seems
>>>>>>>>>>>> to
>>>>>>>>>>>> also affect other analyses such as calculations of RMSD along
>>>>>>>>>>>> the
>>>>>>>>>>>> trajectory. I got some impossible RMSD values like 1.5 nm,
>>>>>>>>>>>> considering
>>>>>>>>>>>> other replicate simulations give only ~0.3 nm.
>>>>>>>>>>>>
>>>>>>>>>>>> Any idea on how to fix this (the effect on RMSD calculations
>>>>>>>>>>>> etc.,
>>>>>>>>>>>> not
>>>>>>>>>>>> visualization)?
>>>>>>>>>>>>
>>>>>>>>>>>>
>>>>>>>>>>>>      trjconv -center -pbc mol -ur compact usually solves all
>>>>>>>>>>>> problems
>>>>>>>>>>>> for
>>>>>>>>>>>>
>>>>>>>>>>>>    such
>>>>>>>>>>>>
>>>>>>>>>>>>  simple systems.  Center on the protein, everything else gets
>>>>>>>>>>> re-wrapped
>>>>>>>>>>> around it. More complex operations would be needed for protein
>>>>>>>>>>> complexes,
>>>>>>>>>>> membranes, etc. But proteins in water (with or without ligands)
>>>>>>>>>>> are
>>>>>>>>>>> generally straightforward.
>>>>>>>>>>>
>>>>>>>>>>>
>>>>>>>>>>>      I did try that, and other combinations of -pbc and -ur
>>>>>>>>>>> options.
>>>>>>>>>>> Sorry
>>>>>>>>>>>
>>>>>>>>>>>    that
>>>>>>>>>>>
>>>>>>>>>>>  I should have been more specific:
>>>>>>>>>>
>>>>>>>>>> My system is a homo dimer of protein, each with two ligands
>>>>>>>>>> (cofactor
>>>>>>>>>> +
>>>>>>>>>> inhibitor) bound symmetrically. Some trjconv options can put the
>>>>>>>>>> dimer
>>>>>>>>>> protein together, but the two ligands of the second monomer are
>>>>>>>>>> always
>>>>>>>>>> in
>>>>>>>>>> the other periodic boxes.
>>>>>>>>>>
>>>>>>>>>>
>>>>>>>>>>     This changes the approach considerably; always provide full
>>>>>>>>>> details
>>>>>>>>>> in
>>>>>>>>>>
>>>>>>>>>>   the
>>>>>>>>>>
>>>>>>>>> first message - you'll arrive at a solution faster!
>>>>>>>>>
>>>>>>>>> The centering in this case should be done with respect to one
>>>>>>>>> protein
>>>>>>>>> or
>>>>>>>>> some chosen custom group of residues that makes a logical center
>>>>>>>>> (e.g.
>>>>>>>>> the
>>>>>>>>> interfacial residues of one subunit).  The same concept applies -
>>>>>>>>> the
>>>>>>>>> remaining elements of the system (the other protein, the ligands,
>>>>>>>>> etc)
>>>>>>>>> will
>>>>>>>>> be wrapped with respect to the centered subunit.
>>>>>>>>>
>>>>>>>>>
>>>>>>>>>     With trjconv -s x.tpr -f x.gro -n x.ndx -center -pbc mol -ur
>>>>>>>>> compact
>>>>>>>>> -o
>>>>>>>>>
>>>>>>>>>   x.gro, I tried centering on monomer A, monomer B, and two
>>>>>>>>>
>>>>>>>> interfacial
>>>>>>>> residues of monomer A. None of them worked, and the two ligands of
>>>>>>>> monomer
>>>>>>>> B are always outside the box.
>>>>>>>>
>>>>>>>> Can I send .gro files to the mailing list?
>>>>>>>>
>>>>>>>>
>>>>>>>>    No. To share files you need to upload them to a file-sharing
>>>>>>>> service
>>>>>>>>
>>>>>>>>  and
>>>>>>> provide a link.  The input files (.tpr, .gro, and .ndx) would be
>>>>>>> useful.
>>>>>>>
>>>>>>>
>>>>>>> Please have a look at these files at
>>>>>>>
>>>>>>>   https://drive.google.com/folderview?id=0B8xuDbuW-
>>>>>>>
>>>>>> ACESnlTcUVrUEdKSGs&usp=sharing
>>>>>>
>>>>>> Please make the .tpr file with the .mdp, .top, and .gro files
>>>>>> provided,
>>>>>> since different versions of gromacs may not be compatible with all
>>>>>> .tpr
>>>>>> file versions.
>>>>>>
>>>>>>
>>>>>>   Well, I have just about every GROMACS version available to me, and
>>>>>> you
>>>>>>
>>>>> could just tell me what version you're using.  There's a reason I want
>>>>> to
>>>>> see your .tpr file; its contents are part of the consistency check I
>>>>> want
>>>>> to do. Moreover, all I'm going to get by using these files is generate
>>>>> a
>>>>> .tpr with the coordinates from nvt.gro; there's no purpose to fitting
>>>>> this
>>>>> structure to itself and there's no way to check anything this way.
>>>>>
>>>>> Please provide the .tpr, as requested, as well as provide the *exact*
>>>>> command you're using (no "x" substitutions; I need to know exactly what
>>>>> you're doing if you want me to try to reproduce or fix it).
>>>>>
>>>>>
>>>>>   I tried many commands that did not work, and the most recent one is:
>>>>>
>>>>
>>>> trjconv -s nvt.tpr -f nvt.gro -pbc mol -center -ur compact -n chain.ndx
>>>> -o
>>>> nvt2.gro, with selected indexes as:
>>>>
>>>> Select group for centering
>>>> ......
>>>> Select a group: 36
>>>> Selected 36: 'enzymeA'
>>>> Select group for output
>>>> ......
>>>> Select a group: 32
>>>> Selected 32: 'dimtprz'
>>>>
>>>>
>>>> Note nvt.tpr is generated with em2.gro, a structural file with no
>>>> periodicity and looks alright, and nvt.gro is the output from running
>>>> nvt.tpr
>>>>
>>>>
>>>>  trjconv -s nvt.tpr -f nvt.gro -o nvt2.gro -pbc nojump
>>> trjconv -s nvt.tpr -f nvt2.gro -o nvt3.gro -n r_90_187.ndx -pbc mol -ur
>>> compact -center
>>>
>>> That works for me, though the index file you provided doesn't have
>>> anything called "enzymeA" so I used r_90_187 as the group for centering.
>>>
>>>
>>>  Alright:\  Is there a reason that you did not use -pbc whole before
>> using
>> -pbc nojump?
>>
>>
> Because the molecules are already whole.  mdrun writes an intact structure
> when it writes the final frame.  Trajectories can be "broken" but the final
> structure is not.


Good to know that! But if I have this nvt3.gro, do I need to use it to make
subsequent trajectories whole? I usually just do a -pbc nojump followed by
-fit rot+trans for trajectories.


>
>
>  But when I used your exact commands (assuming you selected the entire
>> system including solvents for -pbc nojump) for npt, it did not work (files
>> at the same link
>> https://drive.google.com/folderview?id=0B8xuDbuW-
>> ACESnlTcUVrUEdKSGs&usp=sharing
>> ).
>>
>> trjconv -s npt.tpr -f npt.gro -o npt2.gro -pbc nojump
>> trjconv -s npt.tpr -f npt2.gro -o npt3.gro -n r_90_187.ndx -pbc mol -ur
>> compact -center
>>
>>
> Are you doing something to the coordinates before creating your .tpr
> files? Everything you've provided so far has shown the unit cell shifted or
> rotated with respect to the box.  If you're trying to mess with periodicity
> in between runs or before trying to do this fitting, don't.  The input
> coordinates in npt.tpr are unsuitable for removing jumps (as is clear from
> visualizing npt2.gro); for what reason I do not know.
>
>
That's weird, I basically followed the standard gromacs MD workflow (from
old tutorials), and no special manipulation of the coordinates at all. In
other words, only grompp and mdrun commands have been used from energy
minimization to production runs.

BTW, what is the difference between "unit cell" and (pbc?) "box"?



>  Does this mean I have to use nvt.tpr for fitting subsequent structural
>> files (.gro and .trr)?
>>
>
> Use whatever reference makes sense for what you're trying to do.
>
>
> -Justin
>
> --
> ==================================================
>
> Justin A. Lemkul, Ph.D.
> Ruth L. Kirschstein NRSA Postdoctoral Fellow
>
> Department of Pharmaceutical Sciences
> School of Pharmacy
> Health Sciences Facility II, Room 629
> University of Maryland, Baltimore
> 20 Penn St.
> Baltimore, MD 21201
>
> jalemkul at outerbanks.umaryland.edu | (410) 706-7441
> http://mackerell.umaryland.edu/~jalemkul
>
> ==================================================
> --
> Gromacs Users mailing list
>
> * Please search the archive at http://www.gromacs.org/
> Support/Mailing_Lists/GMX-Users_List before posting!
>
> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>
> * For (un)subscribe requests visit
> https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or
> send a mail to gmx-users-request at gromacs.org.
>

More information about the gromacs.org_gmx-users mailing list