VMD-L Mailing List

From: Martin, Erik W (Erik.Martin_at_stjude.org)
Date: Fri Jan 04 2013 - 14:05:14 CST

For the record, since I first sent this a couple hours ago, I fixed most
of the problems thusly:

set filename "contactstest"
puts $filename
set cutoff 0.75
set numfr [molinfo $traj get numframes]
set firstframe 0
set findP [atomselect $traj "protein and name P"]
set Pres [$findP get resid]

set wd40 [atomselect $reference "chain D"]
set refsel [atomselect $reference "chain F and resid 4 and name N P OG1 or
resid 5 and name N CA CG"]
set lastframe [expr {$numfr-1}]
puts "last frame: $lastframe"
set mv [atomselect $traj all]

        set n 0
        foreach j $Pres {
                incr n
                set j2 [expr {$j+1}]
                                set tempsel [atomselect $traj "protein and resid $j and name N P OG or
resid $j2 and name N CA CG"]
                                
                for {set i $firstframe} {$i < $numfr} {incr i} {
                        puts "Measuring P site $j in frame $i"
                        set outfile [open $filename$n.csv a]
                        $mv frame $i
                        $tempsel frame $i
                        set M [measure fit $tempsel $refsel]
                        $mv move $M
                        set tempsel2 [atomselect $traj "protein and noh" frame $i]
                        set tmp [measure contacts $cutoff $wd40 $tempsel2]
                                                
                        set tmp2 [llength [lindex $tmp 0]]
                        puts "$tmp2"
                        puts $outfile "$i, $tmp2"
                        close $outfile
                        }
                        
                        }

But It still runs really slowly, so any other tips are always appreciated!

Thanks,
Erik

On 1/4/13 11:17 AM, "Martin, Erik W" <Erik.Martin_at_stjude.org> wrote:

>Hi, I quickly wrote this script that I was hoping would perform the
>following function: align each phosphorylated residue in a peptide with
>a ligand in a crystal structure and loop through frames in a trajectory
>measuring contacts between the simulated peptide and the protein its
>"bound" to. This seemed simple enough, but I'm running into all kinds
>of problems... The output file only contains the info from the last
>frame, the "tempsel" atomselection does not seem to be updating (or at
>least the alignment I'm doing doesn't work. Finally, it does not seem to
>actually be closing the files that I label as "outfile" and I get this
>error: couldn't open "contactstest4.csv": too many open files
>
>If anyone enjoys writing TCL scripts for trajectory analysis and would
>like to help me with this troubleshooting I'd appreciate it! I thought
>this was really straightforward...
>
>set filename "contactstest"
>puts $filename
>set cutoff 1
>set numfr [molinfo $traj get numframes]
>set firstframe 0
>set findP [atomselect $traj "protein and name P"]
>set Pres [$findP get resid]
>set mv [atomselect $traj all]
>set wd40 [atomselect $reference "chain D"]
>set refsel [atomselect $reference "chain F and resid 4 and name P OG1 or
>resid 5 and name N CA CG"]
>set lastframe [expr {$numfr-1}]
>puts "last frame: $lastframe"
>
> set n 1
> foreach j $Pres {
> incr n
> set j2 [expr {$j+1}]
> set tempsel [atomselect $traj "protein and resid $j and
>name P OG or resid $j2 and name N CA CG"]
> puts "test $j"
> for {set i $firstframe} {$i < $numfr} {incr i} {
> puts "testA $i"
> set outfile [open $filename$n.csv w]
> $tempsel frame $i
> set M [measure fit $tempsel $refsel]
> $mv move $M
> set tempsel2 [atomselect $traj "protein and noh"]
> set tmp [measure contact $cutoff $wd40 $tempsel2]
> set tmp2 [llength [lindex $tmp 0]]
>
> puts $outfile "$i, $tmp2"
>
> }
> close $outfile
> }
>
>Email Disclaimer: www.stjude.org/emaildisclaimer
>Consultation Disclaimer: www.stjude.org/consultationdisclaimer
>
>
>