<html><head><meta name="color-scheme" content="light dark"></head><body><pre style="word-wrap: break-word; white-space: pre-wrap;">#!/usr/bin/env python

"""
Title: CHASA.py
Author: Pat Fleming, pat.fleming@jhu.edu
Date: Fall, 2004
Requirements: Numeric

Updated: Sun Jan 23 11:35:53 EST 2005

Calculate Conditional Hydrophobic Accessible Surface Area (CHASA), backbone
N and O solvation accessiblity and total solvation free energy for polypeptides. 

Also identify unsatisfied hydrogen bond donors (N) and acceptors (O) in backbone.

This script calculates the accessible surface area (ASA) conditional upon prior 
solvation of the backbone N and O atoms by placing an oxygen atom proximate to 
these backbone atoms at acceptable hydrogen bond distance and orientation and 
then including these solvation "waters" in the ASA calculation.

The CHASA method and calculation of total solvation energy are described in:
------------------------------------------------------------------------
Patrick J. Fleming, Nicholas C. Fitzkee, Mihaly Mezei, Rajgopal Srinivasan,
and George D. Rose (2004)

"A novel method reveals that solvent water favors polyproline II over
beta-strand conformation in peptides and unfolded proteins:
Conditional Hydrophobic Accessible Surface Area (CHASA)"

Proteins (2005) 14:111-118
------------------------------------------------------------------------
Input:

A standard PDB file is the input. The current version does not deal well with
hydrogens and they should be stripped. (The code deals only with backbone hydrogens
used in current versions of LINUS).

The N-terminal nitrogen is not solvated and thus the CHASA for some atoms proximate to 
this atom in cartesian space will be incorrect.

Output:

The output is in the form of a PDB file with results in the occupancy and B-factor
columns and header and TER records.

Example:

COMPND numintHbd     numSolv    num_nonHbd      total_bb_polar
COMPND     84          147           3              111
.
.
.
ATOM      9  N   THR     2      13.719  19.413  27.573  5.00  0.00
ATOM     10  CA  THR     2      13.088  19.661  26.283  0.00  0.00
ATOM     11  C   THR     2      13.561  18.631  25.300  0.00  0.00
ATOM     12  O   THR     2      14.763  18.432  25.121  3.00  1.57
ATOM     13  CB  THR     2      13.527  20.980  25.667  0.00  7.80
ATOM     14  OG1 THR     2      13.307  22.020  26.627  0.00  5.84
ATOM     15  CG2 THR     2      12.704  21.284  24.409  0.00 16.10
.
.
.
ATOM    149  N   ALA    20       9.346  17.206  29.144  0.00  0.49
ATOM    150  CA  ALA    20       8.985  15.930  29.750  0.00  2.64
ATOM    151  C   ALA    20      10.067  15.607  30.760  0.00  3.52
ATOM    152  O   ALA    20      11.193  16.119  30.686 -1.00  1.77
ATOM    153  CB  ALA    20       8.856  14.815  28.714  0.00  2.39
.
.
.
ATOM    499  O   HOH   120       7.787   1.168  13.283  0.00  0.00
ATOM    500  O   HOH   121       9.368  -0.535   6.166  0.00  0.00
ATOM    501  O   HOH   122      22.918  12.605  15.164  0.00  0.00
ATOM    502  O   HOH   123      16.522  12.543  26.647  0.00  0.00
ATOM    503  O   HOH   124      22.421   8.292  28.255  0.00  0.00
ATOM    504  O   HOH   125      11.087  12.195  10.474  0.00  0.00
TER 1424.592  -18.062

where:

	numintHbd = number of internally hydrogen bonded backbone N and O

	numSolv = number of solvation "waters" (max = 5 per backbone polar group)

	num_nonHbd = number of backbone polar groups not satisfied by hydrogen bonding

        total_bb_polar = number of backbone polar groups in the molecule that should
                         be hydrogen bonded ([2 x Number of residues] -1)

        occupancy column (for N and O only):
            (if &gt; 0.00) = number of solvation "waters" accessible to that atom (max = 5)
            (if = -1.00) = this atom is not hydrogen bond satisfied

        B factor column = CHASA in square angstroms

        TER record:
            (first number) = total CHASA for molecule
            (second number) = solvation free energy for molecule

        HOH atoms = solvation "waters" in hydrogen bonding proximity to backbone polar groups
                    (these are the conditional atoms added prior to ASA calculation)

------------------------------------------------------------------------
Much of the code here is modified from the LINUS suite of programs originally 
written by Raj Srinivasan (http://roselab.jhu.edu/dist/).

The code for calculation of ASA by the method of Shrake and Rupley was ported from
C code originally written by Nick Fitzkee.

The script requires that Numeric be installed.
------------------------------------------------------------------------
"""

USAGE = """
python chasa.py pdbinfile &gt; [pdboutfile]

For non-conditional (traditional) ASA uncomment the following line
near bottom of script:

#   wat_list = None

"""

import sys, os, time
from math import sqrt, pi, sin, cos, acos, fabs
from string import strip
import gzip
import Numeric

RADIANS_TO_DEGREES = 180.0/pi
DEGREES_TO_RADIANS = pi/180.0
StringType = type('')
IntType = type(0)

class Linus:
    """Class to describe a (modified) LINUS object

    Instantiation

        An instance is created by calling the class constructor with
        one argument - name of the pdb file. 

    Attributes

        o *protein* - instance of LinusMol - the protein that is being
        simulated

        o *hbdpar* - hydrogen bonding energy parameters.

            i.   *use_hbond* - boolean - if true backbone - backbone
                 hydrogen bonds are enabled

            ii.  *use_sidechain_hbond* - boolean - if true sidechain
                 acceptor to backbone donor hydrogen bonds are
                 enabled

            iii. *hbond_distance* - float - optimal distance between
                 backbone donor and backbone acceptor

            iv.  *hbdond_torsion* - float - minium torsion angle between
                 the backbone acceptor and the backbone donor and it's two
                 antecedent atoms( 'C' of previous residue and 'CA' of the
                 residue

            v.   *hbond_probe* - float - distance over which the hydrogen
                 bond energy for backbone donor to backbone acceptor
                 bonds will be scaled to zero from the maximum value

            vi.  *hbond_score_short* - float - backbone donor to backbone
                 acceptor hydrogen bond energy for residues separations of
                 5 or less

            vii. *hbond_score_long* - float - backbone donor to backbone
                 acceptor hydrogen bond energy for residues separations of
                 6 or less

            viii. *sidechain_hbond_distance* - float - optimal distance between
                  backbone donor and sidechain acceptor

            ix.  *sidechain_hbond_torsion* - float - minium torsion angle
                 between the sidechain acceptor and the backbone donor and
                 it's two antecedent atoms ( 'C' of previous residue and 'CA' 
                 of the residue)

            x.   *sidechain_hbond_score* - float - energy of a sidechain
                 to backbone hydrogen bond

            xi.  *hbond_winmin* - integer - minimum separation between
                 two residues that can participate in backbone to backbone
                 hydrogen bonds

            xii. *hbond_winmax* - integer - maximum  separation between
                 two residues that can participate in backbone/sidechain
                 acceptor to backbone donor hydrogen bonds

    """

    def __init__(self, pdbfile):
        """Create a LINUS simulation object

            pdbfile - file containing coordinates of the
                      molecule to be simulated in PDB format


        """
        
        pdbfile = os.path.expanduser(pdbfile)
        self.protein = LinusProtein(pdbfile)
        self.hbdpar = default_hbdpar.copy()
        self.hblist = self.conlist = self.dscnstrlist = None

    def __getattr__(self, name):
        if self.__dict__.has_key(name):
            return self.__dict__[name]
        elif self.hbdpar.has_key(name):
            return self.hbdpar[name]
        else:
            raise AttributeError, 'Undefined attribute %s' % name

    def __setattr__(self, name, value):
        if name in ('pdbfile', 'protein', 'hbdpar', 'hblist'):
            self.__dict__[name] = value
        elif not hasattr(self, name):
            self.__dict__[name] = value
        else:
            raise AttributeError, 'Attribute %s cannot be set directly' % name

    def set_hbond_parameters(self, **kw):
        """Set the hydrogen bond parameters

        Arguments

            series of keyword=value arguments

        Errors

            will raise AttributeError if invalid parameter name is specified
            and TypeError if invalid value for the parameter is supplied

        """    
        
        for key, value in kw.items():
            if self.hbdpar.has_key(key):
                vtype = type(self.hbdpar[key])
                self.hbdpar[key] = converters[vtype](value)
            else:
                raise KeyError, 'Undefined hbond parameter %s' % key

class LinusProtein:
    """Class to describe a protein in LINUS.

    Instantiation

       Instantiation of this class requires one argument - the
       name of the file containing a description of the molecule
       in PDB format.

    Class Attributes

       o *filename* - string - name of file from which the molecule
       information was read

       o *atoms* - list - list of atoms in the molecule. each item
       in this list is an instance of *Atom3d*

       o *num_atoms* - integer - number of atoms in molecule

       o *num_residues* - integer - number of residues in molecule

       o *residue_names* - list - name of each residue in molecule

       o *res_pdb_number* - list - number of each residue in input PDB file 

       o *residue_first_atom_indices* - list - index of the the
       first atom of each residue (this is horrendously named)

       o *phi_atoms* - list - contains a reference to the 'C' atom
       of each residue

       o *psi_atoms* - list - contains a reference to the 'O' atom
       of each residue

       o *omega_atoms* - list - contains a reference to the 'CA' atom
       of each residue

       o *chi_torsions* - list of lists - references to the atoms
       corresponding to the various chi torsions of each residue

    Methods

        o *get_atom* - get a reference to an atom in a molecule given
        the serial number of the residue it belongs to and the atoms'
        name

        o *update_coordinates* - recompute cartesian coordinates for
        a molecule after changes to the internal coordinates

        o *radius_of_gyration* - compute the radius of gyration of
        the molecule

        o *reset_conformation* - restore the conformation of the
        molecule to that specified in the input PDB file from which
        the molecule was first created
    """

    def __init__(self, filename):
        """New instance of LinusMol

        Arguments

            o *filename* - string - name of the PDB format file containing
            a description of the molecule

        Result

            New instace of LinusMol
        """
        
        self.filename = os.path.expanduser(filename)
        protein_from_pdb(self)
        nr = self.num_residues
        fp = self._orig_fp_distance = []
        sp = self._orig_sp_angle = []
        tp = self._orig_tp_torsion = []
        for atom in self.atoms:
            fp.append(atom.fp_distance)
            sp.append(atom.sp_angle)
            tp.append(atom.tp_torsion)

    def get_atom(self, resid, atomname):
        """Get a reference to an atom in a molecule

        Arguments

            o *resid* - integer - serial number of residue in which
            atom is present

            o *atomname* - string - name of atom

        Result

            Returns a reference to the atom.  If an atom with required
            name is not found in the resiude, *ValueError* is raised.
        """    
        i1 = self.residue_first_atom_indices[resid]
        i2 = self.residue_first_atom_indices[resid+1]
        for i in range(i1, i2):
            if self.atoms[i].name == atomname:
                return self.atoms[i]
        else:
            raise ValueError, 'No such atom: %s' % atomname


class Atom3d:
    """Class to describe an atom in a molecule

    Attributes

        o *name* - string - name of an atom in PDB format

        o *radius* - float - hard sphere radius of an atom

        o *x*, *y*, *z* - float - cartesian coordinates of atom

        o *first_parent*, *second_parent*, *third_parent* - instances of
        Atom3d - Z-matrix description of an atom

        o *fp_distance*, *sp_angle*, *tp_torsion* - float - distance,
        angle in degrees and torsion angle in degrees made by atom with its
        first, second and third parents

        o *resnum* - integer - index of residue that atom is a part of

    """
    def __init__(self):
        self.name = ''
        self.radius = 0.0
        self.x = self.y = self.z = 0.0
        self.first_parent = self.second_parent = self.third_parent = None
        self.fp_distance = self.sp_angle = self.tp_torsion = 0.0
        self.fp_distance_tmp = self.sp_angle_tmp = self.tp_torsion_tmp = 0.0
        self.resnum = 0

    def __hash__(self):
        return hash(repr(self.resnum) + '_' + self.name)

    def coords(self):
        """Cartesian coordinates of an atom as a list"""
        return [self.x, self.y, self.z]

    def distance(self, other):
        """Euclidean distance between two atoms"""
        return sqrt((self.x - other.x)**2.0 + (self.y - other.y)**2.0 +
                    (self.z - other.z)**2.0)

    def close(self, other, dist=None):
        """check if two atoms are within a specified distance"""

        if dist is None:
            r1 = 0.90*self.radius
            r2 = 0.90*other.radius
            if r1 &lt; 1.0e-6 or r2 &lt; 1.0e-6: return 0
            dist = r1 + r2
        d2 = dist * dist
        dx = (self.x - other.x)**2.0
        if dx &gt;= d2: return 0
        dy = (self.y - other.y)**2.0
        if dy &gt;= d2: return 0
        dz = (self.z - other.z)**2.0
        if dz &gt;= d2: return 0
        return dx + dy + dz &lt; d2

    def int_to_cart(self):
        """generate cartesian coordinates for an atom from it's
        internal coordinates"""

        
        a2 = self.first_parent
        if a2 is None:
            self.x = self.y = self.z = 0.0
            return

        a3 = self.second_parent
        if a3 is None:
            self.x = self.fp_distance
            self.y = self.z = 0.0
            return


        ang = self.sp_angle * DEGREES_TO_RADIANS
        sina = sin(ang)
        cosa = -cos(ang)

        a4 = self.third_parent
        if a4 is None:
            self.z = 0.0
            self.x = a2.x + self.fp_distance*cosa
            self.y = self.fp_distance * sina
            return

        tor = self.tp_torsion * DEGREES_TO_RADIANS
        sint = sina * sin(tor)
        cost = sina * cos(tor)

        u1x = a3.x - a4.x
        u1y = a3.y - a4.y
        u1z = a3.z - a4.z
        d = 1.0 / sqrt(u1x**2.0 + u1y**2.0 + u1z**2.0)
        u1x = u1x * d
        u1y = u1y * d
        u1z = u1z * d

        u2x = a2.x - a3.x
        u2y = a2.y - a3.y
        u2z = a2.z - a3.z
        d = 1.0 / sqrt(u2x**2.0 + u2y**2.0 + u2z**2.0)
        u2x = u2x * d
        u2y = u2y * d
        u2z = u2z * d

        cosine = u1x*u2x + u1y*u2y + u1z*u2z

        if abs(cosine) &lt; 1.0:
            sine = 1.0/sqrt(1.0 - cosine*cosine)
        else:
            sine = 1.0/sqrt(cosine*cosine - 1.0)
            
        u3x = sine * (u1y*u2z - u1z*u2y)
        u3y = sine * (u1z*u2x - u1x*u2z)
        u3z = sine * (u1x*u2y - u1y*u2x)

        u4x = cost * (u3y*u2z - u3z*u2y)
        u4y = cost * (u3z*u2x - u3x*u2z)
        u4z = cost * (u3x*u2y - u3y*u2x)

        dist = self.fp_distance
        self.x = a2.x + dist*(u2x*cosa + u4x + u3x*sint)
        self.y = a2.y + dist*(u2y*cosa + u4y + u3y*sint)
        self.z = a2.z + dist*(u2z*cosa + u4z + u3z*sint)

    def commit_internal_coords(self):
        """Save internal coordinates following acceptance of a conformation
        change"""
        self.fp_distance_tmp = self.fp_distance
        self.sp_angle_tmp = self.sp_angle
        self.tp_torsion_tmp = self.tp_torsion


# RESIDUES - description of residues known to LINUS.  Each residue is 
# described as a tuple.  The items in tuple are themselves tuples describing
# each atom in the residue.  The fields of the tuple are:

#  0. name of atom
#  1, 2, 3 - first second and third parent of the atom.  The parents 
#            are described as a tuple of two values - offset and name.
#            The legal values for offset are 0 or -1. A value of 0
#            means that the parent for the atom is in the same residue
#            and -1 means previous residue
#  4 - Hard sphere radius of atom
#  5 - Radius of atom for use in contact energy calculations.  A value of
#      0 indicates the atom does not participate in contact energy
#      calculations
#  6 - is 1 if atom is a hydrogen bond donor and 0 otherwise
#  7 - is 1 if atom is a hydrogen bond acceptor and 0 otherwise
#  8 - number of bonds separating atom from its own backbone amide 
#      nitrogen
#  9 - number of bonds separating atom from the backbone amide nitrogen
#      of the succeeding residue

ACE = (
    (' CA ', (-1, ' C  '), (-1, ' CA '), (-1, ' O  '), 1.7, 0.0, 0, 0, 1, 2),
    (' C  ', ( 0, ' CA '), (-1, ' C  '), (-1, ' CA '), 1.50, 0.0, 0, 0, 2, 1),
    (' O  ', ( 0, ' C  '), ( 0, ' CA '), (-1, ' C  '), 1.35, 0.0, 0, 1, 3, 2)
)

NH2 = (
    (' N  ', (-1, ' C  '), (-1, ' CA '), (-1, ' O  '), 1.35, 0.0, 1, 0, 0, 3)
)

NME = (
    (' N  ', (-1, ' C  '), (-1, ' CA '), (-1, ' O  '), 1.35, 0.0, 1, 0, 0, 3),
    (' CA ', ( 0, ' N  '), (-1, ' C  '), (-1, ' CA '), 1.5, 0.0, 0, 0, 1, 2),
    (' H  ', ( 0, ' N  '), (-1, ' C  '), (-1, ' CA '), 0.90, 0.0, 0, 0, 1, 4),
)

ALA = (
    (' N  ', (-1, ' C  '), (-1, ' CA '), (-1, ' O  '), 1.35, 0.0, 1, 0, 0, 3),
    (' CA ', ( 0, ' N  '), (-1, ' C  '), (-1, ' CA '), 1.5, 0.0, 0, 0, 1, 2),
    (' C  ', ( 0, ' CA '), ( 0, ' N  '), (-1, ' C  '), 1.50, 0.0, 0, 0, 2, 1),
    (' O  ', ( 0, ' C  '), ( 0, ' CA '), ( 0, ' N  '), 1.34, 0.0, 0, 1, 3, 2),
    (' CB ', ( 0, ' CA '), ( 0, ' N  '), ( 0, ' C  '), 1.70, 2.0, 0, 0, 2, 3),
    (' H  ', ( 0, ' N  '), (-1, ' C  '), (-1, ' CA '), 0.90, 0.0, 0, 0, 1, 4),
    (' HA ', ( 0, ' CA '), ( 0, ' N  '), ( 0, ' C  '), 1.00, 0.0, 0, 0, 2, 3)
)

ARG = (
    (' N  ', (-1, ' C  '), (-1, ' CA '), (-1, ' O  '), 1.35, 0.0, 1, 0, 0, 3),
    (' CA ', ( 0, ' N  '), (-1, ' C  '), (-1, ' CA '), 1.5, 0.0, 0, 0, 1, 2),
    (' C  ', ( 0, ' CA '), ( 0, ' N  '), (-1, ' C  '), 1.50, 0.0, 0, 0, 2, 1),
    (' O  ', ( 0, ' C  '), ( 0, ' CA '), ( 0, ' N  '), 1.34, 0.0, 0, 1, 3, 2),
    (' CB ', ( 0, ' CA '), ( 0, ' N  '), ( 0, ' C  '), 1.7, 2.0, 0, 0, 2, 3),
    (' H  ', ( 0, ' N  '), (-1, ' C  '), (-1, ' CA '), 0.90, 0.0, 0, 0, 1, 4),
    (' HA ', ( 0, ' CA '), ( 0, ' N  '), ( 0, ' C  '), 1.00, 0.0, 0, 0, 2, 3),
    (' CG ', ( 0, ' CB '), ( 0, ' CA '), ( 0, ' N  '), 1.7, 2.0, 0, 0, 3, 4),
    (' CD ', ( 0, ' CG '), ( 0, ' CB '), ( 0, ' CA '), 1.7, 2.0, 0, 0, 4, 5),
    (' NE ', ( 0, ' CD '), ( 0, ' CG '), ( 0, ' CB '), 1.35, 0.0, 1, 0, 5, 6),
    (' CZ ', ( 0, ' NE '), ( 0, ' CD '), ( 0, ' CG '), 1.50, 0.0, 0, 0, 6, 7),
    (' NH1', ( 0, ' CZ '), ( 0, ' NE '), ( 0, ' CD '), 1.35, 0.0, 1, 0, 7, 8),
    (' NH2', ( 0, ' CZ '), ( 0, ' NE '), ( 0, ' NH1'), 1.35, 0.0, 1, 0, 7, 8)
)

ASN = (
    (' N  ', (-1, ' C  '), (-1, ' CA '), (-1, ' O  '), 1.35, 0.0, 1, 0, 0, 3),
    (' CA ', ( 0, ' N  '), (-1, ' C  '), (-1, ' CA '), 1.5, 0.0, 0, 0, 1, 2),
    (' C  ', ( 0, ' CA '), ( 0, ' N  '), (-1, ' C  '), 1.50, 0.0, 0, 0, 2, 1),
    (' O  ', ( 0, ' C  '), ( 0, ' CA '), ( 0, ' N  '), 1.34, 0.0, 0, 1, 3, 2),
    (' CB ', ( 0, ' CA '), ( 0, ' N  '), ( 0, ' C  '), 1.7, 2.0, 0, 0, 2, 3),
    (' H  ', ( 0, ' N  '), (-1, ' C  '), (-1, ' CA '), 0.90, 0.0, 0, 0, 1, 4),
    (' HA ', ( 0, ' CA '), ( 0, ' N  '), ( 0, ' C  '), 1.00, 0.0, 0, 0, 2, 3),
    (' CG ', ( 0, ' CB '), ( 0, ' CA '), ( 0, ' N  '), 1.7, 0.0, 0, 0, 3, 4),
    (' OD1', ( 0, ' CG '), ( 0, ' CB '), ( 0, ' CA '), 1.35, 0.0, 1, 1, 4, 5),
    (' ND2', ( 0, ' CG '), ( 0, ' CB '), ( 0, ' OD1'), 1.35, 0.0, 1, 1, 4, 5)
)

ASP = (
    (' N  ', (-1, ' C  '), (-1, ' CA '), (-1, ' O  '), 1.35, 0.0, 1, 0, 0, 3),
    (' CA ', ( 0, ' N  '), (-1, ' C  '), (-1, ' CA '), 1.5, 0.0, 0, 0, 1, 2),
    (' C  ', ( 0, ' CA '), ( 0, ' N  '), (-1, ' C  '), 1.50, 0.0, 0, 0, 2, 1),
    (' O  ', ( 0, ' C  '), ( 0, ' CA '), ( 0, ' N  '), 1.34, 0.0, 0, 1, 3, 2),
    (' CB ', ( 0, ' CA '), ( 0, ' N  '), ( 0, ' C  '), 1.7, 2.0, 0, 0, 2, 3),
    (' H  ', ( 0, ' N  '), (-1, ' C  '), (-1, ' CA '), 0.90, 0.0, 0, 0, 1, 4),
    (' HA ', ( 0, ' CA '), ( 0, ' N  '), ( 0, ' C  '), 1.00, 0.0, 0, 0, 2, 3),
    (' CG ', ( 0, ' CB '), ( 0, ' CA '), ( 0, ' N  '), 1.7, 0.0, 0, 0, 3, 4),
    (' OD1', ( 0, ' CG '), ( 0, ' CB '), ( 0, ' CA '), 1.35, 0.0, 1, 1, 4, 5),
    (' OD2', ( 0, ' CG '), ( 0, ' CB '), ( 0, ' OD1'), 1.35, 0.0, 1, 1, 4, 5)
)

CYS = (
    (' N  ', (-1, ' C  '), (-1, ' CA '), (-1, ' O  '), 1.35, 0.0, 1, 0, 0, 3),
    (' CA ', ( 0, ' N  '), (-1, ' C  '), (-1, ' CA '), 1.5, 0.0, 0, 0, 1, 2),
    (' C  ', ( 0, ' CA '), ( 0, ' N  '), (-1, ' C  '), 1.50, 0.0, 0, 0, 2, 1),
    (' O  ', ( 0, ' C  '), ( 0, ' CA '), ( 0, ' N  '), 1.34, 0.0, 0, 1, 3, 2),
    (' CB ', ( 0, ' CA '), ( 0, ' N  '), ( 0, ' C  '), 1.7, 2.0, 0, 0, 2, 3),
    (' H  ', ( 0, ' N  '), (-1, ' C  '), (-1, ' CA '), 0.90, 0.0, 0, 0, 1, 4),
    (' HA ', ( 0, ' CA '), ( 0, ' N  '), ( 0, ' C  '), 1.00, 0.0, 0, 0, 2, 3),
    (' SG ', ( 0, ' CB '), ( 0, ' CA '), ( 0, ' N  '), 1.7, 0.0, 1, 1, 3, 4),
)

