VMD-L Mailing List

From: Axel Kohlmeyer (akohlmey_at_gmail.com)
Date: Thu Feb 18 2021 - 01:13:34 CST

Same problem as before. 

--
Dr. Axel Kohlmeyer akohlmey_at_gmail.com https://urldefense.com/v3/__http://goo.gl/1wk0__;!!DZ3fjg!vUrMzduLli0ZEXu_x7__X4CLWKDfGxQEqXCRjsC_biyldGcvX0h88lP7n4ccISp3OQ$ 
College of Science & Technology, Temple University, Philadelphia PA, USA
International Centre for Theoretical Physics, Trieste, Italy
On Thu, Feb 18, 2021, 01:39 Sunita Patel <sunita.patel_at_cbs.ac.in> wrote:
>  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;
>>> }
>>> }
>>> ====================================================
>>>
>>>