VMD-L Mailing List

From: Myunggi Yi (myunggi_at_gmail.com)
Date: Tue Nov 21 2006 - 09:47:43 CST

Dear vmd users,

I'm running a script using text mode.

The script quit too early.
The following is my script.

vmd -dispdev text -eofexit < closest.scr > vmd.log

====closest.scr========
source set-chain.tcl
source closest.tk
set sel1 [atomselect 0 "resname HMT"]
set sel2 [atomselect 0 "noh and resid 27 to 37"]
closest $sel1 $sel2 closest-hmt.dat
=================

======set-chain.tcl===========
# loading trajectory
mol load parm7 hmt-in.parm7
mol addfile org1-10.trj type crdbox first 250 last 999 step 1
mol addfile org11-15.trj type crdbox

# setup chain ID for helices
[atomselect 0 "resid 1 to 27"] set chain A
[atomselect 0 "resid 28 to 54"] set chain B
[atomselect 0 "resid 55 to 81"] set chain C
[atomselect 0 "resid 82 to 108"] set chain D

# change resid according to the first unit
foreach a_list [[atomselect 0 "chain A"] get resid] {
    set new_list [expr $a_list + 20]
    lappend res_list $new_list
}
[atomselect 0 "chain A"] set resid $res_list
[atomselect 0 "chain B"] set resid $res_list
[atomselect 0 "chain C"] set resid $res_list
foreach a_list [[atomselect 0 "chain D"] get resid] {
    set new_list [expr $a_list - 61]
    lappend b_list $new_list
}
[atomselect 0 "chain D"] set resid $b_list
======================

=========closest.tk=============
# Usage: source closest.tk
# set sel1 [atomselect 0 "resname HMT"]
# set sel2 [atomselect 0 "noh and resid 27 to 37"]
# closest $sel1 $sel2 closest-hmt.dat

proc closest { sel1 sel2 file } {
    set fout [open $file w]
    set nf [molinfo [$sel1 molid] get numframes]
    set list1 [$sel1 list]
    set list2 [$sel2 list]
    set rilist [$sel2 get resid]
    set nlist1 [$sel1 get name]
    set nlist2 [$sel2 get name]

    # find distances between each pair
    for { set i 0 } { $i < $nf } { incr i } {
        $sel1 frame $i
        $sel2 frame $i
        set crd1 [$sel1 get {x y z}]
        set crd2 [$sel2 get {x y z}]
        set min 1000.0

        foreach atom1 $crd1 id1 $list1 {
            foreach atom2 $crd2 id2 $list2 {
                set dist [veclength [vecsub $atom2 $atom1]]
                if {$dist < $min} then {
                    set min $dist
                    set clidx1 $id1
                    set clidx2 $id2
                }
            }
        }
        puts $fout "[expr ($i + 1)*10] [lindex $nlist1 [lsearch -exact
$list1 $clidx1]] [lindex $rilist [lsearch -exact $list2 $clidx2]]
[lindex $nlist2 [lsearch -exact $list2 $clidx2]] $min"
    }
    close $fout
}
==================================

This generated result with only a few snap shots.

10 H13 30 O 2.55249404907
20 H1 35 N 2.5135371685
30 H2 31 CA 2.72125959396
40 H1 35 N 2.57328534126
50 H12 31 O 2.73753976822
60 H9 31 OG 2.11878705025

When I run "closest.tk" with GUI mode after loading trajectory,
It worked and generated results with all snap shots.

Do you have any idea?
What is wrong with my script? 

-- 
Best wishes,
MYUNGGI YI
==================================
KLB 419
Institute of Molecular Biophysics
Florida State University
Tallahassee, FL 32306
Office: (850) 645-1334
http://www.scs.fsu.edu/~myunggi