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

# Copyright 2005-2007 Nanorex, Inc.  See LICENSE file for details. 
"""
pi_bond_sp_chain.py -- geometric info for individual pi bonds, or chains
of them connected by sp atoms.

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

Note, 070414: it turns out a lot of the same concepts and similar code
ought to be useful for keeping track of chains of "directional bonds"
(such as pseudo-atom DNA backbones) in which direction-inference
and consistency checks should occur. Maybe I'll even create a general class
for perceiving bond-chains (or rings), with subclasses which know
about specific ways of recognizing bonds they apply to; or a class
whose instances hold bond-classifiers and act as perceivers of chains
of bonds of that type. (Once the chain is perceived, by generic code,
it does need specific code for its chain-type to figure out what to
do, though the basic fact of being destroyed when involved-atom-structure
changes is probably common. In the bond direction case, unlike in PiBondSpChain,
destroying the perceived strand will leave behind some persistent state,
namely the per-bond directions themselves; the perceived strand's role
is to help change those directions in an organized way.)
"""

import math
from Numeric import dot

from utilities import debug_flags
from model.jigs import Jig
from geometry.VQT import V, Q, A, cross, vlen, norm, twistor_angle

from operations.bond_chains import grow_bond_chain

from model.bond_constants import V_SINGLE
from model.bond_constants import V_DOUBLE
from model.bond_constants import V_TRIPLE
from model.bond_constants import V_AROMATIC
from model.bond_constants import V_GRAPHITE
from model.bond_constants import V_CARBOMERIC

DFLT_OUT = V(0.0, -0.6, 0.8) # these are rotated from standard out = V(0,0,1), up = V(0,1,0) but still orthonormal
DFLT_UP =  V(0.0,  0.8, 0.6)

# _recompile_counter is useful for development, and is harmless otherwise
try:
    _recompile_counter
except:
    _recompile_counter = 1
else:
    # increment the counter only when the module got recompiled, not just when it got reloaded from the same .pyc file
    if __file__.endswith(".py"):
        _recompile_counter += 1
        print "_recompile_counter incremented to", _recompile_counter
    else:
        pass ## print "_recompile_counter unchanged at", _recompile_counter
    pass

def bond_get_pi_info(bond, **kws):
    """
    #doc; **kws include out, up, abs_coords
    """
    obj = bond.pi_bond_obj
    if obj is not None:
        # let the obj update itself and return the answer.
        # (I think this should always work, since it gets invalled and can remove itself as needed.
        #  If not, it must resort to similar code to remainder of this routine,
        #  since it can't return anything here to tell this routine to do that for it.)
        if obj._recompile_counter != _recompile_counter:
            # it would not be good to do this except when the module got recompiled due to being modified
            # (ie not merely when it got reloaded), since that would defeat persistence of these objects, which we need to debug
            print "destroying %r of obsolete class" % obj # happens when we reload this module at runtime, during development
            obj.destroy()
            obj = None
    # special case for triple bonds [bruce 050728 to fix bug 821]: they are always arbitrarily oriented and not twisted
    if bond.v6 == V_TRIPLE:
        return pi_vectors(bond, **kws)
    # not worth optimizing for single bonds, since we're never called for them for drawing,
    # and if we're called for them when depositing sp2-sp-sp2 to position singlets on last sp2,
    # it would not even be correct to exclude them here!
    if obj is not None:
        return obj.get_pi_info(bond, **kws) # no need to pass an index -- that method can find one on bond if it stored one
    # if the obj is not there, several situations are possible; we only store an obj if the pi bond chain length is >1.
    if not bond.potential_pi_bond():
        return None # note, this does *not* happen for single bonds, if they could in principle be changed to higher bond types
    if pi_bond_alone_in_its_sp_chain_Q(bond):
        # optimize this common case, with simpler code and by not storing an object for it
        # (since too many atoms would need to invalidate it)
        return pi_vectors(bond, **kws)
    obj = bond.pi_bond_obj = make_pi_bond_obj(bond) # make an obj that updates itself as needed when atoms/bonds change or move
    return obj.get_pi_info(bond, **kws)

class PerceivedStructureType(Jig): #e rename, not really a Type (eg not a pattern), more like a set or layer or record
    """
    ... Subclasses are used for specific kinds of perceived structures
    (e.g. sp chains, pi systems, diamond or graphite rings).
    """
    pass