GLN = (
    (' N  ', (-1, ' C  '), (-1, ' CA '), (-1, ' O  '), 1.35, 0.0, 1, 0, 0, 3),
    (' CA ', ( 0, ' N  '), (-1, ' C  '), (-1, ' CA '), 1.5, 0.0, 0, 0, 1, 2),
    (' C  ', ( 0, ' CA '), ( 0, ' N  '), (-1, ' C  '), 1.50, 0.0, 0, 0, 2, 1),
    (' O  ', ( 0, ' C  '), ( 0, ' CA '), ( 0, ' N  '), 1.34, 0.0, 0, 1, 3, 2),
    (' CB ', ( 0, ' CA '), ( 0, ' N  '), ( 0, ' C  '), 1.7, 2.0, 0, 0, 2, 3),
    (' H  ', ( 0, ' N  '), (-1, ' C  '), (-1, ' CA '), 0.90, 0.0, 0, 0, 1, 4),
    (' HA ', ( 0, ' CA '), ( 0, ' N  '), ( 0, ' C  '), 1.00, 0.0, 0, 0, 2, 3),
    (' CG ', ( 0, ' CB '), ( 0, ' CA '), ( 0, ' N  '), 1.7, 2.0, 0, 0, 3, 4),
    (' CD ', ( 0, ' CG '), ( 0, ' CB '), ( 0, ' CA '), 1.50, 0.0, 0, 0, 4, 5),
    (' OE1', ( 0, ' CD '), ( 0, ' CG '), ( 0, ' CB '), 1.35, 0.0, 1, 1, 5, 6),
    (' NE2', ( 0, ' CD '), ( 0, ' CG '), ( 0, ' OE1'), 1.35, 0.0, 1, 1, 5, 6)
)

GLU = (
    (' N  ', (-1, ' C  '), (-1, ' CA '), (-1, ' O  '), 1.35, 0.0, 1, 0, 0, 3),
    (' CA ', ( 0, ' N  '), (-1, ' C  '), (-1, ' CA '), 1.5, 0.0, 0, 0, 1, 2),
    (' C  ', ( 0, ' CA '), ( 0, ' N  '), (-1, ' C  '), 1.50, 0.0, 0, 0, 2, 1),
    (' O  ', ( 0, ' C  '), ( 0, ' CA '), ( 0, ' N  '), 1.34, 0.0, 0, 1, 3, 2),
    (' CB ', ( 0, ' CA '), ( 0, ' N  '), ( 0, ' C  '), 1.7, 2.0, 0, 0, 2, 3),
    (' H  ', ( 0, ' N  '), (-1, ' C  '), (-1, ' CA '), 0.90, 0.0, 0, 0, 1, 4),
    (' HA ', ( 0, ' CA '), ( 0, ' N  '), ( 0, ' C  '), 1.00, 0.0, 0, 0, 2, 3),
    (' CG ', ( 0, ' CB '), ( 0, ' CA '), ( 0, ' N  '), 1.7, 2.0, 0, 0, 3, 4),
    (' CD ', ( 0, ' CG '), ( 0, ' CB '), ( 0, ' CA '), 1.50, 0.0, 0, 0, 4, 5),
    (' OE1', ( 0, ' CD '), ( 0, ' CG '), ( 0, ' CB '), 1.35, 0.0, 1, 1, 5, 6),
    (' OE2', ( 0, ' CD '), ( 0, ' CG '), ( 0, ' OE1'), 1.35, 0.0, 1, 1, 5, 6)
)

GLY = (
    (' N  ', (-1, ' C  '), (-1, ' CA '), (-1, ' O  '), 1.35, 0.0, 1, 0, 0, 3),
    (' CA ', ( 0, ' N  '), (-1, ' C  '), (-1, ' CA '), 1.5, 0.0, 0, 0, 1, 2),
    (' C  ', ( 0, ' CA '), ( 0, ' N  '), (-1, ' C  '), 1.50, 0.0, 0, 0, 2, 1),
    (' O  ', ( 0, ' C  '), ( 0, ' CA '), ( 0, ' N  '), 1.34, 0.0, 0, 1, 3, 2),
    (' H  ', ( 0, ' N  '), (-1, ' C  '), (-1, ' CA '), 0.90, 0.0, 0, 0, 1, 4),
    (' HA ', ( 0, ' CA '), ( 0, ' N  '), ( 0, ' C  '), 1.00, 0.0, 0, 0, 2, 3),
    (' HB ', ( 0, ' CA '), ( 0, ' N  '), ( 0, ' C  '), 1.00, 0.0, 0, 0, 2, 3),
)

HIS = (
    (' N  ', (-1, ' C  '), (-1, ' CA '), (-1, ' O  '), 1.35, 0.0, 1, 0, 0, 3),
    (' CA ', ( 0, ' N  '), (-1, ' C  '), (-1, ' CA '), 1.5, 0.0, 0, 0, 1, 2),
    (' C  ', ( 0, ' CA '), ( 0, ' N  '), (-1, ' C  '), 1.50, 0.0, 0, 0, 2, 1),
    (' O  ', ( 0, ' C  '), ( 0, ' CA '), ( 0, ' N  '), 1.34, 0.0, 0, 1, 3, 2),
    (' CB ', ( 0, ' CA '), ( 0, ' N  '), ( 0, ' C  '), 1.7, 2.0, 0, 0, 2, 3),
    (' H  ', ( 0, ' N  '), (-1, ' C  '), (-1, ' CA '), 0.90, 0.0, 0, 0, 1, 4),
    (' HA ', ( 0, ' CA '), ( 0, ' N  '), ( 0, ' C  '), 1.00, 0.0, 0, 0, 2, 3),
    (' CG ', ( 0, ' CB '), ( 0, ' CA '), ( 0, ' N  '), 1.7, 0.0, 0, 0, 3, 4),
    (' ND1', ( 0, ' CG '), ( 0, ' CB '), ( 0, ' CA '), 1.35, 0.0, 1, 1, 4, 5),
    (' CD2', ( 0, ' CG '), ( 0, ' CB '), ( 0, ' ND1'), 1.50, 0.0, 0, 0, 4, 5),
    (' CE1', ( 0, ' ND1'), ( 0, ' CG '), ( 0, ' CB '), 1.50, 0.0, 0, 0, 5, 6),
    (' NE2', ( 0, ' CD2'), ( 0, ' CG '), ( 0, ' CE1'), 1.35, 0.0, 1, 1, 5, 6),
)

ILE = (
    (' N  ', (-1, ' C  '), (-1, ' CA '), (-1, ' O  '), 1.35, 0.0, 1, 0, 0, 3),
    (' CA ', ( 0, ' N  '), (-1, ' C  '), (-1, ' CA '), 1.5, 0.0, 0, 0, 1, 2),
    (' C  ', ( 0, ' CA '), ( 0, ' N  '), (-1, ' C  '), 1.50, 0.0, 0, 0, 2, 1),
    (' O  ', ( 0, ' C  '), ( 0, ' CA '), ( 0, ' N  '), 1.34, 0.0, 0, 1, 3, 2),
    (' CB ', ( 0, ' CA '), ( 0, ' N  '), ( 0, ' C  '), 1.7, 2.0, 0, 0, 2, 3),
    (' H  ', ( 0, ' N  '), (-1, ' C  '), (-1, ' CA '), 0.90, 0.0, 0, 0, 1, 4),
    (' HA ', ( 0, ' CA '), ( 0, ' N  '), ( 0, ' C  '), 1.00, 0.0, 0, 0, 2, 3),
    (' CG1', ( 0, ' CB '), ( 0, ' CA '), ( 0, ' N  '), 1.7, 2.0, 0, 0, 3, 4),
    (' CG2', ( 0, ' CB '), ( 0, ' CA '), ( 0, ' CG1'), 1.7, 2.0, 0, 0, 3, 4),
    (' CD1', ( 0, ' CG1'), ( 0, ' CB '), ( 0, ' CA '), 1.7, 2.0, 0, 0, 4, 5),
)

LEU = (
    (' N  ', (-1, ' C  '), (-1, ' CA '), (-1, ' O  '), 1.35, 0.0, 1, 0, 0, 3),
    (' CA ', ( 0, ' N  '), (-1, ' C  '), (-1, ' CA '), 1.5, 0.0, 0, 0, 1, 2),
    (' C  ', ( 0, ' CA '), ( 0, ' N  '), (-1, ' C  '), 1.50, 0.0, 0, 0, 2, 1),
    (' O  ', ( 0, ' C  '), ( 0, ' CA '), ( 0, ' N  '), 1.34, 0.0, 0, 1, 3, 2),
    (' CB ', ( 0, ' CA '), ( 0, ' N  '), ( 0, ' C  '), 1.7, 2.0, 0, 0, 2, 3),
    (' H  ', ( 0, ' N  '), (-1, ' C  '), (-1, ' CA '), 0.90, 0.0, 0, 0, 1, 4),
    (' HA ', ( 0, ' CA '), ( 0, ' N  '), ( 0, ' C  '), 1.00, 0.0, 0, 0, 2, 3),
    (' CG ', ( 0, ' CB '), ( 0, ' CA '), ( 0, ' N  '), 1.7, 2.0, 0, 0, 3, 4),
    (' CD1', ( 0, ' CG '), ( 0, ' CB '), ( 0, ' CA '), 1.7, 2.0, 0, 0, 4, 5),
    (' CD2', ( 0, ' CG '), ( 0, ' CB '), ( 0, ' CD1'), 1.7, 2.0, 0, 0, 4, 5),
)

LYS = (
    (' N  ', (-1, ' C  '), (-1, ' CA '), (-1, ' O  '), 1.35, 0.0, 1, 0, 0, 3),
    (' CA ', ( 0, ' N  '), (-1, ' C  '), (-1, ' CA '), 1.5, 0.0, 0, 0, 1, 2),
    (' C  ', ( 0, ' CA '), ( 0, ' N  '), (-1, ' C  '), 1.50, 0.0, 0, 0, 2, 1),
    (' O  ', ( 0, ' C  '), ( 0, ' CA '), ( 0, ' N  '), 1.34, 0.0, 0, 1, 3, 2),
    (' CB ', ( 0, ' CA '), ( 0, ' N  '), ( 0, ' C  '), 1.7, 2.0, 0, 0, 2, 3),
    (' H  ', ( 0, ' N  '), (-1, ' C  '), (-1, ' CA '), 0.90, 0.0, 0, 0, 1, 4),
    (' HA ', ( 0, ' CA '), ( 0, ' N  '), ( 0, ' C  '), 1.00, 0.0, 0, 0, 2, 3),
    (' CG ', ( 0, ' CB '), ( 0, ' CA '), ( 0, ' N  '), 1.7, 2.0, 0, 0, 3, 4),
    (' CD ', ( 0, ' CG '), ( 0, ' CB '), ( 0, ' CA '), 1.7, 2.0, 0, 0, 4, 5),
    (' CE ', ( 0, ' CD '), ( 0, ' CG '), ( 0, ' CB '), 1.7, 0.0, 0, 0, 5, 6),
    (' NZ ', ( 0, ' CE '), ( 0, ' CD '), ( 0, ' CG '), 1.35, 0.0, 1, 0, 6, 7),
)

MET = (
    (' N  ', (-1, ' C  '), (-1, ' CA '), (-1, ' O  '), 1.35, 0.0, 1, 0, 0, 3),
    (' CA ', ( 0, ' N  '), (-1, ' C  '), (-1, ' CA '), 1.5, 0.0, 0, 0, 1, 2),
    (' C  ', ( 0, ' CA '), ( 0, ' N  '), (-1, ' C  '), 1.50, 0.0, 0, 0, 2, 1),
    (' O  ', ( 0, ' C  '), ( 0, ' CA '), ( 0, ' N  '), 1.34, 0.0, 0, 1, 3, 2),
    (' CB ', ( 0, ' CA '), ( 0, ' N  '), ( 0, ' C  '), 1.7, 2.0, 0, 0, 2, 3),
    (' H  ', ( 0, ' N  '), (-1, ' C  '), (-1, ' CA '), 0.90, 0.0, 0, 0, 1, 4),
    (' HA ', ( 0, ' CA '), ( 0, ' N  '), ( 0, ' C  '), 1.00, 0.0, 0, 0, 2, 3),
    (' CG ', ( 0, ' CB '), ( 0, ' CA '), ( 0, ' N  '), 1.7, 2.0, 0, 0, 3, 4),
    (' SD ', ( 0, ' CG '), ( 0, ' CB '), ( 0, ' CA '), 1.80, 0.0, 0, 1, 4, 5),
    (' CE ', ( 0, ' SD '), ( 0, ' CG '), ( 0, ' CB '), 1.7, 2.0, 0, 0, 5, 6),
)

PHE = (
    (' N  ', (-1, ' C  '), (-1, ' CA '), (-1, ' O  '), 1.35, 0.0, 1, 0, 0, 3),
    (' CA ', ( 0, ' N  '), (-1, ' C  '), (-1, ' CA '), 1.5, 0.0, 0, 0, 1, 2),
    (' C  ', ( 0, ' CA '), ( 0, ' N  '), (-1, ' C  '), 1.50, 0.0, 0, 0, 2, 1),
    (' O  ', ( 0, ' C  '), ( 0, ' CA '), ( 0, ' N  '), 1.34, 0.0, 0, 1, 3, 2),
    (' CB ', ( 0, ' CA '), ( 0, ' N  '), ( 0, ' C  '), 1.7, 2.0, 0, 0, 2, 3),
    (' H  ', ( 0, ' N  '), (-1, ' C  '), (-1, ' CA '), 0.90, 0.0, 0, 0, 1, 4),
    (' HA ', ( 0, ' CA '), ( 0, ' N  '), ( 0, ' C  '), 1.00, 0.0, 0, 0, 2, 3),
    (' CG ', ( 0, ' CB '), ( 0, ' CA '), ( 0, ' N  '), 1.50, 2.0, 0, 0, 3, 4),
    (' CD1', ( 0, ' CG '), ( 0, ' CB '), ( 0, ' CA '), 1.50, 2.0, 0, 0, 4, 5),
    (' CD2', ( 0, ' CG '), ( 0, ' CB '), ( 0, ' CD1'), 1.50, 2.0, 0, 0, 4, 5),
    (' CE1', ( 0, ' CD1'), ( 0, ' CG '), ( 0, ' CB '), 1.50, 2.0, 0, 0, 5, 6),
    (' CE2', ( 0, ' CD2'), ( 0, ' CG '), ( 0, ' CE1'), 1.50, 2.0, 0, 0, 5, 6),
    (' CZ ', ( 0, ' CE1'), ( 0, ' CD1'), ( 0, ' CG '), 1.50, 2.0, 0, 0, 6, 7),
)

PRO = (
    (' N  ', (-1, ' C  '), (-1, ' CA '), (-1, ' O  '), 1.35, 0.0, 0, 0, 0, 3),
    (' CA ', ( 0, ' N  '), (-1, ' C  '), (-1, ' CA '), 1.5, 0.0, 0, 0, 1, 2),
    (' C  ', ( 0, ' CA '), ( 0, ' N  '), (-1, ' C  '), 1.50, 0.0, 0, 0, 2, 1),
    (' O  ', ( 0, ' C  '), ( 0, ' CA '), ( 0, ' N  '), 1.34, 0.0, 0, 1, 3, 2),
    (' CB ', ( 0, ' CA '), ( 0, ' N  '), ( 0, ' C  '), 1.7, 2.0, 0, 0, 2, 3),
    (' HA ', ( 0, ' CA '), ( 0, ' N  '), ( 0, ' C  '), 1.00, 0.0, 0, 0, 2, 3),
    (' CG ', ( 0, ' CB '), ( 0, ' CA '), ( 0, ' N  '), 1.7, 2.0, 0, 0, 3, 2),
    (' CD ', ( 0, ' CG '), ( 0, ' CB '), ( 0, ' CA '), 1.45, 2.0, 0, 0, 1, 1)
)

SER = (
    (' N  ', (-1, ' C  '), (-1, ' CA '), (-1, ' O  '), 1.35, 0.0, 1, 0, 0, 3),
    (' CA ', ( 0, ' N  '), (-1, ' C  '), (-1, ' CA '), 1.5, 0.0, 0, 0, 1, 2),
    (' C  ', ( 0, ' CA '), ( 0, ' N  '), (-1, ' C  '), 1.50, 0.0, 0, 0, 2, 1),
    (' O  ', ( 0, ' C  '), ( 0, ' CA '), ( 0, ' N  '), 1.34, 0.0, 0, 1, 3, 2),
    (' CB ', ( 0, ' CA '), ( 0, ' N  '), ( 0, ' C  '), 1.7, 2.0, 0, 0, 2, 3),
    (' H  ', ( 0, ' N  '), (-1, ' C  '), (-1, ' CA '), 0.90, 0.0, 0, 0, 1, 4),
    (' HA ', ( 0, ' CA '), ( 0, ' N  '), ( 0, ' C  '), 1.00, 0.0, 0, 0, 2, 3),
    (' OG ', ( 0, ' CB '), ( 0, ' CA '), ( 0, ' N  '), 1.35, 0.0, 1, 1, 3, 4)
)

THR = (
    (' N  ', (-1, ' C  '), (-1, ' CA '), (-1, ' O  '), 1.35, 0.0, 1, 0, 0, 3),
    (' CA ', ( 0, ' N  '), (-1, ' C  '), (-1, ' CA '), 1.5, 0.0, 0, 0, 1, 2),
    (' C  ', ( 0, ' CA '), ( 0, ' N  '), (-1, ' C  '), 1.50, 0.0, 0, 0, 2, 1),
    (' O  ', ( 0, ' C  '), ( 0, ' CA '), ( 0, ' N  '), 1.34, 0.0, 0, 1, 3, 2),
    (' CB ', ( 0, ' CA '), ( 0, ' N  '), ( 0, ' C  '), 1.7, 2.0, 0, 0, 2, 3),
    (' H  ', ( 0, ' N  '), (-1, ' C  '), (-1, ' CA '), 0.90, 0.0, 0, 0, 1, 4),
    (' HA ', ( 0, ' CA '), ( 0, ' N  '), ( 0, ' C  '), 1.00, 0.0, 0, 0, 2, 3),
    (' OG1', ( 0, ' CB '), ( 0, ' CA '), ( 0, ' N  '), 1.35, 0.0, 1, 1, 3, 4),
    (' CG2', ( 0, ' CB '), ( 0, ' CA '), ( 0, ' OG1'), 1.7, 2.0, 0, 0, 3, 4)
)

TRP = (
    (' N  ', (-1, ' C  '), (-1, ' CA '), (-1, ' O  '), 1.35, 0.0, 1, 0, 0, 3),
    (' CA ', ( 0, ' N  '), (-1, ' C  '), (-1, ' CA '), 1.5, 0.0, 0, 0, 1, 2),
    (' C  ', ( 0, ' CA '), ( 0, ' N  '), (-1, ' C  '), 1.50, 0.0, 0, 0, 2, 1),
    (' O  ', ( 0, ' C  '), ( 0, ' CA '), ( 0, ' N  '), 1.34, 0.0, 0, 1, 3, 2),
    (' CB ', ( 0, ' CA '), ( 0, ' N  '), ( 0, ' C  '), 1.7, 2.0, 0, 0, 2, 3),
    (' H  ', ( 0, ' N  '), (-1, ' C  '), (-1, ' CA '), 0.90, 0.0, 0, 0, 1, 4),
    (' HA ', ( 0, ' CA '), ( 0, ' N  '), ( 0, ' C  '), 1.00, 0.0, 0, 0, 2, 3),
    (' CG ', ( 0, ' CB '), ( 0, ' CA '), ( 0, ' N  '), 1.50, 2.0, 0, 0, 3, 4),
    (' CD1', ( 0, ' CG '), ( 0, ' CB '), ( 0, ' CA '), 1.50, 2.0, 0, 0, 4, 5),
    (' CD2', ( 0, ' CG '), ( 0, ' CB '), ( 0, ' CD1'), 1.50, 2.0, 0, 0, 4, 5),
    (' NE1', ( 0, ' CD1'), ( 0, ' CG '), ( 0, ' CB '), 1.35, 0.0, 1, 0, 5, 6),
    (' CE2', ( 0, ' NE1'), ( 0, ' CD1'), ( 0, ' CG '), 1.50, 2.0, 0, 0, 5, 6),
    (' CE3', ( 0, ' CD2'), ( 0, ' CG '), ( 0, ' CE2'), 1.50, 2.0, 0, 0, 5, 6),
    (' CZ2', ( 0, ' CE2'), ( 0, ' NE1'), ( 0, ' CD1'), 1.50, 2.0, 0, 0, 6, 7),
    (' CZ3', ( 0, ' CE3'), ( 0, ' CD2'), ( 0, ' CG '), 1.50, 2.0, 0, 0, 6, 7),
    (' CH2', ( 0, ' CZ3'), ( 0, ' CE3'), ( 0, ' CD2'), 1.50, 2.0, 0, 0, 7, 8),
)

TYR = (
    (' N  ', (-1, ' C  '), (-1, ' CA '), (-1, ' O  '), 1.35, 0.0, 1, 0, 0, 3),
    (' CA ', ( 0, ' N  '), (-1, ' C  '), (-1, ' CA '), 1.5, 0.0, 0, 0, 1, 2),
    (' C  ', ( 0, ' CA '), ( 0, ' N  '), (-1, ' C  '), 1.50, 0.0, 0, 0, 2, 1),
    (' O  ', ( 0, ' C  '), ( 0, ' CA '), ( 0, ' N  '), 1.34, 0.0, 0, 1, 3, 2),
    (' CB ', ( 0, ' CA '), ( 0, ' N  '), ( 0, ' C  '), 1.7, 2.0, 0, 0, 2, 3),
    (' H  ', ( 0, ' N  '), (-1, ' C  '), (-1, ' CA '), 0.90, 0.0, 0, 0, 1, 4),
    (' HA ', ( 0, ' CA '), ( 0, ' N  '), ( 0, ' C  '), 1.00, 0.0, 0, 0, 2, 3),
    (' CG ', ( 0, ' CB '), ( 0, ' CA '), ( 0, ' N  '), 1.50, 2.0, 0, 0, 3, 4),
    (' CD1', ( 0, ' CG '), ( 0, ' CB '), ( 0, ' CA '), 1.50, 2.0, 0, 0, 4, 5),
    (' CD2', ( 0, ' CG '), ( 0, ' CB '), ( 0, ' CD1'), 1.50, 2.0, 0, 0, 4, 5),
    (' CE1', ( 0, ' CD1'), ( 0, ' CG '), ( 0, ' CB '), 1.50, 2.0, 0, 0, 5, 6),
    (' CE2', ( 0, ' CD2'), ( 0, ' CG '), ( 0, ' CE1'), 1.50, 2.0, 0, 0, 5, 6),
    (' CZ ', ( 0, ' CE1'), ( 0, ' CD1'), ( 0, ' CG '), 1.50, 2.0, 0, 0, 6, 7),
    (' OH ', ( 0, ' CZ '), ( 0, ' CE1'), ( 0, ' CE2'), 1.35, 0.0, 1, 1, 7, 8),
)

VAL = (
    (' N  ', (-1, ' C  '), (-1, ' CA '), (-1, ' O  '), 1.35, 0.0, 1, 0, 0, 3),
    (' CA ', ( 0, ' N  '), (-1, ' C  '), (-1, ' CA '), 1.5, 0.0, 0, 0, 1, 2),
    (' C  ', ( 0, ' CA '), ( 0, ' N  '), (-1, ' C  '), 1.50, 0.0, 0, 0, 2, 1),
    (' O  ', ( 0, ' C  '), ( 0, ' CA '), ( 0, ' N  '), 1.34, 0.0, 0, 1, 3, 2),
    (' CB ', ( 0, ' CA '), ( 0, ' N  '), ( 0, ' C  '), 1.7, 2.0, 0, 0, 2, 3),
    (' H  ', ( 0, ' N  '), (-1, ' C  '), (-1, ' CA '), 0.90, 0.0, 0, 0, 1, 4),
    (' HA ', ( 0, ' CA '), ( 0, ' N  '), ( 0, ' C  '), 1.00, 0.0, 0, 0, 2, 3),
    (' CG1', ( 0, ' CB '), ( 0, ' CA '), ( 0, ' N  '), 1.7, 2.0, 0, 0, 3, 4),
    (' CG2', ( 0, ' CB '), ( 0, ' CA '), ( 0, ' CG1'), 1.7, 2.0, 0, 0, 3, 4),
)
RESIDUES = {
    'ACE': ACE,
    'NH2': NH2,
    'NME': NME,
    'ALA': ALA,
    'ARG': ARG,
    'ASN': ASN,
    'ASP': ASP,
    'CYS': CYS,
    'GLN': GLN,
    'GLU': GLU,
    'GLY': GLY,
    'HIS': HIS,
    'ILE': ILE,
    'LEU': LEU,
    'LYS': LYS,
    'MET': MET,
    'PHE': PHE,
    'PRO': PRO,
    'SER': SER,
    'THR': THR,
    'TRP': TRP,
    'TYR': TYR,
    'VAL': VAL,
}

PhiAtoms = {
    'ANY': ' C  '
}

PsiAtoms = {
    'ANY': ' O  '
}

OmeAtoms = {
    'ANY': ' CA '
}

