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