path: root/cad/src/model/bond_constants.py

# Copyright 2005-2008 Nanorex, Inc.  See LICENSE file for details. 
"""
bond_constants.py -- constants and simple functions for use with class Bond
(which can be defined without importing that class).

@author: Bruce
@version: $Id$
@copyright: 2005-2008 Nanorex, Inc.  See LICENSE file for details. 

History:

050429 - Started out as part of bonds.py. Gradually extended.

050707 - Split into separate file, largely to avoid recursive import problems
(since these constants need to be imported by many bond-related modules).
Many of them are still imported via bonds module, by code in other modules.

050920 - Full mmp support for Carbomeric bonds.
[FYI: As of today, sim executable reportedly accepts them and uses same params
as for bond2.]
"""

from math import floor, ceil

from geometry.VQT import Q
from utilities.debug import print_compact_traceback    
from utilities import debug_flags
import foundation.env as env
from simulation.PyrexSimulator import thePyrexSimulator

# ==

# Bond valence constants -- exact ints, 6 times the numeric valence they represent.
# If these need an order, their standard order is the same as the order of their numeric valences
# (as in the constant list BOND_VALENCES).

V_SINGLE     = 6 * 1   #  6
V_GRAPHITE   = 6 * 4/3 #  8 (this can't be written 6 * (1+1/3) or 6 * (1+1/3.0) - first one is wrong, second one is not an exact int)
V_AROMATIC   = 6 * 3/2 #  9
V_DOUBLE     = 6 * 2   # 12
V_CARBOMERIC = 6 * 5/2 # 15 for the bonds in a carbomer of order 2.5 (which alternate with aromatic bonds); saved as bondc as of 050920
V_TRIPLE     = 6 * 3   # 18

V_UNKNOWN    = 6 * 7/6 #  7 not in most tables here, and not yet used; someday might be used internally by bond-type inference code

BOND_VALENCES = [V_SINGLE, V_GRAPHITE, V_AROMATIC, V_DOUBLE, V_CARBOMERIC, V_TRIPLE]
    # when convenient (e.g. after A8), V_GRAPHITE should be renamed to V_GRAPHITIC [bruce 060629]
BOND_MMPRECORDS = ['bond1', 'bondg', 'bonda', 'bond2', 'bondc', 'bond3']
    # (Some code might assume these all start with "bond".)
    # (These mmp record names are also hardcoded into mmp-reading code in files_mmp.py.)
bond_type_names = {V_SINGLE:'single', V_DOUBLE:'double', V_TRIPLE:'triple',
                   V_AROMATIC:'aromatic', V_GRAPHITE:'graphitic', V_CARBOMERIC:'carbomeric'}

BOND_VALENCES_HIGHEST_FIRST = list(BOND_VALENCES)
BOND_VALENCES_HIGHEST_FIRST.reverse()

V_ZERO_VALENCE = 0 # used as a temporary valence by some code

BOND_LETTERS = ['?'] * (V_TRIPLE+1) # modified just below, to become a string; used in initial Bond.draw method via bond_letter_from_v6

for v6, mmprec in zip( BOND_VALENCES, BOND_MMPRECORDS ):
    BOND_LETTERS[v6] = mmprec[4] # '1','g',etc
    # for this it's useful to also have '?' for in-between values but not for negative or too-high values,
    # so a list or string is more useful than a dict

assert BOND_LETTERS[V_CARBOMERIC] == 'c' # not 'a', not 'b'

BOND_LETTERS[0] = '0' # see comment in bond_letter_from_v6

BOND_LETTERS = "".join(BOND_LETTERS)
    ## print "BOND_LETTERS:",BOND_LETTERS # 0?????1?ga??2?????3

BOND_MIN_VALENCES = [ 999.0] * (V_TRIPLE+1) #bruce 050806; will be modified below
BOND_MAX_VALENCES = [-999.0] * (V_TRIPLE+1)

bond_valence_epsilon = 1.0 / 64 # an exact float; arbitrary, but must be less than 1/(2n) where no atom has more than n bonds