#bruce 060223: How will we (PiBondSpChain) deal with Undo? We could just store every attr in this object, but we don't want to bother.
# We want to think of it as entirely derived from other state (as it in fact is). So if things change, we'd rather invalidate it
# and start over. In theory that should be ok, provided Undo (when changing other atoms & bonds) invals them in all the same
# ways their LL change methods would do, including telling this object they changed (which makes it destroy itself).
# Since this is a Jig, it will at least participate in Jigs' scheme for having atom lists and being on atoms' jig lists,
# and Undo ought to properly handle that... how does that interact with our own policy of destroying ourself as soon as some atom
# we're interested in says it changed? Could that happen with atom->jig and jig->atom connections not yet corresponding?
# Solution to that worry: first let Undo do all its setattrs (so atom->jig and jig->atom connections have all been changed),
# then call all object updaters at once. Ours is our superclass's and does nothing (I think); the one on the first atom we
# encounter will destroy us; this will remove us from an atom's jig list but that should not (I think) invalidate the atom,
# since those are really just "back pointers" rather than properties of the atom (if it did inval it, that would be bad, I think --
# maybe not, since _undo_update is really an inval method and those can in general call other inval methods).
#  So if Atom & Bond _undo_updates are correct, we should be ok. If not, some bugs will be noticed at some point.
#
# One other thing: we're a Jig, but we're not in the node tree (which is why we don't get into mmp files). So Undo will not find us
# as a child object. We'll still get an objkey (I think)... hmm, this might be a problem, since when restoring some atom->jig to us,
# it won't restore our pointer to that atom. Can that ever happen (in light of our atomset being fixed)? I don't know. I think this
# might indeed cause bugs, so we might need to get counted as a child, or be some new undo-kind of object that easily gets
# invalidated by any _undo_update on atoms it touches (maybe we introduce "subscribe to another obj's _undo_update" or to their
# specific attrs?). I think it's easiest to wait and see if bugs happen. If they do, that subscription idea sounds best to me.
# CHANGED MY MIND, I'll just make us an S_CHILD of our bonds... then our atomset should be ok, and we could have _undo_update
# notice it changing too, if we wanted to. Let's do that for safety (maybe we'll never know if it was needed).
# ###@@@

