Compute the number of native contacts

Functions
	residueContacts cutoff query target
	residueContactPairs cutoff sel1 sel2
	firstContactTimeMatrix cutoff sel1 sel2
	prepareNativeContacts cutoff sel1 ?sel2?
	measureNativeContacts nclist cutoff sel1 ?sel2?

Detailed Description

Computation of the number of native contacts requires two steps:

1. creation of a "reference" list, i.e. the native state;
2. the actual measurement.

In the current version, both should come from the same structure. Contacts is based on VMD's measure contacts feature and therefore can be performed on one (intra-molecular) or two (inter-molecular) selections.

Note: For native contacts, the "reference" structure must be equivalent to the one to be analyzed (same number of atoms)

For a GUI to compute native contacts, see the RMSD trajectory tool + NC extension.

Example usage:

set chaina [ atomselect top "chain A"];
set ref [ prepareNativeContacts 7 $chaina];
# (go or select the frame of interest)
measureNativeContacts $ref 7 $chaina;

Another (deliberately verbose) example:

## The atom selection
set chaina [ atomselect top "chain A and name CA"]

## Use first trajectory frame as a reference
$chaina frame 0
set ref [ prepareNativeContacts 7 $chaina]

## Get the number of native contacts (for computing their fraction)
set nnc [ llength $ref]
puts "There are $nnc contacts in the native state"

              # Now, for each frame,
forFrames fno $chaina {
              # compute number of native contacts,
        set nc [measureNativeContacts $ref 7 $chaina]  
              # their fraction,
        set qnc [ expr 100.0 * $nc / $nnc ]
              # and print both.
        puts [ format "Frame %d: %f, %.3f%%" $fno $nc $qnc ]
}

Function Documentation

firstContactTimeMatrix()

firstContactTimeMatrix			cutoff sel1 sel2

Return a matrix (as a serialized array) with the time for first contact of each contact pair. Times may be missing if they never do a first contact! Example: set fctm [firstContactTimeMatrix 5 $a $b]

measureNativeContacts()

measureNativeContacts			nclist cutoff sel1 ?sel2?

TBD.

prepareNativeContacts()

prepareNativeContacts			cutoff sel1 ?sel2?

Prepare the "reference" list of native contacts, e.g. from the crystal structure. Return value: a list of native contact pairs (only useful to be passed as an argument to measureNativeContacts, or to get its length).

Note. The function is subject to change. Later can save sizes of atomselections (for safety checking) and/or cutoff.

residueContactPairs()

residueContactPairs			cutoff sel1 sel2

Like measureNativeContacts, but return resid's rather than indices. Return a list of (unique-ified) residue contacts as { R1 R2 } { R1 R3 } ... where the first element of a pair is a resid in sel1 and the second is in sel2. All computations are taken in the selections' frames.

residueContacts()

residueContacts			cutoff query target

Given the frame, compute whether each residue in query is making contact (< cutoff) with target. Return list like { { R1 0 } {R2 1 } ... } where R1, R2 etc belong to query