for v6 in BOND_VALENCES:
    if v6 % V_SINGLE == 0: # exact valence
        BOND_MIN_VALENCES[v6] = BOND_MAX_VALENCES[v6] = v6 / 6.0
    else:
        # non-integral (and inexact) valence
        BOND_MIN_VALENCES[v6] = floor(v6 / 6.0) + bond_valence_epsilon
        BOND_MAX_VALENCES[v6] = ceil(v6 / 6.0)  - bond_valence_epsilon
    pass

BOND_MIN_VALENCES[V_UNKNOWN] = 1.0 # guess, not yet used
BOND_MAX_VALENCES[V_UNKNOWN] = 3.0 # ditto

# constants returned as statuscode by Atom.directional_bond_chain_status()
# [bruce 071016]
DIRBOND_CHAIN_MIDDLE = 'middle'
DIRBOND_CHAIN_END = 'end'
DIRBOND_NONE = None
DIRBOND_ERROR = 'error'

# ==

# [Note: the functions atoms_are_bonded and find_bond, were moved
#  from bonds.py to here (to remove an import cycle) by bruce 071216;
#  I also removed the old alias 'bonded' for atoms_are_bonded,
#  since 'bonded' is too generic to be searched for.]

def atoms_are_bonded(a1, a2):
    """
    Are these atoms (or singlets) already directly bonded?
    [AssertionError if they are the same atom.]
    """
    #bruce 041119 #e optimized by bruce 050502 (which indirectly added "assert a1 is not a2")
    ## return a2 in a1.neighbors()
    return not not find_bond(a1, a2)

def find_bond(a1, a2):
    """
    If a1 and a2 are bonded, return their Bond object; if not, return None.
    [AssertionError if they are the same atom.]
    """
    #bruce 050502; there might be an existing function in some other file, to merge this with
    assert a1 is not a2
    for bond in a1.bonds:
        if bond.atom1 is a2 or bond.atom2 is a2:
            return bond
    return None

def find_Pl_bonds(atom1, atom2):
    #bruce 080409 moved this here to clear up an import cycle
    """
    return the two bonds in atom1-Pl5-atom2,
    or (None, None) if that structure is not found.
    """
    for bond1 in atom1.bonds:
        Pl = bond1.other(atom1)
        if Pl.element.symbol == 'Pl5': # avoid import of model.elements
            bond2 = find_bond(Pl, atom2)
            if bond2:
                return bond1, bond2
    return None, None

# ==

def min_max_valences_from_v6(v6):
    """
    Given a v6 value (an int representing a bond order times 6),
    return a range of acceptable floating point bond orders.  V_SINGLE
    returns 1.0, 1.0.  V_AROMATIC returns 1+delta, 2-delta, for small
    delta.  etc...
    """
    return BOND_MIN_VALENCES[v6], BOND_MAX_VALENCES[v6]

def valence_to_v6(valence): #bruce 051215
    """
    Given a valence (int or float, single bond is 1 or 1.0),
    return it as a v6, being careful about rounding errors.
    """
    return int(valence * V_SINGLE + 0.01) # kluge: 0.01 is based on knowledge of scale of V_SINGLE (must be > 0, < 1/V_SINGLE)

def bond_letter_from_v6(v6): #bruce 050705
    """
    Return a bond letter summarizing the given v6,
    which for legal values is one of 1 2 3 a g b,
    and for illegal values is one of - 0 ? +
    """
    try:
        ltr = BOND_LETTERS[v6]
            # includes special case of '0' for v6 == 0,
            # which should only show up for transient states that are never drawn, except in case of bugs
    except IndexError: # should only show up for transient states...
        if v6 < 0:
            ltr = '-'
        else:
            ltr = '+'
    return ltr

def btype_from_v6(v6): #bruce 050705
    """
    Given a legal v6, return 'single', 'double', etc.
    For illegal values, return 'unknown'.
    For V_CARBOMERIC this returns 'carbomeric', not 'aromatic'.
    """
    try:
        return bond_type_names[v6]
    except KeyError:
        if debug_flags.atom_debug:
            print "atom_debug: illegal bond v6 %r, calling it 'unknown'" % (v6,)
        return 'unknown' #e stub for this error return; should it be an error word like this, or single, or closest legal value??
    pass