class PiBondSpChain(PerceivedStructureType):
    """
    Records one chain (or ring) of potential-pi-bonds connected by -sp- atoms;
    as a Jig, sits on all atoms whose structure or motion matters for the extent of this chain
    and the geometry of its pi orbitals (ie all atoms in the chain, and the immediate neighbor atoms of the ends).
       Main job is to calculate pi-orbital angles for drawing the pi bonds.
        Someday we might extend this to also infer bond types along the chain, etc.
    """
    ringQ = False # class constant; in the Ring subclass, this is True
    have_geom = False # initial value of instance variable
    def _undo_update(self): #bruce 060223; see long comment above
        """
        our atomset (a member of our superclass Jig) must have changed...
        """
        self.destroy()
        PerceivedStructureType._undo_update(self) # might be pointless after that, but you never know...
        return
    def __init__(self, listb, lista): ### not like Jig init, might be a problem eg for copying it ###@@@ review whether that happens
        """
        Make one, from lists of bonds and atoms (one more atom than bond,
        even for a ring, where 1st and last atoms are same)
        """
        self._recompile_counter = _recompile_counter
        self.listb = listb
        self.lista = lista
        assert len(lista) == 1 + len(listb)
        assert len(listb) > 0
        #e could assert each bond's atoms are same as prior and next corresponding atom in lista
        # now add ourselves into each bond, exclusively (for cleaner code, this should be a separate method, called by init caller #e)
        for i, bond in zip( range(len(listb)), listb ):
            if bond.pi_bond_obj is not None:
                if debug_flags.atom_debug:
                    print "atom_debug: bug: obs pi_bond_obj found (and discarded) on %r" % (bond,)
                bond.pi_bond_obj.destroy()
            bond.pi_bond_obj = self
            bond.pi_obj_memo = i
            # now inval bond?? no, that would be wrong -- this func is part of updating it after smth else invalled it earlier.
        # figure out which atoms we need to monitor, for invalidating our geometric info or our structure.
        all_atoms = list(lista) # atoms in the chain, plus neighbor atoms whose posns matter
        assert len(lista) >= 2
        end1 = lista[0]
        end2 = lista[-1]
        if end1 is not end2: # not a ring
            nn = end1.neighbors()
            nn.remove(lista[1])
            all_atoms.extend(nn) # order doesn't matter
            nn = end2.neighbors()
            nn.remove(lista[-2])
            all_atoms.extend(nn)
        assy = all_atoms[0].molecule.assy #e need to zap the need for assy from Node
        PerceivedStructureType.__init__(self, assy, all_atoms)
        return
    super = PerceivedStructureType # needed for this to work across reloads
    listb = () # permit repeated destroy [bruce 060322]
    _recompile_counter = None # ditto (not sure if needed)
    def destroy(self):
        super = self.super
        for bond in self.listb:
            self.listb = () #bruce 060322 to save RAM (might not be needed)
            if bond.pi_bond_obj is self:
                bond.pi_bond_obj = None
                bond.pi_obj_memo = None # only needed to cause an exception if something tries to misuse it
            elif bond.pi_bond_obj is not None:
                if debug_flags.atom_debug:
                    print "atom_debug: bug: wrong pi_bond_obj %r found (not changed) on %r as we destroy %r" % (bond.pi_bond_obj, bond, self)
        super.destroy(self)
    def changed_structure(self, atom):
        # ideally this should react differently to neighbor atoms and others; but even neighbors might be brought in, eg C(sp3)-C#C;
        # so it's simplest for now to just forget about it and rescan for it when next needed.
        self.destroy() ###k need to make sure it's ok for jigs (Nodes) not in the model tree
        # no point in calling super.changed_structure after self.destroy!
    def moved_atom(self, atom):
        """
        some atom I'm on moved; invalidate my geometric info
        """
        # must be fast, so don't bother calling Jig.moved_atom, since it's a noop
        self.have_geom = False
    def changed_bond_type(self, bond): #bruce 050726
        """
        some bond I'm on changed type (perhaps due to user change
        or to bond inference code); invalidate my geometric info
        """
        # there is no superclass method for this -- it's only called on bond.pi_bond_obj for our own bonds
        ###e it might be worth only invalidating if the bond type differs from what we last used when recomputing
        self.have_geom = False        
    def get_pi_info(self, bond, out = DFLT_OUT, up = DFLT_UP, abs_coords = False):
        if len(self.listb) == 1:
            if debug_flags.atom_debug:
                print "atom_debug: should never happen (but should work if it does): optim for len1 PiBondSpChain object %r" % self
            return pi_vectors(self.listb[0], out = out, up = up, abs_coords = abs_coords)
                # optimization; should never happen since it's done instead of creating this object
        if not self.have_geom:
            self.recompute_geom(out = out, up = up) # always in abs coords
            self.have_geom = True
        ind = bond.pi_obj_memo # this was stored again whenever our structure changed (for now, only in init)
        assert type(ind) == type(1)
        assert 0 <= ind < len(self.listb)
        assert self.listb[ind] is bond
        biL_pvec, biR_pvec = self.pvecs_i( ind, abs_coords = True)
            #bruce 050727 -- abs_coords arg was used both here & below -- wrong, using True here is right (so I fixed it),
            # but I didn't yet see any effect from that presumed bug. It might be that the chunk quat is
            # always Q(1,0,0,0) whenever the pi_info gets updated.
        bond_axis = self.axes[ind]
            ## bond_axis = self.lista[ ind+1 ].posn() - self.lista[ ind ].posn() # order matters, i think
        if bond.atom1 is self.lista[ind]:
            assert bond.atom2 is self.lista[ind+1]
        else:
            # flipped
            assert bond.atom1 is self.lista[ind+1]
            assert bond.atom2 is self.lista[ind]
            bond_axis = - bond_axis
            biL_pvec, biR_pvec = biR_pvec, biL_pvec
        return pi_info_from_abs_pvecs( bond, bond_axis, biL_pvec, biR_pvec, abs_coords = abs_coords) 

    def pvecs_i(self, i, abs_coords = False):
        """
        Use last recomputed geom to get the 2 pvecs for bond i...
        in the coordsys of the bond, or always abs if flag is passed.
        WARNING: pvec order matches atom order in lista, maybe not in bond.
        """
        biL_pvec = biR_pvec = self.quats_cum[i].rot(self.b0L_pvec) # not yet twisted
        axis = self.axes[i]
##        if self.twist90 and i % 2 == 1:
##            biL_pvec = norm(cross(biL_pvec, axis))
##            biR_pvec = norm(cross(biR_pvec, axis))
        twist = self.twist # this is the twist angle (in radians), for each bond; or None or 0 or 0.0 for no twist
        if twist:
            # twist that much times i and i+1 for L and R, respectively, around axes[i]
            theta = twist
            quatL = Q(axis, theta * (i))
            quatR = Q(axis, theta * (i+1))
            biL_pvec = quatL.rot(biL_pvec)
            biR_pvec = quatR.rot(biR_pvec)
        if not abs_coords:
            # put into bond coords (which might be the same as abs coords)
            bond = self.listb[i]
            quat = bond.bond_to_abs_coords_quat()
            biL_pvec = quat.unrot(biL_pvec)
            biR_pvec = quat.unrot(biR_pvec)
        return biL_pvec, biR_pvec # warning: these two vectors might be the same object

    def recompute_geom(self, out = DFLT_OUT, up = DFLT_UP):
        """
        Compute and store enough geometric info for pvecs_i to run (using submethods for various special cases):
        self.quats_cum, self.twist, self.b0L_pvec.
        """
        # put coords on each bond axis, then link these (gauge transform)... start with first and propogate.
        # warning: each bond in turn might have a different coordinate system!
        # (if they alternate external/internal, different chunks)
        # and the geom we compute for each bond should be in that bond's system, to help draw it properly.
        # fortunately they're only linked by rotations... so we should be able to compensate for coord changes.
        # This also means -- when a bond is invalidated, it might be from changing its coord sys, even if no atom motion
        # or structure change. This needs to invalidate it here, too, or (more simply) inval this whole thing.
        # or we could let the pi info be in abs coords and rotate the vecs before use?
        # yes, this is effectively what we ended up doing, except even the abs coords get computed each time they're asked for
        # from the quats stored by this routine.
        
        self.twist = None # default value
        self.axes = self.quats_cum = self.chain_quat_cum = None # invalid values (for catching bugs)

        bonds = self.listb
        lista = self.lista
        