ChiAtoms = {
    'ARG': (' CG ', ' CD ', ' NE ', ' CZ '), 
    'ASN': (' CG ', ' OD1'),
    'ASP': (' CG ', ' OD1'),
    'CYS': (' SG ',),
    'GLN': (' CG ', ' CD ', ' OE1'),
    'GLU': (' CG ', ' CD ', ' OE1'),
    'HIS': (' CG ', ' ND1'),
    'ILE': (' CG1', ' CD1'),
    'LEU': (' CG ', ' CD1'),
    'LYS': (' CG ', ' CD ', ' CE ', ' NZ '),
    'MET': (' CG ', ' SD ', ' CE '),
    'PHE': (' CG ', ' CD1'),
    'SER': (' OG ',),
    'THR': (' OG1',),
    'TRP': (' CG ', ' CD1'),
    'TYR': (' CG ', ' CD1'),
    'VAL': (' CG1',),
}


HYDROPHOBIC = {
    ' C  ' : 1,
    ' CA ' : 1,
    ' CB ' : 1,
    ' CD ' : 1,
    ' CD1' : 1,
    ' CD2' : 1,
    ' CE ' : 1,
    ' CE1' : 1,
    ' CE2' : 1,
    ' CE3' : 1,
    ' CG ' : 1,
    ' CG1' : 1,
    ' CG2' : 1,
    ' CH2' : 1,
    ' CH3' : 1,
    ' CZ ' : 1,
    ' CZ2' : 1,
    ' CZ3' : 1,
    ' N  ' : 0,
    ' ND1' : 0,
    ' ND2' : 0,
    ' NE ' : 0,
    ' NE1' : 0,
    ' NE2' : 0,
    ' NH1' : 0,
    ' NH2' : 0,
    ' NZ ' : 0,
    ' O  ' : 0,
    ' OD1' : 0,
    ' OD2' : 0,
    ' OE1' : 0,
    ' OE2' : 0,
    ' OG ' : 0,
    ' OG1' : 0,
    ' OH ' : 0,
    ' OW1' : 0,
    ' OW2' : 0,
    ' OXT' : 0,
    ' SD ' : 0,
    ' SG ' : 0,
    ' H  ' : 1,
    ' HA ' : 1,
    ' HB ' : 1,
    ' HD1' : 1,
    ' HD2' : 1,
    ' HE ' : 1,
    ' HE1' : 1,
    ' HE2' : 1,
    ' HE3' : 1,
    ' HG ' : 1,
    ' HG1' : 1,
    ' HH ' : 1,
    ' HH2' : 1,
    ' HZ ' : 1,
    ' HZ2' : 1,
    ' HZ3' : 1,
    '1H  ' : 1,
    '1HA ' : 1,
    '1HB ' : 1,
    '1HD ' : 1,
    '1HD1' : 1,
    '1HD2' : 1,
    '1HE ' : 1,
    '1HE2' : 1,
    '1HG ' : 1,
    '1HG1' : 1,
    '1HG2' : 1,
    '1HH1' : 1,
    '1HH2' : 1,
    '1HH3' : 1,
    '1HW1' : 1,
    '1HW2' : 1,
    '1HZ ' : 1,
    '2H  ' : 1,
    '2HA ' : 1,
    '2HB ' : 1,
    '2HD ' : 1,
    '2HD1' : 1,
    '2HD2' : 1,
    '2HE ' : 1,
    '2HE2' : 1,
    '2HG ' : 1,
    '2HG1' : 1,
    '2HG2' : 1,
    '2HH1' : 1,
    '2HH2' : 1,
    '2HH3' : 1,
    '2HW1' : 1,
    '2HW2' : 1,
    '2HZ ' : 1,
    '3H  ' : 1,
    '3HB ' : 1,
    '3HD1' : 1,
    '3HD2' : 1,
    '3HE ' : 1,
    '3HG1' : 1,
    '3HG2' : 1,
    '3HH3' : 1,
    '3HZ ' : 1
    }

HYDROGEN = {
    ' H  ' : 1,
    ' HA ' : 1,
    ' HB ' : 1,
    ' HD1' : 1,
    ' HD2' : 1,
    ' HE ' : 1,
    ' HE1' : 1,
    ' HE2' : 1,
    ' HE3' : 1,
    ' HG ' : 1,
    ' HG1' : 1,
    ' HH ' : 1,
    ' HH2' : 1,
    ' HZ ' : 1,
    ' HZ2' : 1,
    ' HZ3' : 1,
    '1H  ' : 1,
    '1HA ' : 1,
    '1HB ' : 1,
    '1HD ' : 1,
    '1HD1' : 1,
    '1HD2' : 1,
    '1HE ' : 1,
    '1HE2' : 1,
    '1HG ' : 1,
    '1HG1' : 1,
    '1HG2' : 1,
    '1HH1' : 1,
    '1HH2' : 1,
    '1HH3' : 1,
    '1HW1' : 1,
    '1HW2' : 1,
    '1HZ ' : 1,
    '2H  ' : 1,
    '2HA ' : 1,
    '2HB ' : 1,
    '2HD ' : 1,
    '2HD1' : 1,
    '2HD2' : 1,
    '2HE ' : 1,
    '2HE2' : 1,
    '2HG ' : 1,
    '2HG1' : 1,
    '2HG2' : 1,
    '2HH1' : 1,
    '2HH2' : 1,
    '2HH3' : 1,
    '2HW1' : 1,
    '2HW2' : 1,
    '2HZ ' : 1,
    '3H  ' : 1,
    '3HB ' : 1,
    '3HD1' : 1,
    '3HD2' : 1,
    '3HE ' : 1,
    '3HG1' : 1,
    '3HG2' : 1,
    '3HH3' : 1,
    '3HZ ' : 1
    }


def get_residue_atoms(res_name):
    res = RESIDUES[res_name]
    atoms = []
    for arec in res:
        a = Atom3d()
        a.name = arec[0]
        a.radius = arec[4]
        a.bsep0 = arec[-2]
        a.bsep1 = arec[-1]
        atoms.append(a)
    return atoms

def get_atom_with_name(atoms, start, end, name):
    for i in range(start, end):
        if atoms[i].name == name: return atoms[i]
    return None

def protein_from_pdb(p):
    """Fill in the attributes of a LinusMol

    Arguments

        o *p* - instance of LinusMol

    Result

        None
    """    
    if p.filename[-3:] == '.gz':
        data = gzip.GzipFile(p.filename, 'rb').readlines()
    else:
        data = open(p.filename, 'rb').readlines()
    na = nr = -1
    a = p.atoms = []
    r = p.residue_names = []
    f = p.residue_first_atom_indices = []
    respdbnum = p.res_pdb_number = []

    nfstr = 'Unable to find atom %s in residue %i, deleting from molecule...'
    current_res = None
    aname_seen = []
    
    for line in data:
        if line[:4] == 'ATOM':
            na = na + 1
            rnam = line[17:20]
            rnum = line[22:27]
            if rnum &lt;&gt; current_res:
                alim = 0
                if len(f): alim = f[-1]
                for aptr in xrange(len(a)-1, alim, -1):
                    if a[aptr].name not in aname_seen:
#                           print (nfstr % (a[aptr].name, int(current_res)))
                        del a[aptr]
                aname_seen = []
                nr = nr + 1
                f.append(na)
                r.append(rnam)
                respdbnum.append(rnum)
                current_res = rnum
                a.extend(get_residue_atoms(rnam))
                start = na
                end = len(a)
            aname = line[12:16]
            if aname == ' OXT': 
                na = na -1
                continue
            aname_seen.append(aname)
            atom = get_atom_with_name(a, start, end, aname)
            if atom is None:
                raise NameError, 'Unknown atom "%s" in resiude %s' % (aname,
                                                                      rnam)
            atom.x = float(line[30:38])
            atom.y = float(line[38:46])
            atom.z = float(line[46:54])
            atom.resnum = nr
        elif line[:4] == 'COMP':
            p.name = strip(line[10:70])
        elif line[:3] in ('TER', 'END'):
            break


    alim = 0
    if len(f): alim = f[-1]
    for aptr in xrange(len(a)-1, alim, -1):
        if a[aptr].name not in aname_seen:
#           print (nfstr % (a[aptr].name, int(current_res)))
            del a[aptr]


    p.num_atoms = na + 1
    nr = p.num_residues = nr + 1
    f.append(p.num_atoms)

    # create and initialize Zmatrix

    for i in range(nr):
        res = RESIDUES[r[i]]
        for atom_record in res:
            atom = get_atom_with_name(a, f[i], f[i+1],
                                      atom_record[0])
            if not atom: continue
            
            fp, fpname = atom_record[1]
            fp = fp + i

            if fp &gt;= 0 and fp &lt; nr:
                atom.first_parent = get_atom_with_name(a, f[fp],
                                                       f[fp+1], fpname)
                atom.fp_distance = atom.distance(atom.first_parent)
            fp, fpname = atom_record[2]
            fp = fp + i
            if fp &gt;= 0 and fp &lt; nr:
                atom.second_parent = get_atom_with_name(a, f[fp],
                                                       f[fp+1], fpname)
                atom.sp_angle = angle(atom, atom.first_parent,
                                             atom.second_parent)
            fp, fpname = atom_record[3]
            fp = fp + i
            if fp &gt;= 0 and fp &lt; nr:
                atom.third_parent = get_atom_with_name(a, f[fp],
                                                       f[fp+1], fpname)
                atom.tp_torsion = torsion(atom, atom.first_parent,
                                                 atom.second_parent,
                                                 atom.third_parent)
#           print i,atom.name,atom.fp_distance,atom.sp_angle,atom.tp_torsion

    commit_internal_coords(a)

    # setup pointers to the phi, psi, omega and chi torsions for each
    # residue

    p.phi_atoms = [None]*nr
    p.psi_atoms = [None]*nr
    p.omega_atoms = [None]*nr
    p.chi_atoms = [None]*nr
    
            
    for i in range(p.num_residues):
        name = r[i]
        try:
            phiatom = PhiAtoms[name]
        except KeyError:
            phiatom = PhiAtoms['ANY']

        try:
            psiatom = PsiAtoms[name]
        except KeyError:
            psiatom = PsiAtoms['ANY']

        try:
            omeatom = OmeAtoms[name]
        except KeyError:
            omeatom = OmeAtoms['ANY']

        chiatom = ChiAtoms.get(name, None)

        start = f[i]
        end = f[i+1]

        atom = get_atom_with_name(a, start, end, phiatom)
        if atom: p.phi_atoms[i] = atom

        atom = get_atom_with_name(a, start, end, psiatom)
        if atom: p.psi_atoms[i] = atom

        atom = get_atom_with_name(a, start, end, omeatom)
        if atom: p.omega_atoms[i] = atom

        if chiatom:
            chis = ()
            for atomname in chiatom:
                chis = chis + (get_atom_with_name(a, start, end, atomname),)
            p.chi_atoms[i] = chis
            
def angle(a1, a2, a3):
    """Angle in degrees between 3 atoms"""
    x1 = a1.x - a2.x
    y1 = a1.y - a2.y
    z1 = a1.z - a2.z

    x2 = a3.x - a2.x
    y2 = a3.y - a2.y
    z2 = a3.z - a2.z

    try:
        ang = (x1*x2 + y1*y2 + z1*z2)/sqrt((x1**2.0 + y1**2.0 + z1**2.0) *
                                           (x2**2.0 + y2**2.0 + z2**2.0))
    except ZeroDivisionError:
        return -999.999
    
    try:
        return acos(ang)*RADIANS_TO_DEGREES
    except ValueError:
        if ang &lt; -1.0:
            return 180.0
        else:
            return 0.0

def torsion(a1, a2, a3, a4):
    """Torsion angle in degrees between 4 atoms"""
    a = a1.x - a2.x; b = a1.y - a2.y; c = a1.z - a2.z;
    d = a3.x - a2.x; e = a3.y - a2.y; f = a3.z - a2.z;
    g = a4.x - a3.x; h = a4.y - a3.y; i = a4.z - a3.z;

    ax = b*f - e*c; ay = c*d - f*a; az = a*e - b*d;
    bx = h*f - e*i; by = i*d - f*g; bz = g*e - h*d;

    try:
        ang = (ax*bx+ay*by+az*bz)/sqrt((ax**2.0 + ay**2.0 + az**2.0) *
                                       (bx**2.0 + by**2.0 + bz**2.0))
    except ZeroDivisionError:
        return -999.999
    else:
        try:
            ang = acos(ang) * RADIANS_TO_DEGREES
        except ValueError:
            if ang &lt; -1.0:
                return 180.0
            else:
                return 0.0
        else:
            if ax*(e*bz-f*by) + ay*(f*bx-d*bz) + az*(d*by-e*bx) &gt;= 0.0:
                return -ang
            else:
                return ang

def commit_internal_coords(atoms, start=None, end=None):
    if start is None: start = 0
    if end is None: end = len(atoms)
    for i in range(start, end):
        atoms[i].commit_internal_coords()

default_hbdpar = {
    'use_hbond': 0,
    'use_sidechain_hbond': 0,
    'hbond_distance': 3.0,
    'hbond_probe': 1.5,
    'hbond_score_short': 0.5,
    'hbond_score_long': 1.0,
    'hbond_torsion': 130.0,
    'sidechain_hbond_distance': 3.0,
    'sidechain_hbond_score': 1.0,
    'sidechain_hbond_torsion': 130.0,
    'hbond_winmin': 2,
    'hbond_winmax': 6
}

converters = {
    type(0.0): float,
    type(0): int,
    type(''): str
}

def print_hbs_chasa(file, mol, flags, probe=1.4, ndiv=3,
                    ext_atoms=None, solv_list=None,ext_radius=1.25, title=''):
    """Prints PDB file containing:
           o number of solvating waters on backbone groups - occupancy column
           o CASA - Bfactor column
           o number of unsatisfied backbone Hbond D and A - in TER record

    Arguments

      o *filename* - filename or open file - location to which ASA
        information will be written

      o *mol* - instance of linusMol - molecule whose accessible surface
        area is to be calculated

      o *flags* - object - the object returned by the *make_asa_list*
        function in the **construct** module.

      o *probe* - float - the water probe size to use for accessible surface

      o *ndiv* - integer - a number from 1 to 5 representing the number of
        sampling points used in the area calculation.  From least (1) to
        most (5) accurate, this is either 60, 240, 960, 3840, or 15,350
        points.

      o *ext_atoms* - list of coordinates - a list of tuples, each
        tuple containing (x, y, z) of an additional atom to include
        in the ASA calculation

      o *ext_radius* - float - the radius of the atoms in ext_atoms.

      o *title* - string - An optional title for the PDB output
    """

    atoms = mol.atoms
    fai = mol.residue_first_atom_indices
    residue_names = mol.residue_names
    minres, maxres = get_res_extents(mol)

  
    use_ext = 0
    ext_coords = None

    use_data = 1
    data = Numeric.zeros(len(atoms), 'd')
  
    if ext_atoms:
        ext_coords = []
        map(lambda x: map(lambda y: ext_coords.append(y), x), ext_atoms)
        ext_coords = Numeric.array(ext_coords, 'd')
        use_ext = len(ext_atoms)

    tot_asa =  asa_evaluate(atoms, data, ext_coords, flags, probe,
                            ext_radius, use_data, use_ext, ndiv)

    ff = file
    if type(file) == type('hi'):
        ff = open(file, 'wb')

    write = ff.write
    pdbfmt = 'ATOM  %5d %-4s %-3s  %4d    %8.3f%8.3f%8.3f%6.2f%6.2f\n'

    p_solv_nrg = 0.0
    ap_solv_nrg = 0.0
    Gamma_p = 3.0/5.0
    Gamma_hb_oxy = 2.0
    Gamma_ap = 0.03
    CHASA = 0.0
    for i in xrange(minres,maxres):
        rname = residue_names[i]
        start = fai[i]
        end = fai[i+1]
        occ = 0.0
        for j in range(start, end):
            atom = atoms[j]
            residue_num = int(mol.res_pdb_number[atom.resnum])
            if atom.name == ' N  ':
                if solv_list[i][0][2] &gt; 0:
                    p_solv_nrg = p_solv_nrg - (Gamma_p *(solv_list[i][0][2]))
                elif solv_list[i][0][2] &lt; 0:
                    p_solv_nrg = p_solv_nrg + 1.0
    
                write(pdbfmt % (j+1, atom.name, rname, residue_num, atom.x, atom.y,
                            atom.z, solv_list[i][0][2],data[j]))
            elif atom.name == ' O  ':
                if solv_list[i][1][2] &gt; 0:
                    if solv_list[i][1][3] == 0:
                        p_solv_nrg = p_solv_nrg - (Gamma_p *(solv_list[i][1][2]))
                    elif solv_list[i][1][3] &gt; 0:
                        p_solv_nrg = p_solv_nrg - (Gamma_hb_oxy)
                elif solv_list[i][1][2] &lt; 0:
                    p_solv_nrg = p_solv_nrg + 1.0
    
                write(pdbfmt % (j+1, atom.name, rname, residue_num, atom.x, atom.y,
                            atom.z, solv_list[i][1][2],data[j]))
            else:
                write(pdbfmt % (j+1, atom.name, rname, residue_num, atom.x, atom.y,
                            atom.z, occ,data[j]))
            if 'C' in atom.name:
                ap_solv_nrg = ap_solv_nrg + (Gamma_ap * data[j])
                CHASA = CHASA + data[j]

    aidx = j+2
    ridx = residue_num+1

    if ext_atoms:
        for x,y,z in ext_atoms:
            write(pdbfmt % (aidx, ' O  ', 'HOH', ridx, x, y, z, 0.0, 0.0))
            aidx = aidx + 1
            ridx = ridx + 1

    tot_solv_nrg = ap_solv_nrg + p_solv_nrg
    write('TER %8.3f %8.3f\n' % (CHASA,tot_solv_nrg))

    if ff &lt;&gt; file:
        write('END\n')
        ff.close()

    return tot_asa,tot_solv_nrg

def find_neighbors(atoms, coords, flags, probe, ext_rad,
         k, numatm, numext,nghlst):
    """
    Given an atom index k, this function calculates which atoms are close
    enough to be relevant in solvent accessible surface area calculations.
    This set of atom indices are stored in nhhlst, and the size of nghlst
    is returned.  Positive indices in nghlst represent atoms in the protein
    atom list atoms.  Negative indices (offset by one) represent atom
    coordinates stored in the external coordinate list.
    """
    numngh = 0

    a = atoms[k]
    arad = a.radius + probe + probe
    
    for i in range(0,k):
        if not (flags[i]): continue
        b = atoms[i]
        d = a.distance(b)
        if d &lt; (arad + b.radius):
#           print 'found neighbor ',atoms[i].name
            nghlst.append(i)

    for i in range(k+1,numatm):
        if not (flags[i]): continue
        b = atoms[i]
        d = a.distance(b)   
        if d &lt; (arad + b.radius):
#           print 'found neighbor ',atoms[i].name
            nghlst.append(i)    

    numngh = len(nghlst)
    if not numext: return numngh

    ax = a.x
    ay = a.y
    az = a.z
    for i in range(numext):
        pos = i*3
        bx = coords[pos] - ax
        by = coords[pos+1] - ay
        bz = coords[pos+2] - az
        d = sqrt(bx*bx + by*by + bz*bz)
        if d &lt; (arad + ext_rad):
            nghlst.append(-(i+1))

    numngh = len(nghlst)
    return numngh

def tri_buried(atoms, coords, tx, ty, tz, probe, ext_rad,
                         nghlst, numngh, nghstrt):
    """
     This function determines whether one of the radial points of an atom 
    is bured by steric contacts with other solvated radii.  tx, ty, tz    
    are points on the solvated radii of the atom in question, i.e.
    r = atom radius + solvent radius.  This function checks whether
    any of the other solvated atoms collide with this particular point.
    If it is buried within the molecule, this function returns true.

    """
    for k in range(nghstrt,numngh):
        i = nghlst[k]
        if i &gt;= 0:
            a = atoms[i]
            dx = a.x - tx
            dy = a.y - ty
            dz = a.z - tz
            d = dx*dx + dy*dy + dz*dz
            r = probe + a.radius
            if d &lt; (r*r):
                nghstrt = k
                return 1,nghstrt

        else:
            pos = 3*(-(i+1))
            dx = coords[pos] - tx
            dy = coords[pos+1] - ty
            dz = coords[pos+2] - tz
            d = dx*dx + dy*dy + dz*dz
            r = ext_rad + probe
            if d &lt; (r*r):
                nghstrt = k
                return 1,nghstrt

    for k in range(nghstrt):
        i = nghlst[k]
        if i &gt;= 0:
            a = atoms[i]
            dx = a.x - tx
            dy = a.y - ty
            dz = a.z - tz
            d = dx*dx + dy*dy + dz*dz
            r = probe + a.radius
            if d &lt; (r*r):
                nghstrt = k
                return 1,nghstrt

        else:
            pos = 3*(-(i+1))
            dx = coords[pos] - tx
            dy = coords[pos+1] - ty
            dz = coords[pos+2] - tz
            d = dx*dx + dy*dy + dz*dz
            r = ext_rad + probe
            if d &lt; (r*r):
                nghstrt = k
                return 1,nghstrt

    return 0,nghstrt


def asa_evaluate(atoms, data, ext_coords, flags, probe,
                        ext_radius, use_data, use_ext, ndiv):
    """
    Given a list of atoms and a Numeric array of flags, for each atom,   
    this function determines the accessible surface area for each atom    
    according to the flags given.  If use_data is true, then the ASA for
    each flagged atom will be stored in data (a Numeric array if it used).
    If use_ext is nonzero, coords is assumed to hold a list of use_ext 
    atom coordinates to include in the calculation. ext_rad is the radius
    of external water atoms in coords.  probe is the probe water radius,
    and triangles and ndiv both specify the precision of the calculation.
    To cite this algorithm, use
   
    Shrake, A., and J. A. Rupley.  "Environment and Exposure to Solvent of
      Protein Atoms.  Lysozyme and Insulin."  J. Mol. Bio. 79 (1973): 351-
      371.

    """
    ntrian = 960
    triarea = 4.0 * pi / ntrian
    numatm = len(atoms)
    if use_ext:
        numext = len(ext_coords)/3
    else:
        numext = 0
    PRESENT = 1
    CONTRIB = 2
    MAXNEIGHBORS = 200
    atot = 0.0
    for k in range(len(atoms)):
        nghlst = []
        if not (flags[k]): 
            continue
        a = atoms[k]
       
        numngh = find_neighbors(atoms, ext_coords, flags, probe, ext_radius,  
                 k, numatm, numext,nghlst)
        nghstrt = 0
        nts = 0
        r = probe + a.radius
        x = a.x
        y = a.y
        z = a.z
        for i in range(ntrian):
            tx = triangles[i][0]*r + x
            ty = triangles[i][1]*r + y
            tz = triangles[i][2]*r + z
            buried,nghstrt = tri_buried(atoms, ext_coords, tx, ty, tz, probe, ext_radius,
                         nghlst, numngh, nghstrt)
            if not buried:
                nts = nts + 1
        area = nts * triarea * r * r 
        if use_data:
            data[k] = area
        
        if (flags[k] &amp; CONTRIB):
            atot += area

    return atot

