It might be sometimes interesting to know the equilibrium distance between two atom types.

For instance: when calulating molecule contacts in a trajectory. Depending of the atom type, the distance of equilibrium is different. It is not the same the contact distance between two waters (that can vary from 2.8 Angstroms if they make a Hydrogen Bond to about 3.5 if they are not hydrogen bonding), that a contact between two CT (a lipophillic interaction that can be around 3.9 angstroms). In this simple case, if we choose a single cutoff to identify all possible contacts in a system, I would take the larger one, but then we might be including second solvating shell of waters!

Thus, to make accurate calculations, one must identify first what are the atom types of the atoms we are treating and then know the equilibrium distance between them to decide if they make or not a contact.

This equilibrium distance for a pairwise interaction can be calculated as follows:

r = (2 * (A / B)) ^ (1/6)

Where A and B are the Lennard Jones coefficients for the pairwise interaction. Then  it is the 6th root of two times A divided by B.

Then basically from the topology file we have to extract A and B coefficients. The problem is the encryptation of this file that makes it very hard to acces certain information. There is two ways of procedure now depending on the information we have:

• If we know the index of the two atoms we are studing.
• If we don't know the index but we want to know the equilibrium distance between two atom types. I.e. 'CT' and 'OW' types.

Let's go step by step:

1. Identify the amber atom types (when we know the types and not the indexes).
   
   Each atom is ordered in the topology file and for each atom we have an entry under the %FLAG AMBER_ATOM_TYPE that tells you exactly the corresponding atom type. So if we want to know the index of the atom type 'CT' for instance, we just count where is it located in this list of values. We can have several of them but any of the index will be the same later on as the equilibrium distance depends on the type and not really on the atom index.
    
2. Get the atom type index
   
   Fine, we have the indexes for our atoms (either you knew it previously or you obtained them by browsing the amber_atom_type flag), then we have to find now the atom type index for them. It is an internal value that Amber uses to refer to the different atom types present in the system. Amber has a list of unique atom types of the system and this next list points to a position of this list. This information is stored in the %FLAG ATOM_TYPE_INDEX and has one entry for atom in the system.
   
   Let's say we want to calculate the equilibrium distance between atom 5 and atom 20 that are a CT and a OW. In this list we obtain index 2 for the atom 5 and a 7 for the atom 20. This means that internally, Amber has the CT type in position 2 of the types list and OW in the position 7 (thus all atoms in the system with type OW will end up in this position).
    
3. %FLAG NONBONDED_PARM_INDEX
   
   Under this flag we have a list of length the total of possible interactions between all the atom types in the system. So if we have 15 atom types in the system, the length of this list will be 15^2 = 225.
   But we only have 120 unique combinations, as order is not important: it is the same distance from CT to OH than from OH to CT.
   
   Exactly, this list provides the key to the unique combination list. So it's a way to transform the 225 possible interactions to the 120 actual possible values we have (also the actual number of A and B coefficients we have).
   
   Ok. Then how I know what is the number I have to look for? What is the corresponding list position for my two atom types interaction?
   
   This index we are looking for is obtained from the indexes we obtained previously:
   2 for CT (Index1) and 7 for OW (Index2). NDA is the Number of Different Atom Types in the system.
   
   NonBondedIndex = NDA * (Index1 - 1 ) + Index 2
   NI =  15 * (2 -1) + 7 = 22
   NI = 15 * (7-1) + 2 = 92
   Then we look for the value in those index and it will be the same for both of them: for instance 4 .
    
4. Finally get A and B coefficients
   
   Finally we know what are the corresponding A and B coefficients to the CT and OW interaction. We go to %FLAG LENNARD_JONES_ACOEF and %FLAG LENNARD_JONES_BCOEF and take the 4th value of each of them.
    
5. END :)
   
   We can now calculate the equilibrium distance with the 2 coefficients and we are done!

It's a bit more tricky when programming as we have to be very carefull with the list indexing stuff (depending on the language you use as some of them start with zero). For the rest, all is now easier.

In the next post I will explain how to apply it inside a code and how to easy calculate the different contacts taking into account this equilibrium distance using Biskit and numpy.