VMD-L Mailing List

From: Axel Kohlmeyer (akohlmey_at_gmail.com)
Date: Sun Oct 21 2012 - 11:06:32 CDT

On Sun, Oct 21, 2012 at 5:35 PM, yi <yihang2008_at_126.com> wrote:
> Dear VMD users
>
> I am using the LAMMPS dump command to dump a file, namely dump 1 all custom
> 500 dump.one id type x y z. Then I load the dump file into VMD, I found that
> water molecule around the pbc box border split, for example, the hydrogen

this is because of the data that you write out to the dump.

> atom at the left side of the box, the other two atoms of the same water
> molecule at the right side of the box. And I search Axel’s program code to
> make the molecule whole again, when this code was executed in VMD, the

when referring to code, please don't just refer to a person, but
provide a link to the mailing list archive or some other webpage
where you picked up this code. i have posted *so* many script
fragments to the VMD mailing list (and added code to VMD)
that "axel's code" is a *very* vague reference.

> molecule around the boundary still split. The attachment is one of my
> current image. Please let me know how to make the separate molecule whole?
> Thanks for your help. Axel’ code as follows:

>
> Introduction: In LAMMPS trajectories are PBC are applied per-atom in output.

this is incorrect. you can have unwrapped trajectories
by specifying xu yu zu or dumping the image flags.
have a closer look at the LAMMPS documentation.

> We can use VMD to make molecules whole again:

this code is overly convoluted and in the light of better
alternatives like the pbctools plugin should not be propagated.

what you should do instead is one of the two following options.

option a)
dump unwrapped coordinates. this will keep molecules
whole, but won't wrap molecules back into the principal
unit cell. this can be easily done using the pbctools
plugin by choosing the "wrap" subcommand and
the fragment option to the -component flag.

option b)
dump wrapped by atom coordinates as before.
now you first need to provide VMD with the information
of what *is* a molecule. your dump doesn't have that.
the easiest way would be to load your data file with
the topotools plugin and then create a suitable .psf file
and then load the .psf file and add the dump files on top
of that. with the residue information, you can use
pbc join on the first frame of your trajectory and then
use pbc unwrap to get an unwrapped trajectory and
then use pbc wrap -compount res to get the per molecule
wrapped coordinates like in option a). option a) is
obviously the simpler path. remember, the more information
you provide VMD, the easier it is to get it to do the
right thing. VMD can only know what you "tell" it.

https://sites.google.com/site/akohlmey/software/topotools/topotools-tutorial---various-tips-tricks#TOC-Reconstruct-a-.pdb-.psf-files-from-a-LAMMPS-data-file-for-use-with-visualization-programs

axel.

> set mol [molinfo top]
>
> set ring [atomselect $mol {type 3}]
>
> set nf [molinfo $mol get numframes]
>
> for {set i 0} {$i < $nf} {incr i} {
>
> make_whole make_whole $mol $ring $i 8
>
> }
>
> proc make_whole {mol sel frame num} {
>
> molinfo $mol set frame $frame
>
> $sel frame $frame
>
> set allcoord [$sel get {x y z}]
>
> set num1 [expr {$num - 1}]
>
> set boxhalf [vecscale 0.5 [molinfo $mol get {a b c}]]
>
> set newcoord {}
>
> while {[llength $allcoord] > 0} {
>
> set coord [lrange $allcoord 0 $num1]
>
> set allcoord [lrange $allcoord $num end]
>
> set ref [lindex $coord 0]
>
> lappend newcoord $ref
>
> foreach atom [lrange $coord 1 end] {
>
> set newatom {}
>
> set dist [vecsub $atom $ref]
>
> foreach d $dist b $boxhalf r $atom {
>
> if {$d < -$b} { set r [expr {$r + 2.0*$b}]}
>
> if {$d > $b} { set r [expr {$r - 2.0*$b}]}
>
> lappend newatom $r }
>
> lappend newcoord $newatom
>
> } }
>
> $sel set {x y z} $newcoord
>
> )
>
>
> 

-- 
Dr. Axel Kohlmeyer  akohlmey_at_gmail.com  http://goo.gl/1wk0
International Centre for Theoretical Physics, Trieste. Italy.