def make_hbond_list(p, wmin, wmax, hbparms):
    """
    Make list of all possible hydrogen bonding atom pairs. 
    """

    fai = p.residue_first_atom_indices
    rn = p.residue_names
    atoms = p.atoms

    minres, maxres = get_res_extents(p)

    hblist = []
    hbdist = hbparms['hbond_distance']
    hbdmax = hbdist + hbparms['hbond_probe']
    hbdtor = hbparms['hbond_torsion']
    hbenes = hbparms['hbond_score_short']
    hbenel = hbparms['hbond_score_long']
    shbdis = hbparms['sidechain_hbond_distance']
    shbdmax = shbdis + hbparms['hbond_probe']
    shbtor = hbparms['sidechain_hbond_torsion']
    shbene = hbparms['sidechain_hbond_score']

    for i in xrange(minres, maxres):
        start = fai[i]
        end = fai[i+1]
        if rn[i] &lt;&gt; 'PRO':
            if not ((i == 0) and (minres == 0)):
                # for NH as the donor
                donor = get_atom_named(atoms, ' N  ', start, end)
                da1 = get_atom_named(atoms, ' CA ', start, end)
                da2 = get_atom_named(atoms, ' C  ', fai[i-1], start)
                acps = []
                for j in range(wmin, wmax):
                    pos = i-j
                    if pos &gt;= minres:
                        if hbparms['use_hbond']:
                            acp = get_atom_named(atoms, ' O  ', fai[pos],
                                                 fai[pos+1])
                            acp1 = get_atom_named(atoms, ' C  ', fai[pos],
                                                 fai[pos+1])
                            if j &lt; 6:
                                acps.append((acp, acp1, hbdist, hbdmax, hbdtor, hbenes))
                            else:
                                acps.append((acp, acp1, hbdist, hbdmax, hbdtor, hbenel))
    
                    pos = i+j
                    if pos &lt; maxres:
                        if hbparms['use_hbond']:
                            acp = get_atom_named(atoms, ' O  ', fai[pos],
                                                fai[pos+1])
                            acp1 = get_atom_named(atoms, ' C  ', fai[pos],
                                                 fai[pos+1])
                            if j &lt; 6:
                                acps.append((acp, acp1, hbdist, hbdmax, hbdtor, hbenes))
                            else:
                                acps.append((acp, acp1, hbdist, hbdmax, hbdtor, hbenel))
    
                # look for sidechain acceptors with wmin = 0, so GLU can Hbond itself,etc
                wmin = 0
                for j in range(wmin, wmax):
                    pos = i-j
                    if pos &gt;= minres:
                        if hbparms['use_sidechain_hbond']:
                            scatoms = get_scacceptor(atoms, fai[pos], fai[pos+1],
                                                     rn[pos])
                            if scatoms:
                                for acp,acp1, in scatoms:
                                    acps.append((acp,acp1, shbdis, shbdmax, shbtor,
                                                 shbene))
    
                    pos = i+j
                    if pos &lt; maxres:
                        if hbparms['use_sidechain_hbond']:
                            scatoms = get_scacceptor(atoms, fai[pos], fai[pos+1],
                                                     rn[pos])
                            if scatoms:
                                for acp,acp1, in scatoms:
                                    acps.append((acp, acp1,shbdis, shbdmax, shbtor,
                                                 shbene))
    
                hblist.append(((donor, da1, da2), tuple(acps)))

            # for ARG NH1 or NH2 as the donor
            if rn[i] == 'ARG':
                donor = get_atom_named(atoms, ' NH1', start, end)
                da1 = get_atom_named(atoms, ' CZ ', start, end)
                da2 = get_atom_named(atoms, ' NE ', start, end)
                acps = []
                for j in range(wmin, wmax):
                    pos = i-j
                    if pos &gt;= minres:
                        if hbparms['use_sidechain_hbond']:
                            acp = get_atom_named(atoms, ' O  ', fai[pos],
                                                 fai[pos+1])
                            acp1 = get_atom_named(atoms, ' C  ', fai[pos],
                                             fai[pos+1])
                            acps.append((acp, acp1, shbdis, shbdmax, shbtor,shbene))

                    pos = i+j
                    if pos &lt; maxres:
                        if hbparms['use_sidechain_hbond']:
                            acp = get_atom_named(atoms, ' O  ', fai[pos],
                                                fai[pos+1])
                            acp1 = get_atom_named(atoms, ' C  ', fai[pos],
                                             fai[pos+1])
                            acps.append((acp, acp1, shbdis, shbdmax, shbtor,shbene))
                hblist.append(((donor, da1, da2), tuple(acps)))

                donor = get_atom_named(atoms, ' NH2', start, end)
                da1 = get_atom_named(atoms, ' CZ ', start, end)
                da2 = get_atom_named(atoms, ' NE ', start, end)
                acps = []
                for j in range(wmin, wmax):
                    pos = i-j
                    if pos &gt;= minres:
                        if hbparms['use_sidechain_hbond']:
                            acp = get_atom_named(atoms, ' O  ', fai[pos],
                                                 fai[pos+1])
                            acp1 = get_atom_named(atoms, ' C  ', fai[pos],
                                             fai[pos+1])
                            acps.append((acp, acp1, shbdis, shbdmax, shbtor,shbene))

                    pos = i+j
                    if pos &lt; maxres:
                        if hbparms['use_sidechain_hbond']:
                            acp = get_atom_named(atoms, ' O  ', fai[pos],
                                                fai[pos+1])
                            acp1 = get_atom_named(atoms, ' C  ', fai[pos],
                                             fai[pos+1])
                            acps.append((acp, acp1, shbdis, shbdmax, shbtor,shbene))
                hblist.append(((donor, da1, da2), tuple(acps)))

            # for ASN OD1 as the donor(we don't know about the flip state)
            if rn[i] == 'ASN':
                donor = get_atom_named(atoms, ' OD1', start, end)
                da1 = get_atom_named(atoms, ' CG ', start, end)
                da2 = get_atom_named(atoms, ' CB ', start, end)
                acps = []
                for j in range(wmin, wmax):
                    pos = i-j
                    if pos &gt;= minres:
                        if hbparms['use_sidechain_hbond']:
                            acp = get_atom_named(atoms, ' O  ', fai[pos],
                                                 fai[pos+1])
                            acp1 = get_atom_named(atoms, ' C  ', fai[pos],
                                             fai[pos+1])
                            acps.append((acp, acp1, shbdis, shbdmax, shbtor,shbene))

                    pos = i+j
                    if pos &lt; maxres:
                        if hbparms['use_sidechain_hbond']:
                            acp = get_atom_named(atoms, ' O  ', fai[pos],
                                                fai[pos+1])
                            acp1 = get_atom_named(atoms, ' C  ', fai[pos],
                                             fai[pos+1])
                            acps.append((acp, acp1, shbdis, shbdmax, shbtor,shbene))
                hblist.append(((donor, da1, da2), tuple(acps)))

            # for ASN ND1 as the donor
            if rn[i] == 'ASN':
                donor = get_atom_named(atoms, ' ND2', start, end)
                da1 = get_atom_named(atoms, ' CG ', start, end)
                da2 = get_atom_named(atoms, ' CB ', start, end)
                acps = []
                for j in range(wmin, wmax):
                    pos = i-j
                    if pos &gt;= minres:
                        if hbparms['use_sidechain_hbond']:
                            acp = get_atom_named(atoms, ' O  ', fai[pos],
                                                 fai[pos+1])
                            acp1 = get_atom_named(atoms, ' C  ', fai[pos],
                                             fai[pos+1])
                            acps.append((acp, acp1, shbdis, shbdmax, shbtor,shbene))

                    pos = i+j
                    if pos &lt; maxres:
                        if hbparms['use_sidechain_hbond']:
                            acp = get_atom_named(atoms, ' O  ', fai[pos],
                                                fai[pos+1])
                            acp1 = get_atom_named(atoms, ' C  ', fai[pos],
                                             fai[pos+1])
                            acps.append((acp, acp1, shbdis, shbdmax, shbtor,shbene))
                hblist.append(((donor, da1, da2), tuple(acps)))

            # for ASP OD1 as the donor
            if rn[i] == 'ASN':
                donor = get_atom_named(atoms, ' OD1', start, end)
                da1 = get_atom_named(atoms, ' CG ', start, end)
                da2 = get_atom_named(atoms, ' CB ', start, end)
                acps = []
                for j in range(wmin, wmax):
                    pos = i-j
                    if pos &gt;= minres:
                        if hbparms['use_sidechain_hbond']:
                            acp = get_atom_named(atoms, ' O  ', fai[pos],
                                                 fai[pos+1])
                            acp1 = get_atom_named(atoms, ' C  ', fai[pos],
                                             fai[pos+1])
                            acps.append((acp, acp1, shbdis, shbdmax, shbtor,shbene))

                    pos = i+j
                    if pos &lt; maxres:
                        if hbparms['use_sidechain_hbond']:
                            acp = get_atom_named(atoms, ' O  ', fai[pos],
                                                fai[pos+1])
                            acp1 = get_atom_named(atoms, ' C  ', fai[pos],
                                             fai[pos+1])
                            acps.append((acp, acp1, shbdis, shbdmax, shbtor,shbene))
                hblist.append(((donor, da1, da2), tuple(acps)))

            # for ASP OD2 as the donor
            if rn[i] == 'ASP':
                donor = get_atom_named(atoms, ' OD2', start, end)
                da1 = get_atom_named(atoms, ' CG ', start, end)
                da2 = get_atom_named(atoms, ' CB ', start, end)
                acps = []
                for j in range(wmin, wmax):
                    pos = i-j
                    if pos &gt;= minres:
                        if hbparms['use_sidechain_hbond']:
                            acp = get_atom_named(atoms, ' O  ', fai[pos],
                                                 fai[pos+1])
                            acp1 = get_atom_named(atoms, ' C  ', fai[pos],
                                             fai[pos+1])
                            acps.append((acp, acp1, shbdis, shbdmax, shbtor,shbene))

                    pos = i+j
                    if pos &lt; maxres:
                        if hbparms['use_sidechain_hbond']:
                            acp = get_atom_named(atoms, ' O  ', fai[pos],
                                                fai[pos+1])
                            acp1 = get_atom_named(atoms, ' C  ', fai[pos],
                                             fai[pos+1])
                            acps.append((acp, acp1, shbdis, shbdmax, shbtor,shbene))
                hblist.append(((donor, da1, da2), tuple(acps)))

            # for GLN OE1 as the donor
            if rn[i] == 'GLN':  
                donor = get_atom_named(atoms, ' OE1', start, end)
                da1 = get_atom_named(atoms, ' CD ', start, end)
                da2 = get_atom_named(atoms, ' CG ', start, end)
                acps = []
                for j in range(wmin, wmax):
                    pos = i-j
                    if pos &gt;= minres:
                        if hbparms['use_sidechain_hbond']:
                            acp = get_atom_named(atoms, ' O  ', fai[pos],
                                                 fai[pos+1])
                            acp1 = get_atom_named(atoms, ' C  ', fai[pos],
                                             fai[pos+1])
                            acps.append((acp, acp1, shbdis, shbdmax, shbtor,shbene))

                    pos = i+j
                    if pos &lt; maxres:
                        if hbparms['use_sidechain_hbond']:
                            acp = get_atom_named(atoms, ' O  ', fai[pos],
                                                fai[pos+1])    
                            acp1 = get_atom_named(atoms, ' C  ', fai[pos],
                                             fai[pos+1])
                            acps.append((acp, acp1, shbdis, shbdmax, shbtor,shbene))
                hblist.append(((donor, da1, da2), tuple(acps)))

            # for GLN NE2 as the donor
            if rn[i] == 'GLN':  
                donor = get_atom_named(atoms, ' NE2', start, end)
                da1 = get_atom_named(atoms, ' CD ', start, end)
                da2 = get_atom_named(atoms, ' CG ', start, end)
                acps = []
                for j in range(wmin, wmax):
                    pos = i-j
                    if pos &gt;= minres:
                        if hbparms['use_sidechain_hbond']:
                            acp = get_atom_named(atoms, ' O  ', fai[pos],
                                                 fai[pos+1])
                            acp1 = get_atom_named(atoms, ' C  ', fai[pos],
                                             fai[pos+1])
                            acps.append((acp, acp1, shbdis, shbdmax, shbtor,shbene))

                    pos = i+j
                    if pos &lt; maxres:
                        if hbparms['use_sidechain_hbond']:
                            acp = get_atom_named(atoms, ' O  ', fai[pos],
                                                fai[pos+1])
                            acp1 = get_atom_named(atoms, ' C  ', fai[pos],
                                             fai[pos+1])
                            acps.append((acp, acp1, shbdis, shbdmax, shbtor,shbene))
                hblist.append(((donor, da1, da2), tuple(acps)))

            # for GLU OE1 as the donor
            if rn[i] == 'GLU':
                donor = get_atom_named(atoms, ' OE1', start, end)
                da1 = get_atom_named(atoms, ' CD ', start, end)
                da2 = get_atom_named(atoms, ' CG ', start, end)
                acps = []
                for j in range(wmin, wmax):
                    pos = i-j
                    if pos &gt;= minres:
                        if hbparms['use_sidechain_hbond']:
                            acp = get_atom_named(atoms, ' O  ', fai[pos],
                                                 fai[pos+1])
                            acp1 = get_atom_named(atoms, ' C  ', fai[pos],
                                             fai[pos+1])
                            acps.append((acp, acp1, shbdis, shbdmax, shbtor,shbene))

                    pos = i+j
                    if pos &lt; maxres:
                        if hbparms['use_sidechain_hbond']:
                            acp = get_atom_named(atoms, ' O  ', fai[pos],
                                                fai[pos+1])
                            acp1 = get_atom_named(atoms, ' C  ', fai[pos],
                                             fai[pos+1])
                            acps.append((acp, acp1, shbdis, shbdmax, shbtor,shbene))
                hblist.append(((donor, da1, da2), tuple(acps)))

            # for GLU OE2 as the donor
            if rn[i] == 'GLU':
                donor = get_atom_named(atoms, ' OE2', start, end)
                da1 = get_atom_named(atoms, ' CD ', start, end)
                da2 = get_atom_named(atoms, ' CG ', start, end)
                acps = []
                for j in range(wmin, wmax):
                    pos = i-j
                    if pos &gt;= minres:
                        if hbparms['use_sidechain_hbond']:
                            acp = get_atom_named(atoms, ' O  ', fai[pos],
                                                 fai[pos+1])
                            acp1 = get_atom_named(atoms, ' C  ', fai[pos],
                                             fai[pos+1])
                            acps.append((acp, acp1, shbdis, shbdmax, shbtor,shbene))

                    pos = i+j
                    if pos &lt; maxres:
                        if hbparms['use_sidechain_hbond']:
                            acp = get_atom_named(atoms, ' O  ', fai[pos],
                                                fai[pos+1])
                            acp1 = get_atom_named(atoms, ' C  ', fai[pos],
                                             fai[pos+1])
                            acps.append((acp, acp1, shbdis, shbdmax, shbtor,shbene))
                hblist.append(((donor, da1, da2), tuple(acps)))

            # for HIS NE2 as the donor
            if rn[i] == 'HIS':
                donor = get_atom_named(atoms, ' NE2', start, end)
                da1 = get_atom_named(atoms, ' CD2', start, end)
                da2 = get_atom_named(atoms, ' CG ', start, end)
                acps = []
                for j in range(wmin, wmax):
                    pos = i-j
                    if pos &gt;= minres:
                        if hbparms['use_sidechain_hbond']:
                            acp = get_atom_named(atoms, ' O  ', fai[pos],
                                                 fai[pos+1])
                            acp1 = get_atom_named(atoms, ' C  ', fai[pos],
                                             fai[pos+1])
                            acps.append((acp, acp1, shbdis, shbdmax, shbtor,shbene))

                    pos = i+j
                    if pos &lt; maxres:
                        if hbparms['use_sidechain_hbond']:
                            acp = get_atom_named(atoms, ' O  ', fai[pos],
                                                fai[pos+1])
                            acp1 = get_atom_named(atoms, ' C  ', fai[pos],
                                             fai[pos+1])
                            acps.append((acp, acp1, shbdis, shbdmax, shbtor,shbene))
                hblist.append(((donor, da1, da2), tuple(acps)))

            # for HIS ND1 as the donor (see 5cyt, HIS 18 ND1 to PRO 30 O)
            if rn[i] == 'HIS':
                donor = get_atom_named(atoms, ' ND1', start, end)
                da1 = get_atom_named(atoms, ' CG ', start, end)
                da2 = get_atom_named(atoms, ' CB ', start, end)
                acps = []
                for j in range(wmin, wmax):
                    pos = i-j
                    if pos &gt;= minres:
                        if hbparms['use_sidechain_hbond']:
                            acp = get_atom_named(atoms, ' O  ', fai[pos],
                                                 fai[pos+1])
                            acp1 = get_atom_named(atoms, ' C  ', fai[pos],
                                             fai[pos+1])
                            acps.append((acp, acp1, shbdis, shbdmax, shbtor,shbene))

                    pos = i+j
                    if pos &lt; maxres:
                        if hbparms['use_sidechain_hbond']:
                            acp = get_atom_named(atoms, ' O  ', fai[pos],
                                                fai[pos+1])
                            acp1 = get_atom_named(atoms, ' C  ', fai[pos],
                                             fai[pos+1])
                            acps.append((acp, acp1, shbdis, shbdmax, shbtor,shbene))
                hblist.append(((donor, da1, da2), tuple(acps)))

            # for LYS NZ as the donor
            if rn[i] == 'LYS':
                donor = get_atom_named(atoms, ' NZ ', start, end)
                da1 = get_atom_named(atoms, ' CE ', start, end)
                da2 = get_atom_named(atoms, ' CD ', start, end)
                acps = []
                for j in range(wmin, wmax):
                    pos = i-j
                    if pos &gt;= minres:
                        if hbparms['use_sidechain_hbond']:
                            acp = get_atom_named(atoms, ' O  ', fai[pos],
                                                 fai[pos+1])
                            acp1 = get_atom_named(atoms, ' C  ', fai[pos],
                                             fai[pos+1])
                            acps.append((acp, acp1, shbdis, shbdmax, shbtor,shbene))

                    pos = i+j
                    if pos &lt; maxres:
                        if hbparms['use_sidechain_hbond']:
                            acp = get_atom_named(atoms, ' O  ', fai[pos],
                                                fai[pos+1])
                            acp1 = get_atom_named(atoms, ' C  ', fai[pos],
                                             fai[pos+1])
                            acps.append((acp, acp1, shbdis, shbdmax, shbtor,shbene))
                hblist.append(((donor, da1, da2), tuple(acps)))

            # for SER OG as the donor  
            if rn[i] == 'SER':
                donor = get_atom_named(atoms, ' OG ', start, end)
                da1 = get_atom_named(atoms, ' CB ', start, end)
                da2 = get_atom_named(atoms, ' CA ', start, end)
                acps = []
                for j in range(wmin, wmax):
                    pos = i-j
                    if pos &gt;= minres:
                        if hbparms['use_sidechain_hbond']:
                            acp = get_atom_named(atoms, ' O  ', fai[pos],
                                                 fai[pos+1])
                            acp1 = get_atom_named(atoms, ' C  ', fai[pos],
                                             fai[pos+1])
                            acps.append((acp, acp1, shbdis, shbdmax, shbtor,shbene))

                    pos = i+j      
                    if pos &lt; maxres:
                        if hbparms['use_sidechain_hbond']:
                            acp = get_atom_named(atoms, ' O  ', fai[pos],
                                                fai[pos+1])
                            acp1 = get_atom_named(atoms, ' C  ', fai[pos],
                                             fai[pos+1])
                            acps.append((acp, acp1, shbdis, shbdmax, shbtor,shbene))
                hblist.append(((donor, da1, da2), tuple(acps)))

            # for THR OG1 as the donor
            if rn[i] == 'THR':
                donor = get_atom_named(atoms, ' OG1', start, end)
                da1 = get_atom_named(atoms, ' CB ', start, end)
                da2 = get_atom_named(atoms, ' CA ', start, end)
                acps = []
                for j in range(wmin, wmax):
                    pos = i-j
                    if pos &gt;= minres:
                        if hbparms['use_sidechain_hbond']:
                            acp = get_atom_named(atoms, ' O  ', fai[pos],
                                                 fai[pos+1])
                            acp1 = get_atom_named(atoms, ' C  ', fai[pos],
                                             fai[pos+1])
                            acps.append((acp, acp1, shbdis, shbdmax, shbtor,shbene))

                    pos = i+j
                    if pos &lt; maxres:
                        if hbparms['use_sidechain_hbond']:
                            acp = get_atom_named(atoms, ' O  ', fai[pos],
                                                fai[pos+1])
                            acp1 = get_atom_named(atoms, ' C  ', fai[pos],
                                             fai[pos+1])
                            acps.append((acp, acp1, shbdis, shbdmax, shbtor,shbene))
                hblist.append(((donor, da1, da2), tuple(acps)))      

            # for TYR OH as the donor
            if rn[i] == 'TYR':
                donor = get_atom_named(atoms, ' OH ', start, end)
                da1 = get_atom_named(atoms, ' CZ ', start, end)
                da2 = get_atom_named(atoms, ' CG ', start, end)
                acps = []
                for j in range(wmin, wmax):
                    pos = i-j
                    if pos &gt;= minres:
                        if hbparms['use_sidechain_hbond']:
                            acp = get_atom_named(atoms, ' O  ', fai[pos],
                                                 fai[pos+1])
                            acp1 = get_atom_named(atoms, ' C  ', fai[pos],
                                             fai[pos+1])
                            acps.append((acp, acp1, shbdis, shbdmax, shbtor,shbene))

                    pos = i+j
                    if pos &lt; maxres:
                        if hbparms['use_sidechain_hbond']:
                            acp = get_atom_named(atoms, ' O  ', fai[pos],
                                                fai[pos+1])
                            acp1 = get_atom_named(atoms, ' C  ', fai[pos],
                                             fai[pos+1])
                            acps.append((acp, acp1, shbdis, shbdmax, shbtor,shbene))
                hblist.append(((donor, da1, da2), tuple(acps)))

            # for ARG NH1 or NH2 as the donor
            if rn[i] == 'ARG':
                donor = get_atom_named(atoms, ' NH1', start, end)
                da1 = get_atom_named(atoms, ' CZ ', start, end)
                da2 = get_atom_named(atoms, ' NE ', start, end)
                acps = []
                for j in range(wmin, wmax):
                    pos = i-j
                    if pos &gt;= minres:
                        if hbparms['use_sidechain_hbond']:
                            acp = get_atom_named(atoms, ' O  ', fai[pos],
                                                 fai[pos+1])
                            acp1 = get_atom_named(atoms, ' C  ', fai[pos],
                                             fai[pos+1])
                            acps.append((acp, acp1, shbdis, shbdmax, shbtor,shbene))

                    pos = i+j
                    if pos &lt; maxres:
                        if hbparms['use_sidechain_hbond']:
                            acp = get_atom_named(atoms, ' O  ', fai[pos],
                                                fai[pos+1])
                            acp1 = get_atom_named(atoms, ' C  ', fai[pos],
                                             fai[pos+1])
                            acps.append((acp, acp1, shbdis, shbdmax, shbtor,shbene))
                hblist.append(((donor, da1, da2), tuple(acps)))

                donor = get_atom_named(atoms, ' NH2', start, end)
                da1 = get_atom_named(atoms, ' CZ ', start, end)
                da2 = get_atom_named(atoms, ' NE ', start, end)
                acps = []
                for j in range(wmin, wmax):
                    pos = i-j
                    if pos &gt;= minres:
                        if hbparms['use_sidechain_hbond']:
                            acp = get_atom_named(atoms, ' O  ', fai[pos],
                                                 fai[pos+1])
                            acp1 = get_atom_named(atoms, ' C  ', fai[pos],
                                             fai[pos+1])
                            acps.append((acp, acp1, shbdis, shbdmax, shbtor,shbene))

                    pos = i+j
                    if pos &lt; maxres:
                        if hbparms['use_sidechain_hbond']:
                            acp = get_atom_named(atoms, ' O  ', fai[pos],
                                                fai[pos+1])
                            acp1 = get_atom_named(atoms, ' C  ', fai[pos],
                                             fai[pos+1])
                            acps.append((acp, acp1, shbdis, shbdmax, shbtor,shbene))
                hblist.append(((donor, da1, da2), tuple(acps)))



    return hblist

def is_hydrophobic(res, atm):
    """returns true if an atom in a residue is hydrophobic"""
    return HYDROPHOBIC[atm]

def is_hydrogen(res, atm):
    """returns true if an atom in a residue is a hydrogen"""
    return HYDROGEN.has_key(atm)

def make_asa_list(mol, hphob=1, hphil=0, hydrogens=0):
    """Construct parameters for accessible surface area scoring

    Arguments

        o *mol* - instance of linusMol - molecule whose accessible surface
          area is to be calculated

        o *hphob* - boolean - whether to include hydrophobic atoms in
          the evaluation of surface area (by default, all atoms are used
          in steric calculations except ACE, NME)

        o *hphil* - boolean - whether to include hydrophillic atoms in
          the evaluation of surface area (by default, all atoms are used
          in steric calculations except ACE, NME)

        o *hydrogens* boolean - whether to include hydrogens in the
          scoring function (as contributing to the return value or
          as steric hinderances to solvation).  If included, hydrogens
          are hydrophobic, according to *HYDROPHILLIC* in **linusRes**.

    Result

        Returns an object that can be passed to *asa_score* or
        *print_asa_score* functions in the **linusScore** module.
    """

    # first  bit: whether or not the atom is sterically present
    # second bit: whether or not the atom contributes to total ASA

    PRESENT = 1
    CONTRIB = 2

    atoms = mol.atoms
    fai = mol.residue_first_atom_indices
    flags = [0]*len(atoms)
   
    minres, maxres = get_res_extents(mol)

    for r in xrange(minres, maxres):
        rnam = mol.residue_names[r]
        first, last = fai[r], fai[r+1]
        for at in xrange(first, last):
            if not hydrogens and is_hydrogen(rnam, atoms[at].name):
                continue

            hydrophobic = is_hydrophobic(rnam, atoms[at].name)

            flags[at] = PRESENT

            if hphob and hydrophobic:
                flags[at] = flags[at] | CONTRIB

            if hphil and not hydrophobic:
                flags[at] = flags[at] | CONTRIB

    return Numeric.array(flags, 'i')

def get_atom_named(atoms, atom_name, start, end):
    for i in range(start, end):
        atom = atoms[i]
        if atom.name == atom_name:
            return atom

def get_res_extents(p):
    """get_res_extents (linusMol p)

    Returns the residue index extents suitable for a range statement that
    scans the protein.  This information is returned as a tuple of minres,
    maxres.  Capping residues, such as ACE and NME, are ignored.
    """
    rnam = p.residue_names

    minres = 0
    maxres = p.num_residues

    if len(rnam) and rnam[minres] == 'ACE' : minres = 1
    if len(rnam) and rnam[-1] == 'NME' : maxres = maxres-1

    return minres, maxres

def angle(a1, a2, a3):
    """Angle in degrees between 3 atoms"""
    x1 = a1.x - a2.x
    y1 = a1.y - a2.y
    z1 = a1.z - a2.z      

    x2 = a3.x - a2.x
    y2 = a3.y - a2.y
    z2 = a3.z - a2.z

    try:
        ang = (x1*x2 + y1*y2 + z1*z2)/sqrt((x1**2.0 + y1**2.0 + z1**2.0) *
                                           (x2**2.0 + y2**2.0 + z2**2.0))
    except ZeroDivisionError:
        return -999.999

    try:
        return acos(ang)*RADIANS_TO_DEGREES
    except ValueError:
        if ang &lt; -1.0:
            return 180.0
        else:
            return 0.0

