VMD-L Mailing List

From: Sunita Patel (sunita.patel_at_cbs.ac.in)
Date: Thu Feb 18 2021 - 00:36:29 CST

 Dear VMD-users,

I have rewritten the vmd script considering your suggestions. After each
frame comparison I am deleting the list variable (unset var) and
re-assigning new values to the same variable for comparison with the next
frame. It is working fine up to 80 frames. For hundreds of frames it is
again crashing in GUI mode. I tried running it in vmdtext mode which is
working fine but it is very slow.

Could anybody see my code and recommend where I need to improve. The tcl
script is attached and also appended in this email. Should I send test
input files?

Thank you so much for your help.

Sincerely,
Sunita

==========
# Usage "source vmd_native_contacts.tcl" source vmd_native_contact 0.
# First frame is considered as reference native structure.
proc vmd_native_contacts { mol } {
puts "Enter output file name:"
set outputfilename [gets stdin]
set output [open $outputfilename w]
set allCa [atomselect $mol "name CA"]
#get the index of CA atoms
set rid [$allCa get index]
# No. of C-alphas
set ridLength [llength $rid]
set num [molinfo $mol get numframes]
set k 0 ;
# Iteration over frames
while {$k < $num} {
set allpairs {};
   set NC 0;
#compare the CA atoms pair which are separated by 6.5 Ang or less
for {set i 0} {$i < $ridLength} {incr i} {
for {set j [expr $i+4]} {$j < $ridLength} {incr j} {
set fir [lindex $rid $i];
set sec [lindex $rid $j];
set firNdx [atomselect $mol "index $fir" frame $k];
set secNdx [atomselect $mol "index $sec" frame $k];
set firCen [measure center $firNdx];
set secCen [measure center $secNdx];
set diff [vecsub $firCen $secCen];
set dist [veclength $diff];
  if {$dist < 6.5} {
if {$k == 0} {
    lappend refContacts "$fir $sec" ; }
             if {$k != 0} {
lappend allpairs "$fir $sec" ; }
}
}
}
        set LrefNC [llength $refContacts]
for {set n 0} {$n < $LrefNC} {incr n} {
        set refNC [lindex $refContacts $n];
set search [lsearch $allpairs "$refNC"]
if {$search != -1} {
set NC [expr $NC + 1]
}
}
unset allpairs;
# It gives frame versus number of native contacts.
set pcgNC [expr ([format %8.1f $NC] / [format %8.1f $LrefNC]) * 100]
puts $output "$k [format %8.1f $pcgNC]" ;
flush $output;
set k [expr $k + 1]
}
}
================================

On Tue, Feb 16, 2021 at 4:57 PM Ashar Malik <asharjm_at_gmail.com> wrote:

> Haven't had a detailed look through your code.
>
> Just a quick point though, that if you make a selection, do delete it when
> you are done with using it.
>
> This should be a fix if memory is the only issue with your code.
>
>
>
> On Tue, Feb 16, 2021 at 7:05 PM Sunita Patel <sunita.patel_at_cbs.ac.in>
> wrote:
>
>> Dear VMD-users,
>>
>> I wrote a tcl script to find the number of native contacts in each frame.
>> The reference frame (frame 0) whose native contacts have been matched with
>> other frames. This program is working fine with a small number of frames
>> (frames 30). However it crashed when I ran the program for 400 frames.
>>
>> I need your help to rectify the program so that it will work for a larger
>> number of frames. Please find the program appended below or also as an
>> attachment.
>>
>> Your help will be highly appreciated.
>>
>> Sincerely,
>> Sunita
>> ============================================================
>> # Usage "source vmd_native_contacts.tcl" source vmd_native_contact 0.
>> # First frame is considered as reference native structure.
>> proc vmd_native_contacts { mol } {
>> puts "Enter output file name:"
>> set outputfilename [gets stdin]
>> set output [open $outputfilename w]
>>
>> set allCa [atomselect $mol "name CA"]
>> #get the index of CA atoms
>> set rid [$allCa get index]
>> # No. of C-alphas
>> set ridLength [llength $rid]
>> set num [molinfo $mol get numframes]
>> set k 0 ;
>> set allcontacts {}
>> # Iteration over frames
>> while {$k < $num} {
>> set ctr 0
>> set contactCa {}
>> #compare the CA atoms pair which are separated by 6.5 Ang or less
>> for {set i 0} {$i < $ridLength} {incr i} {
>> for {set j [expr $i+4]} {$j < $ridLength} {incr j} {
>> set fir [lindex $rid $i];
>> set sec [lindex $rid $j];
>> set firNdx [atomselect $mol "index $fir" frame $k];
>> set secNdx [atomselect $mol "index $sec" frame $k];
>> set firCen [measure center $firNdx];
>> set secCen [measure center $secNdx];
>> set diff [vecsub $firCen $secCen];
>> set dist [veclength $diff];
>> if {$dist < 6.5} {
>> set ctr [expr $ctr + 1]
>> lappend contactCa "$fir $sec" ;
>> }
>> }
>> }
>> # CA-atomd distance less than 6.5 are put in the list 'allcontacts'.
>> lappend allcontacts "$contactCa" ;
>> set k [expr $k + 1]
>> }
>>
>> # No. of native contacts in the reference structure that is 'frame 0'.
>> set nrefNc [llength [lindex $allcontacts 0]];
>> # No. of frames 'fr'.
>> set fr [llength [lindex $allcontacts]];
>> # Searching of native contacts from frame 1 onwards (frame 0 is the
>> reference)
>> for {set m 1} {$m < $fr} {incr m} {
>> set Ncontact 0
>> for {set n 0} {$n < $nrefNc} {incr n} {
>> set refNc [lindex $allcontacts 0 $n];
>> set search [lsearch [lindex $allcontacts $m] "$refNc"]
>> if {$search != -1} {
>> set Ncontact [expr $Ncontact + 1]
>> }
>> }
>> set pcgNC [expr ([format %8.1f $Ncontact] / [format %8.1f $nrefNc]) * 100]
>> puts $output "$m [format %8.1f $pcgNC]";
>> flush $output;
>> }
>> }
>> ====================================================
>>
>>