VMD-L Mailing List

From: pl (ladam_at_smbh.univ-paris13.fr)
Date: Fri Nov 12 2004 - 07:59:48 CST

Hi John,

Well after imerging myself into tcl I finally find it's a cool
language and could do the job with the following macro I post
at the end.
But now I have to face another problem:
My macro extracts from each frame the solute plus H2O that is
Hbonded to it (just the donor or acceptot in fact), and saves a pdb
file.
At the end, I just cat all pdb files into one and read it in
vmd. This works nicely if I have a set of pdb files
that contains same amount of atoms and are separated by a END
command. But of course here, from frame to frame the amount of waters
fluctuates and the animation stops as soon as there is a change in
amount of atom.

I include here my macro (I am so proud of...) for other people that
would need it, this was built for AMBER pdb files with wrapped waters
in PBC (with either iwrap=1 or 'image' with ptraj ):

------------------------------------------------------------------------------
# DNA HBond donors and acceptors, we also need heteroatoms(G-N2,N1 . C-
N4 . T-N3 . A-N6) ->DNA
set AHD [atomselect 0 "(resname DG DG3 DG5 and name H1 H2 O6 N7 N3 N1
N2) or (resname DT DT3 DT5 and name O4 O2 N3) \
                    or (resname DA DA3 DA5 and name N1 N7 N3 N6) or
(resname DC DC3 DC5 and name N3 O2 N4) \
    or (name 1H2 2H2 H3 1H4 2H4 1H6 2H6)"]
set WAT [atomselect 0 water]
set DNA [atomselect 0 {all not resname WAT "Na+"}]

set cutoff 3.0
set angle 50.0
#measure hbonds $cutoff $angle $DNA $WAT

# -------------- Main trajectory Loop --------------

animate goto start
set num_steps [molinfo top get numframes]
puts "Processing $num_steps frames..."

for {set frame 0} {$frame < $num_steps} {incr frame} {

puts " Frame... $frame"
animate goto $frame
display update
# Measure HBonds between selections
set idx [measure hbonds $cutoff $angle $AHD $WAT]
# Now we extract the list of atom index from HBonds
set sel1 [atomselect top \
"index [lindex $idx 0] [lindex $idx 1] [lindex $idx 2]"]
set sel1list [$sel1 list]
set n [llength $sel1list]
#puts " Lenght of sel1list is $n "

# Now we extract the list of atom index from DNA
set sel2 [atomselect top {all not resname WAT "Na+"}]
set sel2list [$sel2 list]
set m [llength $sel2list]
#puts " Lenght of sel2list is $m "

# Erase duplicate atoms in lists (for writepdb later on)
for {set i 0} {$i < $n} {incr i} {
    #puts " Reading index $i "
    set atom1 [lindex $sel1list $i]
    foreach atom2 $sel2list {
    if {$atom1 == $atom2} {
#puts "Atom1: $atom1 is the same as Atom2: $atom2"
set sel1list [lreplace $sel1list $i $i ]
}
}
}
# Output length of HBonds list after cleaning...
set n [llength $sel1list]
#puts " Lenght of sel1list is now $n "

# Merge AHD list at the end of DNA atom list
for {set j 0} {$j < $n} {incr j} {
lappend sel2list [lindex $sel1list $j]
}
# Output length of Merged list...
set q [llength $sel2list]
#puts " Lenght of sel2list is now $q "

# Save as pdb
set final [atomselect top "index [lindex $sel2list]"]
$final writepdb final_$frame.pdb
#puts " EOLoop... $frame"

}
    
--------------------------------------------------------------------------
Bye

Le jeudi 11 novembre 2004 à 10:39 -0600, John Stone a écrit :
> Hi,
> There are some inconsistencies in the script you provided.
> Why are you doing atom selectionsin the foreach loop?
> You're not using the results from the atom selection for anything?
> $sel1, $sel2, $sel3 aren't used anywhere?
>
> The "writepdb" feature is a function of an atom selection, so you'd
> need to use something like $sel1 writepdb foo.pdb or something like that.
>
> If you want a list of integers from an atom selection, you'd do
> $sel1 list
>
> It isn't clear to me what you're actually trying to get out of all
> of this, are you trying to get a PDB file, or just a text file containing
> a list of integers?
>
> John Stone
> vmd_at_ks.uiuc.edu
>
> On Tue, Nov 09, 2004 at 12:52:51PM +0100, pl wrote:
> > Hi all,
> > I want to use the writepdb function to save a list of selected
> > atoms. To my understanding writepdb needs a list of integers and
> > my list is text. This list was created using the 'measure hbond'
> > command like this:
> >
> > set hblist [measure hbonds $cutoff $angle $DNA $WAT]
> > set atom1 {}
> > set atom2 {}
> > set atom3 {}
> > foreach atom1 [lindex $hblist 0] atom2 [lindex $hblist 1] atom3 [lindex
> > $hblist 2] {
> > set sel1 [atomselect top "index $atom1"]
> > set sel2 [atomselect top "index $atom2"]
> > set sel3 [atomselect top "index $atom3"]
> > lappend pdblist $atom1
> > lappend pdblist $atom2
> > lappend pdblist $atom3
> > }
> >
> > Now I want to throw this pdblist into a file using writepdb:
> >
> > $pdblist writepdb my_file.pdb
> >
> > As I said, this does not work as $pdblist cantains text and not
> > a list of integer.
> > How do I turn my list into integers. This is a dumb tcl question, boy
> > how I hate this language....
> >
> > Bye
> >
>