VMD-L Mailing List

From: Bob Johnson (robertjo_at_physics.upenn.edu)
Date: Fri Mar 16 2007 - 09:05:18 CDT

Hi Olaf,
Thanks for your messages. In response to your first email, no the new new
PBCtools unwrap feature doesn't work. I installed the vmd beta version in
order to try it. The unwrap tool basically gives the same result as the
Gromacs trjconv -pbc nojump tool: a highly distorted trajectory due to
erroneously detected jumps.

I'm sure you are familiar with the long bonds that are drawn whenever a
segment crosses into the next box and is wrapped back around into the
primary one. Thus, these long bonds are an indication of a "real jump". This
is why I thought a script that would detect these long bonds and correct for
it would be useful and error proof against MC (or in my case REMD swaps)
moves. I think that my need for a corrected REMD trajectory is pretty large.
The fact that my ssDNA molecule doesn't remain whole makes it difficult to
proceed with analysis. I'll take a look at your scripts and give this issue
more thought.
Thanks,
Bob

----- Original Message -----
From: "Olaf Lenz" <olenz_at_fias.uni-frankfurt.de>
To: "Bob Johnson" <robertjo_at_physics.upenn.edu>
Cc: "VMD Mailing List" <vmd-l_at_ks.uiuc.edu>
Sent: Friday, March 16, 2007 4:57 AM
Subject: Re: vmd-l: Scripts to unwrap PBC from REMD trajectory

> Hello again!
>
> After rereading your description of the problem, I'm not completely sure
> anymore whether my script solves your problem. What the new version of
> the plugin does is to compensate jumps in the whole coordinate system,
> as opposed to jumps in coordinates of single atoms.
>
>> However, in REMD there are large displacements between successive
>> frames due to the instantaenous coordinate swaps. Thus, attempting to
>> use the aforementioned tools on a REMD trajectory leads to detection
>> of erroneous "jumps".
>
> To see whether I have understood the problem, I will try to restate it:
> you're employing some kind of hybrid MC-like algorithm, where the
> coordinates of atoms are swapped? If so, also the new unwrapping
> procedure will not work correctly.
>
>> I had the idea of writing a script that would remove PBC wrapping based
>> on bond
>> length. If the bond in VMD gets larger than a certain amount, the script
>> would
>> detect this as a jump.
>> Then, the script would translate the appropriate atoms
>> in the appropriate direction to make the molecule whole again.
>
> In fact, I have already written the documentation and some pseudo code
> of a procedure "pbcjoin" that should join wrapped
> residues/chains/segments, however, I have not implemented it so far. I
> have attached the Tcl file containing the documentation so that you can
> have a look at it.
>
>> I'm assuming VMD
>> stores information about what atoms share a "chemical bond" (i.e. are
>> connected
>> by a line). Is there a way to access this information? My idea is still
>> somewhat
>> rough so any suggestions or insights on this would be appreciated!
>
> Usually, I would think that all atoms that are bonded belong either to a
> segment, a chain or a residue, so the actual bonding information is not
> really required. If this is not so, the bond information can be obtained
> from selections done via the "atomselect" command ("$sel getbonds"). It
> does make sense to extend the pbcjoin procedure so that it can also join
> any bonded structure, so I have added it to the file.
>
> I'm not sure whether I will manage to implement the procedure during the
> next days, but it may very well happen. How strong is your need?
>
> Best regards
> Olaf
>

--------------------------------------------------------------------------------

> ############################################################
> #
> # This file contains procedures to join compounds of atoms that are
> # wrapped around unit cell boundaries.
> #
> # $Id$
> #
>
> package provide pbctools 2.0
>
> namespace eval ::PBCTools:: {
> namespace export pbc*
> ############################################################
> #
> # pbcjoin $compound [OPTIONS...]
> #
> # Joins compounds of type $compound of atoms that have been
> # split due to wrapping around the unit cell boundaries, so that
> # they are not split anymore. $compound must be one of the values
> # "residue", "chain", "segment" or "bonded".
> #
> # OPTIONS:
> # -molid $molid|top the molecule to use (default: top)
> # -first $first|first|now
> # the first frame to use (default: now)
> # -last $last|last|now
> # the last frame to use (default: now)
> # -all|allframes equivalent to "-first first -last last"
> # -now equivalent to "-first now -last now"
> # -sel $sel the selection of atoms to be joined
> # (default: all)
> # -noref|-ref $sel the reference atom in the compound:
> #
> # AUTHOR: Olaf
> #
> proc pbcjoin { compound args } {
> # loop over all frames
> # loop over all compounds
> # get the coordinates of all atoms in the compound
> # get the coordinates of the reference atom in the compound
> # wrap the coordinates of the compound into the box
> # centered around the reference atom
> }
>
> }
>