def invert_dict(dict1): #bruce 050705
    res = {}
    for key, val in dict1.items():
        res[val] = key
    return res

bond_type_names_inverted = invert_dict(bond_type_names)

def v6_from_btype(btype): #bruce 050705
    """
    Return the v6 corresponding to the given bond-type name
    ('single', 'double', etc). Exception if name not legal.
    """
    return bond_type_names_inverted[btype]

_bond_arrows = {
    0: "<- ->".split(),
    1: "-- ->".split(), # from atom1 to atom2, i.e. to the right
   -1: "<- --".split(),
}
    
def bonded_atoms_summary(bond, quat = Q(1,0,0,0)): #bruce 050705; direction feature, bruce 070414. ###e SHOULD CALL bond_left_atom
    """
    Given a bond, and an optional quat describing the orientation it's shown in,
    order the atoms left to right based on that quat,
    and return a text string summarizing the bond
    in the form C26(sp2) <-2-> C34(sp3) or so,
    leaving out the < or > if the bond has a direction.
    """
    a1 = bond.atom1
    a2 = bond.atom2
    direction = bond._direction
    vec = a2.posn() - a1.posn()
    vec = quat.rot(vec)
    if vec[0] < 0.0:
        a1, a2 = a2, a1
        direction = - direction
    a1s = describe_atom_and_atomtype(a1)
    a2s = describe_atom_and_atomtype(a2)
    bondletter = bond_letter_from_v6(bond.v6)
    if bondletter == '1':
        bondletter = ''
    arrows = _bond_arrows.get(direction, ("<-", " (invalid direction) ->"))
    return "%s %s%s%s %s" % (a1s, arrows[0], bondletter, arrows[1], a2s)
    
def bond_left_atom(bond, quat = Q(1,0,0,0)): #bruce 070415, modified from bonded_atoms_summary, which ought to call this now ##e
    # TODO: make this method name clearer: bond_leftmost_atom? bond_get_atom_on_left? [bruce 080807 comment]
    """
    Given a bond, and an optional quat describing the orientation it's shown in,
    order the atoms left to right based on that quat
    (i.e. as the bond would be shown on the screen using it),
    and return the leftmost atom.
    """
    a1 = bond.atom1
    a2 = bond.atom2
    vec = a2.posn() - a1.posn()
    vec = quat.rot(vec)
    if vec[0] < 0.0:
        a1, a2 = a2, a1
    return a1

def describe_atom_and_atomtype(atom): #bruce 050705, revised 050727 #e refile?
    """
    If a standard atom, return a string like C26(sp2) with atom name and
    atom hybridization type, but only include the type if more than one is 
    possible for the atom's element and the atom's type is not the default 
    type for that element.
    
    If a PAM Ss or Sj atom, returns a string like Ss28(A) with atom name
    and dna base name. 
    
    @deprecated: For some purposes use L{Atom.getInformationString()} instead.
                 (But for others, that might not be suitable. Needs review.)
    """
    res = str(atom)
    if atom.atomtype is not atom.element.atomtypes[0]:
        res += "(%s)" % atom.atomtype.name
    if atom.getDnaBaseName():
        res += "(%s)" % atom.getDnaBaseName()
    
    return res

# ==

_bond_params = {} # maps triple of atomtype codes and v6 to (rcov1, rcov2) pairs

def bond_params(atomtype1, atomtype2, v6 = V_SINGLE):
    #bruce 060324 for bug 900; made v6 optional, 080405
    """
    Given two atomtypes and an optional bond order encoded as v6,
    look up or compute the parameters for that kind of bond.
    For now, the return value is just a pair of numbers, rcov1 and rcov2,
    for use as the covalent radii for atom1 and atom2 respectively for this kind of bond
    (with their sum adjusted to the equilibrium bond length if this is known).
    """
    atcode1 = id(atomtype1)
        # maybe: add a small int code attr to atomtype, for efficiency
    atcode2 = id(atomtype2)
    try:
        return _bond_params[(atcode1, atcode2, v6)]
    except KeyError:
        res = _bond_params[(atcode1, atcode2, v6)] = \
              _compute_bond_params(atomtype1, atomtype2, v6)
        return res
    pass

