tclBC incorrect output

From: P.-L. Chau (pc104_at_pasteur.fr)
Date: Wed May 02 2012 - 22:18:21 CDT

I have adapted a TclBC script from the examples in the user-defined forces
tutorial. I find that the coordinates output by the script are incorrect,
and a search of the tutorial does not help me to understand why. Could I
ask for some help, please?

The tclBC script I have is as follows:

#!/usr/bin/tcl
set pdbSource nachr_popcwi_open.pdb
set aalist {} ;# a list of atom IDs of the 5 key amino acids
# chain A: Leu250
# chain B: Leu256
# chain C: Leu264
# chain D: Leu250
# chain E: Leu259
set atomID 1
set inStream [open $pdbSource r]
foreach line [split [read $inStream] \n] {
  set string1 [string range $line 13 14]
  set string2 [string range $line 17 19]
  set string2 [string trim $string2]
  set string3 [string range $line 21 21]
  set string4 [string range $line 23 25]
if { [string equal $string2 "LEU"] } \
{
if { [string equal $string3 "A"] && [string equal $string4 "250"] || \
     [string equal $string3 "B"] && [string equal $string4 "256"] || \
     [string equal $string3 "C"] && [string equal $string4 "264"] || \
     [string equal $string3 "D"] && [string equal $string4 "250"] || \
     [string equal $string3 "E"] && [string equal $string4 "259"] } \
{
if { [string equal $string1 "CA"] } {
puts "$string2 $string3 $string4 $string1"
lappend aalist $atomID
}
}
}
incr atomID
}
close $inStream
puts $aalist
#####################################################################
wrapmode cell
proc calcforces {step unique Rrate Rtarget K} {
    global aalist
    while {[nextatom]} {
        set atomid [getid]
        if { [lsearch $aalist $atomid] >= 0 } {
            set rvec [getcoord]
            puts $rvec
            foreach { x y z } $rvec {break}
            puts "$x $y $z"
        }
    }
}

and I have included the following in the NAMD conf file:

tclBC on
tclBCScript {
  set tclBCScript opentest4.tcl
  set pdbSource nachr_popcwi_open.pdb
  source $tclBCScript
}
tclBCArgs {15. 0.01 5.}

to drive the tclBC files. Examining the files manually, I know that the
first atom to be picked out should be the LEU A250 CA, which has its
coordinates as (-7.644,8.145,-4.409) but on execution by NAMD 2.8b2 the
ouptut is:

[...]
LEU A 250 CA
LEU B 256 CA
LEU C 264 CA
LEU D 250 CA
LEU E 259 CA
71170 77036 83081 89113 95134
Info: Startup phase 8 took 1.28413 s, 271.677 MB of memory in use
Info: Startup phase 9 took 0.000862837 s, 345.993 MB of memory in use
Info: Finished startup at 4.43505 s, 345.993 MB of memory in use

TCL: Running for 1 steps
1.4202225083 -10.3271917571 -3.90947041257
[...]

The coordinates are clearly incorrect. I have tried "wrapmode input" but
that does not change the output.

Could I ask for some advice on where I go wrong, please? Thank you very
much indeed.

P-L Chau

email: pc104_at_pasteur.fr
Bioinformatique Structurale
CNRS URA 2185
Institut Pasteur
75724 Paris
France
tel: +33 1 45688546
fax: +33 1 45688719

This archive was generated by hypermail 2.1.6 : Mon Dec 31 2012 - 23:21:29 CST