##        # first, do we want that 90-deg offset? Only if double bonds, I think.
##        # (It's just a kluge to avoid returning pi orders of (1,0) for some of them and (0,1) for the others.
##        #  We use it so we can return pi orders of (1,0) for all double bonds.)
##        #
##        # Let's assume bond inference was done, so if one bond is double, they all are.
##        # BUT WAIT, are all the middle atoms C(sp), not some other element?? I think so... ###@@@  check this somewhere,
##        #  and if another is poss, probably don't even make it form a chain like this
##        self.twist90 = ( bonds[0].v6 == V_DOUBLE )

        posns = A( map( lambda atom: atom.posn(), lista ) )
        self.axes = axes = posns[1:] - posns[:-1] # axes along bonds, all in consistent directions & abs coords
        quats_incr = map( lambda (axis1, axis2): Q(axis1, axis2), zip( axes[:-1], axes[1:] ) )
            # how to transform vectors perp to one bond, to those perp to the next
        # cumulative quats (there might be a more compact python notation for this, and/or a way to do more of it in Numeric)
        qq = Q(1,0,0,0)
        quats_cum = [qq] # how to get from bonds[0] to each bonds[i] starting from i == 0, using only bend-projections at atoms
        for qi, i in zip( quats_incr, range(len(quats_incr)) ):
            # qi tells us how to map (perps to axes[i]) to (perps to axes[i+1])
            qq = qq + qi # make a new quat, don't use +=
            ###e ideally we'd now enforce qq turning axes[0] into axes[i], to compensate for cumulative small errors [nim]
            if self.adjacent_double_bonds( i, i+1 ):
                axis = axes[i+1]
                theta = math.pi / 2 # 90 degrees
                qq += Q(axis, theta)
            quats_cum.append(qq)
        assert len(quats_cum) == len(bonds)
        self.quats_cum = quats_cum
        self.chain_quat_cum = quats_cum[-1]
        if self.ringQ:
            # for rings we need to know how to get all the way around the ring, ie from bonds[-1] to bonds[0]
            # (but we won't add these quats to our lists)
            self.ringquat_incr = Q( axes[-1], axes[0] )
            if self.adjacent_double_bonds( -1, 0 ):
                self.ringquat_incr += Q(axes[0], math.pi / 2)
            ## self.ringquat_cum = self.chain_quat_cum + self.ringquat_incr #k probably not used
        # now we know total twist, so we can use code similar to   pi_vectors   function to decide what to do at the ends.
        # BTW all this only mattered if we weren't alternating single-triple bonds... ideally we'd rule that out first.
        # but no harm if we don't.
        # BTW2, at what point should this code do any bond-type inference? and how -- should it take most recently user-changed
        # bond in it, or most recently inferred end-bond, and propogate that, in some cases? btw3, remember it might be a ring.

        # call subclass-specific funcs to do the rest -- compute the twist and start/end p orbital vectors, used later by pvecs_i...
        self.out = out
        self.up = up
        self.recompute_geom_from_quats()
        return # from recompute_geom

    def adjacent_double_bonds( self, i1, i2 ): # replaces self.twist90 code, 050726
        """
        #doc;
        i1 might be -1
        """
        bonds = self.listb
        v1 = bonds[i1].v6 # ok for negative indices i1
        v2 = bonds[i2].v6
        #bruce 050728 kluge bugfix: treat V_SINGLE as V_DOUBLE so this gives user benefit of the doubt if they
        # just haven't yet bothered to put in the correct bondtypes, and so this works from depositmode growing an sp chain
        # when it's making new singlets on an sp2 atom it deposits at the end of the sp chain.
        if v1 == V_SINGLE:
            v1 = V_DOUBLE
        if v2 == V_SINGLE:
            v2 = V_DOUBLE
        # if both are double, or one is double and one is aromatic or graphite (but not carbomeric)
        if v1 == V_DOUBLE and v2 in (V_DOUBLE, V_AROMATIC, V_GRAPHITE):
            return True
        if v2 == V_DOUBLE and v1 in (V_AROMATIC, V_GRAPHITE):
            return True
        return False
    
    def recompute_geom_from_quats(self): # non-ring class
        """
        [not a ring; various possible end situations]
        """
        out = self.out
        up = self.up
        atoms = self.lista
        bonds = self.listb
            # (I wish I could rename these attrs to "atoms" and "bonds",
            #  but self.atoms conflicts with the Jig attr, and self.chain_atoms is too long.)
        # default values:
        self.bm1R_pvec = "not needed" # and bug if used, we hope, unless we later change this value
        self.twist = None
        # Figure out p vectors for atom[0]-end (left end) of bond[0], and atom[-1]-end (right end) of bond[-1].
        # If both are determined from outside this chain (i.e. if the subrs here don't return None),
        # or if we are a ring (so they are forced to be the same, except for the projection due to bond bending at the ring-joining atom),
        # then [using self.recompute_geom_both_ends_constrained()] they must be made to match (up to an arbitrary sign),
        # using a combination of some twist (to be computed here) along each bond,
        # plus a projection and (in some cases, depending on bond types i think -- see self.twist90)
        # 90-degree turn at each bond-bond connection;
        # the projection part for all the bond-bond connections is already accumulated in self.chain_quat_cum.
        
        # note: the following p-vec retvals are in abs coords, as they should be
        #e rename the funcs, since they are not only for sp2, but for any atom that ends our chain of pi bonds
        pvec1 = p_vector_from_sp2_atom(atoms[0], bonds[0], out = out, up = up) # might be None
        pvec2 = p_vector_from_sp2_atom(atoms[-1], bonds[-1], out = out, up = up) # ideally, close to negative or positive of pvec1 ###@@@ handle neg
        # handle one being None (use other one to determine the twist) or both being None (use arb vectors)
        if pvec1 is None:
            if pvec2 is None:
                # use arbitrary vectors on left end of bonds[0], perp to bond and to out; compute differently if bond axis ~= out
                axis = self.axes[0]
                pvec = cross(out, axis)
                lenpvec = vlen(pvec)
                if lenpvec < 0.01:
                    # bond axis is approx parallel to out
                    pvec = cross(up, axis)
                    lenpvec = vlen(pvec) # won't be too small -- bond can't be parallel to both up and out
                pvec /= lenpvec
                self.b0L_pvec = pvec
            else:
                # pvec2 is defined, pvec1 is not. Need to transport pvec2 back to coords of pvec1
                # so our standard code (pvec_i, which wants pvec1, i.e. self.b0L_pvec) can be used.
                self.b0L_pvec = self.chain_quat_cum.unrot(pvec2)
        else:
            if pvec2 is None:
                self.b0L_pvec = pvec1
            else:
                # both vectors not None -- use recompute_geom_both_ends_constrained
                self.b0L_pvec = pvec1
                self.bm1R_pvec = pvec2
                self.recompute_geom_both_ends_constrained()
        return # from non-ring recompute_geom_from_quats

    def recompute_geom_both_ends_constrained(self):
        """
        Using self.b0L_pvec and self.chain_quat_cum and self.bm1R_pvec,
        figure out self.twist...
        [used for rings, and for linear chains with both ends constrained]
        """
        bonds = self.listb
        # what pvec would we have on the right end of the chain, if there were no per-bond twist?
        bm1R_pvec_no_twist = self.chain_quat_cum.rot( self.b0L_pvec ) # note, this is a vector perp to bond[-1]
        axism1 = self.axes[-1]
        # what twist is needed to put that vector over the actual one (or its neg), perhaps with 90-deg offset?
        i = len(bonds) - 1
