VMD-L Mailing List

From: Olaf Lenz (olenz_at_icp.uni-stuttgart.de)
Date: Thu May 22 2014 - 01:17:51 CDT

Hi!

2014-05-21 23:05 GMT+02:00 Hamid <hamidhasanzade_at_yahoo.com>:

> When I use VMD to visualize the traject.vtf file, the bonds do not cross
> the periodic boundaries. They cross the whole simulation box to connect two
> beads of a molecule
>

I am posting this reply also to the VMD mailing list, so that it is
archived and others can also profit from the reply. I have been asked this
question quite a few times, maybe it would make sense to include it into
some FAQ?

The problem with bonds crossing the whole box is not a problem of any
specific trajectory format, but a more generic problem. Ultimately, VMD
cannot *know* whether a bond between two atoms that are far apart is
supposed to cross the whole box or via the periodic boundaries. At most,
VMD could make the guess that a bond usually should connect an atom to its
minimal image in a periodic system.

There is no generic correct solution to avoid these overlong bonds, but
there are a few following possible ways around it:

1. If you have a simulation trajectory that starts without overlong bonds
but starts to develop them only during the simulation, you can use the
function "pbc unwrap" from the pbctools package. This function will prevent
that atoms are wrapped periodically. Like this, some of the atoms will walk
outside of the central image, but the bonds will be correct. If some
molecules walk too far, you can use "pbc wrap" with the correct "-compound"
argument to wrap the back into the central image.
Note that this method will not work when atoms can move further than half a
box size between two frames, e.g. when using parallel tempering (aka REMD),
or in some MC simulations.

2. The function "pbc join" can join compounds (i.e. structures of bonded
atoms) over periodic boundaries, and thus prevent overlong bonds. The
function is rather slow and cannot solve all possible problems, e.g. when a
compound stretches more than half a box.

3. You can simply remove those bonds that stretch the whole box. To do
that, you can use the following Tcl function. Copy it into your ~/.vmdrc
file, then you can call "remove_long_bonds 10" in the VMD console to remove
all bonds that are longer than 10. The number depends on your specific
system, it just needs to be longer than any real bond in the system. A
reasonable value would be half the size of the simulation box.
Beware that this has a few consequences. After the operation, the bonds are
gone from the system, so when you have a simulation over several timesteps,
the bonds will be gone even after the atoms have moved. Also, if you would
like to write the molecule's data to another file afterwards, the bonds
will not be included. Furthermore,

proc remove_long_bonds { max_length } {
    set n [ molinfo top get numatoms ]
    for { set i 0 } { $i < $n } { incr i } {
        set bead [ atomselect top "index $i" ]
        set bonds [ lindex [$bead getbonds] 0 ]

        if { $i % 1000 == 0 || $i == $n-1 } then {
            vmdcon -info "$i / $n"
        }

        if { [ llength bonds ] > 0 } {
            set bonds_new {}
            set xyz [ lindex [$bead get {x y z}] 0 ]

            foreach j $bonds {
                set bead_to [ atomselect top "index $j" ]
                set xyz_to [ lindex [$bead_to get {x y z}] 0 ]
                $bead_to delete
                if { [ vecdist $xyz $xyz_to ] < $max_length } {
                    lappend bonds_new $j
                }
            }
            $bead setbonds [ list $bonds_new ]
        }
        $bead delete
    }
    vmdcon -info "$n / $n"
}

Olaf 

-- 
Dr. rer. nat. Olaf Lenz
Institut für Computerphysik, Allmandring 3, D-70569 Stuttgart
Phone: +49-711-685-63607