def norm(x1, y1, z1, x2, y2, z2, sqrt=sqrt):
    dx = x1 - x2; dy = y1 - y2; dz = z1 - z2;
    d = 1.0 / sqrt(dx*dx + dy*dy + dz*dz)
    return dx*d, dy*d, dz*d

def ztox(coordsfp, coordssp, coordstp, zmvalu, zmindx):
    """Convert internal coordinates to cartesian coordinates

    Arguments

        o *coords* - list of type 'f' and 3 elements

        o *zmvalu* - list of type 'f' and 3 elements, where each
        element is respecively the distance, angle and torsion angle of
        the atom to it's parent, second parent and third parent.

        o *zmindx* - list of type 'i' and 3 elements, where each
        element is respectively the integer index of the atoms' first,
        second and third parents.

    Returns

        atcoor - the converted x,y,z coords
    """

    fp, sp, tp = zmindx
    fpd, spa, tpt = zmvalu

    spa = spa * DEGREES_TO_RADIANS
    tpt = tpt * DEGREES_TO_RADIANS

    sina = sin(spa)
    cosa = -cos(spa)

    sint = sina * sin(tpt)
    cost = sina * cos(tpt)

    fpx, fpy, fpz = coordsfp
    spx, spy, spz = coordssp
    tpx, tpy, tpz = coordstp

    u2x, u2y, u2z = norm(fpx, fpy, fpz, spx, spy, spz)
    u1x, u1y, u1z = norm(spx, spy, spz, tpx, tpy, tpz)

    cosine = u1x*u2x + u1y*u2y + u1z*u2z
    if abs(cosine) &lt; 1.0:
        sine = 1.0/sqrt(1.0 - cosine*cosine)
    else:
        sine = 1.0/sqrt(cosine*cosine - 1.0)

    u3x = sine * (u1y*u2z - u1z*u2y)
    u3y = sine * (u1z*u2x - u1x*u2z)
    u3z = sine * (u1x*u2y - u1y*u2x)

    u4x = cost * (u3y*u2z - u3z*u2y)
    u4y = cost * (u3z*u2x - u3x*u2z)
    u4z = cost * (u3x*u2y - u3y*u2x)

    dummy = Numeric.zeros((1, 3), 'f')
    atcoor = dummy[0]

    atcoor[0] = fpx + fpd*(u2x*cosa + u3x*sint + u4x)
    atcoor[1] = fpy + fpd*(u2y*cosa + u3y*sint + u4y)
    atcoor[2] = fpz + fpd*(u2z*cosa + u3z*sint + u4z)

    return atcoor

def get_atom_named(atoms, atom_name, start, end):
    for i in range(start, end):
        atom = atoms[i]
        if atom.name == atom_name:
            return atom

def get_scacceptor(atoms, start, end, res): 
    scatoms = []
    for arec in RESIDUES[res]:
        if arec[7] and arec[0] &lt;&gt; ' O  ':
            acp = (get_atom_named(atoms, arec[0], start, end))
            acp1 = acp.first_parent
            scatoms.append((acp,acp1))
    return scatoms

def make_hbtab_loos(p, hblist, ANGLE=angle,
             ABS=abs, TORSION=torsion):
    hbtab_loos = []
    for dlist, alist in hblist:
        d1, da1, da2 = dlist
        Close = d1.close

        for acp, acp1, hbd, hbdmax, hbtors, hbene in alist:
            # for sidechains relax orientation of donor
            if (d1.name == ' OG ' or d1.name == ' OG1' or d1.name == ' NZ ' or
                d1.name == ' OH ' or d1.name == ' NH1' or d1.name == ' NH2' or
                d1.name == ' ND1' or d1.name == ' ND2' or d1.name == ' OD2' or
                d1.name == ' SG ' or d1.name == ' NE2' or d1.name == ' OE2' or
                d1.name == ' NE1' or d1.name == ' OE1' or d1.name == ' OD1'):
                if Close(acp, hbdmax):
                    if ANGLE(d1, acp, acp1) &gt; 90.0:
                        dist = d1.distance(acp)
#                       print dist, d1.name, p.res_pdb_number[d1.resnum], acp.name, \
#                             p.res_pdb_number[acp.resnum]
                        hbtab_loos.append((d1, acp, dist))
            elif (acp.name == ' OG ' or acp.name == ' OG1' or 
                acp.name == ' OH ' or 
                acp.name == ' ND1' or acp.name == ' ND2' or acp.name == ' OD2' or
                acp.name == ' SG ' or acp.name == ' NE2' or acp.name == ' OE2' or
                acp.name == ' NE1' or acp.name == ' OE1' or acp.name == ' OD1'):
                if Close(acp, hbdmax):
                    if ANGLE(d1, acp, acp1) &gt; 90.0:
                        dist = d1.distance(acp)   
#                       print dist, d1.name, p.res_pdb_number[d1.resnum], acp.name, \
#                             p.res_pdb_number[acp.resnum]
                        hbtab_loos.append((d1, acp, dist))
            else:
                if Close(acp, hbdmax):
                    # first three keep oxygen in approximate +/- 70 deg cone of NH
                    # (it's actually a square)
                    # Equivalent to theta angle &gt;= 110 in Kortemme et al. JMB, 2003
                    # Last angle keeps N in +/- 90 deg cone of CO
                    # Equivalent to psi angle &gt;= 90 in Kortemme
                    if ABS(TORSION(acp, d1, da1, da2)) &gt; 130.0 and \
                    ANGLE(da1, d1, acp) &gt; 69.0 and \
                    ANGLE(da2, d1, acp) &gt; 69.0 and \
                    ANGLE(d1, acp, acp1) &gt; 90.0:
                        dist = d1.distance(acp)
#                       print dist, d1.name, p.res_pdb_number[d1.resnum], acp.name, \
#                             p.res_pdb_number[acp.resnum]
                        hbtab_loos.append((d1, acp, dist))

    return hbtab_loos
def _ishbonded_loos(p,atom, hbtab_loos,numint_loos):
    """
    Tests for hydrogen bond between atom and other internal atom
    """
    rn = p.residue_names
    for rec in hbtab_loos:
        d1, acp, dist = rec
        if (atom.resnum==d1.resnum and \
            rn[atom.resnum]==rn[d1.resnum] and \
            atom.name==d1.name):
            numint_loos = numint_loos + 1
            return 1,numint_loos
        elif (atom.resnum==acp.resnum and \
            rn[atom.resnum]==rn[acp.resnum] and \
            atom.name==acp.name):
            numint_loos = numint_loos + 1
            return 1,numint_loos
    return 0,numint_loos


def _close(coords1, coords2, hsd):
    """ Local version that deals with coordinates, not atom objects
    """
    x1,y1,z1 = coords1
    x2,y2,z2 = coords2
    dist = sqrt((x1 - x2)**2 + (y1 - y2)**2 + (z1 - z2)**2)
    if (dist &lt; hsd):
        return 1

    return 0

def mk_o2_solv(prot,atoms,start,end,acp,sp,tp,numvirt,
               tryall=0, watpdb=None, ext_rad=1.25):
    wat_coord = []
    #define a new single list zmatrix index
    ozmindx = [0,0,0]

    #use the O as the first parent     
    ozmindx[0] = acp      

    #use the C (same res) as the second parent
    ozmindx[1] = sp

    #use the CA (same res) as the third parent
    ozmindx[2] = tp

    #coords of the first parent, second and third    
    coordsfp = acp.coords()
    coordssp = sp.coords()
    coordstp = tp.coords()

    #define the distance, angle, and torsion for the 180 deg "on top"
    #sol1
    ozmvalu = [2.95, 180.0, 0.0]

    owat2_crds = []
    owat2_crds = ztox(coordsfp, coordssp, coordstp, ozmvalu, ozmindx)

    # check for bumps of virt water with other atoms
    isbump = 0
    gotone = 0
    scale = 0.9
    for atom in atoms:
        if (prot.residue_names[atom.resnum] != 'NME') and (prot.residue_names[atom.resnum] != 'ACE'): 
            atm_crds = atom.coords()
            atm_rad = scale * atom.radius
            if atom.name == ' H  ':
                # NH radius is 0.9, and H--O dist is 1.95
                # and H won't bump its own water of hydration
                hsd = 1.94
            else:      
                hsd = (ext_rad + atm_rad)

            if _close(owat2_crds, atm_crds, hsd):
                #here only need the first instance of a bump
                isbump = 1
                break

    # if no bump write out virtual water
    if not isbump:
        numvirt = numvirt + 1
        xx2,yy2,zz2 = owat2_crds
        format16 = 'ATOM    900  O   HOH   921    %8.3f%8.3f%8.3f%6.2f%6.2f\n'
        if watpdb and (tryall or not gotone):
            watpdb.write(format16 % (xx2,yy2,zz2, 0.0, 50.0))

        if gotone == 0:
            wat_coord.append((xx2,yy2,zz2))
            gotone = 1
    elif watpdb and tryall:
        xx2,yy2,zz2 = owat2_crds
        format16 = 'ATOM    900  O   HOH   921    %8.3f%8.3f%8.3f%6.2f%6.2f\n'
        watpdb.write(format16 % (xx2,yy2,zz2, 0.0, 0.0))
        
    #If the 180 worked bail out and go to NH
#   if isbump &lt; 1 :
#       return numvirt

    #define the distance, angle, and torsion for the "lower side"
    #sol2
    ozmvalu = [2.95, 130.0, 90.0]

    owat1_crds = []
    owat1_crds = ztox(coordsfp, coordssp, coordstp, ozmvalu, ozmindx)

    # check for bumps of virt water with other atoms
    isbump = 0
    for atom in atoms:
        if (prot.residue_names[atom.resnum] != 'NME') and (prot.residue_names[atom.resnum] != 'ACE'): 
            atm_crds = atom.coords()
            atm_rad = scale * atom.radius
            if atom.name == ' H  ':
                # NH radius is 0.9, and H--O dist is 1.95        
                # and H won't bump its own water of hydration        
                hsd = 1.94
            else:        
                hsd = (ext_rad + atm_rad)        

            if _close(owat1_crds, atm_crds, hsd):
                #here only need the first instance of a bump
                isbump = 1
                break

    # if no bump write out virtual water
    if not isbump:
        numvirt = numvirt + 1
        xx1,yy1,zz1 = owat1_crds
        format15 = 'ATOM    900  O   HOH   922    %8.3f%8.3f%8.3f%6.2f%6.2f\n'
        if watpdb and (tryall or not gotone):
            watpdb.write(format15 % (xx1,yy1,zz1, 0.0, 50.0))
            
        if gotone == 0:
            wat_coord.append((xx1,yy1,zz1))
            gotone = 1
    elif watpdb and tryall:
        xx1,yy1,zz1 = owat1_crds
        format15 = 'ATOM    900  O   HOH   922    %8.3f%8.3f%8.3f%6.2f%6.2f\n'
        watpdb.write(format15 % (xx1,yy1,zz1, 0.0, 0.0))

    #If this worked bail out and go to NH
#   if isbump &lt; 1 :
#       return numvirt

    #define the distance, angle, and torsion for the "upper side"
    #sol3
    ozmvalu = [2.95, 130.0, -90.0]

    owat3_crds = []
    owat3_crds = ztox(coordsfp, coordssp, coordstp, ozmvalu, ozmindx)

    # check for bumps of virt water with other atoms
    isbump = 0
    for atom in atoms:
        if (prot.residue_names[atom.resnum] != 'NME') and (prot.residue_names[atom.resnum] != 'ACE'): 
            atm_crds = atom.coords()
            atm_rad = scale * atom.radius
            if atom.name == ' H  ':
                # NH radius is 0.9, and H--O dist is 1.95        
                # and H won't bump its own water of hydration        
                hsd = 1.94
            else:        
                hsd = (ext_rad + atm_rad)        

            if _close(owat3_crds, atm_crds, hsd):
                #here only need the first instance of a bump
                isbump = 1
                break

    # if no bump write out virtual water
    if not isbump:
        numvirt = numvirt + 1
        xx3,yy3,zz3 = owat3_crds
        format17 = 'ATOM    900  O   HOH   923    %8.3f%8.3f%8.3f%6.2f%6.2f\n'
        if watpdb and (tryall or not gotone):
            watpdb.write(format17 % (xx3,yy3,zz3, 0.0, 50.0))
        
        if gotone == 0:
            wat_coord.append((xx3,yy3,zz3))
            gotone = 1
    elif watpdb and tryall:
        xx3,yy3,zz3 = owat3_crds
        format17 = 'ATOM    900  O   HOH   923    %8.3f%8.3f%8.3f%6.2f%6.2f\n'
        watpdb.write(format17 % (xx3,yy3,zz3, 0.0, 0.0))


    #If this worked bail out and go to NH
#   if isbump &lt; 1 :
#       return numvirt

    #define the distance, angle, and torsion for the near lone pair
    #sol4
    ozmvalu = [2.95, 130.0, 0.0]

    owat4_crds = []
    owat4_crds = ztox(coordsfp, coordssp, coordstp, ozmvalu, ozmindx)

    # check for bumps of virt water with other atoms
    isbump = 0
    for atom in atoms:
        if (prot.residue_names[atom.resnum] != 'NME') and (prot.residue_names[atom.resnum] != 'ACE'): 
            atm_crds = atom.coords()
            atm_rad = scale * atom.radius
            if atom.name == ' H  ':
                # NH radius is 0.9, and H--O dist is 1.95        
                # and H won't bump its own water of hydration        
                hsd = 1.94
            else:        
                hsd = (ext_rad + atm_rad)        

            if _close(owat4_crds, atm_crds, hsd):
                #here only need the first instance of a bump
                isbump = 1
                break

    # if no bump write out virtual water
    if not isbump:
        numvirt = numvirt + 1
        xx4,yy4,zz4 = owat4_crds
        format17 = 'ATOM    900  O   HOH   924    %8.3f%8.3f%8.3f%6.2f%6.2f\n'
        if watpdb and (tryall or not gotone):
            watpdb.write(format17 % (xx4,yy4,zz4, 0.0, 50.0))

        if gotone == 0:
            wat_coord.append((xx4,yy4,zz4))
            gotone = 1
    elif watpdb and tryall:
        xx4,yy4,zz4 = owat4_crds
        format17 = 'ATOM    900  O   HOH   924    %8.3f%8.3f%8.3f%6.2f%6.2f\n'
        watpdb.write(format17 % (xx4,yy4,zz4, 0.0, 0.0))


    #If this worked bail out and go to NH
#   if isbump &lt; 1 :
#       return numvirt

    #define the distance, angle, and torsion for the far lone pair
    #sol5
    ozmvalu = [2.95, -130.0, 0.0]

    owat5_crds = []
    owat5_crds = ztox(coordsfp, coordssp, coordstp, ozmvalu, ozmindx)

    # check for bumps of virt water with other atoms
    isbump = 0
    for atom in atoms:
        if (prot.residue_names[atom.resnum] != 'NME') and (prot.residue_names[atom.resnum] != 'ACE'): 
            atm_crds = atom.coords()
            atm_rad = scale * atom.radius
            if atom.name == ' H  ':
                # NH radius is 0.9, and H--O dist is 1.95        
                # and H won't bump its own water of hydration        
                hsd = 1.94
            else:        
                hsd = (ext_rad + atm_rad)        

            if _close(owat5_crds, atm_crds, hsd):
                #here only need the first instance of a bump
                isbump = 1
                break

    # if no bump write out virtual water
    if not isbump:
        numvirt = numvirt + 1
        xx5,yy5,zz5 = owat5_crds
        format17 = 'ATOM    900  O   HOH   925    %8.3f%8.3f%8.3f%6.2f%6.2f\n'
        if watpdb and (tryall or not gotone):
            watpdb.write(format17 % (xx5,yy5,zz5, 0.0, 50.0))
        
        if gotone == 0:
            wat_coord.append((xx5,yy5,zz5))
            gotone = 1
    elif watpdb and tryall:
        xx5,yy5,zz5 = owat5_crds
        format17 = 'ATOM    900  O   HOH   925    %8.3f%8.3f%8.3f%6.2f%6.2f\n'
        watpdb.write(format17 % (xx5,yy5,zz5, 0.0, 0.0))
        
    return numvirt,wat_coord

def mk_nh_solv(prot,atoms,start,end,donor,sp,tp,numvirt,
               tryall=0, watpdb=None, ext_rad=1.25):
    wat_coord = []
    #define a new single list zmatrix index
    nzmindx = [0,0,0]        

    #use the N as the first parent
    nzmindx[0] = donor      

    #use the CA as the second parent
    nzmindx[1] = sp

    #use the C (prev res) as the third parent
    nzmindx[2] = tp

    #define the distance, angle, and torsion for 180 deg
    nzmvalu = [2.95, 119.0, 180.0]

    #coords of the first parent, second and third
    coordsfp = donor.coords()
    coordssp = sp.coords()
    coordstp = tp.coords()

    nwat_crds = []
    nwat_crds = ztox(coordsfp, coordssp, coordstp, nzmvalu, nzmindx)

    #check for bumps of N hbonded virtual O with all other atoms
    isbump = 0
    gotone = 0
    scale = 0.9
    for atom in atoms:
        if (prot.residue_names[atom.resnum] != 'NME') and (prot.residue_names[atom.resnum] != 'ACE'): 
            atm_crds = atom.coords()
            atm_rad = scale * atom.radius
            if atom.name == ' H  ':
                # NH radius is 0.9, and H--O dist is 1.95        
                # and H won't bump its own water of hydration        
                hsd = 1.94
            else:        
                hsd = (ext_rad + atm_rad)        

            if _close(nwat_crds, atm_crds, hsd):
                #here only need the first instance of a bump
                isbump = 1
                break

        # if no bump write out virtual water
    if not isbump:
        numvirt =  numvirt + 1
        xx,yy,zz = nwat_crds
        format14 = 'ATOM    900  O   HOH   910    %8.3f%8.3f%8.3f%6.2f%6.2f\n'
        if watpdb and (tryall or not gotone):
            watpdb.write(format14 % (xx,yy,zz, 0.0, 50.0))
            
        if gotone == 0:
            wat_coord.append((xx,yy,zz))
            gotone = 1
    elif tryall and watpdb:
        xx,yy,zz = nwat_crds
        format14 = 'ATOM    900  O   HOH   910    %8.3f%8.3f%8.3f%6.2f%6.2f\n'
        watpdb.write(format14 % (xx,yy,zz, 0.0, 0.0))

#       return numvirt

    #define the distance, angle, and torsion for CA side
    nzmvalu = [2.95, 109.0, 180.0]

    #coords of the first parent, second and third
    coordsfp = donor.coords()
    coordssp = sp.coords()
    coordstp = tp.coords()

    nwat_crds = []
    nwat_crds = ztox(coordsfp, coordssp, coordstp, nzmvalu, nzmindx)

    #check for bumps of N hbonded virtual O with all other atoms
    isbump = 0
    for atom in atoms:
        if (prot.residue_names[atom.resnum] != 'NME') and (prot.residue_names[atom.resnum] != 'ACE'):
            atm_crds = atom.coords()
            atm_rad = scale * atom.radius
            if atom.name == ' H  ':
                # NH radius is 0.9, and H--O dist is 1.95
                # and H won't bump its own water of hydration
                hsd = 1.94
            else:
                hsd = (ext_rad + atm_rad)

            if _close(nwat_crds, atm_crds, hsd):
                #here only need the first instance of a bump
                isbump = 1
                break

        # if no bump write out virtual water
    if not isbump:
        numvirt =  numvirt + 1
        xx,yy,zz = nwat_crds
        format14 = 'ATOM    900  O   HOH   910    %8.3f%8.3f%8.3f%6.2f%6.2f\n'
        if watpdb and (tryall or not gotone):
            watpdb.write(format14 % (xx,yy,zz, 0.0, 50.0))

        if gotone == 0:
            wat_coord.append((xx,yy,zz))
            gotone = 1
    elif tryall and watpdb:
        xx,yy,zz = nwat_crds
        format14 = 'ATOM    900  O   HOH   910    %8.3f%8.3f%8.3f%6.2f%6.2f\n'
        watpdb.write(format14 % (xx,yy,zz, 0.0, 0.0))
#       return numvirt

    #define the distance, angle, and torsion for C side
    nzmvalu = [2.95, 129.0, 180.0]

    #coords of the first parent, second and third
    coordsfp = donor.coords()
    coordssp = sp.coords()
    coordstp = tp.coords()

    nwat_crds = []
    nwat_crds = ztox(coordsfp, coordssp, coordstp, nzmvalu, nzmindx)

    #check for bumps of N hbonded virtual O with all other atoms
    isbump = 0
    for atom in atoms:
        if (prot.residue_names[atom.resnum] != 'NME') and (prot.residue_names[atom.resnum] != 'ACE'):
            atm_crds = atom.coords()
            atm_rad = scale * atom.radius
            if atom.name == ' H  ':
                # NH radius is 0.9, and H--O dist is 1.95
                # and H won't bump its own water of hydration
                hsd = 1.94
            else:
                hsd = (ext_rad + atm_rad)

            if _close(nwat_crds, atm_crds, hsd):
                #here only need the first instance of a bump
                isbump = 1
                break

        # if no bump write out virtual water
    if not isbump:
        numvirt =  numvirt + 1
        xx,yy,zz = nwat_crds
        format14 = 'ATOM    900  O   HOH   910    %8.3f%8.3f%8.3f%6.2f%6.2f\n'
        if watpdb and (tryall or not gotone):
            watpdb.write(format14 % (xx,yy,zz, 0.0, 50.0))

        if gotone == 0:
            wat_coord.append((xx,yy,zz))
            gotone = 1
    elif watpdb and tryall:
        xx,yy,zz = nwat_crds
        format14 = 'ATOM    900  O   HOH   910    %8.3f%8.3f%8.3f%6.2f%6.2f\n'
        watpdb.write(format14 % (xx,yy,zz, 0.0, 0.0))
        
#       return numvirt

    #define the distance, angle, and torsion for right side
    nzmvalu = [2.95, 118.5, 170.0]

    #coords of the first parent, second and third
    coordsfp = donor.coords()
    coordssp = sp.coords()
    coordstp = tp.coords()

    nwat_crds = []
    nwat_crds = ztox(coordsfp, coordssp, coordstp, nzmvalu, nzmindx)

    #check for bumps of N hbonded virtual O with all other atoms
    isbump = 0
    for atom in atoms:
        if (prot.residue_names[atom.resnum] != 'NME') and (prot.residue_names[atom.resnum] != 'ACE'):
            atm_crds = atom.coords()
            atm_rad = scale * atom.radius
            if atom.name == ' H  ':
                # NH radius is 0.9, and H--O dist is 1.95
                # and H won't bump its own water of hydration
                hsd = 1.94
            else:
                hsd = (ext_rad + atm_rad)

            if _close(nwat_crds, atm_crds, hsd):
                #here only need the first instance of a bump
                isbump = 1
                break

        # if no bump write out virtual water
    if not isbump:
        numvirt =  numvirt + 1
        xx,yy,zz = nwat_crds
        format14 = 'ATOM    900  O   HOH   910    %8.3f%8.3f%8.3f%6.2f%6.2f\n'
        if watpdb and (tryall or not gotone):
            watpdb.write(format14 % (xx,yy,zz, 0.0, 50.0))

        if gotone == 0:
            wat_coord.append((xx,yy,zz))
            gotone = 1
    elif watpdb and tryall:
        xx,yy,zz = nwat_crds
        format14 = 'ATOM    900  O   HOH   910    %8.3f%8.3f%8.3f%6.2f%6.2f\n'
        watpdb.write(format14 % (xx,yy,zz, 0.0, 0.0))
#       return numvirt

    #define the distance, angle, and torsion for left side
    nzmvalu = [2.95, 118.5, -170.0]

    #coords of the first parent, second and third
    coordsfp = donor.coords()
    coordssp = sp.coords()
    coordstp = tp.coords()

    nwat_crds = []
    nwat_crds = ztox(coordsfp, coordssp, coordstp, nzmvalu, nzmindx)

    #check for bumps of N hbonded virtual O with all other atoms
    isbump = 0
    for atom in atoms:
        if (prot.residue_names[atom.resnum] != 'NME') and (prot.residue_names[atom.resnum] != 'ACE'):
            atm_crds = atom.coords()
            atm_rad = scale * atom.radius
            if atom.name == ' H  ':
                # NH radius is 0.9, and H--O dist is 1.95
                # and H won't bump its own water of hydration
                hsd = 1.94
            else:
                hsd = (ext_rad + atm_rad)

            if _close(nwat_crds, atm_crds, hsd):
                #here only need the first instance of a bump
                isbump = 1
                break

        # if no bump write out virtual water
    if not isbump:
        numvirt =  numvirt + 1
        xx,yy,zz = nwat_crds
        format14 = 'ATOM    900  O   HOH   910    %8.3f%8.3f%8.3f%6.2f%6.2f\n'
        if watpdb and (tryall or not gotone):
            watpdb.write(format14 % (xx,yy,zz, 0.0, 50.0))

        if gotone == 0:
            wat_coord.append((xx,yy,zz))
            gotone = 1
    elif watpdb and tryall:
        xx,yy,zz = nwat_crds
        format14 = 'ATOM    900  O   HOH   910    %8.3f%8.3f%8.3f%6.2f%6.2f\n'
        watpdb.write(format14 % (xx,yy,zz, 0.0, 0.0))
        