##        if self.twist90 and i % 2 == 1:
##            bm1R_pvec_no_twist = norm(cross(axism1, bm1R_pvec_no_twist))
        # use negative if that fits better
        if dot(bm1R_pvec_no_twist, self.bm1R_pvec) < 0:
            bm1R_pvec_no_twist = - bm1R_pvec_no_twist
        # total twist is shared over every bond
        # note: we care which direction we rotate around axism1
        total_twist = twistor_angle( axism1, bm1R_pvec_no_twist, self.bm1R_pvec)
            # in radians; angular range is undefined, for now
            # [it's actually +- 2pi, twice what it would need to be in general],
            # so we need to coerce it into the range we want, +- pi/2 (90 degrees); here too we use the fact
            # that in this case, a twist of pi (180 degrees) is the same as 0 twist.
        while total_twist > math.pi/2:
            total_twist -= math.pi
        while total_twist <= - math.pi/2:
            total_twist += math.pi
        self.twist = total_twist / len(bonds)
        return

    pass # end of class PiBondSpChain

class RingPiBondSpChain( PiBondSpChain):
    """
    Class for a perceived ring of =C= or -C# atoms,
    or other 2-bond sp atoms if there are any
    (not sure if there can be)
    """
    ringQ = True
    # Note: I've assumed a ring means no constraints on angles... what if other bonds? but sp2 atoms count as endpoints,
    # so it's ok, this is just ring of sp, so it's true.
    def recompute_geom_from_quats(self): # ring class
        # get arb p-vector for use at bond[-1] right end
        self.bm1R_pvec = arb_perp_unit_vector(self.axes[-1])
        # then project that through ring-joining-atom to a pvec for use at bond[0] left end
        self.b0L_pvec = self.ringquat_incr.rot(self.bm1R_pvec)
        # then use same code as non-ring with constrained pvecs
        self.recompute_geom_both_ends_constrained( ) # fyi: one of several methods to finish up; chains can use this one or others
        return
    pass # end of class RingPiBondSpChain

