VMD-L Mailing List

From: Sunita Patel (sunita.patel_at_cbs.ac.in)
Date: Tue Feb 16 2021 - 05:01:30 CST

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;
}
}
====================================================