#       return numvirt

    return numvirt,wat_coord

def mk_amid_solv(prot,atoms,start,end,donor,sp,tp,numvirt,
                 tryall=0,watpdb=None, ext_rad=1.25):
    wat_coord = []
    #define a new single list zmatrix index
    nzmindx = [0,0,0]

    #use the N as the first parent
    nzmindx[0] = donor

    #use the C (prev res) as the second parent
    nzmindx[1] = sp

    #use the O (prev res) as the third parent          
    nzmindx[2] = tp

    #define the distance, angle, and torsion
    nzmvalu = [2.95, 180.0, 180.0]

    #coords of the first parent, second and third
    coordsfp = donor.coords()
    coordssp = sp.coords()
    coordstp = tp.coords() 
   
    nwat_crds = []
    nwat_crds = ztox(coordsfp, coordssp, coordstp, nzmvalu, nzmindx)

    #check for bumps of N hbonded virtual O with all other atoms
    isbump = 0
    scale = 0.9
    for atom in atoms:
        if (prot.residue_names[atom.resnum] != 'NME') and (prot.residue_names[atom.resnum] != 'ACE'):
            atm_crds = atom.coords()
            atm_rad = scale * atom.radius
            if atom.name == ' H  ':
                # NH radius is 0.9, and H--O dist is 1.95        
                # and H won't bump its own water of hydration        
                hsd = 1.94
            else:        
                hsd = (ext_rad + atm_rad)        

            if _close(nwat_crds, atm_crds, hsd):
                #here only need the first instance of a bump
                isbump = 1
                break
        # if no bump write out virtual water
    if not isbump:
        numvirt =  numvirt + 1
        xx,yy,zz = nwat_crds
        format14 = 'ATOM    900  O   HOH   911    %8.3f%8.3f%8.3f%6.2f%6.2f\n'
        if watpdb: watpdb.write(format14 % (xx,yy,zz, 0.0, 50.0))
        wat_coord.append((xx,yy,zz))
    elif tryall and watpdb:
        xx,yy,zz = nwat_crds
        format14 = 'ATOM    900  O   HOH   911    %8.3f%8.3f%8.3f%6.2f%6.2f\n'
        watpdb.write(format14 % (xx,yy,zz, 0.0, 0.0))
        

    return numvirt,wat_coord

def mk_oh_solv(prot,atoms,start,end,acp,sp,tp,numvirt,tryall=0,watpdb=None,ext_rad=1.25):
    wat_coord = []
    #define a new single list zmatrix index
    ozmindx = [0,0,0]

    #use the O as the first parent 
    ozmindx[0] = acp

    #use the C (same res) as the second parent
    ozmindx[1] = sp

    #use the CA (same res) as the third parent
    ozmindx[2] = tp

    #coords of the first parent, second and third
    coordsfp = acp.coords()
    coordssp = sp.coords()
    coordstp = tp.coords()

    #define the distance, angle, and torsion for the "lower side"
    #sol2
    ozmvalu = [2.95, 130.0, 90.0]

    owat1_crds = []
    owat1_crds = ztox(coordsfp, coordssp, coordstp, ozmvalu, ozmindx)

    # check for bumps of virt water with other atoms
    isbump = 0
    scale = 0.9
    for atom in atoms:
        if (prot.residue_names[atom.resnum] != 'NME') and (prot.residue_names[atom.resnum] != 'ACE'): 
            atm_crds = atom.coords()
            atm_rad = scale * atom.radius
            if atom.name == ' H  ':
                # NH radius is 0.9, and H--O dist is 1.95        
                # and H won't bump its own water of hydration        
                hsd = 1.94
            else:        
                hsd = (ext_rad + atm_rad)        

            if _close(owat1_crds, atm_crds, hsd):
                #here only need the first instance of a bump
                isbump = 1
                break

    # if no bump write out virtual water
    if not isbump:
        numvirt = numvirt + 1
        xx1,yy1,zz1 = owat1_crds
        format15 = 'ATOM    900  O   HOH   922    %8.3f%8.3f%8.3f%6.2f%6.2f\n'
        if watpdb: watpdb.write(format15 % (xx1,yy1,zz1, 0.0, 50.0))
        wat_coord.append((xx1,yy1,zz1))
    elif tryall and watpdb:
        xx1,yy1,zz1 = owat1_crds
        format15 = 'ATOM    900  O   HOH   922    %8.3f%8.3f%8.3f%6.2f%6.2f\n'
        watpdb.write(format15 % (xx1,yy1,zz1, 0.0, 0.0))
        
    #If this worked bail out and go to NH
    if isbump &lt; 1 and not tryall:
        return numvirt,wat_coord

    #define the distance, angle, and torsion for the "upper side"
    #sol3
    ozmvalu = [2.95, 130.0, -90.0]

    owat3_crds = []
    owat3_crds = ztox(coordsfp, coordssp, coordstp, ozmvalu, ozmindx)

    # check for bumps of virt water with other atoms
    isbump = 0
    for atom in atoms:
        if (prot.residue_names[atom.resnum] != 'NME') and (prot.residue_names[atom.resnum] != 'ACE'): 
            atm_crds = atom.coords()
            atm_rad = scale * atom.radius
            if atom.name == ' H  ':
                # NH radius is 0.9, and H--O dist is 1.95        
                # and H won't bump its own water of hydration        
                hsd = 1.94
            else:        
                hsd = (ext_rad + atm_rad)        

            if _close(owat3_crds, atm_crds, hsd):
                #here only need the first instance of a bump
                isbump = 1
                break

    # if no bump write out virtual water   
    if not isbump:  
        numvirt = numvirt + 1
        xx3,yy3,zz3 = owat3_crds
        format17 = 'ATOM    900  O   HOH   923    %8.3f%8.3f%8.3f%6.2f%6.2f\n'
        if watpdb: watpdb.write(format17 % (xx3,yy3,zz3, 0.0, 50.0))
        wat_coord.append((xx3,yy3,zz3))
    elif tryall and watpdb:
        xx3,yy3,zz3 = owat3_crds
        format17 = 'ATOM    900  O   HOH   923    %8.3f%8.3f%8.3f%6.2f%6.2f\n'
        watpdb.write(format17 % (xx3,yy3,zz3, 0.0, 0.0))
 
    #If this worked bail out and go to NH
    if isbump &lt; 1 and not tryall:
        return numvirt,wat_coord

    #define the distance, angle, and torsion for the near lone pair
    #sol4
    ozmvalu = [2.95, 130.0, 0.0]

    owat4_crds = []
    owat4_crds = ztox(coordsfp, coordssp, coordstp, ozmvalu, ozmindx)

    # check for bumps of virt water with other atoms
    isbump = 0
    for atom in atoms:
        if (prot.residue_names[atom.resnum] != 'NME') and (prot.residue_names[atom.resnum] != 'ACE'): 
            atm_crds = atom.coords()
            atm_rad = scale * atom.radius
            if atom.name == ' H  ':
                # NH radius is 0.9, and H--O dist is 1.95        
                # and H won't bump its own water of hydration        
                hsd = 1.94
            else:        
                hsd = (ext_rad + atm_rad)        

            if _close(owat4_crds, atm_crds, hsd):
                #here only need the first instance of a bump
                isbump = 1
                break

    # if no bump write out virtual water   
    if not isbump:  
        numvirt = numvirt + 1
        xx4,yy4,zz4 = owat4_crds
        format17 = 'ATOM    900  O   HOH   924    %8.3f%8.3f%8.3f%6.2f%6.2f\n'
        if watpdb: watpdb.write(format17 % (xx4,yy4,zz4, 0.0, 50.0))
        wat_coord.append((xx4,yy4,zz4))
    elif tryall and watpdb:
        xx4,yy4,zz4 = owat4_crds
        format17 = 'ATOM    900  O   HOH   924    %8.3f%8.3f%8.3f%6.2f%6.2f\n'
        watpdb.write(format17 % (xx4,yy4,zz4, 0.0, 0.0))
        
    #If this worked bail out and go to NH
    if isbump &lt; 1 and not tryall :
        return numvirt,wat_coord

    #define the distance, angle, and torsion for the far lone pair
    #sol5
    ozmvalu = [2.95, -130.0, 0.0]

    owat5_crds = []
    owat5_crds = ztox(coordsfp, coordssp, coordstp, ozmvalu, ozmindx)

    # check for bumps of virt water with other atoms
    isbump = 0
    for atom in atoms:
        if (prot.residue_names[atom.resnum] != 'NME') and (prot.residue_names[atom.resnum] != 'ACE'): 
            atm_crds = atom.coords()
            atm_rad = scale * atom.radius
            if atom.name == ' H  ':
                # NH radius is 0.9, and H--O dist is 1.95        
                # and H won't bump its own water of hydration        
                hsd = 1.94
            else:        
                hsd = (ext_rad + atm_rad)        

            if _close(owat5_crds, atm_crds, hsd):
                #here only need the first instance of a bump
                isbump = 1
                break

    # if no bump write out virtual water   
    if not isbump:  
        numvirt = numvirt + 1
        xx5,yy5,zz5 = owat5_crds
        format17 = 'ATOM    900  O   HOH   925    %8.3f%8.3f%8.3f%6.2f%6.2f\n'
        if watpdb: watpdb.write(format17 % (xx5,yy5,zz5, 0.0, 50.0))
        wat_coord.append((xx5,yy5,zz5))
    elif tryall and watpdb:
        xx5,yy5,zz5 = owat5_crds
        format17 = 'ATOM    900  O   HOH   925    %8.3f%8.3f%8.3f%6.2f%6.2f\n'
        watpdb.write(format17 % (xx5,yy5,zz5, 0.0, 0.0))

    return numvirt,wat_coord

def mk_virt_bb_cntmult_loosHbd(prot,hblist,tryall=0,watpdb=None,ext_rad=1.25):
    """
    Places oxygen in hydrogen bonding orientation to all donors and acceptors.
    Goes thru polypeptide and calls,
		mk_o2_solv, mk_nh_solv, mk_oh_solv, mk_amid_solv
    as appropriate.
    """

    numres = prot.num_residues
    atoms = prot.atoms
    rn = prot.residue_names
    fai = prot.residue_first_atom_indices

    minres,maxres = get_res_extents(prot)

    hbtab_loos = make_hbtab_loos(prot, hblist)

    wat_list = []
    cntmult_list = []
    for i in range(minres,maxres+1):
        cntmult_list.append(([i,1,0,0],[i,2,0,0]))
    numbb = 0
    numposs = 0
    bbtot = 0
    numint_loos=0
    numint_sc_loos=0
    numvirt_loos=0
    numvirt_sc_loos=0
    #backbone first
    for i in range(minres,maxres):
        start = fai[i]
        end = fai[i+1]
        acp = get_atom_named(atoms, ' O  ', start, end)
        sp = get_atom_named(atoms, ' C  ', start, end)
        tp = get_atom_named(atoms, ' CA ', start, end )

        bbtot = bbtot + 1
        numposs = numposs + 1

        hbond, numint_loos = _ishbonded_loos(prot, acp, hbtab_loos,numint_loos)

        old_numvirt_loos = numvirt_loos
        numvirt_loos,wcoord = mk_o2_solv(prot,atoms,start,end,acp,sp,tp, \
             numvirt_loos,tryall=tryall, watpdb=watpdb,ext_rad=ext_rad)
        if len(wcoord) &gt; 0:
            wat_list.append(wcoord[0])

        if hbond:
            numbb = numbb + 1
            if numvirt_loos &gt; old_numvirt_loos:
                # if hbonded and solvent accessible just add one water
                cntmult_list[i][1][2] = 1
                # but flag this water as special for linusScore
                cntmult_list[i][1][3] = 1

        if not hbond:
            if numvirt_loos &gt; old_numvirt_loos:
                numbb = numbb + 1
                num_solv = numvirt_loos - old_numvirt_loos
                cntmult_list[i][1][2] = num_solv
            else:
                cntmult_list[i][1][2] = -1

#               print prot.res_pdb_number[acp.resnum], prot.residue_names[acp.resnum] \
#                    , acp.name
 
    #construct a virtual water along the NH axis at 2.95 from the N
    for i in range(minres,maxres):
        if rn[i] &lt;&gt; 'PRO':
            if not ((i == 0) and (minres == 0)):
                start = fai[i]
                end = fai[i+1]
                prev = fai[i-1]      
                donor = get_atom_named(atoms, ' N  ', start, end)
                sp = get_atom_named(atoms, ' CA ', start, end)
                tp = get_atom_named(atoms, ' C  ', prev, start)
    
                bbtot = bbtot + 1
                numposs = numposs + 1
    
                hbond, numint_loos = _ishbonded_loos(prot, donor, hbtab_loos,numint_loos)
                if hbond:
                    numbb = numbb + 1
    
                if not hbond:
                    old_numvirt_loos = numvirt_loos
                    numvirt_loos,wcoord = mk_nh_solv(prot,atoms,start,end,donor,sp,tp, \
                     numvirt_loos,tryall=tryall, watpdb=watpdb,ext_rad=ext_rad)
                    if len(wcoord) &gt; 0:
                        wat_list.append(wcoord[0])
                    if numvirt_loos &gt; old_numvirt_loos:
                        numbb = numbb + 1
                        num_solv = numvirt_loos - old_numvirt_loos
                        cntmult_list[i][0][2] = num_solv
    
                    else:
                        cntmult_list[i][0][2] = -1
#                   print prot.res_pdb_number[donor.resnum], prot.residue_names[donor.resnum] \
#                     , donor.name

    # now for sidechains
    for i in range(minres,maxres):
        if rn[i] == 'ARG':
            start = fai[i]
            end = fai[i+1]

            #construct virtual waters at 2.95 from the N
            acp = get_atom_named(atoms, ' NH1', start, end)
            sp = get_atom_named(atoms, ' CZ ', start, end)
            tp = get_atom_named(atoms, ' NE ', start, end )

            numposs = numposs + 1
            hbond, numint_sc_loos = _ishbonded_loos(prot, acp, hbtab_loos,numint_sc_loos)

            if not hbond:
                numvirt_sc_loos,wcoord = mk_oh_solv(prot,atoms,start,end,acp,sp,tp, \
                   numvirt_sc_loos, tryall=tryall, watpdb=watpdb,ext_rad=ext_rad)
                if len(wcoord) &gt; 0:
                    wat_list.append(wcoord[0])

            #construct virtual waters at 2.95 from the N
            acp = get_atom_named(atoms, ' NH2', start, end)
            sp = get_atom_named(atoms, ' CZ ', start, end)
            tp = get_atom_named(atoms, ' NE ', start, end )

            numposs = numposs + 1
            hbond, numint_sc_loos = _ishbonded_loos(prot, acp, hbtab_loos,numint_sc_loos)

            # try all 4 positions
            if not hbond:
                numvirt_sc_loos,wcoord = mk_oh_solv(prot,atoms,start,end,acp,sp,tp, \
                    numvirt_sc_loos, tryall=tryall, watpdb=watpdb,ext_rad=ext_rad)
                if len(wcoord) &gt; 0:
                    wat_list.append(wcoord[0])

        if rn[i] == 'ASN':
            start = fai[i]
            end = fai[i+1]

            #construct virtual waters at 2.95 from the O
            # one along the C=O axis
            acp = get_atom_named(atoms, ' OD1', start, end)
            sp = get_atom_named(atoms, ' CG ', start, end)
            tp = get_atom_named(atoms, ' CB ', start, end )

            numposs = numposs + 1
            hbond, numint_sc_loos = _ishbonded_loos(prot, acp, hbtab_loos,numint_sc_loos)
            # if C=O not hbonded try all 5 positions
            if not hbond:
                numvirt_sc_loos,wcoord = mk_o2_solv(prot,atoms,start,end,acp,sp,tp, \
                   numvirt_sc_loos, tryall=tryall, watpdb=watpdb, ext_rad=ext_rad)
                if len(wcoord) &gt; 0:
                    wat_list.append(wcoord[0])

            donor = get_atom_named(atoms, ' ND2', start, end)
            sp = get_atom_named(atoms, ' CG ', start, end)
            tp = get_atom_named(atoms, ' CB ', start, end )

            numposs = numposs + 1
            hbond, numint_sc_loos = _ishbonded_loos(prot, donor, hbtab_loos,numint_sc_loos)

            if not hbond:
                numvirt_sc_loos,wcoord = mk_amid_solv(prot,atoms,start,end,donor,sp,tp, \
                    numvirt_sc_loos,tryall=tryall, watpdb=watpdb,ext_rad=ext_rad)
                if len(wcoord) &gt; 0:
                    wat_list.append(wcoord[0])

        if rn[i] == 'ASP':   
            start = fai[i] 
            end = fai[i+1]

            #construct virtual waters at 2.95 from the O
            # one along the C=O axis
            acp = get_atom_named(atoms, ' OD1', start, end)
            sp = get_atom_named(atoms, ' CG ', start, end)
            tp = get_atom_named(atoms, ' CB ', start, end )

            numposs = numposs + 1
            hbond, numint_sc_loos = _ishbonded_loos(prot, acp, hbtab_loos,numint_sc_loos)

            # if C=O not hbonded try all 5 positions
            if not hbond:
                numvirt_sc_loos,wcoord = mk_o2_solv(prot,atoms,start,end,acp,sp,tp, \
                    numvirt_sc_loos,tryall=tryall, watpdb=watpdb, ext_rad=ext_rad)
                if len(wcoord) &gt; 0:
                    wat_list.append(wcoord[0])

            #construct virtual waters at 2.95 from the O
            # one along the C=O axis
            acp = get_atom_named(atoms, ' OD2', start, end)
            sp = get_atom_named(atoms, ' CG ', start, end)
            tp = get_atom_named(atoms, ' CB ', start, end )

            numposs = numposs + 1
            hbond, numint_sc_loos = _ishbonded_loos(prot, acp, hbtab_loos,numint_sc_loos)

            # if C=O not hbonded try all 5 positions
            if not hbond:
                numvirt_sc_loos,wcoord = mk_o2_solv(prot,atoms,start,end,acp,sp,tp, \
                    numvirt_sc_loos, tryall=tryall, watpdb=watpdb, ext_rad=ext_rad)
                if len(wcoord) &gt; 0:
                    wat_list.append(wcoord[0])

        if rn[i] == 'GLN':   
            start = fai[i] 
            end = fai[i+1]

            #construct virtual waters at 2.95 from the O
            # one along the C=O axis
            acp = get_atom_named(atoms, ' OE1', start, end)
            sp = get_atom_named(atoms, ' CD ', start, end)
            tp = get_atom_named(atoms, ' CG ', start, end )

            numposs = numposs + 1
            hbond, numint_sc_loos = _ishbonded_loos(prot, acp, hbtab_loos,numint_sc_loos)

            # if C=O not hbonded try all 5 positions
            if not hbond:
                numvirt_sc_loos,wcoord = mk_o2_solv(prot,atoms,start,end,acp,sp,tp, \
                   numvirt_sc_loos, tryall=tryall, watpdb=watpdb, ext_rad=ext_rad)
                if len(wcoord) &gt; 0:
                    wat_list.append(wcoord[0])

            donor = get_atom_named(atoms, ' NE2', start, end)
            sp = get_atom_named(atoms, ' CD ', start, end)
            tp = get_atom_named(atoms, ' CG ', start, end )

            numposs = numposs + 1
            hbond, numint_sc_loos = _ishbonded_loos(prot, donor, hbtab_loos,numint_sc_loos)

            if not hbond:
                numvirt_sc_loos,wcoord = mk_amid_solv(prot,atoms,start,end,donor,sp,tp, \
                      numvirt_sc_loos, tryall=tryall, watpdb=watpdb,ext_rad=ext_rad)
                if len(wcoord) &gt; 0:
                    wat_list.append(wcoord[0])

        if rn[i] == 'GLU':
            start = fai[i]
            end = fai[i+1]

            #construct virtual waters at 2.95 from the O
            # one along the C=O axis
            acp = get_atom_named(atoms, ' OE1', start, end)
            sp = get_atom_named(atoms, ' CD ', start, end)
            tp = get_atom_named(atoms, ' CG ', start, end )

            numposs = numposs + 1
            hbond, numint_sc_loos = _ishbonded_loos(prot, acp, hbtab_loos,numint_sc_loos)

            # if C=O not hbonded try all 5 positions
            if not hbond:
                numvirt_sc_loos,wcoord = mk_o2_solv(prot,atoms,start,end,acp,sp,tp, \
                   numvirt_sc_loos, tryall=tryall, watpdb=watpdb, ext_rad=ext_rad)
                if len(wcoord) &gt; 0:
                    wat_list.append(wcoord[0])

            #construct virtual waters at 2.95 from the O
            # one along the C=O axis
            acp = get_atom_named(atoms, ' OE2', start, end)
            sp = get_atom_named(atoms, ' CD ', start, end)
            tp = get_atom_named(atoms, ' CG ', start, end )

            numposs = numposs + 1
            hbond, numint_sc_loos = _ishbonded_loos(prot, acp, hbtab_loos,numint_sc_loos)

            # if C=O not hbonded try all 5 positions
            if not hbond:
                numvirt_sc_loos,wcoord = mk_o2_solv(prot,atoms,start,end,acp,sp,tp, \
                  numvirt_sc_loos, tryall=tryall, watpdb=watpdb, ext_rad=ext_rad)
                if len(wcoord) &gt; 0:
                    wat_list.append(wcoord[0])

        if rn[i] == 'HIS':  
            start = fai[i]
            end = fai[i+1]

            donor = get_atom_named(atoms, ' NE2', start, end)
            sp = get_atom_named(atoms, ' CD2', start, end)
            tp = get_atom_named(atoms, ' CG ', start, end )

            numposs = numposs + 1
            hbond, numint_sc_loos = _ishbonded_loos(prot, donor, hbtab_loos,numint_sc_loos)
            if not hbond:
                numvirt_sc_loos,wcoord = mk_nh_solv(prot,atoms,start,end,donor,sp,tp, \
                    numvirt_sc_loos, tryall=tryall, watpdb=watpdb,ext_rad=ext_rad)
                if len(wcoord) &gt; 0:
                    wat_list.append(wcoord[0])

            donor = get_atom_named(atoms, ' ND1', start, end)
            sp = get_atom_named(atoms, ' CG ', start, end)
            tp = get_atom_named(atoms, ' CB ', start, end )

            numposs = numposs + 1
            hbond, numint_sc_loos = _ishbonded_loos(prot, donor, hbtab_loos,numint_sc_loos)
            if not hbond:
                numvirt_sc_loos,wcoord = mk_nh_solv(prot,atoms,start,end,donor,sp,tp, \
                    numvirt_sc_loos, tryall=tryall, watpdb=watpdb,ext_rad=ext_rad)
                if len(wcoord) &gt; 0:
                    wat_list.append(wcoord[0])

        if rn[i] == 'LYS':
            start = fai[i]
            end = fai[i+1]

            #construct virtual waters at 2.95 from the N
            acp = get_atom_named(atoms, ' NZ ', start, end)
            sp = get_atom_named(atoms, ' CE ', start, end)
            tp = get_atom_named(atoms, ' CD ', start, end )

            numposs = numposs + 1
            hbond, numint_sc_loos = _ishbonded_loos(prot, acp, hbtab_loos,numint_sc_loos)
            # try all 4 positions
            if not hbond:
                numvirt_sc_loos,wcoord = mk_oh_solv(prot,atoms,start,end,acp,sp,tp, \
                    numvirt_sc_loos, tryall=tryall, watpdb=watpdb,ext_rad=ext_rad)
                if len(wcoord) &gt; 0:
                    wat_list.append(wcoord[0])

        if rn[i] == 'SER':
            start = fai[i]
            end = fai[i+1]

            #construct virtual waters at 2.95 from the O
            # one along the C=O axis
            acp = get_atom_named(atoms, ' OG ', start, end)
            sp = get_atom_named(atoms, ' CB ', start, end)
            tp = get_atom_named(atoms, ' CA ', start, end )

            numposs = numposs + 1
            hbond, numint_sc_loos = _ishbonded_loos(prot, acp, hbtab_loos,numint_sc_loos)
            # try all 4 positions
            if not hbond:
                numvirt_sc_loos,wcoord = mk_oh_solv(prot,atoms,start,end,acp,sp,tp, \
                   numvirt_sc_loos, tryall=tryall, watpdb=watpdb,ext_rad=ext_rad)
                if len(wcoord) &gt; 0:
                    wat_list.append(wcoord[0])

        if rn[i] == 'THR':
            start = fai[i]
            end = fai[i+1]

            #construct virtual waters at 2.95 from the O
            # one along the C=O axis
            acp = get_atom_named(atoms, ' OG1', start, end)
            sp = get_atom_named(atoms, ' CB ', start, end)
            tp = get_atom_named(atoms, ' CA ', start, end )

            numposs = numposs + 1
            hbond, numint_sc_loos = _ishbonded_loos(prot, acp, hbtab_loos,numint_sc_loos)
            # try all 4 positions
            if not hbond:
                numvirt_sc_loos,wcoord = mk_oh_solv(prot,atoms,start,end,acp,sp,tp, \
                     numvirt_sc_loos, tryall=tryall, watpdb=watpdb,ext_rad=ext_rad)
                if len(wcoord) &gt; 0:
                    wat_list.append(wcoord[0])

        if rn[i] == 'TYR':
            start = fai[i]
            end = fai[i+1]

            #construct virtual waters at 2.95 from the O
            # one along the C=O axis
            acp = get_atom_named(atoms, ' OH ', start, end)
            sp = get_atom_named(atoms, ' CZ ', start, end)
            tp = get_atom_named(atoms, ' CE1', start, end )

            numposs = numposs + 1
            hbond, numint_sc_loos = _ishbonded_loos(prot, acp, hbtab_loos,numint_sc_loos)
            # try all 4 positions
            if not hbond:
                numvirt_sc_loos,wcoord = mk_oh_solv(prot,atoms,start,end,acp,sp,tp, \
                    numvirt_sc_loos, tryall=tryall, watpdb=watpdb,ext_rad=ext_rad)
                if len(wcoord) &gt; 0:
                    wat_list.append(wcoord[0])

    total_loos = (numint_loos) + numvirt_loos
    bb_nonbd_loos = bbtot - numbb
    return numint_loos, numvirt_loos, bbtot, numbb,wat_list,cntmult_list