def make_pi_bond_obj(bond): # see also find_chain_or_ring_from_bond, which generalizes this [bruce 071126]
    """
    #doc ...
    make one for a pi bond, extended over sp atoms in the middle...
    what if sp3-sp-sp3, no pi bonds? then no need for one!
    """
    atom1 = bond.atom1
    atom2 = bond.atom2
    # more atoms to the right? what we need to grow: sp atom, exactly one other bond, it's a pi bond. (and no ring formed)
    ringQ, listb1, lista1 = grow_pi_sp_chain(bond, atom1)
    if ringQ:
        assert atom2 is lista1[-1]
        return RingPiBondSpChain( [bond] + listb1 , [atom2, atom1] + lista1 )
    ringQ, listb2, lista2 = grow_pi_sp_chain(bond, atom2)
    assert not ringQ
    listb1.reverse()
    lista1.reverse()
    return PiBondSpChain( listb1 + [bond] + listb2, lista1 + [atom1, atom2] + lista2 ) # one more atom than bond

def sp_atom_2bonds(atom):
    """
    """
    ## warning: this being true is *not* enough to know that a pi bond containing atom has an sp-chain length > 1.
    return atom.atomtype.is_linear() and len(atom.bonds) == 2 and atom.atomtype.numbonds == 2

def next_bond_in_sp_chain(bond, atom):
    """
    Given a pi bond and one of its atoms, try to grow the chain
    beyond that atom... return next bond, or None.

    @note: This function not returning None (for either atom on the end
    of a potential pi bond) is the definitive condition for that bond
    not being the only one in its sp-chain.
    """
    assert bond.potential_pi_bond()
    if not sp_atom_2bonds(atom):
        return None #k retval
    bonds = atom.bonds
    # len(bonds) == 2, known from subr conds above
    if bond is bonds[0]:
        obond = bonds[1]
    else:
        assert bond is bonds[1]
        obond = bonds[0]
    if obond.potential_pi_bond():
        return obond
    return None

def pi_bond_alone_in_its_sp_chain_Q(bond):
    return next_bond_in_sp_chain(bond, bond.atom1) is None and next_bond_in_sp_chain(bond, bond.atom2) is None

def grow_pi_sp_chain(bond, atom): # WARNING: superceded by grow_pi_sp_chain_NEWER_BETTER, except that one is untested. see its comment.
    """
    Given a potential pi bond and one of its atoms,
    grow the pi-sp-chain containing bond in the direction of atom,
    adding newly found bonds and atoms to respective lists (listb, lista) which we'll return,
    until you can't or until you notice that it came back to bond and formed a ring
    (in which case return as much as possible, but not another ref to bond or atom).
       Return value is the tuple (ringQ, listb, lista) where ringQ says whether a ring was detected
    and len(listb) == len(lista) == number of new (bond, atom) pairs found.
    """
    listb, lista = [], []
    origbond = bond # for detecting a ring
    while 1:
        nextbond = next_bond_in_sp_chain(bond, atom) # this is the main difference from grow_bond_chain
        if nextbond is None:
            return False, listb, lista
        if nextbond is origbond:
            return True, listb, lista
        nextatom = nextbond.other(atom)
        listb.append(nextbond)
        lista.append(nextatom)
        bond, atom = nextbond, nextatom
    pass

def grow_pi_sp_chain_NEWER_BETTER(bond, atom):
    #bruce 070415 -- newer & better implem of grow_pi_sp_chain (equiv,
    # but uses general helper func), but untested. Switch to this one when there's time to test it.
    return grow_bond_chain(bond, atom, next_bond_in_sp_chain)

# ==