def ideal_bond_length(atom1, atom2, v6 = V_SINGLE): #bruce 080404
    """
    Return the ideal bond length between atom1 and atom2 (using their
    current atomtypes) for a bond of the given order (default V_SINGLE).
    The atoms need not be bonded; if they are, that bond's current order
    is ignored. [review: revise that rule, to use it if they are bonded??]
    
    (Use bond_params, which uses getEquilibriumDistanceForBond if ND-1
     is available.)
    """
    rcov1, rcov2 = bond_params(atom1.atomtype, atom2.atomtype, v6)
    return rcov1 + rcov2

def _compute_bond_params(atomtype1, atomtype2, v6): #bruce 080405 revised this
    """
    [private helper function for bond_params]
    """
    # note: this needn't be fast, since its results for given arguments
    # are cached for the entire session by the caller.
    # (note: as of 041217 rcovalent is always a number; it's 0.0 for Helium,
    #  etc, so for nonsense bonds like He-He the entire bond is drawn as if
    #  "too long". [review: what should we do for a nonsense bond like C-He?])
    rcov1 = atomtype1.rcovalent
    rcov2 = atomtype2.rcovalent
    rcovsum = rcov1 + rcov2 # ideal length according to our .rcovalent tables,
        # used as a fallback if we can't get a better length
    if not rcovsum:
        print "error: _compute_bond_params for nonsense bond:", \
              atomtype1, atomtype2, v6
        rcov1 = rcov2 = 0.5 # arbitrary
        rcovsum = rcov1 + rcov2
    # now adjust rcov1 and rcov2 to make their sum the equilibrium bond length
    # [bruce 060324 re bug 900]
    eltnum1 = atomtype1.element.eltnum
        # note: both atoms and atomtypes have .element
    eltnum2 = atomtype2.element.eltnum
    ltr = bond_letter_from_v6(v6)
    
    # determine ideal bond length (special case for one element being Singlet)
    assert eltnum1 or eltnum2, "can't bond bondpoints to each other"
    if eltnum1 == 0 or eltnum2 == 0:
        # one of them is Singlet (bondpoint); getEquilibriumDistanceForBond
        # doesn't know about those, so work around this somehow
        nicelen = None
        if atomtype1.element.pam or atomtype2.element.pam:
            # in this case we prefer the native formula (see if this fixes bug 2944) [bruce 081210]
            nicelen = rcovsum
        if nicelen is None:
            # use half the distance for a bond from the non-Singlet atomtype to itself
            eltnum = eltnum1 or eltnum2
            nicelen = _safely_call_getEquilibriumDistanceForBond( eltnum, eltnum, ltr)
            if nicelen:
                nicelen = nicelen / 2.0
            pass
        pass
    else:
        nicelen = _safely_call_getEquilibriumDistanceForBond( eltnum1, eltnum2, ltr)
    if nicelen is None:
        # the call failed, use our best guess
        nicelen = rcovsum
    assert nicelen > 0.0
    # now adjust rcov1 and rcov2 so their sum equals nicelen
    # (this works even if one of them is 0, presumably for a bondpoint)
    ratio = nicelen / float(rcovsum)
    rcov1 *= ratio
    rcov2 *= ratio
    
    return rcov1, rcov2

def _safely_call_getEquilibriumDistanceForBond( eltnum1, eltnum2, ltr): #bruce 080405 split this out
    """
    #doc
    eg: args for C-C are (6, 6, '1')

    @return: ideal length in Angstroms (as a positive float),
             or None if the call of getEquilibriumDistanceForBond failed
    """
    try:
        pm = thePyrexSimulator().getEquilibriumDistanceForBond(eltnum1,
                                                               eltnum2,
                                                               ltr)
        assert pm > 2.0, "too-low pm %r for getEquilibriumDistanceForBond%r" % \
               (pm, (eltnum1, eltnum2, ltr))
            # 1.0 means an error occurred; 2.0 is still ridiculously low
            # [not as of 070410]; btw what will happen for He-He??
            # update 070410: it's 1.8 for (202, 0, '1').
            # -- but a new sim-params.txt today changes the above to 170
        nicelen = pm / 100.0 # convert picometers to Angstroms
        return nicelen
    except:
        # be fast when this happens a lot (not important now that our retval
        # is cached, actually; even so, don't print too much)
        if not env.seen_before("error in getEquilibriumDistanceForBond"):
            #e include env.redraw_counter to print more often? no.
            msg = "bug: ignoring exceptions when using " \
                  "getEquilibriumDistanceForBond, like this one: "
            print_compact_traceback(msg)
        return None
    pass
    
