VMD-L Mailing List

From: Ashar Malik (asharjm_at_gmail.com)
Date: Thu Feb 18 2021 - 01:02:58 CST

As a note all "atomselect" selections should be removed.

for example when done with allCa

use

$allCa delete

Another point is that you can make the nested loop better:

Currently its

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];

But there is no need for it to be like this - you make many many
firNdx selections unnecessarily. I don't understand the full logic of your
code from just a quick screen but I think it will give you the same result
if you restructure it like this. (Do test it).

for {set i 0} {$i < $ridLength} {incr i} {
    set fir [lindex $rid $i];
    set firNdx [atomselect $mol "index $fir" frame $k];
    for {set j [expr $i+4]} {$j < $ridLength} {incr j} {
        set sec [lindex $rid $j];
        set secNdx [atomselect $mol "index $sec" frame $k];
        ...
        ...
        $secNdx delete
    }
    $firNdx delete
}

This way you make far less selections and you remove when done.

Hopefully this will be a fix.

p.s. you should always run the script in text mode - so instead of writing
it as a proc - you can write a simple script which loads a
structure/trajectory file without any gui elements and runs the rest of
your code.

p.s.2. you can remove waters if you have them since you want only the
within protein contacts.

p.s.3 you can skip frames - depending on your simulation time step - if its
too small and you have too many frames processing consecutive frames might
give you unnecessary results.

On Thu, Feb 18, 2021 at 2:39 PM 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;
>>> }
>>> }
>>> ====================================================
>>>
>>>