# Given a unit sphere, the array below is the coordinates of points
# lying on the sphere that can be used to estimate accessible surface area.
# There are 960 points.

triangles = [
	[ 0.142065,  0.276550,  0.950441], [-0.010435,  0.180640,  0.983494],
	[ 0.093860,  0.192683,  0.976762], [ 0.148649,  0.117947,  0.981831],
	[-0.164403,  0.075573,  0.983494], [-0.309310, -0.031464,  0.950441],
	[-0.213641, -0.017153,  0.976762], [-0.164023, -0.095418,  0.981831],
	[-0.061830,  0.090608,  0.993965], [-0.061717, -0.080548,  0.994838],
	[-0.008704,  0.012755,  0.999881], [ 0.097504,  0.028104,  0.994838],
	[-0.010177, -0.156740,  0.987587], [ 0.144419, -0.211636,  0.966620],
	[ 0.095021, -0.139247,  0.985688], [ 0.149663, -0.047667,  0.987587],
	[ 0.659784, -0.031464,  0.750796], [ 0.539741,  0.075573,  0.838432],
	[ 0.582308, -0.017153,  0.812787], [ 0.538734, -0.095418,  0.837055],
	[ 0.398310,  0.180640,  0.899288], [ 0.245163,  0.276550,  0.929201],
	[ 0.299846,  0.192683,  0.934326], [ 0.251522,  0.117947,  0.960638],
	[ 0.449659,  0.090608,  0.888593], [ 0.303643,  0.028104,  0.952371],
	[ 0.403196,  0.012755,  0.915025], [ 0.449900, -0.080548,  0.889439],
	[ 0.252866, -0.047667,  0.966326], [ 0.249395, -0.211636,  0.944993],
	[ 0.302308, -0.139247,  0.942985], [ 0.399691, -0.156740,  0.903150],
	[-0.348009,  0.509982,  0.786643], [-0.364005,  0.361769,  0.858268],
	[-0.306882,  0.449713,  0.838798], [-0.204165,  0.470843,  0.858268],
	[-0.368013,  0.203516,  0.907275], [-0.360619,  0.043725,  0.931688],
	[-0.316346,  0.133354,  0.939224], [-0.215792,  0.150880,  0.964712],
	[-0.317683,  0.294553,  0.901286], [-0.164646,  0.241277,  0.956387],
	[-0.214764,  0.314722,  0.924568], [-0.158462,  0.403205,  0.901286],
	[-0.061824,  0.255946,  0.964712], [ 0.090756,  0.351740,  0.931688],
	[-0.008845,  0.343190,  0.939224], [-0.055341,  0.416881,  0.907275],
	[-0.211240, -0.625845,  0.750796], [-0.107557, -0.536442,  0.837055],
	[-0.196321, -0.548484,  0.812787], [-0.267147, -0.475041,  0.838432],
	[ 0.000236, -0.429325,  0.903150], [ 0.106147, -0.309387,  0.944993],
	[ 0.019450, -0.332267,  0.942985], [-0.047802, -0.252840,  0.966326],
	[-0.089017, -0.448301,  0.889439], [-0.136870, -0.272499,  0.952371],
	[-0.158871, -0.370796,  0.915025], [-0.248305, -0.385678,  0.888593],
	[-0.201528, -0.191210,  0.960638], [-0.346897, -0.127467,  0.929201],
	[-0.288738, -0.208962,  0.934326], [-0.313420, -0.305039,  0.899288],
	[ 0.699503,  0.043725,  0.713291], [ 0.696646,  0.203516,  0.687943],
	[ 0.661814,  0.133354,  0.737712], [ 0.579521,  0.150880,  0.800868],
	[ 0.673594,  0.361769,  0.644511], [ 0.630591,  0.509983,  0.585041],
	[ 0.613427,  0.449713,  0.649204], [ 0.526769,  0.470843,  0.707687],
	[ 0.648047,  0.294553,  0.702334], [ 0.501791,  0.403205,  0.765266],
	[ 0.562711,  0.314722,  0.764400], [ 0.529250,  0.241277,  0.813437],
	[ 0.409434,  0.416881,  0.811526], [ 0.284882,  0.351740,  0.891696],
	[ 0.379351,  0.343190,  0.859252], [ 0.438090,  0.255946,  0.861724],
	[ 0.276003, -0.309387,  0.910001], [ 0.356751, -0.429325,  0.829704],
	[ 0.354847, -0.332267,  0.873889], [ 0.425849, -0.252840,  0.868749],
	[ 0.429643, -0.536442,  0.726386], [ 0.490790, -0.625845,  0.606170],
	[ 0.501588, -0.548484,  0.669010], [ 0.576783, -0.475041,  0.664573],
	[ 0.433318, -0.448301,  0.781832], [ 0.579301, -0.385678,  0.718097],
	[ 0.507597, -0.370796,  0.777725], [ 0.502148, -0.272499,  0.820727],
	[ 0.643341, -0.305039,  0.702184], [ 0.685915, -0.127466,  0.716431],
	[ 0.634517, -0.208961,  0.744126], [ 0.564808, -0.191210,  0.802764],
	[-0.445452,  0.498942,  0.743391], [-0.584618,  0.427737,  0.689393],
	[-0.499292,  0.427913,  0.753391], [-0.459801,  0.350916,  0.815746],
	[-0.707664,  0.342972,  0.617723], [-0.808789,  0.249821,  0.532399],
	[-0.747671,  0.259481,  0.611275], [-0.715434,  0.181892,  0.674589],
	[-0.633362,  0.349398,  0.690488], [-0.639427,  0.187484,  0.745642],
	[-0.597104,  0.271403,  0.754856], [-0.509028,  0.272874,  0.816352],
	[-0.595202,  0.107892,  0.796300], [-0.456318,  0.032883,  0.889209],
	[-0.507548,  0.111690,  0.854353], [-0.463503,  0.192697,  0.864889],
	[ 0.135742,  0.436900,  0.889209], [ 0.116528,  0.593572,  0.796300],
	[ 0.081035,  0.513335,  0.854354], [-0.010454,  0.501854,  0.864889],
	[ 0.091454,  0.732507,  0.674589], [ 0.062235,  0.844202,  0.532399],
	[ 0.030958,  0.790812,  0.611275], [-0.061373,  0.783997,  0.617723],
	[ 0.058537,  0.663770,  0.745643], [-0.094445,  0.717151,  0.690488],
	[-0.035036,  0.654954,  0.754856], [-0.068516,  0.573476,  0.816352],
	[-0.185164,  0.700322,  0.689393], [-0.302204,  0.596693,  0.743391],
	[-0.216434,  0.620933,  0.753391], [-0.159133,  0.556089,  0.815746],
	[-0.131372, -0.685178,  0.716431], [ 0.049500, -0.710273,  0.702184],
	[-0.036748, -0.667028,  0.744125], [-0.027863, -0.595645,  0.802764],
	[ 0.232067, -0.710273,  0.664573], [ 0.403843, -0.685178,  0.606170],
	[ 0.327870, -0.667028,  0.669010], [ 0.342885, -0.595645,  0.726386],
	[ 0.147936, -0.680038,  0.718097], [ 0.259470, -0.566933,  0.781832],
	[ 0.160220, -0.607844,  0.777725], [ 0.070675, -0.566933,  0.820726],
	[ 0.269715, -0.488718,  0.829704], [ 0.187471, -0.369801,  0.910001],
	[ 0.180031, -0.451560,  0.873889], [ 0.080186, -0.488718,  0.868749],
	[-0.431392, -0.137040,  0.891696], [-0.537458, -0.229271,  0.811526],
	[-0.457871, -0.228124,  0.859251], [-0.398047, -0.314627,  0.861724],
	[-0.630483, -0.318858,  0.707687], [-0.704780, -0.401265,  0.585040],
	[-0.642400, -0.407254,  0.649204], [-0.582445, -0.495343,  0.644511],
	[-0.558394, -0.320257,  0.765266], [-0.510541, -0.496059,  0.702334],
	[-0.498210, -0.409243,  0.764399], [-0.417621, -0.404861,  0.813437],
	[-0.443487, -0.574503,  0.687943], [-0.295735, -0.635419,  0.713291],
	[-0.365455, -0.567647,  0.737712], [-0.351773, -0.484629,  0.800869],
	[ 0.953363,  0.249821,  0.169376], [ 0.894196,  0.342972,  0.287722],
	[ 0.928397,  0.259481,  0.265986], [ 0.923809,  0.181892,  0.336886],
	[ 0.809496,  0.427738,  0.402190], [ 0.703004,  0.498942,  0.506796],
	[ 0.756413,  0.427913,  0.494702], [ 0.744783,  0.350916,  0.567588],
	[ 0.854704,  0.349398,  0.383930], [ 0.790242,  0.272874,  0.548687],
	[ 0.846840,  0.271403,  0.457387], [ 0.882075,  0.187484,  0.432196],
	[ 0.767608,  0.192697,  0.611266], [ 0.770620,  0.032883,  0.636446],
	[ 0.803902,  0.111690,  0.584180], [ 0.861473,  0.107892,  0.496208],
	[ 0.153263,  0.844202,  0.513647], [ 0.182623,  0.732507,  0.655807],
	[ 0.213168,  0.790812,  0.573738], [ 0.300529,  0.783997,  0.543167],
	[ 0.207697,  0.593572,  0.777518], [ 0.226770,  0.436900,  0.870456],
	[ 0.263245,  0.513335,  0.816816], [ 0.351449,  0.501854,  0.790333],
	[ 0.240944,  0.663770,  0.708065], [ 0.385598,  0.573477,  0.722799],
	[ 0.330539,  0.654954,  0.679543], [ 0.359669,  0.717152,  0.596935],
	[ 0.468598,  0.556089,  0.686426], [ 0.571420,  0.596694,  0.563414],
	[ 0.496587,  0.620933,  0.606501], [ 0.442567,  0.700323,  0.560073],
	[ 0.878628, -0.401265,  0.258841], [ 0.858857, -0.318858,  0.400866],
	[ 0.846689, -0.407254,  0.342435], [ 0.789760, -0.495343,  0.361821],
	[ 0.814449, -0.229271,  0.533018], [ 0.748706, -0.137040,  0.648582],
	[ 0.760206, -0.228124,  0.608314], [ 0.706230, -0.314627,  0.634231],
	[ 0.815396, -0.320257,  0.482250], [ 0.705126, -0.404861,  0.582139],
	[ 0.759771, -0.409243,  0.505241], [ 0.746566, -0.496059,  0.443356],
	[ 0.639671, -0.484629,  0.596620], [ 0.553582, -0.635419,  0.538322],
	[ 0.627277, -0.567647,  0.533199], [ 0.679284, -0.574503,  0.456640],
	[-0.589414,  0.770371,  0.243147], [-0.691910,  0.643752,  0.326871],
	[-0.616865,  0.717414,  0.323721], [-0.560020,  0.727500,  0.396384],
	[-0.771571,  0.492950,  0.402092], [-0.822949,  0.328280,  0.463667],
	[-0.775962,  0.416237,  0.473952], [-0.721793,  0.421259,  0.549141],
	[-0.709675,  0.576731,  0.404650], [-0.661673,  0.506272,  0.553062],
	[-0.653675,  0.584862,  0.480256], [-0.579294,  0.662218,  0.475274],
	[-0.598792,  0.506276,  0.620591], [-0.459869,  0.578830,  0.673406],
	[-0.527761,  0.585660,  0.615200], [-0.516093,  0.662828,  0.542501],
	[-0.763951, -0.314554,  0.563414], [-0.688655, -0.233611,  0.686426],
	[-0.759240, -0.236034,  0.606500], [-0.813471, -0.156790,  0.560073],
	[-0.595443, -0.144298,  0.790333], [-0.489504, -0.051880,  0.870456],
	[-0.573977, -0.057979,  0.816816], [-0.628440,  0.022998,  0.777518],
	[-0.674586, -0.149986,  0.722799], [-0.705928,  0.017632,  0.708064],
	[-0.730382, -0.069011,  0.679543], [-0.798919, -0.073461,  0.596935],
	[-0.748672,  0.096998,  0.655807], [-0.841976,  0.165058,  0.513646],
	[-0.814100,  0.089811,  0.573738], [-0.839603,  0.005978,  0.543167],
	[-0.502466,  0.829703,  0.243147], [-0.473263,  0.786703,  0.396384],
	[-0.443147,  0.835958,  0.323721], [-0.347195,  0.878983,  0.326871],
	[-0.429056,  0.722221,  0.542501], [-0.371337,  0.639244,  0.673406],
	[-0.352945,  0.704953,  0.615200], [-0.253129,  0.742154,  0.620591],
	[-0.405447,  0.780851,  0.475274], [-0.230200,  0.800706,  0.553062],
	[-0.306298,  0.821910,  0.480256], [-0.278310,  0.871092,  0.404650],
	[-0.129121,  0.825694,  0.549142], [-0.005662,  0.885991,  0.463667],
	[-0.104696,  0.874305,  0.473952], [-0.177730,  0.898184,  0.402092],
	[ 0.053324, -0.964447,  0.258841], [ 0.173326, -0.915993,  0.361821],
	[ 0.070546, -0.936889,  0.342435], [-0.016203, -0.915993,  0.400866],
	[ 0.287315, -0.841980,  0.456640], [ 0.389865, -0.747138,  0.538322],
	[ 0.299891, -0.791053,  0.533199], [ 0.218068, -0.772328,  0.596620],
	[ 0.189740, -0.876033,  0.443356], [ 0.119927, -0.804197,  0.582139],
	[ 0.104085, -0.856678,  0.505241], [ 0.000945, -0.876033,  0.482250],
	[ 0.035501, -0.772328,  0.634231], [-0.145350, -0.747137,  0.648582],
	[-0.064727, -0.791053,  0.608314], [-0.083434, -0.841980,  0.533018],
	[-0.580197, -0.796670,  0.169376], [-0.506169, -0.793915,  0.336886],
	[-0.580090, -0.769900,  0.265986], [-0.645367, -0.707614,  0.287722],
	[-0.414536, -0.762848,  0.496208], [-0.311566, -0.705594,  0.636446],
	[-0.397084, -0.707854,  0.584180], [-0.459283, -0.644525,  0.611266],
	[-0.496161, -0.753015,  0.432196], [-0.542194, -0.636371,  0.548687],
	[-0.561458, -0.689610,  0.457387], [-0.636953, -0.668498,  0.383929],
	[-0.598292, -0.565589,  0.567588], [-0.720899, -0.472719,  0.506796],
	[-0.674230, -0.548348,  0.494701], [-0.693419, -0.597841,  0.402189],
	[ 0.939203,  0.328280,  0.100644], [ 0.867672,  0.492950,  0.064390],
	[ 0.900107,  0.416238,  0.128663], [ 0.880067,  0.421259,  0.219140],
	[ 0.764766,  0.643752,  0.026779], [ 0.637523,  0.770371, -0.009616],
	[ 0.694586,  0.717414,  0.053548], [ 0.671090,  0.727501,  0.142762],
	[ 0.811827,  0.576732,  0.091203], [ 0.719976,  0.662218,  0.207610],
	[ 0.790269,  0.584862,  0.182787], [ 0.826393,  0.506272,  0.246504],
	[ 0.688492,  0.662828,  0.294343], [ 0.688586,  0.578830,  0.436812],
	[ 0.727944,  0.585660,  0.356510], [ 0.795322,  0.506276,  0.333388],
	[ 0.793692, -0.051880,  0.606104], [ 0.859336, -0.144298,  0.490632],
	[ 0.850086, -0.057979,  0.523443], [ 0.884581,  0.022998,  0.465819],
	[ 0.903889, -0.233611,  0.358344], [ 0.924434, -0.314554,  0.215588],
	[ 0.937136, -0.236034,  0.257028], [ 0.968602, -0.156790,  0.192946],
	[ 0.905342, -0.149985,  0.397316], [ 0.969804, -0.073461,  0.232558],
	[ 0.939499, -0.069011,  0.335529], [ 0.928308,  0.017632,  0.371394],
	[ 0.985924,  0.005978,  0.167088], [ 0.976435,  0.165058,  0.139033],
	[ 0.974580,  0.089811,  0.205249], [ 0.946917,  0.096998,  0.306496],
	[ 0.607263,  0.639244,  0.471804], [ 0.608542,  0.722222,  0.328744],
	[ 0.567363,  0.704953,  0.425606], [ 0.477805,  0.742154,  0.470010],
	[ 0.591397,  0.786703,  0.177053], [ 0.557655,  0.829704,  0.024750],
	[ 0.535013,  0.835958,  0.122209], [ 0.448119,  0.878983,  0.163028],
	[ 0.560284,  0.780851,  0.276323], [ 0.415585,  0.871092,  0.261700],
	[ 0.471177,  0.821910,  0.320087], [ 0.430052,  0.800706,  0.417043],
	[ 0.322184,  0.898184,  0.299105], [ 0.188464,  0.885991,  0.423675],
	[ 0.283499,  0.874305,  0.393979], [ 0.335654,  0.825695,  0.453393],
	[ 0.537751, -0.705594,  0.461477], [ 0.576908, -0.762848,  0.291960],
	[ 0.595647, -0.707854,  0.379666], [ 0.663487, -0.644525,  0.379963],
	[ 0.598108, -0.793915,  0.109393], [ 0.599900, -0.796670, -0.073737],
	[ 0.637987, -0.769900,  0.015049], [ 0.706539, -0.707614,  0.009214],
	[ 0.626586, -0.753014,  0.200898], [ 0.736836, -0.668498,  0.100914],
	[ 0.696523, -0.689610,  0.198229], [ 0.714913, -0.636371,  0.289709],
	[ 0.795922, -0.597841,  0.095369], [ 0.862509, -0.472719,  0.180596],
	[ 0.814860, -0.548347,  0.187932], [ 0.773913, -0.565589,  0.284898],
	[-0.635102,  0.754881,  0.163708], [-0.741664,  0.665917,  0.080560],
	[-0.708318,  0.686408,  0.164710], [-0.737669,  0.628238,  0.247309],
	[-0.831159,  0.556009, -0.005306], [-0.898352,  0.430169, -0.088985],
	[-0.885636,  0.464328, -0.007038], [-0.913858,  0.399457,  0.072784],
	[-0.807220,  0.584943,  0.078979], [-0.889599,  0.428997,  0.156767],
	[-0.837161,  0.522654,  0.161227], [-0.801228,  0.545692,  0.245467],
	[-0.903436,  0.359676,  0.233317], [-0.868637,  0.312790,  0.384229],
	[-0.867415,  0.385231,  0.314941], [-0.817330,  0.477436,  0.322531],
	[-0.893722,  0.147515,  0.423675], [-0.953826,  0.027443,  0.299105],
	[-0.917487,  0.054759,  0.393979], [-0.891236, -0.011527,  0.453392],
	[-0.981859, -0.096824,  0.163028], [-0.975905, -0.216788,  0.024750],
	[-0.973474, -0.193423,  0.122209], [-0.948166, -0.263883,  0.177053],
	[-0.962650, -0.069407,  0.261700], [-0.931373, -0.237045,  0.276323],
	[-0.937120, -0.139103,  0.320087], [-0.902383, -0.108539,  0.417042],
	[-0.894372, -0.303357,  0.328744], [-0.816640, -0.332417,  0.471803],
	[-0.863280, -0.271307,  0.425606], [-0.865270, -0.174351,  0.470010],
	[ 0.025419,  0.922888,  0.384229], [-0.005554,  0.972385,  0.233317],
	[-0.042483,  0.948160,  0.314941], [-0.146601,  0.935137,  0.322531],
	[-0.038797,  0.996593,  0.072784], [-0.073048,  0.993351, -0.088985],
	[-0.109492,  0.993963, -0.007038], [-0.214725,  0.976660, -0.005306],
	[-0.075148,  0.984772,  0.156767], [-0.250395,  0.964917,  0.078979],
	[-0.181476,  0.970089,  0.161227], [-0.216030,  0.945028,  0.245467],
	[-0.349695,  0.933393,  0.080560], [-0.471386,  0.866600,  0.163708],
	[-0.380933,  0.909814,  0.164710], [-0.316066,  0.915937,  0.247309],
	[ 0.125737, -0.975487,  0.180596], [ 0.266524, -0.959098,  0.095369],
	[ 0.213532, -0.958689,  0.187932], [ 0.244515, -0.926847,  0.284898],
	[ 0.401328, -0.915888,  0.009214], [ 0.523132, -0.849056, -0.073737],
	[ 0.484320, -0.874762,  0.015049], [ 0.521220, -0.846382,  0.109393],
	[ 0.353858, -0.929839,  0.100914], [ 0.472752, -0.857989,  0.200898],
	[ 0.388214, -0.899997,  0.198229], [ 0.331935, -0.897712,  0.289709],
	[ 0.500020, -0.815316,  0.291960], [ 0.460982, -0.757980,  0.461477],
	[ 0.441980, -0.812716,  0.379666], [ 0.358276, -0.852799,  0.379963],
	[-0.241049, -0.757980,  0.606104], [-0.343909, -0.815316,  0.465819],
	[-0.255929, -0.812716,  0.523443], [-0.178924, -0.852799,  0.490632],
	[-0.435542, -0.846382,  0.306496], [-0.509680, -0.849056,  0.139033],
	[-0.438935, -0.874762,  0.205249], [-0.365007, -0.915888,  0.167088],
	[-0.354854, -0.857989,  0.371394], [-0.285159, -0.929839,  0.232558],
	[-0.278254, -0.899997,  0.335529], [-0.190400, -0.897712,  0.397316],
	[-0.207127, -0.959099,  0.192946], [-0.044118, -0.975487,  0.215588],
	[-0.121865, -0.958690,  0.257028], [-0.112000, -0.926847,  0.358344],
	[-0.790032, -0.430169,  0.436811], [-0.868215, -0.399457,  0.294343],
	[-0.810741, -0.464328,  0.356510], [-0.761384, -0.556010,  0.333388],
	[-0.922091, -0.359676,  0.142762], [-0.949774, -0.312790, -0.009616],
	[-0.921265, -0.385231,  0.053548], [-0.878258, -0.477436,  0.026779],
	[-0.879125, -0.428997,  0.207610], [-0.833008, -0.545692,  0.091203],
	[-0.832720, -0.522654,  0.182787], [-0.772708, -0.584943,  0.246504],
	[-0.775352, -0.628238,  0.064390], [-0.648094, -0.754881,  0.100644],
	[-0.715744, -0.686408,  0.128662], [-0.713115, -0.665917,  0.219140],
	[ 0.790032,  0.430169, -0.436811], [ 0.761384,  0.556010, -0.333388],
	[ 0.810741,  0.464328, -0.356510], [ 0.868215,  0.399457, -0.294343],
	[ 0.713115,  0.665917, -0.219140], [ 0.648094,  0.754881, -0.100644],
	[ 0.715744,  0.686408, -0.128662], [ 0.775352,  0.628238, -0.064390],
	[ 0.772708,  0.584943, -0.246504], [ 0.833008,  0.545692, -0.091203],
	[ 0.832720,  0.522654, -0.182787], [ 0.879125,  0.428997, -0.207610],
	[ 0.878258,  0.477436, -0.026779], [ 0.949774,  0.312790,  0.009616],
	[ 0.921265,  0.385231, -0.053548], [ 0.922091,  0.359676, -0.142762],
	[ 0.906224, -0.216788, -0.362990], [ 0.966347, -0.096824, -0.238325],
	[ 0.942512, -0.193423, -0.272506], [ 0.940940, -0.263883, -0.212124],
	[ 0.994380,  0.027443, -0.102248], [ 0.988407,  0.147515,  0.035935],
	[ 0.998499,  0.054760, -0.000736], [ 0.997870, -0.011527,  0.064215],
	[ 0.987702, -0.069407, -0.140095], [ 0.993741, -0.108539,  0.026419],
	[ 0.987329, -0.139103, -0.076371], [ 0.964751, -0.237045, -0.114300],
	[ 0.980585, -0.174351,  0.089743], [ 0.936624, -0.332417,  0.110612],
	[ 0.961207, -0.271307,  0.049741], [ 0.951483, -0.303357, -0.051523],
	[ 0.031928,  0.993351, -0.110612], [ 0.064406,  0.996593,  0.051523],
	[ 0.097795,  0.993963, -0.049741], [ 0.195143,  0.976660, -0.089743],
	[ 0.097319,  0.972385,  0.212125], [ 0.128517,  0.922888,  0.362990],
	[ 0.163503,  0.948160,  0.272506], [ 0.262143,  0.935137,  0.238325],
	[ 0.130991,  0.984773,  0.114300], [ 0.295459,  0.945028,  0.140095],
	[ 0.230424,  0.970089,  0.076371], [ 0.261222,  0.964917, -0.026419],
	[ 0.388078,  0.915937,  0.102248], [ 0.497708,  0.866600, -0.035935],
	[ 0.415016,  0.909814,  0.000736], [ 0.353062,  0.933394, -0.064215],
	[ 0.868637, -0.312790, -0.384229], [ 0.903436, -0.359676, -0.233317],
	[ 0.867415, -0.385231, -0.314941], [ 0.817330, -0.477436, -0.322531],
	[ 0.913858, -0.399457, -0.072784], [ 0.898352, -0.430169,  0.088985],
	[ 0.885636, -0.464328,  0.007038], [ 0.831159, -0.556009,  0.005306],
	[ 0.889599, -0.428997, -0.156767], [ 0.807220, -0.584943, -0.078979],
	[ 0.837161, -0.522654, -0.161227], [ 0.801228, -0.545692, -0.245467],
	[ 0.741664, -0.665917, -0.080560], [ 0.635102, -0.754881, -0.163708],
	[ 0.708318, -0.686408, -0.164710], [ 0.737669, -0.628238, -0.247309],
	[-0.537751,  0.705594, -0.461477], [-0.663487,  0.644525, -0.379963],
	[-0.595647,  0.707854, -0.379666], [-0.576908,  0.762848, -0.291960],
	[-0.773913,  0.565589, -0.284898], [-0.862509,  0.472719, -0.180596],
	[-0.814860,  0.548347, -0.187932], [-0.795922,  0.597841, -0.095369],
	[-0.714913,  0.636371, -0.289709], [-0.736836,  0.668498, -0.100914],
	[-0.696523,  0.689610, -0.198229], [-0.626586,  0.753014, -0.200898],
	[-0.706539,  0.707614, -0.009214], [-0.599900,  0.796670,  0.073737],
	[-0.637987,  0.769900, -0.015049], [-0.598108,  0.793915, -0.109393],
	[-0.988407, -0.147515, -0.035935], [-0.994380, -0.027443,  0.102248],
	[-0.998499, -0.054760,  0.000736], [-0.997870,  0.011527, -0.064215],
	[-0.966347,  0.096824,  0.238325], [-0.906224,  0.216788,  0.362990],
	[-0.942512,  0.193423,  0.272506], [-0.940940,  0.263883,  0.212124],
	[-0.987702,  0.069407,  0.140095], [-0.964751,  0.237045,  0.114300],
	[-0.987329,  0.139103,  0.076371], [-0.993741,  0.108539, -0.026419],
	[-0.951483,  0.303357,  0.051523], [-0.936624,  0.332417, -0.110612],
	[-0.961207,  0.271307, -0.049741], [-0.980585,  0.174351, -0.089743],
	[-0.460982,  0.757980, -0.461477], [-0.500020,  0.815316, -0.291960],
	[-0.441980,  0.812716, -0.379666], [-0.358276,  0.852799, -0.379963],
	[-0.521220,  0.846382, -0.109393], [-0.523132,  0.849056,  0.073737],
	[-0.484320,  0.874762, -0.015049], [-0.401328,  0.915888, -0.009214],
	[-0.472752,  0.857989, -0.200898], [-0.353858,  0.929839, -0.100914],
	[-0.388214,  0.899997, -0.198229], [-0.331935,  0.897712, -0.289709],
	[-0.266524,  0.959098, -0.095369], [-0.125737,  0.975487, -0.180596],
	[-0.213532,  0.958689, -0.187932], [-0.244515,  0.926847, -0.284898],
	[-0.025419, -0.922888, -0.384229], [ 0.146601, -0.935137, -0.322531],
	[ 0.042483, -0.948160, -0.314941], [ 0.005554, -0.972385, -0.233317],
	[ 0.316066, -0.915937, -0.247309], [ 0.471386, -0.866600, -0.163708],
	[ 0.380933, -0.909814, -0.164710], [ 0.349695, -0.933393, -0.080560],
	[ 0.216030, -0.945028, -0.245467], [ 0.250395, -0.964917, -0.078979],
	[ 0.181476, -0.970089, -0.161227], [ 0.075148, -0.984772, -0.156767],
	[ 0.214725, -0.976660,  0.005306], [ 0.073048, -0.993351,  0.088985],
	[ 0.109492, -0.993963,  0.007038], [ 0.038797, -0.996593, -0.072784],
	[-0.128517, -0.922888, -0.362990], [-0.097319, -0.972385, -0.212125],
	[-0.163503, -0.948160, -0.272506], [-0.262143, -0.935137, -0.238325],
	[-0.064406, -0.996593, -0.051523], [-0.031928, -0.993351,  0.110612],
	[-0.097795, -0.993963,  0.049741], [-0.195143, -0.976660,  0.089743],
	[-0.130991, -0.984773, -0.114300], [-0.261222, -0.964917,  0.026419],
	[-0.230424, -0.970089, -0.076371], [-0.295459, -0.945028, -0.140095],
	[-0.353062, -0.933394,  0.064215], [-0.497708, -0.866600,  0.035935],
	[-0.415016, -0.909814, -0.000736], [-0.388078, -0.915937, -0.102248],
	[-0.688586, -0.578830, -0.436812], [-0.688492, -0.662828, -0.294343],
	[-0.727944, -0.585660, -0.356510], [-0.795322, -0.506276, -0.333388],
	[-0.671090, -0.727501, -0.142762], [-0.637523, -0.770371,  0.009616],
	[-0.694586, -0.717414, -0.053548], [-0.764766, -0.643752, -0.026779],
	[-0.719976, -0.662218, -0.207610], [-0.811827, -0.576732, -0.091203],
	[-0.790269, -0.584862, -0.182787], [-0.826393, -0.506272, -0.246504],
	[-0.867672, -0.492950, -0.064390], [-0.939203, -0.328280, -0.100644],
	[-0.900107, -0.416238, -0.128663], [-0.880067, -0.421259, -0.219140],
	[ 0.720899,  0.472719, -0.506796], [ 0.598292,  0.565589, -0.567588],
	[ 0.674230,  0.548348, -0.494701], [ 0.693419,  0.597841, -0.402189],
	[ 0.459283,  0.644525, -0.611266], [ 0.311566,  0.705594, -0.636446],
	[ 0.397084,  0.707854, -0.584180], [ 0.414536,  0.762848, -0.496208],
	[ 0.542194,  0.636371, -0.548687], [ 0.496161,  0.753015, -0.432196],
	[ 0.561458,  0.689610, -0.457387], [ 0.636953,  0.668498, -0.383929],
	[ 0.506169,  0.793915, -0.336886], [ 0.580197,  0.796670, -0.169376],
	[ 0.580090,  0.769900, -0.265986], [ 0.645367,  0.707614, -0.287722],
	[ 0.975905,  0.216788, -0.024750], [ 0.981859,  0.096824, -0.163028],
	[ 0.973474,  0.193423, -0.122209], [ 0.948166,  0.263883, -0.177053],
	[ 0.953826, -0.027443, -0.299105], [ 0.893722, -0.147515, -0.423675],
	[ 0.917487, -0.054759, -0.393979], [ 0.891236,  0.011527, -0.453392],
	[ 0.962650,  0.069407, -0.261700], [ 0.902383,  0.108539, -0.417042],
	[ 0.937120,  0.139103, -0.320087], [ 0.931373,  0.237045, -0.276323],
	[ 0.865270,  0.174351, -0.470010], [ 0.816640,  0.332417, -0.471803],
	[ 0.863280,  0.271307, -0.425606], [ 0.894372,  0.303357, -0.328744],
	[ 0.509680,  0.849056, -0.139033], [ 0.435542,  0.846382, -0.306496],
	[ 0.438935,  0.874762, -0.205249], [ 0.365007,  0.915888, -0.167088],
	[ 0.343909,  0.815316, -0.465819], [ 0.241049,  0.757980, -0.606104],
	[ 0.255929,  0.812716, -0.523443], [ 0.178924,  0.852799, -0.490632],
	[ 0.354854,  0.857989, -0.371394], [ 0.190400,  0.897712, -0.397316],
	[ 0.278254,  0.899997, -0.335529], [ 0.285159,  0.929839, -0.232558],
	[ 0.112000,  0.926847, -0.358344], [ 0.044118,  0.975487, -0.215588],
	[ 0.121865,  0.958690, -0.257028], [ 0.207127,  0.959099, -0.192946],
	[ 0.589414, -0.770371, -0.243147], [ 0.560020, -0.727500, -0.396384],
	[ 0.616865, -0.717414, -0.323721], [ 0.691910, -0.643752, -0.326871],
	[ 0.516093, -0.662828, -0.542501], [ 0.459869, -0.578830, -0.673406],
	[ 0.527761, -0.585660, -0.615200], [ 0.598792, -0.506276, -0.620591],
	[ 0.579294, -0.662218, -0.475274], [ 0.661673, -0.506272, -0.553062],
	[ 0.653675, -0.584862, -0.480256], [ 0.709675, -0.576731, -0.404650],
	[ 0.721793, -0.421259, -0.549141], [ 0.822949, -0.328280, -0.463667],
	[ 0.775962, -0.416237, -0.473952], [ 0.771571, -0.492950, -0.402092],
	[-0.553582,  0.635419, -0.538322], [-0.639671,  0.484629, -0.596620],
	[-0.627277,  0.567647, -0.533199], [-0.679284,  0.574503, -0.456640],
	[-0.706230,  0.314627, -0.634231], [-0.748706,  0.137040, -0.648582],
	[-0.760206,  0.228124, -0.608314], [-0.814449,  0.229271, -0.533018],
	[-0.705126,  0.404861, -0.582139], [-0.815396,  0.320257, -0.482250],
	[-0.759771,  0.409243, -0.505241], [-0.746566,  0.496059, -0.443356],
	[-0.858857,  0.318858, -0.400866], [-0.878628,  0.401265, -0.258841],
	[-0.846689,  0.407254, -0.342435], [-0.789760,  0.495343, -0.361821],
	[-0.924434,  0.314554, -0.215588], [-0.903889,  0.233611, -0.358344],
	[-0.937136,  0.236034, -0.257028], [-0.968602,  0.156790, -0.192946],
	[-0.859336,  0.144298, -0.490632], [-0.793692,  0.051880, -0.606104],
	[-0.850086,  0.057979, -0.523443], [-0.884581, -0.022998, -0.465819],
	[-0.905342,  0.149985, -0.397316], [-0.928308, -0.017632, -0.371394],
	[-0.939499,  0.069011, -0.335529], [-0.969804,  0.073461, -0.232558],
	[-0.946917, -0.096998, -0.306496], [-0.976435, -0.165058, -0.139033],
	[-0.974580, -0.089811, -0.205249], [-0.985924, -0.005978, -0.167088],
	[-0.053324,  0.964447, -0.258841], [ 0.016203,  0.915993, -0.400866],
	[-0.070546,  0.936889, -0.342435], [-0.173326,  0.915993, -0.361821],
	[ 0.083434,  0.841980, -0.533018], [ 0.145350,  0.747137, -0.648582],
	[ 0.064727,  0.791053, -0.608314], [-0.035501,  0.772328, -0.634231],
	[-0.000945,  0.876033, -0.482250], [-0.119927,  0.804197, -0.582139],
	[-0.104085,  0.856678, -0.505241], [-0.189740,  0.876033, -0.443356],
	[-0.218068,  0.772328, -0.596620], [-0.389865,  0.747138, -0.538322],
	[-0.299891,  0.791053, -0.533199], [-0.287315,  0.841980, -0.456640],
	[ 0.005662, -0.885991, -0.463667], [ 0.129121, -0.825694, -0.549142],
	[ 0.104696, -0.874305, -0.473952], [ 0.177730, -0.898184, -0.402092],
	[ 0.253129, -0.742154, -0.620591], [ 0.371337, -0.639244, -0.673406],
	[ 0.352945, -0.704953, -0.615200], [ 0.429056, -0.722221, -0.542501],
	[ 0.230200, -0.800706, -0.553062], [ 0.405447, -0.780851, -0.475274],
	[ 0.306298, -0.821910, -0.480256], [ 0.278310, -0.871092, -0.404650],
	[ 0.473263, -0.786703, -0.396384], [ 0.502466, -0.829703, -0.243147],
	[ 0.443147, -0.835958, -0.323721], [ 0.347195, -0.878983, -0.326871],
	[-0.557655, -0.829704, -0.024750], [-0.591397, -0.786703, -0.177053],
	[-0.535013, -0.835958, -0.122209], [-0.448119, -0.878983, -0.163028],
	[-0.608542, -0.722222, -0.328744], [-0.607263, -0.639244, -0.471804],
	[-0.567363, -0.704953, -0.425606], [-0.477805, -0.742154, -0.470010],
	[-0.560284, -0.780851, -0.276323], [-0.430052, -0.800706, -0.417043],
	[-0.471177, -0.821910, -0.320087], [-0.415585, -0.871092, -0.261700],
	[-0.335654, -0.825695, -0.453393], [-0.188464, -0.885991, -0.423675],
	[-0.283499, -0.874305, -0.393979], [-0.322184, -0.898184, -0.299105],
	[-0.953363, -0.249821, -0.169376], [-0.923809, -0.181892, -0.336886],
	[-0.928397, -0.259481, -0.265986], [-0.894196, -0.342972, -0.287722],
	[-0.861473, -0.107892, -0.496208], [-0.770620, -0.032883, -0.636446],
	[-0.803902, -0.111690, -0.584180], [-0.767608, -0.192697, -0.611266],
	[-0.882075, -0.187484, -0.432196], [-0.790242, -0.272874, -0.548687],
	[-0.846840, -0.271403, -0.457387], [-0.854704, -0.349398, -0.383930],
	[-0.744783, -0.350916, -0.567588], [-0.703004, -0.498942, -0.506796],
	[-0.756413, -0.427913, -0.494702], [-0.809496, -0.427738, -0.402190],
	[ 0.431392,  0.137040, -0.891696], [ 0.398047,  0.314627, -0.861724],
	[ 0.457871,  0.228124, -0.859251], [ 0.537458,  0.229271, -0.811526],
	[ 0.351773,  0.484629, -0.800869], [ 0.295735,  0.635419, -0.713291],
	[ 0.365455,  0.567647, -0.737712], [ 0.443487,  0.574503, -0.687943],
	[ 0.417621,  0.404861, -0.813437], [ 0.510541,  0.496059, -0.702334],
	[ 0.498210,  0.409243, -0.764399], [ 0.558394,  0.320257, -0.765266],
	[ 0.582445,  0.495343, -0.644511], [ 0.704780,  0.401265, -0.585040],
	[ 0.642400,  0.407254, -0.649204], [ 0.630483,  0.318858, -0.707687],
	[ 0.489504,  0.051880, -0.870456], [ 0.595443,  0.144298, -0.790333],
	[ 0.573977,  0.057979, -0.816816], [ 0.628440, -0.022998, -0.777518],
	[ 0.688655,  0.233611, -0.686426], [ 0.763951,  0.314554, -0.563414],
	[ 0.759240,  0.236034, -0.606500], [ 0.813471,  0.156790, -0.560073],
	[ 0.674586,  0.149986, -0.722799], [ 0.798919,  0.073461, -0.596935],
	[ 0.730382,  0.069011, -0.679543], [ 0.705928, -0.017632, -0.708064],
	[ 0.839603, -0.005978, -0.543167], [ 0.841976, -0.165058, -0.513646],
	[ 0.814100, -0.089811, -0.573738], [ 0.748672, -0.096998, -0.655807],
	[ 0.456318, -0.032883, -0.889209], [ 0.595202, -0.107892, -0.796300],
	[ 0.507548, -0.111690, -0.854353], [ 0.463503, -0.192697, -0.864889],
	[ 0.715434, -0.181892, -0.674589], [ 0.808789, -0.249821, -0.532399],
	[ 0.747671, -0.259481, -0.611275], [ 0.707664, -0.342972, -0.617723],
	[ 0.639427, -0.187484, -0.745642], [ 0.633362, -0.349398, -0.690488],
	[ 0.597104, -0.271403, -0.754856], [ 0.509028, -0.272874, -0.816352],
	[ 0.584618, -0.427737, -0.689393], [ 0.445452, -0.498942, -0.743391],
	[ 0.499292, -0.427913, -0.753391], [ 0.459801, -0.350916, -0.815746],
	[-0.276003,  0.309387, -0.910001], [-0.425849,  0.252840, -0.868749],
	[-0.354847,  0.332267, -0.873889], [-0.356751,  0.429325, -0.829704],
	[-0.564808,  0.191210, -0.802764], [-0.685915,  0.127466, -0.716431],
	[-0.634517,  0.208961, -0.744126], [-0.643341,  0.305039, -0.702184],
	[-0.502148,  0.272499, -0.820727], [-0.579301,  0.385678, -0.718097],
	[-0.507597,  0.370796, -0.777725], [-0.433318,  0.448301, -0.781832],
	[-0.576783,  0.475041, -0.664573], [-0.490790,  0.625845, -0.606170],
	[-0.501588,  0.548484, -0.669010], [-0.429643,  0.536442, -0.726386],
	[-0.187471,  0.369801, -0.910001], [-0.269715,  0.488718, -0.829704],
	[-0.180031,  0.451560, -0.873889], [-0.080186,  0.488718, -0.868749],
	[-0.342885,  0.595645, -0.726386], [-0.403843,  0.685178, -0.606170],
	[-0.327870,  0.667028, -0.669010], [-0.232067,  0.710273, -0.664573],
	[-0.259470,  0.566933, -0.781832], [-0.147936,  0.680038, -0.718097],
	[-0.160220,  0.607844, -0.777725], [-0.070675,  0.566933, -0.820726],
	[-0.049500,  0.710273, -0.702184], [ 0.131372,  0.685178, -0.716431],
	[ 0.036748,  0.667028, -0.744125], [ 0.027863,  0.595645, -0.802764],
	[-0.135742, -0.436900, -0.889209], [ 0.010454, -0.501854, -0.864889],
	[-0.081035, -0.513335, -0.854354], [-0.116528, -0.593572, -0.796300],
	[ 0.159133, -0.556089, -0.815746], [ 0.302204, -0.596693, -0.743391],
	[ 0.216434, -0.620933, -0.753391], [ 0.185164, -0.700322, -0.689393],
	[ 0.068516, -0.573476, -0.816352], [ 0.094445, -0.717151, -0.690488],
	[ 0.035036, -0.654954, -0.754856], [-0.058537, -0.663770, -0.745643],
	[ 0.061373, -0.783997, -0.617723], [-0.062235, -0.844202, -0.532399],
	[-0.030958, -0.790812, -0.611275], [-0.091454, -0.732507, -0.674589],
	[-0.226770, -0.436900, -0.870456], [-0.207697, -0.593572, -0.777518],
	[-0.263245, -0.513335, -0.816816], [-0.351449, -0.501854, -0.790333],
	[-0.182623, -0.732507, -0.655807], [-0.153263, -0.844202, -0.513647],
	[-0.213168, -0.790812, -0.573738], [-0.300529, -0.783997, -0.543167],
	[-0.240944, -0.663770, -0.708065], [-0.359669, -0.717152, -0.596935],
	[-0.330539, -0.654954, -0.679543], [-0.385598, -0.573477, -0.722799],
	[-0.442567, -0.700323, -0.560073], [-0.571420, -0.596694, -0.563414],
	[-0.496587, -0.620933, -0.606501], [-0.468598, -0.556089, -0.686426],
	[-0.284882, -0.351740, -0.891696], [-0.409434, -0.416881, -0.811526],
	[-0.379351, -0.343190, -0.859252], [-0.438090, -0.255946, -0.861724],
	[-0.526769, -0.470843, -0.707687], [-0.630591, -0.509983, -0.585041],
	[-0.613427, -0.449713, -0.649204], [-0.673594, -0.361769, -0.644511],
	[-0.501791, -0.403205, -0.765266], [-0.648047, -0.294553, -0.702334],
	[-0.562711, -0.314722, -0.764400], [-0.529250, -0.241277, -0.813437],
	[-0.696646, -0.203516, -0.687943], [-0.699503, -0.043725, -0.713291],
	[-0.661814, -0.133354, -0.737712], [-0.579521, -0.150880, -0.800868],
	[ 0.346897,  0.127467, -0.929201], [ 0.201528,  0.191210, -0.960638],
	[ 0.288738,  0.208962, -0.934326], [ 0.313420,  0.305039, -0.899288],
	[ 0.047802,  0.252840, -0.966326], [-0.106147,  0.309387, -0.944993],
	[-0.019450,  0.332267, -0.942985], [-0.000236,  0.429325, -0.903150],
	[ 0.136870,  0.272499, -0.952371], [ 0.089017,  0.448301, -0.889439],
	[ 0.158871,  0.370796, -0.915025], [ 0.248305,  0.385678, -0.888593],
	[ 0.107557,  0.536442, -0.837055], [ 0.211240,  0.625845, -0.750796],
	[ 0.196321,  0.548484, -0.812787], [ 0.267147,  0.475041, -0.838432],
	[ 0.348009, -0.509982, -0.786643], [ 0.204165, -0.470843, -0.858268],
	[ 0.306882, -0.449713, -0.838798], [ 0.364005, -0.361769, -0.858268],
	[ 0.055341, -0.416881, -0.907275], [-0.090756, -0.351740, -0.931688],
	[ 0.008845, -0.343190, -0.939224], [ 0.061824, -0.255946, -0.964712],
	[ 0.158462, -0.403205, -0.901286], [ 0.164646, -0.241277, -0.956387],
	[ 0.214764, -0.314722, -0.924568], [ 0.317683, -0.294553, -0.901286],
	[ 0.215792, -0.150880, -0.964712], [ 0.360619, -0.043725, -0.931688],
	[ 0.316346, -0.133354, -0.939224], [ 0.368013, -0.203516, -0.907275],
	[-0.249395,  0.211636, -0.944993], [-0.252866,  0.047667, -0.966326],
	[-0.302308,  0.139247, -0.942985], [-0.399691,  0.156740, -0.903150],
	[-0.251522, -0.117947, -0.960638], [-0.245163, -0.276550, -0.929201],
	[-0.299846, -0.192683, -0.934326], [-0.398310, -0.180640, -0.899288],
	[-0.303643, -0.028104, -0.952371], [-0.449659, -0.090608, -0.888593],
	[-0.403196, -0.012755, -0.915025], [-0.449900,  0.080548, -0.889439],
	[-0.539741, -0.075573, -0.838432], [-0.659784,  0.031464, -0.750796],
	[-0.582308,  0.017153, -0.812787], [-0.538734,  0.095418, -0.837055],
	[-0.142065, -0.276550, -0.950441], [-0.148649, -0.117947, -0.981831],
	[-0.093860, -0.192683, -0.976762], [ 0.010435, -0.180640, -0.983494],
	[-0.149663,  0.047667, -0.987587], [-0.144419,  0.211636, -0.966620],
	[-0.095021,  0.139247, -0.985688], [ 0.010177,  0.156740, -0.987587],
	[-0.097504, -0.028104, -0.994838], [ 0.061717,  0.080548, -0.994838],
	[ 0.008704, -0.012755, -0.999881], [ 0.061830, -0.090608, -0.993965],
	[ 0.164023,  0.095418, -0.981831], [ 0.309310,  0.031464, -0.950441],
	[ 0.213641,  0.017153, -0.976762], [ 0.164403, -0.075573, -0.983494]
]