# ==

# Here's an old long comment which is semi-obsolete now [050707], but which motivates the term "v6".
# Note that I'm gradually replacing the term "bond valence" with whichever of "bond order" or "bond type"
# (related but distinct concepts) is appropriate. Note also that all the bond orders we deal with in this code
# are "structural bond orders" (used by chemists to talk about bonding structure), not "physical bond orders"
# (real numbers related to estimates of occupancy of molecular orbitals by electrons).

#bruce 050429: preliminary plan for higher-valence bonds (might need a better term for that):
#
# - Bond objects continue to compare equal when on same pair of atoms (even if they have a
# different valence), and (partly by means of this -- probably it's a kluge) they continue
# to allow only one Bond between any two atoms (two real atoms, or one real atom and one singlet).
#
# - I don't think we need to change anything basic about "internal vs external bonds",
# coordinates, basic inval/draw schemes (except to properly draw new kinds of bonds),
# etc. (Well, not due to bond valence -- we might change those things for other reasons.)
#
# - Each Bond object has a valence. Atoms often sum the valences of their bonds
# and worry about this, but they no longer "count their bonds" -- at least not as a
# substitute for summing the valences. (To prevent this from being done by accident,
# we might even decide that their list of bonds is not really a list, at least temporarily
# while this is being debugged. #?)
#
# This is the first time bonds have any state that needs to be saved,
# except for their existence between their two atoms. This will affect mmpfile read/write,
# copying of molecules (which needs rewriting anyway, to copy jigs/groups/atomsets too),
# lots of things about depositMode, maybe more.
#
# - Any bond object can have its valence change over time (just as the coords,
# elements, or even identities of its atoms can also change). This makes it a lot
# easier to write code which modifies chemical structures in ways which preserve (some)
# bonding but with altered valence on some bonds.
#
# - Atoms might decide they fit some "bonding pattern" and reorder
# their list of bonds into a definite order to match that pattern (this is undecided #?).
# This might mean that code which replaces one bond with a same-valence bond should do it
# in the same place in the list of bonds (no idea if we even have any such code #k).
#
# - We might also need to "invalidate an atom's bonding pattern" when we change anything
# it might care about, about its bonds or even its neighboring elements (two different flags). #?
#
# - We might need to permit atoms to have valence errors, either temporarily or permanently,
# and keep track of this. We might distinguish between "user-permitted" or even "user-intended"
# valence errors, vs "transient undesired" valence errors which we intend to automatically
# quickly get rid of. If valence errors can be long-lasting, we'll want to draw them somehow.
# 
# - Singlets still require exactly one bond (unless they've been killed), but it can have
# any valence. This might affect how they're drawn, how they consider forming new bonds
# (in extrude, fuse chunks, depositMode, etc), and how they're written into sim-input mmp files.
#
# - We represent the bond valence as an integer (6 times the actual valence), since we don't
# want to worry about roundoff errors when summing and comparing valences. (Nor to pay the speed
# penalty for using exactly summable python objects that pretend to have the correct numeric value.)
#
# An example of what we don't want to have to worry about:
#
#   >>> 1/2.0 + 1/3.0 + 1/6.0
#   0.99999999999999989
#   >>> _ >= 1.0
#   False
#
# We do guarantee to all code using these bond-valence constants that they can be subtracted
# and compared as numbers -- i.e. that they are "proportional" to the numeric valence.
# Some operations transiently create bonds with unsupported values of valence, especially bonds
# to singlets, and this is later cleaned up by the involved atoms when they update their bonding
# patterns, before those bonds are ever drawn. Except for bugs or perhaps during debugging,
# only standard-valence bonds will ever be drawn, or saved in files, or seen by most code.

# ==

# end