# 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]
	}
}