def pi_vectors(bond, out = DFLT_OUT, up = DFLT_UP, abs_coords = False): # rename -- pi_info for pi bond in degen sp chain
    # see also PiBondSpChain.get_pi_info
    """
    Given a bond involving some pi orbitals, return the 4-tuple ((a1py, a1pz), (a2py, a2pz), ord_pi_y, ord_pi_z),
    where a1py and a1pz are orthogonal vectors from atom1 giving the direction of its p orbitals for use in drawing
    this bond (for twisted bonds, the details of these might be determined more by graphic design issues than by the
    shapes of the real pi orbitals of the bond, though the goal is to approximate those);
    where a2py and a2pz are the same for atom2 (note that a2py might not be parallel to a1py for twisted bonds);
    and where ord_pi_y and ord_pi_z are order estimates (between 0.0 and 1.0 inclusive) for use in drawing the bond,
    for the pi_y and pi_z orbitals respectively. Note that the choice of how to name the two pi orbitals (pi_y, pi_z)
    is arbitrary and up to this function.
       All returned vectors are in the coordinate system of bond, unless abs_coords is true.
       This should not be called for pi bonds which are part of an "sp chain",
    but it should be equivalent to the code that would be used for a 1-bond sp-chain
    (i.e. it's an optimization of that code, PiBondSpChain.get_pi_info, for that case).
    """
    atom1 = bond.atom1
    atom2 = bond.atom2
    bond_axis = atom2.posn() - atom1.posn() #k not always needed i think
    # following subrs should be renamed, since we routinely call them for sp atoms,
    # e.g. in -C#C- where outer bonds are not potential pi bonds
    pvec1 = p_vector_from_sp2_atom(atom1, bond, out = out, up = up) # might be None
    pvec2 = p_vector_from_sp2_atom(atom2, bond, out = out, up = up) # ideally, close to negative or positive of pvec1
    # handle one being None (use other one in place of it) or both being None (use arb vectors)
    if pvec1 is None:
        if pvec2 is None:
            # use arbitrary vectors perp to the bond; compute differently if bond_axis ~= out [###@@@ make this a subr? dup code?]
            pvec = cross(out, bond_axis)
            lenpvec = vlen(pvec)
            if lenpvec < 0.01:
                pvec = up
            else:
                pvec /= lenpvec
            pvec1 = pvec2 = pvec
        else:
            pvec1 = pvec2
    else:
        if pvec2 is None:
            pvec2 = pvec1
        else:
            # both vectors not None -- use them, but negate pvec2 if this makes them more aligned
            if dot(pvec1, pvec2) < 0:
                pvec2 = - pvec2
    return pi_info_from_abs_pvecs( bond, bond_axis, pvec1, pvec2, abs_coords = abs_coords) 

def pi_info_from_abs_pvecs( bond, bond_axis, pvec1, pvec2, abs_coords = False):
    """
    #doc;
    bond_axis (passed only as an optim) and pvecs are in abs coords; retval is in bond coords unless abs_coords is true
    """
    a1py = pvec1
    a2py = pvec2
    a1pz = norm(cross(bond_axis, pvec1))
    a2pz = norm(cross(bond_axis, pvec2))
    ord_pi_y, ord_pi_z = pi_orders(bond)
    if not abs_coords:
        # put into bond coords (which might be the same as abs coords)
        quat = bond.bond_to_abs_coords_quat()
        a1py = quat.unrot(a1py)
        a2py = quat.unrot(a2py)
        a1pz = quat.unrot(a1pz)
        a2pz = quat.unrot(a2pz)
    return ((a1py, a1pz), (a2py, a2pz), ord_pi_y, ord_pi_z)

def pi_orders(bond):
    try:
        ord_pi_y, ord_pi_z = pi_order_table[bond.v6]
    except:
        if debug_flags.atom_debug:
            print "atom_debug: bug: pi_order_table[bond.v6] for unknown v6 %r in %r" % (bond.v6, bond)
        ord_pi_y, ord_pi_z = 0.5, 0.5 # this combo is not otherwise possible
    return ord_pi_y, ord_pi_z

pi_order_table = {
V_SINGLE:   (0,    0),
V_DOUBLE:   (1,    0), # choice of 1,0 rather than 0,1 is just a convention, but we always use it even for =C=C=C= chains
V_TRIPLE:   (1,    1),
V_AROMATIC: (0.5,  0),
V_GRAPHITE: (0.33, 0),
V_CARBOMERIC: (0.5, 1), # I think the 0.5 should actually be 1-x for adjacent bonds of order x, so it's wrong for graphite ###e
}

# ==
        
def p_vector_from_sp2_atom(atom, bond, out = DFLT_OUT, up = DFLT_UP):
    """
    Given an sp2 atom and a possibly-pi bond to it,
    return its p orbital vector for use in thinking about that pi bond,
    or None if this atom (considering its other bonds) doesn't constrain that direction.
    """
    if not atom.atomtype.is_planar():
        return None # helps make initial stub code simpler
    nbonds = len(atom.bonds)
    if nbonds == 3:
        # we can determine atom's p vector (there is only one) from those bonds
        # (no need to check atomtype -- it can't be sp, and if it was sp3 we should not have been called;
        #  and in case of errors or bugs, maybe the type is not sp2 as we hope, but this code will still "work".)
        return p_vector_from_3_bonds(atom, bond, out = out, up = up)
    elif nbonds == 2:
        return p_vector_from_sp2_2_bonds(atom, bond, out = out, up = up) ### revise to return None if it doesn't know?
    elif nbonds == 1:
        # might be O(sp2); if not for knowing it was sp2, might be open bond or N(sp)... but it won't tell us orientation
        return None
    else:
        assert bond in atom.bonds
        assert len(atom.bonds) <= 3 # will always fail
    pass

