VMD-L Mailing List

From: Jianhui Tian (jianhuitian_at_gmail.com)
Date: Fri Mar 30 2007 - 14:49:45 CDT

Hi all,

I write a tcl script to calculate the helix content of the peptide.
But there is a small bug I can't find out myself to make it run
properly. Can someone help? Thanks a lot.

Jianhui
##########################
set dcdFile n_aka3rm14w6_run10.dcd ;# name of your psf file

set parm aka3rm14w6.prmtop ;# name of parm file

set incrm 1
##########################
mol load parm7 $parm dcd $dcdFile

set num_steps [molinfo top get numframes]
puts $num_steps
set out [open helixcontent.dat w]

for {set frme 1} {$frme <= $num_steps} {incr frme $incrm} {
    puts $frme

    set hlxflag 0
    set hlxcont 0

###### From resid 125 to 148 is amino acid.
    for {set i 125} {$i <= 148 } { incr i} {
        set atom [atomselect top "protein and name CA and resid $i"]
        $atom frame $frme

        set phi [$atom get phi]
        set psi [$atom get psi] if { $phi <= -30 && $phi >= -90
&& $psi <= -17 && $psi >= -77 } {
            if {$hlxflag == 0} {
                 set num 1
                 set hlxflag 1
            } elseif {$hlxflag == 1} {
                 set num [expr $num + 1]
            }
        } else {
            if {$hlxflag == 1} {
               set hlxflag 0
               if { $num > 2 } { set hlxcont [expr $num - 2]
            }
        }

        $atom delete
    }
         set cont [expr $hlxcont / 24]
         puts $out "$frme $cont"
}
exit