def run(file):
#   sys.stdout.write('Simulation starting time: %s\n' % time.ctime(time.time()))
    Simobj = Linus(sys.argv[1])
    prot = Simobj.protein
    numres = prot.num_residues
    atoms = prot.atoms
    fai = prot.residue_first_atom_indices

    hbparms = Simobj.hbdpar
    hbparms['use_hbond'] = 1
    hbparms['use_sidechain_hbond'] = 1
    wmin = 2
    wmax = numres-1       
    hblist = make_hbond_list(prot, wmin, wmax, hbparms)

    numint, numvirt, bbtot, numbb,wat_list,cntmult_list = mk_virt_bb_cntmult_loosHbd(prot,hblist)
    bb_nonbd = bbtot - numbb

    sys.stdout.write('COMPND numintHbd     numSolv    num_nonHbd      total_bb_polar\n')
    sys.stdout.write('COMPND   %4d         %4d        %4d             %4d\n' % \
                    (numint, numvirt,bb_nonbd,bbtot))

    asalist  = make_asa_list(prot)
#   wat_list = None
    atot,tot_solv_nrg = print_hbs_chasa(sys.stdout, prot, 
          asalist, ext_atoms=wat_list,solv_list=cntmult_list)
    int_solv_nrg = tot_solv_nrg - (2.5 * numint)   
    print 'Internal Hbond + Solvation Energy = ', int_solv_nrg

#   sys.stdout.write('Simulation Finished at: %s\n' % time.ctime(time.time()))

if __name__ == "__main__":
    import sys    

    if len(sys.argv) &lt; 2:
        print USAGE 
        sys.exit()

    run(file)

</pre></body></html>