def p_vector_from_3_bonds(atom, bond, out = DFLT_OUT, up = DFLT_UP):
    """
    Given an sp2 atom with 3 bonds, and one of those bonds which we assume has pi orbitals in it,
    return a unit vector from atom along its p orbital, guaranteed perpendicular to bond,
    for purposes of drawing the pi orbital component of bond.
    Note that it's arbitrary whether we return a given vector or its opposite.
    [##e should we fix that, using out and up? I don't think we can in a continuous way, so don't bother.]
       We don't verify the atom is sp2, since we don't need to for this code to work,
    though our result would probably not make sense otherwise.
    """
    others = map( lambda bond: bond.other(atom), atom.bonds)
    assert len(others) == 3
    other1 = bond.other(atom)
    others.remove(other1)
    other2, other3 = others
    apos = atom.posn()
    v1 = other1.posn() - apos
    # if v1 has 0 length, we should return some default value here; this might sometimes happen so I better handle it.
    # actually i'm not sure the remaining code would fail in this case! If not, I might revise this.
    if vlen(v1) < 0.01: # in angstroms
        return + up
    v2 = other2.posn() - apos
    v3 = other3.posn() - apos
    # projecting along v1, we hope v2 and v3 are opposite, and then we return something perpendicular to them.
    # if one is zero, just be perp. to the other one alone.
    # (If both zero? Present code returns near-0. Should never happen, but fix. #e)
    # otherwise if they are not opposite, use perps to each one, "averaged",
    # which means (for normalized vectors), normalize the larger of the sum or difference
    # (equivalent to clustering them in the way (of choice of sign for each) that spans the smallest angle).
    # Optim: no need to project them before taking cross products to get the perps to use.
    ## v2 -= v1 * dot(v2,v1)
    ## v3 -= v1 * dot(v3,v1)
    v2p = cross(v2,v1)
    v3p = cross(v3,v1)
    lenv2p = vlen(v2p)
    if lenv2p < 0.01:
        return norm(v3p)
    v2p /= lenv2p
    lenv3p = vlen(v3p)
    if lenv3p < 0.01:
        return v2p # normalized above
    v3p /= lenv3p
    r1 = v2p + v3p
    r2 = v2p - v3p
    lenr1 = vlen(r1)
    lenr2 = vlen(r2)
    if lenr1 > lenr2:
        return r1/lenr1
    else:
        return r2/lenr2
    pass

def p_vector_from_sp2_2_bonds(atom, bond, out = DFLT_OUT, up = DFLT_UP):
    others = map( lambda bond: bond.other(atom), atom.bonds)
    assert len(others) == 2
    other1 = bond.other(atom)
    others.remove(other1)
    other2, = others
    apos = atom.posn()
    v1 = other1.posn() - apos
    if vlen(v1) < 0.01: # in angstroms
        return + up
    v2 = other2.posn() - apos
    
    v2p = cross(v2,v1)
    lenv2p = vlen(v2p)
    if lenv2p < 0.01:
        # this entire case happens rarely or hopefully never (only when 2 sp2 bonds are parallel)
        res = up - v1 * dot(up,v1)
        if vlen(res) >= 0.01:
            return norm(res)
        res = out - v1 * dot(out,v1)
        return norm(out)
    return v2p / lenv2p

# ==

# pure geometry routines, should be refiled:

def arb_non_parallel_vector(vec):
    """
    Given a nonzero vector, return an arbitrary vector not close to
    being parallel to it.
    """
    x,y,z = vec
    if abs(z) < abs(x) > abs(y):
        # x is biggest, use y
        return V(0,1,0)
    else:
        return V(1,0,0)
    pass

def arb_perp_unit_vector(vec):
    """
    Given a nonzero vector, return an arbitrary unit vector
    perpendicular to it.
    """
    vec2 = arb_non_parallel_vector(vec)
    return norm(cross(vec, vec2))

def arb_ortho_pair(vec): #k not presently used # see also some related code in pi_vectors()
    """
    Given a nonzero vector, return an arbitrary pair of unit vectors
    perpendicular to it and to each other.
    """
    res1 = arb_perp_unit_vector(vec)
    res2 = norm(cross(vec, res1))
    return res1, res2

# end