# *********
#  Analyze RNA sencond structure information from 
#             RNA sequence and duplex Sequences  
# 

import os
from string import lstrip,strip, rstrip
from re import split
from math import  log, sin, cos, sqrt, pi, fabs, radians,degrees,acos, atan2

import numpy as Npy
from numpy import dot, transpose,shape, array, argsort
from numpy.linalg import det, svd  
_array = array 

from libmath import Vector,cross, std,dihedral_vect
from libpdb import PDB, RESIDUE, Atom
from libaction import Rotation
from libntatoms import * 

def getRefBase():
    #fn = 'db/ref_base.bin'
    #if os.path.exists(fn):
        #rf_bases = loadObj(fn)
    #else:
    robj = RefBase()
    rf_bases = robj.getAllRefBase() #getARefBase()
    #    dumpObj(rf_bases, fn)
    return rf_bases


class RefBase:
    ref_bases = { "U": { "C1'":(-2.481,  5.354,   0.000), 
                         "N1": (-1.284,  4.500,   0.000),
                         "C2": (-1.462,  3.131,   0.000),
                         "O2": (-2.563,  2.608,   0.000),
                         "N3": (-0.302,  2.397,   0.000),
                         "C4": ( 0.989,  2.884,   0.000),
                         "O4": ( 1.935,  2.094,  -0.001),
                         "C5": ( 1.089,  4.311,   0.000),
                         "C6": (-0.024,  5.053,   0.000)},
                
                   "A": { "C1'":(-2.479,   5.346,   0.000),
                         "N9": (-1.291,   4.498,   0.000),
                         "C8": ( 0.024,   4.897,   0.000),
                         "N7": ( 0.877,   3.902,   0.000),
                         "C5": ( 0.071,   2.771,   0.000),
                         "C6": ( 0.369,   1.398,   0.000),
                         "N6": ( 1.611,   0.909,   0.000),
                         "N1": (-0.668,   0.532,   0.000),
                         "C2": (-1.912,   1.023,   0.000),
                         "N3": (-2.320,   2.290,   0.000),
                         "C4": (-1.267,   3.124,   0.000)},
                
                  "C": { "C1'":(-2.477,   5.402,   0.000),
                         "N1": (-1.285,   4.542,   0.000),
                         "C2": (-1.472,   3.158,   0.000),
                         "O2": (-2.628,   2.709,   0.001),
                         "N3": (-0.391,   2.344,   0.000),
                         "C4": ( 0.837,   2.868,   0.000),
                         "N4": ( 1.875,   2.027,   0.001),
                         "C5": ( 1.056,   4.275,   0.000),
                         "C6": (-0.023,   5.068,   0.000)},
                 
                 "G": { "C1'": ( -2.477,   5.399,   0.000),
                        "N9": ( -1.289,   4.551,   0.000),
                        "C8": (  0.023,   4.962,   0.000),
                        "N7": (  0.870,   3.969,   0.000),
                        "C5": (  0.071,   2.833,   0.000),
                        "C6": (  0.424,   1.460,   0.000),
                        "O6": (  1.554,   0.955,   0.000),
                        "N1": ( -0.700,   0.641,   0.000),
                        "C2": ( -1.999,   1.087,   0.000),
                        "N2": ( -2.949,   0.139,  -0.001),
                        "N3": ( -2.342,   2.364,   0.001),
                        "C4": ( -1.265,   3.177,   0.000)}}

    def getRingAnames(self, base_name):
        """ return base ring atom names only for RNA (AGCU) """
        if base_name in ["G", "A", "GUA", "ADE"]:
            ring_anames = ["C4", "N3", "C2", "N1", "C6", "C5", "N7", "C8", "N9"]
        elif base_name in ['C','U','CYT','URA']:
            ring_anames = ["C4", "N3", "C2", "N1", "C6", "C5"]
        else:
            self.err_exit("# base name not in GUA(G), ADE(A), CYT(C), and URA(U) ")
        return ring_anames
           
    
    def getRefBase(self, base_name):
        """ return reference base ring atoms' coordinates """
        ring_anames = self.getRingAnames(base_name)
        ref_base = self.ref_bases[base_name[0]]
        base_xyz = []
        for atom_name in ring_anames:
            base_xyz.append(ref_base[atom_name])
        base_xyz = Npy.array(base_xyz)
        return base_xyz
    
    
    def RefBaseCenter(self, base_name, base_xyz):
        """ return base center(only consider ring atoms) and centered ring atoms' coordinates"""
        n_vec, vec_size = base_xyz.shape
        base_center = Npy.sum(base_xyz, 0)/float(n_vec)
        ct_rf_base_xyz_ct = base_xyz.copy() - base_center
        return base_center, ct_rf_base_xyz_ct
    
    
    def prnBase(self,base_name, base_xyz, cmt=""):
        lns = ["# base_name= %s  %s" % (base_name,cmt)]
        num_atom, dim = base_xyz.shape 
        for i in range(num_atom):
            lns +=[ "# %8.3f %8.3f %8.3f " % (base_xyz[i][0], base_xyz[i][1], base_xyz[i][2])]
        self.prnLns(lns)
        
            
    def getAllRefBase(self):
        """ return all reference bases ring atoms coordinates """
        aref_base = {}
        for base_name in ['U', "A", "G", "C"]:
            rbase_xyz = self.getRefBase(base_name)
            base_ct, rbase_xyz_ct = self.RefBaseCenter(base_name, rbase_xyz)
            aref_base[base_name[0]] = (base_ct, Npy.transpose(rbase_xyz_ct), Npy.sum(rbase_xyz_ct * rbase_xyz_ct), rbase_xyz_ct)
        return aref_base




rf_bases = getRefBase()
class RNA2D:
    """ Extract RNA secondary structure from WC basepairs and RNA sequence 
        duplex's functions on secondary structure level, no Cartesian coordinates """
    sec_inf = {}
    def __init__(self):
        #self.sec_inf = {}
        pass 
 
    def getBpSeq(self, bp):
        return "%s%s" % (self.seq_dict[bp[0]], self.seq_dict[bp[1]])
    
    def addSS2inf(self, stype, _ss, vb=0):
        if vb: print "# add %s to 2nd inf: " % stype,  _ss 
        if stype not in self.sec_inf.keys():
            self.sec_inf[stype] = [_ss]
        else:
            self.sec_inf[stype].append(_ss)
   
        
    def prnStr1_2(self):
        """ print out sequence and base pairs of RNA """
        lns = ["# %3d nucleotides " % len(self.resids)]
        lns += ['start_resid= %3d ' % self.resids[0]]
        lns += ['end_resid= %3d ' % self.resids[-1]]
        lns += [" ", "# Seqence ", " "]
        seq_str = ""
        cnt = 0
        for resid in self.resids:
            seq_str += self.seq_dict[resid]
            cnt += 1
            if cnt == 40:
                cnt = 0
                lns.append(seq_str)
                seq_str = ""
        lns.append(seq_str)

        num_helix = len(self.helices)
        lns += [" ",  "# num of helices: %3d " % num_helix ]
        for i in range(num_helix):
            helix = self.helices[i]
            lns += ["", "# helix: %3d " % (i+1) ]
            for j in range(len(helix)):
                _bp = helix[j]
                lns += [ "%3d %s%s %3d " % ( _bp[0],self.seq_dict[_bp[0]], self.seq_dict[_bp[1]], _bp[1])]
        self.prnLines(lns)

        
    def resid2BPresid(self, bp_resids,vb=0):
        """ input a duplex resids  and return basepair resids """
        num_bp = len(bp_resids)/2 
        bps = []
        for i in range(num_bp):
            bps.append((bp_resids[i], bp_resids[-(i+1)]))
        if vb: print "# bps: ", bps 
        return bps 
    

    def bp2BP_Resid(self, bps):
        """ return only base pair resids in terms of (bp_resid1, bp_resid2) from basepairs """
        bp_resids = []
        for _bp in bps:
            bp_resids.append((_bp[0], _bp[1]))
        return bp_resids

    
    def bp2resids(self, bps):
        """ transform basepair resid(bp_left_resid, bp_right_resid) to a resid list: resids 
            a reverse function of resid2BPresid
        """
        resids = []
        for _bp in bps: resids += [_bp[0], _bp[1]]
        return resids 


    def getHelixEnd(self, helices):
        """ return helix's  two end pair residue IDs """
        helix_end = []
        helix_end_resids = []
        #print "helics: ",  helices
        if type(helices) == type({}):
            for h_key in helices.keys():
                helix = helices[h_key] #print "# helix: ", helix 
                helix_end +=  [(helix[0][0],helix[0][1]), (helix[-1][0],helix[-1][1])]
        else:
            for helix in helices:
                helix_end +=  [(helix[0][0],helix[0][1]), (helix[-1][0],helix[-1][1])]
        return helix_end


    def locateHelix(self):
        """ find connected basepairs -- only WC basepairs and GU wobble basepair """
        #hex_bp_names = ("CG","GC","AU","UA","GU","UG")
        #bp_dseq = {"GC":'0',"CG":'1',"AU":'2',"UA":'3',"GU":'4',"UG":'5',
        #           "AC":'6',"CA":'7',"AG":'8',"GA":'9', "CU":'a',"UC":'b',
        #           "UU":'c' }	
        hex_bp_names = ("CG","GC","AU","UA","GU","UG","AC","CA","AG","GA","CU","UC","UU",'GG',"CC","AA")
        _bps = self.BP_resids
        _bp_rkeys = self.BP_resids_keys
        num_bp = len(_bp_rkeys)
        helices = []
        helix_resids = []
        helix = []
        empty_duplex = 1
        for resid in _bp_rkeys:
            if len(_bps[resid]) == 1:  # No triple or more interactions among residues
                _bp = _bps[resid][0]
                if _bp[3] in hex_bp_names:
                    if empty_duplex:
                        helix.append(_bp)
                        empty_duplex = 0
                        continue 
                    else:
                        pre_bp = helix[-1]
                        if pre_bp[0]+1 == _bp[0] and pre_bp[1]-1 == _bp[1]:
                            helix.append(_bp)
                            continue
                        else:
                            if len(helix) > 1:
                                helices.append(helix)
                                helix_resids += self.bp2resids(helix)
                            helix = [_bp]
        if len(helix) > 1:
            helices.append(helix)
            helix_resids += self.bp2resids(helix)
        self.helices = helices
        self.helicx_resids = helix_resids
     
    
    def locateSS(self, helix_resids=None, vb=0):
        """ locate single strand in RNA molecules """
        sstrands = []
        ss_resids = []
        if helix_resids is None: helix_resids = self.helices_resids
        #print "# self.resids: ", self.resids 
        #print "# helix_resids: ", helix_resids
        for resid in self.resids:
            if resid not in helix_resids:
                ss_resids.append(resid)
        
        # find single single strands by residue ID connection 
        _ss = []
        pre_resid = -1 
        for resid in ss_resids:
            if resid == pre_resid+1:
                _ss.append(resid)
            else:
                if _ss: sstrands.append(_ss)
                _ss = [resid]
            pre_resid = resid 
        sstrands.append(_ss)
        if vb: 
            print "# self.sstrands: ", sstrands
            print "# single strand resid: ",ss_resids
        self.sstrands = sstrands
   
        
    def parseSS(self, ss_check=1,  vb=0):
        self.locateSS(vb=vb)
        if ss_check: self.chkSSConnect(vb=vb)
        self.helice2duplex()
        self.setAllDupBP()
        self.parseSingleStrands()
     
    
    def helice2duplex(self):
        """ return duplexs dictionary format, key: duplex_id from helices list """
        if hasattr(self, 'helices'):
            assert type(self.helices) == type([])
            self.duplexs = {}
            cnt = 1
            for _helix in self.helices:
                duplex = []
                for _bp in _helix:
                    duplex.append((_bp[0], _bp[1]))
                self.duplexs[cnt] = duplex
                cnt += 1
            self.sec_inf['dup_resids'] = self.helices_resids
        elif hasattr(self, 'duplexs'):
            self.sec_inf['dup_resids'] = self.helices_resids
        else:
            self.err_exit("Could not find  helices ")
        self.duplex_ids = self.duplexs.keys()
        self.duplex_ids.sort()
        self.sec_inf['duplexs'] = self.duplexs
        
    
    
    def setAllDupBP(self):
        """ set all basepair dictionary """
        dup_bps  = {}
        for dup_id in self.duplex_ids:
            duplex = self.duplexs[dup_id]
            for _bp in duplex:
                dup_bps[_bp[0]] = _bp
                dup_bps[_bp[1]] = _bp
        self.dup_bps = dup_bps
    
    
    def parseSingleStrands(self):
        """ parse single strands into hairpin loop, 
            linker(junction loop which links two duplexes), 
            internal loop(two single trands link two duplexes),
            dangling strand(only linking to one duplex, the other terminal residue is free )
            bulge(linking one chain of two duplexes, the other chains of two dupexes are directly linker together
            
            Note: RNA residue numbers should be continuous, no break up in any two residues  ***
            """
        sstrands = self.sstrands
        bulges = []
        oss  = []
        bp_resids = self.dup_bps.keys()
        for _ss in sstrands:
            if len(_ss) == 0: continue 
            _s5 = _ss[0]-1    # previous resid
            _s3 = _ss[-1] + 1 # next resid 
            
            if _s5 in bp_resids:
                end_bp1 = self.dup_bps[_s5]
            else:
                end_bp1 = None 
            
            if _s3 in bp_resids:
                end_bp2 = self.dup_bps[_s3]
            else:
                end_bp2 = None
            
            if end_bp1 is None:
                if end_bp2: 
                    self.addSS2inf('dss', (end_bp2, _ss))
                    continue
                else:
                    print "# free ss: ", _ss 
                    #self.err_exit(" free single strand ")
                    self.addSS2inf('fss', _ss)
            else:
                if end_bp2 is None:
                    self.addSS2inf('dss',(end_bp1, _ss))
                    continue
                if end_bp1 == end_bp2:
                    self.addSS2inf('tloop', (end_bp1, _ss))
                    continue 
                ibg = self.isBulge(end_bp1, end_bp2)
                if ibg > -1:
                    self.addSS2inf('bulge', (end_bp1, end_bp2, _ss, ibg))
                else:
                    oss.append((end_bp1, end_bp2, _ss))
                    #print "# end: bp1: %s   bp2: %s  :: ss: %s" % (end_bp1, end_bp2, _ss) 
        #print "# oss: ", oss 
        self.__parseIloopLinker(oss)
        #self.prnDict(self.sec_inf)
        
    
    def __parseIloopLinker(self, oss):
        """ parse internal loops and linkers"""
        num_oss = len(oss)
        rm_ssids = []
        for i in range(num_oss):
            if i in rm_ssids: continue 
            ebp1, ebp2, ss = oss[i]
            ILOOP = 0
            for j in range(i+1, num_oss):
                if j in rm_ssids: continue 
                ebpj1, ebpj2, ssj = oss[j]
                if (ebp1 == ebpj2 and ebp2== ebpj1) or (ebp1 == ebpj2 and ebp2== ebpj1):
                    self.addSS2inf('iloop', (ebp1, ebp2, ss, ssj))
                    rm_ssids.append(j)
                    ILOOP = 1
                    break 
            if not ILOOP:
                self.addSS2inf('linker', (ebp1,ebp2,ss))
        #print self.sec_inf 
            


    def isBulge(self, bp1, bp2):
        d5 = fabs(bp1[0] - bp2[0])
        d3 = fabs(bp1[1] - bp2[1])
        if d5 > 1 and d3 == 1:
            return  0  # bulge in 5' chain
        elif d5== 1 and d3 > 1:
            return  1  # bulge in 3' chain 
        else:
            return -1
            
class ConfRMSD:
    """ calculate rmsd between two conformers """
    def calculate_rmsd(self, conformation1, conformation2,conformation3,       
                       is_align = False, is_center = True, vb = False, pdb3p=1):      
        """      
	Calculates the RMSD between two conformations.      

	  * conformation1: array of dimensions [N,3]       
	  * conformation2: array of dimensions [N,3]       
	  * is_align: True = modify conformation2 to be 
			     aligned to conformation1      
	  * is_center: True = remove average center differences      
			      before calculating RMSD       
	"""      
        #nrow, ncol = conformation1.shape
        #print "# cf1: nrow= %3d ncol= %3d " % (nrow, ncol)
        #nrow, ncol = conformation2.shape
        #print "# cf2: nrow= %3d ncol= %3d " % (nrow, ncol)
        #print shape(conformation1) == shape(conformation2)
        if not shape(conformation1) == shape(conformation2):      
            print "# conformation1 shape: ", conformation1.shape
            print "# cofnormation2 shape: ", conformation2.shape
            raise "calculate_rmsd: conformations not the same size."      
        elif not len(shape(conformation1)) == 2:      
            raise "calculate_rmsd: conformations not the correct rank."      

        vecs1 = conformation1.copy()      
        vecs2 = conformation2.copy()      
        vecs3 = conformation3.copy()
        n_vec, vec_size = shape(vecs1)
        n_vec3, vec3_size = shape(vecs3)      
        if is_center:      
            center1 = sum(vecs1,0) / float(n_vec)      
            center2 = sum(vecs2,0) / float(n_vec)
            center3 = sum(vecs3,0) / float(n_vec3)      
            vecs1 -= center1      
            vecs2 -= center2      
            vecs3 -= center3

            E0 = sum(sum(vecs1 * vecs1)) + sum(sum(vecs2 * vecs2))      

            correlation_matrix = dot(transpose(vecs2), vecs1)      

            V, S, W_trans = svd(correlation_matrix)      

            is_reflection = (det(V) * det(W_trans)) < 0.0      
            if is_reflection:      
                # reflect along smallest principal axis      
                S[-1] = -S[-1]      

            rmsd_sq = (E0 - 2. * sum(S)) / float(n_vec)      
            rmsd_sq = max([rmsd_sq, 0.])      

            if is_align:      
                if is_reflection:       
                    V[-1,:] = V[-1,:] * (-1.0)      
                optimal_rotation = dot(V, W_trans)      
                #conformation2[:,:] = dot(vecs2, optimal_rotation) + center2
                conformation2[:,:] = dot(vecs2, optimal_rotation) + center1 
                conformation3[:,:] = dot(vecs3, optimal_rotation) + center1   
            rmsd = sqrt(rmsd_sq)
            if vb: print "#  rmsd= %8.4f " %  rmsd 
            #print conformation2.shape
            return rmsd 


        
class Duplex(Vector, RNA2D, ConfRMSD, Rotation):
    zaxis = array([0,0,1.0])
    fn_bpseg = "bp_seg.bin"
    bp_dseq = {"GC":'0',"CG":'1',"AU":'2',"UA":'3',"GU":'4',"UG":'5',
                   "AC":'6',"CA":'7',"AG":'8',"GA":'9', "CU":'a',"UC":'b',
                   "UU":'c',"GG":'d',"AA":'e',"CC":'f'}	
    def __init__(self):
        RNA2D.__init__(self)
        #print self.sec_inf 
        
    def thePhi2Vector(self, the, phi):
        """ return vector in Cartesian  space from the, phi(in units of degree) in polar space """
        the = radians(the)
        phi = radians(phi)
        x = cos(phi) * sin(the)
        y = sin(phi) * sin(the)
        z = cos(the)
        return array([x,y,z])

    def vector2thePhi(self, axis, vb=1):
        """ transfer vector in Cartesian space to the, phi(in units of degree) in Polar space """
        rAB = self.normVect(axis)
        if vb: self.prnVec(rAB, cmt="# normal axis: ")
        _theta  = acos(rAB[2])
        if fabs(rAB[0]) < 0.001 and fabs(rAB[1]) < 0.001:
            _phi = 0
        else:
            _phi = atan2(rAB[1],rAB[0])
        if vb: print "# bond: theta= %6.1f  phi= %6.1f " % (degrees(_theta), degrees(_phi))
        return degrees(_theta), degrees(_phi)

    def getDupOrient(self, dup_axis, phase_vec, vb=0):
        the, phi = self.axis2ThePhi(dup_axis, vb=vb)
        phase_vec = self.new_arb_rot(phase_vec, self.ZAXIS, -phi)
        phase_vec = self.new_arb_rot(phase_vec, self.YAXIS, -the) #, hex_axis)
        delta, rho = self.vector2thePhi(phase_vec, vb=vb)
        if rho < 0.0: rho += 360.0
        return the, phi, rho, delta 
        
    def chk4_Orient(self, helix, tpt_o4):
        """ find correct orientation from 4 orientations with crystal structure """
        NH_vecs = self.bp2NH(self.bps, self.bp_ids)
        hex_key = 'h%s' % helix
        print "# duple: %s " % hex_key
        std_results = []
        o4 = tpt_o4[hex_key]
        for the, phi, rho in o4:
            print "\n# ---------------------"
            print "# orient %s:  the= %6.1f phi= %6.1f  " % (hex_key, the, phi )
            helix_axis = self.thePhi2Axis(the, phi)
            helix_axis = self.getHexDirection(helix_axis)
            self.NH_Project(NH_vecs, helix_axis)
            deltas = self.getIminoSlantAngle(NH_vecs, helix_axis)
            std_results.append(std(deltas, vb=1))
        print "\n# duplex %s " % hex_key
        for i in range(4):
            the, phi, rho = o4[i]
            av, sigma, dmax, dmin = std_results[i] #std
            print "# the= %6.1f phi= %6.1f rho= %6.1f  delta: av= %8.3f sig= %8.3f max= %8.3f min= %8.3f " % \
                  (the, phi, rho, av, sigma, dmax, dmin)
    
                                  
    
    def getPPvecs(self, bps):
        """ return vectors P --->  P(basepair)"""
        xyz = []
        for _bp in bps:
            dp = self.getBPDp(_bp)
            if dp is not None:
                xyz.append(dp)
            #print "# resid= %3d P --> resid= %3d  P  dp_len= %8.3f " % (resid1, resid2, self.getVecLen(dp))
        #print "#bp_xyz: ", xyz
        return xyz 


    def getBPDp(self, bp_resid, vb=0):
        """ return P-->P vector of basepair """
        #resid1, resid2 = bp_resid
        resid1 = bp_resid[0]
        resid2 = bp_resid[1]
        res1 = self.getRes_by_resid(resid1)
        res2 = self.getRes_by_resid(resid2)
        p1 = res1.getAtomVect('P')
        p2 = res2.getAtomVect('P')
        #print "# p1: ", p1 , "p2: ", p2 
        if p1 is None or  p1[0] < -998:
            return None
        elif p2 is None or p2[0] < -998:
            return None 
        else:
            dp = self.normVect(p2 -p1)
            if vb: print "# dp: %3d->%3d :  %8.3f %8.3f %8.3f " % (resid2, resid1, dp[0], dp[1], dp[2])
            return dp 
    
    
    def dupxyz2pp(self, bps, dup_xyz):
        xyz = []
        for _bp in bps:
            p1 = dup_xyz[_bp[0]][0] 
            p2 = dup_xyz[_bp[1]][0] 
            dp = self.normVect(p2-p1)
            xyz.append(dp)
        return xyz 
    
    
    def vec2planeNorm(self, vecs):
        #print "vecs: ", vecs 
        plane_vecs = []
        for i in range(len(vecs)-1):
            plane_vecs.append(vecs[i+1] - vecs[i])
        plane_norm = array([0., 0., 0.])
        for i in range(len(plane_vecs)-1):
            plane_norm += self.cross(self.normVect(plane_vecs[i]), self.normVect(plane_vecs[i+1]))
        #print "# plane_norm: ", plane_norm
        #print "# plane vecs: ", plane_vecs
        return self.normVect(plane_norm)
    
                
    def pp2axis(self, bp_resids, dup_xyz=None, vb=0):
        """ calculate helix axis by vector P of one strand to the P of base pair """
        #print "bp_resids: ", bp_resids
        num_bp = len(bp_resids)
        if num_bp < 3: 
            #print "# number of basepair should be >= 3 "
            err_exit("# number of basepair should be >= 3: Num(bps)= %3d " % num_bp)
        if dup_xyz is None:
            bp_vecs = self.getPPvecs(bp_resids)
        else:
            bp_vecs = self.dupxyz2pp(bp_resids, dup_xyz)
        #if vb: print "#bp_vecs: \n",  bp_vecs
        #duplex_axis = self.vec2planeNorm(self.getPPvecs(bp_resids))
        duplex_axis = self.vec2planeNorm(bp_vecs)
        self.bp_pp_vecs = bp_vecs
        #print "# before check sign: "
        #self.prnVec(duplex_axis)
        return self.getDuplexAxis(duplex_axis, bp_resids, dup_xyz=dup_xyz, vb=vb)
    
    
    
    def getDuplexAxis(self, duplex_axis, bp_resids, dup_xyz=None, vb=0):
        """ return direction from 5' -> 3' """
        if dup_xyz is None:
            dp = self.getBPDp((bp_resids[0][0], bp_resids[3][0]), vb=0)
            if dp is None:
                dp = self.getBPDp((bp_resids[1][0], bp_resids[4][0]), vb=0)
        else:
            dp = dup_xyz[bp_resids[3][0]][0] - dup_xyz[bp_resids[0][0]][0]
            if self.getVecLen(dp) > 50: 
                dp = dup_xyz[bp_resids[4][0]][0] - dup_xyz[bp_resids[1][0]][0]
            #print dp 
        #print "# dp: ", dp 
        #print "# dup_axis: ", duplex_axis 
        cosp = dot(dp, duplex_axis)
        #print "# cosp: ", cosp 
        if cosp < -0.3: 
            #print "# cosp: ", cosp 
            #print "# acosp: ", degrees(acos(cosp))
            #print "# duplex axis: %8.3f %8.3f %8.3f " % (duplex_axis[0], duplex_axis[1], duplex_axis[2])
            duplex_axis = - duplex_axis 
            #print "# *** reverse Duplex Axis Z direction ***"
            #print "# cosp= %7.3f  cosp= %6.1f " % (cosp, degrees(acos(cosp)))
            #print "# ********     keep axis ponting from 5'--> 3' in lower resid strand of duplex "
        if vb: print "# duplex axis: %8.3f %8.3f %8.3f " % (duplex_axis[0], duplex_axis[1], duplex_axis[2])
        return duplex_axis
    
    
    def axis2ThePhi(self, axis, vb=0):
        #print "# test axis:"
        the, phi = self.vector2thePhi(axis, vb=vb)
        if vb:  print "# Duplex orientation: the= %6.1f phi= %6.1f "  % (the, phi)
        return the, phi 
        
    
    def thephi2rho(self, the, phi, phase_resid,  bps, duplex_xyz, vb=0):
        """ calculate phase angle by rotate duplex to public frame(the, phi in units of degree) """
        #self.prnDict(self.seq_dict)
        nt_name = self.seq_dict[phase_resid]
        xyz = duplex_xyz[phase_resid]
        rNH = self.resXYZ2Imino(Npy.array(xyz), nt_name)
        #print "# rNH: ", rNH
        rNH = self.new_arb_rot(rNH, self.ZAXIS, -phi)
        rNH = self.new_arb_rot(rNH, self.YAXIS, -the) #, hex_axis)
        phase_the, rho = self.vector2thePhi(rNH, vb=vb)
        # first imino proton' projector on xy plane. 
        if vb: print "# duplex  the= %6.1f phi= %6.1f  rho= %6.1f phase_the= %6.1f" % (the, phi, rho, phase_the)
        return rho 
    
    
    def getIminoPhases(self, bps):
        """ directly use imino bond vectors without imino proton addition"""
        drho = 360.0/11.0
        trho = 180.
        for _bp in bps:
            for i in range(2):
                resid = _bp[i]
                _res = self.getRes_by_resid(resid)
                #rNH = _res.getImino()
                rNH = _res.addImino()
                if rNH is None: continue 
                nt_name = _res.nt_name 
                #self.prnVec(rNH, cmt="#  imino rNH: ")
                phase_the, rho = self.vector2thePhi(rNH, vb=0)
                orho = rho 
                if phase_the < 90.:
                    rho += 180.0
                if rho < 0.:
                    rho += 360.0
                dr = trho - rho 
                print "# duplex:   resid= %3d(%s)    orho= %6.1f   rho= %6.1f t_rho= %6.1f  dr= %6.1f   pthe= %6.1f" % (resid, nt_name, orho, rho, trho, dr, phase_the)
                trho += drho
                if trho > 360.: trho -= 360.0
                
                
    def calIminoPhases(self, bps, duplex_xyz):
        """ add imino proton by geometry calculation """
        drho = 360.0/11.0
        trho = 180.
        for _bp in bps:
            for i in range(2):
                resid = _bp[i]
                nt_name = self.seq_dict[resid]
                xyz = duplex_xyz[resid]
                rNH = self.resXYZ2Imino(Npy.array(xyz), nt_name)
                if rNH is None: continue 
                #self.prnVec(rNH, cmt="#  imino rNH: ")
                phase_the, rho = self.vector2thePhi(rNH, vb=0)
                orho = rho 
                if phase_the < 90.:
                    rho += 180.0
                if rho < 0.:
                    rho += 360.0
                dr = trho - rho 
                print "# duplex:   resid= %3d(%s)    orho= %6.1f   rho= %6.1f t_rho= %6.1f  dr= %6.1f   pthe= %6.1f" % (resid, nt_name, orho, rho, trho, dr, phase_the)
                
                trho += drho
                if trho > 360.: trho -= 360.0
    
    def analyzSlantAngle(self, duplex_axis, duplex_resids, vb=0):
        """ NH slant angles"""
        dup_init_resid = duplex_resids[0]
        num_bp = len(duplex_resids)/2
        chain_3_thresh = dup_init_resid + num_bp
        #print "chain3_trhesh: ", chain_3_thresh
        deltas = []
        for i in duplex_resids:
            _residue = self.getRes_by_resid(i)
            NH_imino  = _residue.getImino(vb=vb)
            if NH_imino is None:
                NH_imino  = _residue.addImino()
            #print "NH_imino: ", NH_imino    
            if NH_imino is not None:
                #print "# get NH: ",  NH_imino
                if i >= chain_3_thresh: 
                    reverse = 1
                    NH_imino = -NH_imino  # reverse NH bond vectors in 3' Chain
                else:
                    reverse = 0
                delta = self.getVecAngle(duplex_axis, NH_imino)
                print "# reverse= %3d  %s:  res_sid= %3d  resid= %3d delta= %6.1f   "  % (reverse, _residue.nt_name, i, _residue.resid, degrees(delta))
                deltas.append(degrees(delta))
        return std(deltas, vb=1,name='delta')
   
            
    def res2ChiBonds(self, resids):
        """ return glycosyl bond vectors """
        bonds  = []
        for resid in resids:
            _residue = self.residues[resid]
            #if _residue.nt_name in ['G','A']:
            bonds.append( _residue.getXBond())
        return bonds
    
    def getRingNorms(self, resids):
        """ return ring plane norms """
        bonds  = []
        for resid in resids:
            _residue = self.residues[resid]
            #if _residue.nt_name in ['G','A']:
            bonds.append( _residue.getRingNorm())
        return bonds
    
        
    def prnDuplex(self, duplexes=None):
        """ print duplex """
        if duplexes is None: duplexes = self.helices
        num_helix = len(duplexes)
        print "# num of helices: %3d " % num_helix
        for i in range(num_helix):
            helix = duplexes[i]
            print "\n# helix: %3d " % (i+1)
            print "# bp_id  res_sids  bp_seq   resids" #  resids: residue number in pdb file
            for j in range(len(helix)):
                _bp = helix[j]
                print "%8d  %3d--%3d   %s    %3d--%3d  " % (j+1, _bp[0], _bp[1], _bp[3], _bp[-2],_bp[-1])#,_bp[3], _bp[4], _bp[5],_bp[6])
                
    def chkTriPair(self, bps):
        if len(bps) == 1:
            return bps[0]
        else:
            new_bp = [] # remove intercalated basepair
            # (126, 144, 68, 'UA', '-/- cis ', 0.0027042772300864826, 12.690883623233105, 228, 246)
            #    0    1   2    3        4          5                       5               7     8
            bp_vd = []
            for _bp in bps:
                #if _bp[5]< 1.0: new_bp.append(_bp)
                bp_vd.append(_bp[5])
            min_bp_vd = min(bp_vd)
            min_bp_id = bp_vd.index(min_bp_vd)
            return bps[min_bp_id]
                
            #if len(new_bp) == 1: 
            #    return new_bp[0]
            #else:
            #    return None 
            
    
    def locateHelix3D(self):
        """ find connected basepairs: check the residue connectivity by cooridnates not buy residue resid 
            remove the basepair without P atom(s) 
        """
        #hex_bp_names = ("CG","GC","AU","UA","GU","UG")
        #bp_dseq = {"GC":'0',"CG":'1',"AU":'2',"UA":'3',"GU":'4',"UG":'5',
        #           "AC":'6',"CA":'7',"AG":'8',"GA":'9', "CU":'a',"UC":'b',
        #           "UU":'c' }	
        hex_bp_names = ("CG","GC","AU","UA","GU","UG","AC","CA","AG","GA","CU","UC","UU",'GG',"CC","AA")
        
        _bps = self.BP_resids
        _bp_rkeys = self.BP_resids_keys
        num_bp = len(_bp_rkeys)
        helices = []
        #self.prnBP()
        helix = []
        empty_duplex = 1
        for resid in _bp_rkeys:
            _bp = self.chkTriPair(_bps[resid])
            if _bp is None: 
                print "# No bp for resid: " % _bps[resid]
                continue 
            else:
                if _bp[3] in hex_bp_names:
                    if empty_duplex:
                        helix.append(_bp)
                        empty_duplex = 0
                        continue 
                    else:
                        pre_bp = helix[-1]
                        if self.isConnect(pre_bp[0], _bp[0], vb=0) and self.isConnect(pre_bp[1],_bp[1], vb=0):
                            helix.append(_bp)
                            continue
                        else:
                            if len(helix) > 1:
                                helices.append(helix)
                            helix = [_bp]
                else:
                    print "# bp_name: %s  not in hex_bp_names: %s " % (_bp[3], hex_bp_names)
                    
        if len(helix) > 1:
            helices.append(helix)
        self.helices = helices    
        #self.rmPseudoKnots(helices)
        self.update_duplex_resids(helices)
        
    
    def update_duplex_resids(self, helices):
        """ regenerate resids """
        
        helix_resids = []
        # (202, 210, 101, 'AA', 'S/S cis ', 1.1065983480567945, 47.976257158293805, 225, 233)
        #  iresid1, iresid2, bp_id, bp_seq, bp_type, dv, ang_norm, resid1, resid2
        for _helix in helices:
            print "\n# helix: "
            for _bp in _helix:
                #print _bp 
                print "%3d  %3d %s  %7.2f  %6.1f  %3d--%3d" % (_bp[0], _bp[1], _bp[3], _bp[5], _bp[6], _bp[-2], _bp[-1]) 
            helix_resids += self.bp2resids(_helix)
            #print "# helix_reisds: ", helix_resids
        helix_resids.sort()    
        self.helices_resids = helix_resids
        #print "# helieces_resids: ", self.helices_resids
        #self.addH2inf(helices, helix_resids)
    
    
    def rmPseudoKnots(self, helices): # unfished
        """ remove basepairs in terminal loop regions """
        num_helices = len(helices)
        self.pseduoknots = []
        for i in range(1, num_helices):
            pre_duplex = helices[i-1]
            duplex = helices[i] 
            pre_bp = pre_duplex[-1]
            bp = duplex[0]
            if bp[0] > pre_bp[0] and bp[0] < pre_bp[1] and bp[1] > pre_bp[1]:
                print "# pre_bp: ", pre_bp
                print "# bp: ", bp 
                self.pseduoknots.append(helices[i])
            #elif:
            #    pass 
            
            
        
        
        

    def isBP_has_P(self, bp):
        """ check if both both residues of the basepair have the P:
            if both have P, return 1, else, return 0
        """
        output = 0
        if self.residues[bp[0]].isHasP_O3P() and self.residues[bp[1]].isHasP_O3P():
            output = 1
        #print "# bp: %3d %3d has_p: %d" % (bp[0], bp[1], output)
        return output 
        
            
            
    def removeBP_wo_PO3(self, helix, vb=0):
        """ remove basepair without P or O3' in there nucleotides """

        first_bp = helix[0]
        if self.residues[first_bp[0]].isHasP_O3P() and self.residues[first_bp[1]].isHasP_O3P():
            pass
        else:
            helix = helix[1:]
            if vb: print "# remove bp: ", helix[0]
        _bp = helix[-1]

        if self.residues[ _bp[0]].isHasP_O3P() and self.residues[_bp[1]].isHasP_O3P():
            pass
        else:
            helix = helix[:-1]    
            if vb: print "# remove bp: ", helix[-1]
        return helix


    def getDuplexDSeq(self):
        """ return helices'  digital seqence:
	    GC 0,  CG  1,  AU 2 , UA 3, GU 4,  UG 5
	    input base pair: (13, 83, 'CG')
	    """
        
        helices_dseq = {}
        #print self.duplexs
        duplex_ids = self.duplexs.keys()
        duplex_ids.sort()
        for duplex_id  in duplex_ids:
            #rint "helix: ", self.duplexs[duplex_id]
            helix_dseq = ""
            helix = self.duplexs[duplex_id]
            for _bp in helix:
                helix_dseq += self.bp_dseq[_bp[2]]
            helices_dseq[duplex_id] = (helix_dseq,helix)
        #print helices_dseq
        self.duplex_ids = duplex_ids
        self.duplex_ids.sort()
        return helices_dseq

    
    
    def getHBP_DigitalSeq(self, helix):
        """ return helix digital sequence by basepairs in the helix """
        helix_dseq = ""
        for _bp in helix: helix_dseq += self.bp_dseq[_bp[3]]
        return helix_dseq
    
    def getHBP_ResXYZ(self, helix_resids):
        """ return helix residues' coordinates """
        helices_res_xyz = {}
        for resid in helix_resids:
            #print "resid: ", resid, "iresid: ", self.residues[resid].resid
            res_xyz = self.residues[resid].getResXYZ_Wo_H()
            helices_res_xyz[resid] = res_xyz
        return helices_res_xyz
    
    
    def getDuplexNBP(self, helices, vb=0):
        """ return resids next to duplexes """
        all_nbp_names = []
        for helix in helices:
            first_bp = helix[0]
            lresid = first_bp[0]
            rresid = first_bp[1]
            nlresid = lresid - 1
            nrresid = rresid + 1
            if nlresid in self.resids and nrresid in self.resids:
                nbp_name5 = "%s%s" % (self.residues[nlresid].nt_name, self.residues[nrresid].nt_name)
            else:
                nbp_name5 = "NA"
            last_bp = helix[-1]
            lresid = last_bp[0]
            rresid = last_bp[1]
            nlresid = lresid + 1
            nrresid = rresid - 1
            if nlresid in self.resids and nrresid in self.resids:
                nbp_name3 = "%s%s" % (self.residues[nlresid].nt_name, self.residues[nrresid].nt_name)
            else:
                nbp_name3 = "NA"
            all_nbp_names.append((nbp_name5, nbp_name3))    
        if vb: 
            #print "# pdb: %s " % self.fn_pdb
            print all_nbp_names     
        
            
    def getHelixDSeq(self):
        """ return helices'  digital seqence:
	    GC 0,  CG  1,  AU 2 , UA 3, GU 4,  UG 5
	    bp: (13, 83, 0, 'CG', '+/+ cis ')
	    """
        #self.bp_dseq = {"GC":'0',"CG":'1',"AU":'2',"UA":'3',"GU":'4',"UG":'5'}	
        helices_dseq = []
        all_helices_resids = []
        for helix in self.helices:
            #print "# helix: ", helix
            helix_resids = self.getHBP_Resids(helix)
            #print "# helix_resids: ", helix_resids
            helices_dseq.append(self.getHBP_DigitalSeq(helix))
            all_helices_resids += helix_resids
        #print "# all helices dseq:   ", helices_dseq
        #print "# all helices resids: ", all_helices_resids 
        return helices_dseq, all_helices_resids


    def getBPLnk(self, pre_dseg, next_dseg):
        if pre_dseg[-2:] == next_dseg[:2]:
            #lnk_id = pre_dseg + next_dseg[2:]
            lnk_bp = 2
        elif pre_dseg[-1] == next_dseg[0]:
            #lnk_id = pre_dseg + next_dseg[1:]
            lnk_bp = 1
        return lnk_bp#, lnk_id 

    
    def getLnkBps(self, _bp_seg, pre=1, lnk_bp=2):
        if lnk_bp == 2:
            if pre: 
                lnk_resids = [_bp_seg[-3] , _bp_seg[-2]]
            else:
                lnk_resids = [_bp_seg[0], _bp_seg[1]]
        elif lnk_bp== 1 :
            if pre: 
                lnk_resids = [_bp_seg[-2]]
            else:
                lnk_resids = [_bp_seg[0]]
        #print "# _bpseg: %s  pre= %d  lnk_resid= %s " % (_bp_seg, pre, lnk_resids)	
        return lnk_resids


    def __seek_dseg(self, _dseg):
        db_bp_keys = self.bpseg.keys()
        if _dseg not in db_bp_keys: 
            found = self.seekHseg(_dseg)
            if found: 
                self.bpseg[_dseg] = self.backHseq( _dseg, found )
                self.new_bpseg = 1
                return 1
            else:
                return 0
        return 2 


    def getHelixTemplates(self, helix_id=0):
        """ find helix basepair templates from known strutures """
        seek_hseq, hbp_resids = self.rna_hdeq[helix_id]
        #print "# seek_hhseq: ", seek_hseq
        #print "# hbp_resids: ", hbp_resids 
        print "# helix_id= %s   seek_hseq: %s    hbp_resids= %s " % (helix_id, seek_hseq, hbp_resids)
        bp_seg_ids = []
        bp_segs = []
        num_bp = len(seek_hseq)
        i = 0
        for i in range(num_bp-1): #while i < num_bp:
            _end_id = i+3 
            _dseg = seek_hseq[i:_end_id]
            if self.__seek_dseg(_dseg):
                bp_seg_ids.append(_dseg)
                bp_segs.append(hbp_resids[i:_end_id])
                if _end_id == num_bp: break 
            else:
                _end_id = i+2
                _dseg = seek_hseq[i:_end_id]
                if _dseg == bp_seg_ids[-1][-2:]:continue
                if self.__seek_dseg(_dseg):
                    bp_seg_ids.append(_dseg)
                    bp_segs.append(hbp_resids[i:_end_id])
                    if _end_id == num_bp: break 
                else:
                    print "# could not seg %s in the helix basepair database " % _dseg
                    raise 
        self.saveBP_seg()
        #print "# bp_segs: ", bp_segs 
        #print "# bp_seg_ids: ", bp_seg_ids
        return  bp_seg_ids ,  bp_segs
    
    
    def seekHseg(self, seek_hseq):
        """ find matched helix segment """
        num_sseq = len(seek_hseq)
        found = []
        for t_len in self.temp_lengths:
            if t_len < num_sseq: break 
            templates = self.hdseq_db[t_len]
            for _tt in templates:
                tt_seq = _tt[0]
                if seek_hseq in tt_seq:
                    found.append(_tt)
        return found

    
    def backHseq(self, _seg, _tts):
        """ get matched segment resids 
	    _tt:  ('33232323232322', 0, '1rna')
	"""
        all_bp_segs = []
        for _tt in _tts:
            _rna_str = self.db_2str[_tt[-1]]
            _helix = _rna_str['helices'][_tt[1]]
            start =  _tt[0].find(_seg)
            _bk_segs = []
            for i in range(len(_seg)):
                _bp = _helix[start + i]
                _bk_segs.append( (_bp[0], _bp[1], _bp[3]))
            _bk_segs.append(_tt[-1])
            all_bp_segs.append(_bk_segs)
        return all_bp_segs


    def findBestBPFitTemplates(self,pre_dseg, next_dseg):
        if pre_dseg[-2:] == next_dseg[:2]:
            lnk_id = pre_dseg + next_dseg[2:]
            lnk_bp = 2
        elif pre_dseg[-1] == next_dseg[0]:
            lnk_id = pre_dseg + next_dseg[1:]
            lnk_bp = 1
        else:
            raise "# no link between:  %s -- %s " % (pre_dseg, next_dseg)
        #print "# link_id", link_id
        if lnk_id not in self.LNK_TEMPT.keys():
            lnk_seg = self.bestHbpLnks(lnk_id, lnk_bp, pre_dseg, next_dseg)
            self.LNK_TEMPT[lnk_id] = lnk_seg
            self.new_lnk = 1
        else:
            lnk_seg = self.LNK_TEMPT[lnk_id]
        return lnk_seg, lnk_bp


    def bestHbpLnks(self, lnk_id, lnk_bp, pre_dseg, next_dseg):
        """ find best links between two helix basepairs """
        hbps = self.bpseg
        pre_segs =  hbps[pre_dseg]
        next_segs = hbps[next_dseg]
        num_next_seg = len(next_segs)
        rmsds = []
        rmsd_pairs = []
        cnt = 0
        for i in range(len(pre_segs)):
            _pre_seg = pre_segs[i] #in pre_segs:
            pre_lnk_xyz = self.getBpSegRefXYZ(_pre_seg, pre=1, lnk_bp=lnk_bp)
            for j in range(num_next_seg): # _next_seg in next_segs:
                _next_seg = next_segs[j]
                next_lnk_xyz = self.getBpSegRefXYZ(_next_seg, pre=0, lnk_bp=lnk_bp)
                rmsd = self.calculate_rmsd(pre_lnk_xyz, next_lnk_xyz, vb=0)
                rmsds.append(rmsd)
                rmsd_pairs.append((i,j))
                if rmsd < 0.2: cnt += 1
                if cnt > 5: break 
            if cnt > 5: break  # only allow 20 links with 
        return self.getBest_by_RMSD(rmsds, rmsd_pairs, pre_segs, next_segs)    
        #return best_pairs


    def getBest_by_RMSD(self, rmsds, _pairs, pre_segs, next_segs):
        """ find links with lowest rmsds between the link residues

        """
        rmsds = array(rmsds)
        sort_ids = argsort(rmsds)
        if rmsds[sort_ids[4]] < 0.2:
            num_best_pairs = 5
        else:
            num_best_pairs = 3
        best_pairs = []
        for i in range(5):
            pre_seg_id, next_seg_id = _pairs[sort_ids[i]]
            best_pairs.append((pre_segs[pre_seg_id], next_segs[next_seg_id]))
        return best_pairs
    
    def getBpSegRefXYZ(self, _bp_seg, pre=1, lnk_bp=2):
        #pdb_code = _bp_seg[-1]
        lnk_bps = self.getLnkBps( _bp_seg, pre=pre, lnk_bp=lnk_bp)
        #print "# link_resids: ", lnk_bps 
        return self.__getBpSegRefXYZ(_bp_seg, lnk_bps)
    
    def __getBpSegRefXYZ(self, _bp_seg, lnk_bps):
        """# nt_name= A  atom_id=   7 ref_atom= C1' back_atom= C1' 
           # nt_name= A  atom_id=   9 ref_atom=  C4 back_atom=  C4 
           # nt_name= A  atom_id=   8 ref_atom=  N9 back_atom=  N9 
           # nt_name= C  atom_id=   7 ref_atom= C1' back_atom= C1' 
           # nt_name= C  atom_id=   5 ref_atom= C4' back_atom= C4' 
           # nt_name= C  atom_id=   8 ref_atom=  N1 back_atom=  N1 
        """
        pdb_code = _bp_seg[-1]
        rna_inf = self.db_2str[pdb_code]
        #self.prnDict(rna_inf)
        hres_xyz = rna_inf['hres_xyz']
        ref_xyz = []
        for _bp in lnk_bps:
            left_resid = _bp[0]
            left_nt_name = _bp[2][0]
            ref_xyz += self.__getResRefXYZ(left_nt_name, hres_xyz[left_resid])
            right_resid = _bp[1]
            right_nt_name = _bp[2][1]
            ref_xyz += self.__getResRefXYZ(right_nt_name, hres_xyz[right_resid])
            #print hres_xyz[left_resid]
            #print "# left resid= %3d nt_name= %s right resid= %3d right_nt_name= %s " % (left_resid, left_nt_name, right_resid, right_nt_name)
        return array(ref_xyz)  



    def __getResRefXYZ(self, res_name, res_xyz): 
        _res_ref_ids = self.ref_atom_ids[res_name]
        ref_xyz = []
        for atom_id in _res_ref_ids:
            ref_xyz.append(res_xyz[atom_id,:].copy())
        return ref_xyz


    def getHBP_Resids(self, _hbp):
        """ Return helix basepair resids 
	(80, 144, 'UA'), (81, 143, 'CG'), (82, 142, 'CG'), '2h0x']
	"""
        resids  = []
        for _bp in _hbp:
            resids += [_bp[0], _bp[1]]
        return resids 


    def setHelixAnchorSeg(self, _bp_tempts, _bp_seg):
        self.new_helix = {}
        self.new_helix['bps'] = []
        print "# helix anchor seg: ", _bp_seg 
        #print "# templates: ", _bp_tempts 
        num_bps = len(_bp_seg)
        assert len(_bp_tempts)-1 == num_bps
        hres_xyz = self.bp2ResXYZ(_bp_tempts)
        for i in range(num_bps):
            _bpt = _bp_tempts[i]
            _bp = _bp_seg[i] 
            self.new_helix[_bp[0]] = hres_xyz[_bpt[0]]
            self.new_helix[_bp[1]] = hres_xyz[_bpt[1]]
            self.new_helix['bps'].append(_bp)
            #print "# _bpt= %s _bp= %s " % (_bpt, _bp )
        #print self.helix.keys()
        

    def bp2ResXYZ(self, _bpseg):
        pdb_code = _bpseg[-1]
        rna_inf = self.db_2str[pdb_code]
        hres_xyz= rna_inf['hres_xyz']
        return hres_xyz 


    def getBpSegResXYZ_WoH(self, _bp_seg, list=0):
        pdb_code = _bp_seg[-1]
        fn = self.db_2str[pdb_code]['fn']
        seg_resids = self.getHBP_Resids(_bp_seg)
        seg_xyz = pdb_obj.getSegXYZ_Wo_H(seg_resids,list=list)
        return  seg_xyz


    def getHBP_AnchorRef(self, lnk_bp):
        """ return helix basepair anchor reference frame, previous segment last lnk basepairs' reference frame """
        _bps = self.new_helix['bps']
        if lnk_bp == 2:
            pre_lnk_resids = [_bps[-2][0],_bps[-2][1], _bps[-1][0], _bps[-1][1]]
            nt_names = _bps[-2][2]+_bps[-1][2]
        else:
            pre_lnk_resids = [ _bps[-1][0], _bps[-1][1]]
            nt_names = _bps[-1][2]
        #print "# lnk_bp= %3d ref_resid= %s " % (lnk_bp, pre_lnk_resids)    
        ref_xyz = []
        for i in range(len(pre_lnk_resids)):
            resid = pre_lnk_resids[i]
            nt_name = nt_names[i]
            #print "# ref: resid= %3d nt_name= %s  (bp)nt_name= %s " % (resid, self.seq_dict[resid], nt_name)
            ref_xyz += self.__getResRefXYZ(nt_name, self.new_helix[resid])
        #print "# ref_xyz: ", ref_xyz 
        return array(ref_xyz)
    
    

    def addBpSeg2Helix(self, _bptseg, bpseg,  hm, lnk_bp):
        #print "# before next seg: self.new_helix: ", self.new_helix.keys()
        #print "# before next seg: new_helix_bps: ", self.new_helix['bps']
        #print "bptseg: ", _bptseg, "bpseg: ", bpseg 
        hres_xyz = self.bp2ResXYZ(_bptseg)
        if len(bpseg) == 2:
            lnk_res_ids = [0,1]
        else:
            if lnk_bp == 2:
                lnk_res_ids = [1,2]
            else:
                lnk_res_ids = [0,1,2]
        print "# lnk_resids= %s lnk_bp= %s  bpseg= %s " % (lnk_res_ids, lnk_bp, bpseg)
        first = 1
        for i in lnk_res_ids: #range(1, len(bpseg)):
            #print "# bpt: %s  bp= %s "  % (_bptseg[i], bpseg[i])
            #print "##** bpseg: ", bpseg 
            self.new_helix[bpseg[i][0]] = self.homoMove(hres_xyz[_bptseg[i][0]], hm)
            #self.new_helix[bpseg[i][1]] = self.homoMove(hres_xyz[_bptseg[i][1]], RotM, center1,  center2)
            self.new_helix[bpseg[i][1]] = self.homoMove(hres_xyz[_bptseg[i][1]], hm)
            if first: #i == 1:
                self.new_helix['bps'][-1] = bpseg[i]
                first = 0
            else:
                self.new_helix['bps'].append(bpseg[i])
        #print "# add next seg: self.new_helix: ", self.new_helix.keys()
        #print "# add next seg: new_helix_bps: ", self.new_helix['bps']

    def bestLnkSeg_By_Ref(self, lnk_bp, next_dseg, ref_xyz):
        hbps = self.bpseg
        next_segs = hbps[next_dseg]
        num_next_seg = len(next_segs)
        rmsds = []
        rmsd_pairs = []
        cnt = 0
        for i in range(num_next_seg):
            _next_seg = next_segs[i]
            next_lnk_xyz = self.getBpSegRefXYZ(_next_seg, pre=0, lnk_bp=lnk_bp)
            rmsd = self.calculate_rmsd(ref_xyz, next_lnk_xyz, vb=0)
            rmsds.append(rmsd)
            rmsd_pairs.append(i)
            if rmsd < 0.25: cnt += 1
            if cnt > 5: break 

        rmsds = array(rmsds)
        sort_ids = argsort(rmsds)
        if len(sort_ids) > 5 and rmsds[sort_ids[4]] < 0.2:
            num_best_pairs = 5
        else:
            num_best_pairs = min(3,len(sort_ids))
        best_pairs = []
        for i in range(num_best_pairs):
            next_seg_id = rmsd_pairs[sort_ids[i]]
            best_pairs.append(next_segs[next_seg_id])
        return best_pairs

    def helix2pdb(self, helix):
        bps = helix['bps']
        helix_res = {}
        for _bp in bps:
            helix_res[_bp[0]] = _bp[2][0]
            helix_res[_bp[1]] = _bp[2][1]
        helix_resids = helix_res.keys()
        helix_resids.sort()
        fmt_pdb = '%-6s%5d %-4s %-4s%1s%4d%1s%3s%8.3f%8.3f%8.3f%18s%-4s'
        lns = []
        atom_id = 0
        for resid in helix_resids:
            print "# resid= %3d nt_name= %s " % (resid, helix_res[resid])
            nt_name = helix_res[resid]
            #nt_atoms_woh[nt_name]
            nt_atom_names = nt_atoms_woh[nt_name]
            res_xyz = helix[resid]
            for i in range(len(nt_atom_names)):
                atom_name = nt_atom_names[i]
                atom_xyz = res_xyz[i,:]
                atom_id += 1
                ln =  fmt_pdb  % ('ATOM', atom_id, atom_name, nt_name," ",  resid, "","", atom_xyz[0],atom_xyz[1], atom_xyz[2], "", "") #_atom.x, _atom.y, _atom.z, _atom.misc, _atom.segid)
                lns.append(ln)
                #print ln 
        return lns 
    


    
class Nucleotide(RESIDUE, Rotation):
    def __init__(self, **rest): #name="", chainID="", resid="", segid=""):
        RESIDUE.__init__(self, **rest)
        #self.nt_name = self.name[0].upper()
        self.sugar_dih_atoms = {'mu0':("C4'", "O4'", "C1'", "C2'"),
                                "mu1":("O4'", "C1'", "C2'", "C3'"),
                                "mu2":("C1'", "C2'", "C3'", "C4'"),
                                "mu3":("C2'", "C3'", "C4'", "O4'"),
                                "mu4":("C3'", "C4'", "O4'", "C1'")}
        self.default_atom_xyz = _array([-999., -999., -999.])
        self.setNT_name()
        
    def setNT_name(self):
        name = self.name.upper()
        if len(name) == 1:
            self.nt_name = name[0]
        elif name in ['ADE','GUA', 'CYT','URI']:
            self.nt_name = name[0]
        else:
            self.nt_name = name[-1] # mutated nucleotide
       #print "# nt: resname= %4s  nt_name= %4s " % (self.name, self.nt_name)

       
    def isRNANT(self):
        """ check if it is a valid RNA nucleotide: nt_name in ['A',G','C','U']"""
        if self.nt_name in ['A','G','C','U']:
            return True 
        else:
            return False 
        
            
    def Syn_Anti(self):
        """ chi(R): O4'-C1'-N9-C4; chi(Y): O4'-C1'-N1-C2  """
        self.sugar_syn = 0
        chi_angle = self.getNtXhi()
        if chi_angle >= -90 and chi_angle <= 90:
            self.sugar_syn = 1
        else:
            self.sugar_syn = 0
        print "# sugar_syn: %s" % self.sugar_syn


    def getNtXhi(self):
        """ O4' - C1' - N9' - C4
        """
        ag_atom_names = ("O4'", "C1'", "N9","C4")
        cu_atom_names = ("O4'", "C1'", "N1","C2")
        if self.name in ["GUA", "ADE"]:
            ans = ag_atom_names
        else:
            ans = cu_atom_names
        chi =  dihedral_vect(self.getAtomVect(ans[0]),self.getAtomVect(ans[1]), 
                             self.getAtomVect(ans[2]), self.getAtomVect(ans[3]))
        chi = degrees(chi)
        print "# resid= %4d resname= %4s chi= %7.1f " % (self.resid, self.name, chi)
        return chi 

    def getXBond(self):
        """ return X bond:  C1' - N9'/N1 """
        ag_atom_names = ("C1'", "N9")
        cu_atom_names = ("C1'", "N1")
        if self.name in ["GUA", "ADE"]:
            ans = ag_atom_names
        else:
            ans = cu_atom_names
        chi_bond = self.getAtomVect(ans[1]) - self.getAtomVect(ans[0]) 
        return chi_bond



    def getRingAnames(self):
        if self.nt_name in ["G", "A"]:
            ring_anames = ["C4", "N3", "C2", "N1", "C6", "C5", "N7", "C8", "N9"]
            self.Nsugar = 'N9'
        else:
            ring_anames = ["C4", "N3", "C2", "N1", "C6", "C5"]
            self.Nsugar = 'N1'
        return ring_anames

    def analyzeImino(self):
        if self.nt_name in ['A','C']: return 
        NH = self.getImino()
        C2 = self.getAtomVect('C2')
        if self.nt_name == 'G':
            Natom = 'N1'
            vatom = "C6"
        elif self.nt_name == 'U':
            Natom = 'N3'
            vatom = "C4"
        N = self.getAtomVect(Natom)
        C4_6 = self.getAtomVect(vatom)
        u = self.normVect(C2 - N )
        v = self.normVect(C4_6 - N )
        norm = cross(u, v)
        NH_n = self.normVect(NH)
        n_ang = self.getVecAngle(NH_n, norm)
        c2_ang = self.getVecAngle(NH_n,u)
        c4_6_ang = self.getVecAngle(NH_n, v)
        vn_ang = self.getVecAngle(v, norm)
        un_ang = self.getVecAngle(u, norm)
        uv_ang = self.getVecAngle(u,v)
        #uv_ave = C2 - N + C4_6 - N #u+v
        uv_ave = -(u+v)
        uv_nh = self.getVecAngle(self.normVect(uv_ave), NH_n)
        print "# %3d %s NH= %6.3f C2-NH=%6.1f C4/6-NH=%6.1f u-N-v=%6.1f n-N-H=%6.1f n-N-C2=%6.1f n-N-C4/6=%6.1f u+v_NH= %6.1f" % \
              (self.resid, self.nt_name, self.getVecLen(NH),  degrees(c2_ang), degrees(c4_6_ang), degrees(uv_ang), degrees(n_ang), degrees(un_ang), degrees(vn_ang), degrees(uv_nh) ),
        self.prnVec(NH, "#pdb")

    def addImino(self):
        """ add H atoms to crystall structur by xplor """
        if self.nt_name in ['A','C']: return None
        C2 = self.getAtomVect('C2')
        if self.nt_name == 'G':
            Natom = 'N1'
            vatom = "C6"
            Hatom = 'H1'
        elif self.nt_name == 'U':
            Natom = 'N3'
            vatom = "C4"
            Hatom = 'H3'
        N = self.getAtomVect(Natom)
        C4_6 = self.getAtomVect(vatom)
        u = self.normVect(C2 - N )
        v = self.normVect(C4_6 - N )
        norm = cross(u, v)
        NH_n = - self.normVect(u+v)
        if self.nt_name == 'G': 
            NH_n = self.new_arb_rot(NH_n, norm, -1.4)
        NH = NH_n.copy() * 1.010 
        H_atom = NH + N
        #H_x = self.getAtomVect(Hatom)
        #self.prnVec(H_x, cmt='H_x   ')
        #self.prnVec(H_atom, cmt='H_atom')
        #self.__analyzIminoNH(u,v, norm, NH_n, NH)
        return NH 


    def __analyzIminoNH(self, u, v, norm, NH_n, NH, cmt="addNH"):        
        n_ang = self.getVecAngle(NH_n, norm)
        c2_ang = self.getVecAngle(NH_n,u)
        c4_6_ang = self.getVecAngle(NH_n, v)
        vn_ang = self.getVecAngle(v, norm)
        un_ang = self.getVecAngle(u, norm)
        uv_ang = self.getVecAngle(u,v)
        #uv_ave = C2 - N + C4_6 - N #u+v
        uv_ave = -(u+v)
        uv_nh = self.getVecAngle(self.normVect(uv_ave), NH_n)
        print "# %3d %s NH= %6.3f C2-NH=%6.1f C4/6-NH=%6.1f u-N-v=%6.1f n-N-H=%6.1f n-N-C2=%6.1f n-N-C4/6=%6.1f u+v_NH= %6.1f" % \
              (self.resid, self.nt_name, self.getVecLen(NH),  degrees(c2_ang), degrees(c4_6_ang), degrees(uv_ang), degrees(n_ang), degrees(un_ang), degrees(vn_ang), degrees(uv_nh) ), 
        self.prnVec(NH, cmt)

    def getRingNorm(self):
        """ return ring norm by u x v 
        the base reference frame frame x,y,z  of the reference bases from u, v vector.
	u = N1(Y)-->C2,   N9(R)--> C4
	v = N1(Y)-->C6,   N9(R)--> C8
	y = u + v 
	z = u x v
	x = y x z
	A: N9,C4,C8, G: N9, C4, C8,  C:N1,C2,C6, U:N1,C2,C6 """
        if self.nt_name in ['A','G']:
            N9 = self.getAtomVect('N9')
            C4 = self.getAtomVect('C4')
            C8 = self.getAtomVect('C8')
            u =  C4 - N9 
            v =  C8 - N9 
        else:
            N1 = self.getAtomVect('N1')
            C2 = self.getAtomVect('C2')
            C6 = self.getAtomVect('C6')
            u = C2 - N1 
            v = C6 - N1 
        u  = self.normVect(u)
        v = self.normVect(v)
        y = self.normVect( u + v )
        z = self.normVect( cross(u,v))
        #x = self.normVect( cross(y,z))
        return z
    
    
    def getC6C8(self):
        if self.nt_name in ['A','G']:
            return self.getAtomVect('C8')
        else:
            return self.getAtomVect('C6')
        

    def getBaseNorm(self):
        return self.getRingNorm()

    def getNsugar_XYZ(self):
        """ return Glycosidic Nitrogen coordinates(x,y,z) """
        return getattr(self, self.Nsugar)


    def getRingCenter(self):
        return self.ring_center

    def getFitBaseCenter(self):
        """ return the center of fitted base """
        return self.fitted_base_center


    def lsfitBase(self):#, vecs1, vecs2):
        """ least quare fitting of reference base to the base in PDB """
        rbase = rf_bases[self.nt_name]
        E0 = Npy.sum(self.ring_xyz_ct * self.ring_xyz_ct) + rbase[2] 
        correlation_matrix = Npy.dot(rbase[1], self.ring_xyz_ct)      
        V, S, W_trans = svd(correlation_matrix)      
        is_reflection = (det(V) * det(W_trans)) < 0.0      
        if is_reflection: S[-1] = -S[-1]                   # reflect along smallest principal axis      
        rmsd_sq = (E0 - 2. * Npy.sum(S)) / float(self.num_atom_in_base)      
        rmsd_sq = max([rmsd_sq, 0.])      
        if is_reflection:  V[-1,:] = V[-1,:] * (-1.0)      
        optimal_rotation = Npy.dot(V, W_trans)      
        #fitted_base_xyz = Npy.dot(rbase[3], optimal_rotation) + self.base_center #rbase[0]
        
        self.base_rmsd =  sqrt(rmsd_sq)
        
        #print "# fitting_base_xyz: ", fitted_base_xyz
        #print "# ring_xyz: ", self.base_xyz
        #print fitted_base_xyz - self.base_xyz
        #print "Optimal_rotation: ", optimal_rotation
        
        self.ORot = optimal_rotation
        self.fitted_base_center = self.ring_center - dot(rbase[0], optimal_rotation)  # colum vectors 
        
        #print "E fit bc: ", self.fitted_base_center
        #print " dot: ",  rbase[0] * transpose(optimal_rotation)


    def setBase(self):
        self.base_anames = base_anames[self.nt_name]
        self.setRingCenter(vb=0)
        self.lsfitBase()
        self.HA = H_Acceptor[self.nt_name]
        self.HD = H_Donor[self.nt_name]

        
    def setRingCenter(self, vb=0):
        """  /* get the coordinate of ring center*/"""
        ring_anames = self.getRingAnames()
        ring_xyz = []
        for atom_name in ring_anames:
            ring_xyz.append(self.getAtomVect(atom_name))
        ring_xyz = Npy.array(ring_xyz)    
        n_vec, vec_size = ring_xyz.shape
        ring_center = Npy.sum(ring_xyz, 0)/float(n_vec)
        self.ring_center = ring_center 
        self.ring_xyz = ring_xyz
        self.ring_xyz_ct = ring_xyz.copy() - ring_center 
        self.num_atom_in_base = n_vec
        if vb: print "# base center: ",  self.ring_xyz_ct


    def getNC(self):
        """    get the vector of N1(or N9)-->C1' for the given base.*/
	"""
        C_atom = self.getAtomVect("C1'")
        if self.nt_name in ["A", "G"]:       
            N_atom = self.getAtomVect("N9")
        elif self.nt_name in ["C","U"]:
            N_atom = self.getAtomVect("N1")
        NC_vect = C_atom - N_atom  # get vector of N1(or N9)->C1' for the first base.
        return NC_vect 


    def getResXYZ_Wo_H(self, list=0):
        """ return residue atom  cartesian coordinates  without H"""
        res_xyz = []
        #res_xyz_dt = {}
        for _atom_name in nt_atoms_woh[self.nt_name]: #self.atom_names:
            if _atom_name == "OP1":
                try: 
                    _atom = getattr(self, _atom_name)
                except:
                    _atom = getattr(self, "O1P", self.default_atom_xyz.copy())
            elif _atom_name == "OP2": 
                try: 
                    _atom = getattr(self, _atom_name)
                except:
                    _atom = getattr(self, "O2P", self.default_atom_xyz.copy())
            else:
                try:
                    _atom = getattr(self, _atom_name)
                except:
                    if _atom_name in ['P', 'O1P', 'O2P', 'OP1', 'OP2']:
                        print "# resid= %3d has no  < P > atom "  % self.resid 
                        _atom = _array([-999., -999., -999.])
                        self.prnVec(_atom, cmt="# use defaut P value ")
                    else:
                        print "# could not find  <%s> in residue %3d-%s " % (_atom_name, self.resid, self.name)
                        raise " Atom name error "
            res_xyz.append(_atom)
        if list: 
            return res_xyz
        else:
            return _array(res_xyz).copy()


    def getResXYZ_WoH_atom(self, list=0):
        """ return residue atom  cartesian coordinates  without H"""
        res_xyz = {}
        #res_xyz_dt = {}
        for _atom_name in nt_atoms_woh[self.name[0]]: #self.atom_names:
            if _atom_name[0] != 'H':
                res_xyz[_atom_name] =getattr(self, _atom_name)
        return res_xyz
    
    def resXYZ2Atom(self, res_xyz):
        nt_atom_names = nt_atoms_woh[self.nt_name]
        #print "nt_atoms: ", nt_atom_names
        for i in range(len(nt_atom_names)): 
            _atom_name = nt_atom_names[i]
            setattr(self, _atom_name, res_xyz[i,:])
            #print "# _atom_name: ", _atom_name
        
    def getResRefXYZ(self, list=0):
        """ return reference frame coordinates (C1', C4', N1(Y)/N9(R)"""
        if self.name[0] in ['G','A']:
            ref_atoms = ["C1'","C4","N9"]
        else:
            ref_atoms = ["C1'","C4'","N1"]
        #print "#***  reference resid: %3d   %s " % (self.resid,  self.name )
        res_xyz = []   
        for _atom in ref_atoms:
            res_xyz.append(getattr(self, _atom))
        if list: 
            return res_xyz
        else:
            return _array(res_xyz)

    
    def isHasP_O3P(self):
        """ if residue has not P or O3', this residue is not included in the database """
        if "P" not in self.atom_names or "O3'" not in self.atom_names:
            return False
        elif "O1P" not in self.atom_names and "OP1" not in self.atom_names:
            return False
        else:
            return  True

    
    def getImino(self, vb=0):
        """ return imino NH bond vector """
        #print "# getting imino from  N-->H bond vector"""
        if self.nt_name == "G":
            H_atom_name = 'H1'
            N_atom_name = 'N1'
        elif self.nt_name == 'U':
            H_atom_name = 'H3'
            N_atom_name = 'N3'
        else:
            return None 
        N_atom = self.getAtomVect(N_atom_name, vb=vb)
        H_atom = self.getAtomVect(H_atom_name, vb=vb)
        if N_atom_name is None or H_atom is None:  return None 
        NH = H_atom - N_atom
        return NH


    def getBeta(self):
        """ i P - O5' - C5' - C4' : return beta in degrees
        """
        return self.getDihedraAngle_by_Atoms(("P", "O5'", "C5'", "C4'"))
    

    def getGamma(self):
        """  gamma: O5' - C5' - C4' - C3' : return gamma in degrees 
        """
        return self.getDihedraAngle_by_Atoms(("O5'", "C5'", "C4'", "C3'"))

    def getDelta(self):
        """ C5' - C4' - C3'- O3', return delta angle in degrees
        """
        return self.getDihedraAngle_by_Atoms(("C5'", "C4'", "C3'", "O3'"))


    def getAtomXYZ(self,atom_name):
        #if not self.prime_format and atom_name[-1] == "'":
        #   _atom_name = "%s*" % atom_name[:-1]
        #else:
        #    _atom_name = atom_name
        #return getattr(residue, _atom_name)
        return getattr(self, atom_name)

   
    def getXhi(self):
        """ O4' - C1' - N9' - C4
        """
        if self.nt_name in ['G','A']:
            atom_names = ["O4'", "C1'", "N9", "C4"]
        elif self.nt_name in ['U','C']:
            atom_names = ["O4'", "C1'", "N1", "C2"]
        atoms_xyz = []
        for atom_name in atom_names: 
            atoms_xyz.append(self.getAtomVect(atom_name))
        xhi =  dihedral_vect(atoms_xyz[0],atoms_xyz[1], atoms_xyz[2],atoms_xyz[3])
        return degrees(xhi)


    def getMu1(self):
        """ i O4' - C1' - C2' - C3' 
        """
        #atom_names = [(0,"O4'"),(0,"C1'"), (0, "C2'"),(0,"C3'")]
        
        #self.prnDihAngles(angles, cmt="mu1:  O4' - C1' - C2' - C3' ")
        #self.mu1 = angles 
        
        return self.getDihedraAngle_by_Atoms(("O4'", "C1'", "C2'", "C3'"))


    def getMu2(self):
        """ i C1' - C2' - C3' - C4' 
        """
        #atom_names = [(0,"C1'"),(0,"C2'"), (0, "C3'"),(0,"C4'")]
        return self.getDihedraAngle_by_Atoms(("C1'", "C2'", "C3'", "C4'"))
        


    def toPdb(self, atom_id=0, vb=0):
        """ write residue's coordinates in PDB format """
        fmt_pdb = '%-6s%5d %-4s %-4s%1s%4d%1s%3s%8.3f%8.3f%8.3f%18s%-4s'
        lns = []
        atom_id = atom_id
        nt_name = self.nt_name 
        if vb: print "# resid= %3d nt_name= %s " % (self.resid, self.nt_name)
        for atom in self.atom_list:
            atom_name = atom.name4Pdb() #
            atom_id += 1
            ln =  fmt_pdb  % ('ATOM', atom_id, atom_name, nt_name," ",  self.resid, "","", atom.x,atom.y, atom.z, "", "") #_atom.x, _atom.y, _atom.z, _atom.misc, _atom.segid)
            lns.append(ln)
                #print ln 
        return lns, atom_id  



    def updateAtomXYZ_by_aname(self, aname, _xyz):
        aid = self.atom_names.index(aname)
        atom_obj = self.atom_list[aid]
        assert atom_obj.name == aname 
        atom_obj.setXYZ(xyz=_xyz)
        setattr(self, aname, Npy.array(_xyz))
        #  self.atom_list.append(atom)
        #setattr(self, atom.name, atom.axyz)

    def getBB_BondLengths(self):
        """ return bond lengths of RNA backbone"""
        if "P" in self.atom_names:
            PO5p =  self.getBondLen("O5'", "P")
        else:
            PO5p = -1.0
        O5C5p = self.getBondLen("C5'", "O5'") 
        C5C4p = self.getBondLen("C4'", "C5'")
        C4C3p = self.getBondLen("C3'", "C4'")
        C3O3p = self.getBondLen("O3'", "C3'")
        return [PO5p, O5C5p, C5C4p, C4C3p, C3O3p]

    def getBB_BondAngles(self):
        """ return bond angles of RNA backbone"""
        if "P" in self.atom_names:
            PO5p =  self.getBondVect("O5'", "P")
            has_P = 1
        else:
            has_P = 0
        O5C5p = self.getBondVect("C5'", "O5'") 
        C5C4p = self.getBondVect("C4'", "C5'")
        C4C3p = self.getBondVect("C3'", "C4'")
        C3O3p = self.getBondVect("O3'", "C3'")
        if has_P: 
            #print "# O5C5P: ", O5C5p, "PO5P", PO5p
            P_O5_C5 = self.getVecAngle(-PO5p, O5C5p)
        else:
            P_O5_C5 = -999.0
        O5_C5_C4 = self.getVecAngle(-O5C5p, C5C4p)
        C5_C4_C3 = self.getVecAngle(-C5C4p, C4C3p)
        C4_C3_O3 = self.getVecAngle(-C4C3p, C3O3p)
        return [P_O5_C5, O5_C5_C4, C5_C4_C3, C4_C3_O3]

    
    def getBondAngles_By_atoms(self, ba_atoms):
        """ return bond angles by ba_atoms, 3 successive atoms for a bond angle """
        bond_angles = []
        for  _batoms in ba_atoms:
            v1 = self.getBondVect(_batoms[1], _batoms[0])
            v2 = self.getBondVect(_batoms[1], _batoms[2])
            bond_angles.append(degrees(self.getVecAngle(v1, v2)))
        return bond_angles     

    
    def getTetra(self, tetra_atoms):
        """ analyze tetra center: tetra_atoms: [tetra_center atoms, atom1, atom2, atom3, atom4 ]
	    known_atom[1,2,3, 4] are bonded to tetra_center_atom 
	"""
        v1 = self.getBondVect(tetra_atoms[0], tetra_atoms[1])
        v2 = self.getBondVect(tetra_atoms[0], tetra_atoms[2])
        av = (v1 + v2)/2.0 
        v3 = self.getBondVect(tetra_atoms[0], tetra_atoms[3])
        v4 = self.getBondVect(tetra_atoms[0], tetra_atoms[4])
        #ang = degrees(self.getVecAngle(av, v3))
        av2 = (v3 + v4)/2.0 
        ang = degrees(self.getVecAngle(av, av2))
        return ang 


    def getSugarBondAngles(self):
        """ return bond angles of sugar ring atoms """
        O4_C4_C3 = self.getBondAngle_By_Atoms("O4'", "C4'", "C3'")
        O4_C4_C5 = self.getBondAngle_By_Atoms("O4'", "C4'", "C5'")

        C2_C3_C4 = self.getBondAngle_By_Atoms("C2'", "C3'", "C4'")
        C2_C3_O3 = self.getBondAngle_By_Atoms("C2'", "C3'", "O3'")
        C3_C2_C1 = self.getBondAngle_By_Atoms("C3'", "C2'", "C1'")
        C2_C1_O4 = self.getBondAngle_By_Atoms("C2'", "C1'", "O4'")
        C1_O4_C4 = self.getBondAngle_By_Atoms("C1'", "O4'", "C4'")
        return O4_C4_C3, O4_C4_C5, C2_C3_C4, C2_C3_O3, C3_C2_C1, C2_C1_O4, C1_O4_C4

    
    def getSugarBondLengths(self):
        sugar_atoms = ["C4'", "C3'", "C2'", "C1'", "O4'"]
        sugar_bonds = [("C4'", "C3'"), ("C3'", "C2'"), ("C2'", "C1"), ("C1'", "O4'")]


    def getSurgarDihAngles(self):
        sugar_dih_keys = ["mu0", "mu1","mu2", "mu3", "mu4"]
        #sugar_dih_angles = {}
        sugar_dih_angles = []
        for _sdk in sugar_dih_keys:
            dih_atoms = self.sugar_dih_atoms[_sdk]
            sugar_dih_angles.append( self.getDihedraAngle_by_Atoms(dih_atoms))
        return sugar_dih_angles

    
    def chkMissAtom(self):
        """ check missing heavy atoms """
        nt_atoms = nt_atoms_woh[self.nt_name]
        missing_atoms = []
        for _atom in nt_atoms:
            if _atom not in self.atom_names:
                missing_atoms.append(_atom)
        if len(missing_atoms) > 0:
            self.miss_atom = 1
            self.miss_atoms = missing_atoms
        else:
            self.miss_atom = 0
    
            


            
    
    def getBBXYZ(self):
        """ return residue backbone atoms' cartesian coordinates """
        res_xyz = []
        for _atom_name in nt_bb_atoms:
            atom_xyz = getattr(self, _atom_name, None)
            if atom_xyz: res_xyz.append(atom_xyz)
        return res_xyz
    
        

class RNAPDB(PDB, Rotation): #, MVector):
    def __init__(self, *arg, **argv):
        PDB.__init__(self, *arg, **argv)
        
    def isShortContact(self, res1, res2):
        short_contact = 0
        #base1_xyz = res1.ring_xyz
        #base2_xyz = res2.ring_xyz
        for  atom_name1 in res1.base_anames:
            if atom_name1[0] not in ["O","N"]: continue 
            w_H1, wo_H1 = self.HCatalog(res1.nt_name, atom_name1)
            atom1_xyz = res1.getAtomVect(atom_name1)
            for atom_name2 in res2.base_anames: #j in range(res2.num_atom_in_base):
                if atom_name2[0] not in ["O","N"]: continue 

                w_H2, wo_H2 = self.HCatalog(res2.nt_name, atom_name2)
                if wo_H1 == 1 and wo_H2: continue 

                dxyz = self.getVecLen(atom1_xyz - res2.getAtomVect(atom_name2))
                #print "# dxzy= %8.3f " % dxyz 
                if dxyz <= self.MAX_HB:   # H-bond upper limit 
                    short_contact = 1
                    break 
        return short_contact         


    def calculate_rmsd(self, conformation1, conformation2,       
                       is_align = False, is_center = True, vb = False, pdb3p=1):      
        """      
	Calculates the RMSD between two conformations.      

	  * conformation1: array of dimensions [N,3]       
	  * conformation2: array of dimensions [N,3]       
	  * is_align: True = modify conformation2 to be 
			     aligned to conformation1      
	  * is_center: True = remove average center differences      
			      before calculating RMSD       
	"""      
        #print type(conformation1)
        nrow, ncol = conformation1.shape
        if nrow % 3 == 2:
            conformation1 =conformation1[2:,:]
        nrow, ncol = conformation2.shape
        if nrow % 3 == 2:
            conformation2 = conformation2[2:, :]

        if not shape(conformation1) == shape(conformation2):      
            print "# conformation1 shape: ", conformation1.shape
            print "# cofnormation2 shape: ", conformation2.shape
            raise "calculate_rmsd: conformations not the same size."      
        elif not len(shape(conformation1)) == 2:      
            raise "calculate_rmsd: conformations not the correct rank."      

        vecs1 = conformation1.copy()      
        vecs2 = conformation2.copy()      

        n_vec, vec_size = shape(vecs1)      
        if is_center:      
            center1 = sum(vecs1,0) / float(n_vec)      
            center2 = sum(vecs2,0) / float(n_vec)      
            vecs1 -= center1      
            vecs2 -= center2      

            E0 = sum(sum(vecs1 * vecs1)) + sum(sum(vecs2 * vecs2))      

            correlation_matrix = dot(transpose(vecs2), vecs1)      

            V, S, W_trans = svd(correlation_matrix)      

            is_reflection = (det(V) * det(W_trans)) < 0.0      
            if is_reflection:      
                # reflect along smallest principal axis      
                S[-1] = -S[-1]      

            rmsd_sq = (E0 - 2. * sum(S)) / float(n_vec)      
            rmsd_sq = max([rmsd_sq, 0.])      

            if is_align:      
                if is_reflection:       
                    V[-1,:] = V[-1,:] * (-1.0)      
                optimal_rotation = dot(V, W_trans)      
                conformation2[:,:] = dot(vecs2, optimal_rotation) + center2    
            rmsd = sqrt(rmsd_sq)
            #if vb: print "#  rmsd= %8.4f " %  rmsd 
            print "conformation2.shape"
            return rmsd 


    def ntRMSD(self, resid1, resid2):
        """ least quare fitting of reference base to the base in PDB """
        res1 = self.residues[resid1]
        res2 = self.residues[resid2]
        res1_xyz = res1.getResXYZ()
        res2_xyz = res2.getResXYZ()
        #print type(res1_xyz)
        rmsd = self.calculate_rmsd(res1_xyz, res2_xyz)
        print "# %3d --%3d  %s -- %s  rmsd= %7.3f " % (resid1, resid2, res1.nt_name, res2.nt_name, rmsd)

    def getSegXYZ(self, seg):
        """ return segment's cartesian coordinates """
        seg_xyz = []
        seg_seq = ""
        for resid in seg:
            seg_xyz += self.residues[resid].getResXYZ(list=1)
            seg_seq += self.residues[resid].nt_name

        return _array(seg_xyz), seg_seq

    def prnDict(self, _dict):
        dkeys = _dict.keys()
        dkeys.sort()
        for _key in dkeys:
            print "# %s - > %s " % (_key, _dict[_key])


    def getSegXYZ_Wo_H(self, seg, list=1):
        seg_xyz = []
        for resid in seg:
            _seg_xyz = self.residues[resid].getResXYZ_Wo_H(list=1)
            seg_xyz += _seg_xyz
        if list:
            return seg_xyz
        else:
            return _array(seg_xyz)


    def getSegResXYZ_WoH_atom(self, resids=None, ressids=None, list=1):
        seg_res_xyz = {}
        if resids:
            for resid in resids:
                res = self.getRes_by_resid(resid)
                seg_res_xyz[resid] = res.getResXYZ_Wo_H(list=list)
        else:
            for ressid in ressids:
                res = self.getRes_by_ressid(ressid)
                print res.name 
                seg_res_xyz[res.resid] = res.getResXYZ_Wo_H(list=list)
                #except:
            #    print "# resid= %3d " % resid 
            #    raise 
            #seg_res_xyz[resid] = self.residues[resid].getResXYZ_WoH_atom()
        return seg_res_xyz
    
    def getRes_by_ressid(self, ressid):
        return self.residues[ressid]
    
    
    def getPDBBBXYZ(self, resids=None):
        """ return rna pdb backbone atoms' xyz """
        pdb_bb_xyz = []
        if resids is None:
            resids = self.resids 
        for res_sid in resids:
            res = self.residues[res_sid]
            pdb_bb_xyz += res.getBBXYZ()
        return Npy.array(pdb_bb_xyz)
    
        
    def getSegResXYZ_WoH_dict(self, seg):
        seg_res_xyz = {}
        for resid in seg:
            #seg_res_xyz[resid] = self.residues[resid].getResXYZ_Wo_H(list=1)
            seg_res_xyz[resid] = self.residues[resid].getResXYZ_Wo_H()
        return seg_res_xyz

    def getResSeq(self, resid):
        return self.seq_dict[resid]

    
    def getSegSeq(self, _seg_resids):
        """ input: RNA segment resids
	    return RNA segment sequence """
        _seg_seq = ""
        for resid in _seg_resids:
            _seg_seq += self.getResSeq(resid)
        return _seg_seq

    def getSegRefXYZ(self, seg, list=0):
        """ return reference frame coordinates of nucleotides, C1', C4', N1/N9"""
        seg_xyz = []
        for resid in seg:
            #print "# **********"
            seg_xyz += self.residues[resid].getResRefXYZ(list=1)
        if list:
            return seg_xyz
        else:
            return _array(seg_xyz)



    def getRes_by_resid(self, resid):
        return self.residues[resid]
    
    def getDihedralAngles(self,  atom_names, start_shift=0):
        dih_angles = {}
        _residues = self.residues
        if atom_names[0][0] == -1:
            start_resid = self.pdb_first_resid + 1
        else:
            start_resid = self.pdb_first_resid
        start_resid += start_shift

        if atom_names[-1][0] == 1:
            end_resid = self.pdb_last_resid
        else:
            end_resid = self.pdb_last_resid + 1 
        for resid in range(start_resid, end_resid):
        #for resid in range(2, 5):
            residue = self.residues[resid]
            atoms_xyz = []
            for resid_shift, atom_name in atom_names: 
                residue = _residues[resid+resid_shift]
                atoms_xyz.append(self.getAtomXYZ(residue, atom_name))
                #print resid, residue.name,  atom_name, atoms_xyz[-1]
            dih_angles[resid] = dihedral_vect(atoms_xyz[0],atoms_xyz[1], atoms_xyz[2],atoms_xyz[3])
        return dih_angles


    def segRMSD(self, seg1, seg2):
        """ compare rmsd between segments with same squence"""
        seg1_xyz, seg1_seq = self.getSegXYZ(seg1)
        seg2_xyz, seg2_seq = self.getSegXYZ(seg2)
        rmsd = self.calculate_rmsd(seg1_xyz, seg2_xyz)
        print "# rmsd= %7.3f seg1= %s : %s seg2= %s: %s " % (rmsd, seg1, seg1_seq,  seg2, seg2_seq)


    def isConnect(self, resid1, resid2, vb=0):
        """ check if the two nucleotides are connected:
	    0 : no connect
	    1 : 5' ---> 3' connection
	    2 : 3' ---> 5' connection 
	"""
        res1 = self.residues[resid1]
        res2 = self.residues[resid2]
        o3_p1 = res1.getAtomVect("O3'")
        o5_p1 = res1.getAtomVect("O5'")
        o3_p2 = res2.getAtomVect("O3'")
        o5_p2 = res2.getAtomVect("O5'")
        o35 = self.getAtomDist(o3_p1, o5_p2)
        o53 = self.getAtomDist(o5_p1, o3_p2)
        #print "#  resid %3d -- resid %3d : o35= %7.3f  o53= %7.3f " % (resid1, resid2, o35, o53)
        #self.getDistance
        
        if o35 < 2.7:
            out_put = 1
        elif o53 < 2.7:
            out_put =  2 
        else:
            out_put = 0
        if vb: print "# resid %3d -- %3d: O3/5 -- O3/5: %d " % (resid1, resid2, out_put)
        return out_put


    def setResBase(self):
        _resids = []
        for resid in self.resids:
            #try:
            self.residues[resid].setBase()
            _resids.append(resid)
            #except: 
                ## remove residue with missing atoms in base 
                #print "# remove residue %s with missing atoms in base " % resid 
                #del self.residues[resid]  
        self.resids = _resids 
        
        
    def setSeqDict(self, vb=0):
        """ set the sequence dictionary  resid -- residue name """
        seq_dict = {}
        backid_dt = {}
        ressid_dt = {}
        for resid in self.resids:
            seq_dict[resid] = self.residues[resid].nt_name 
            backid_dt[self.residues[resid].resid] = resid
            ressid_dt[resid] = self.residues[resid].resid
        self.seq_dict = seq_dict
        self.backid_dt = backid_dt
        self.ressid_dt = ressid_dt 
        if vb: 
            print "# backid_dt:\n", backid_dt
            print "# seq_dict: ", self.seq_dict

    def getSequence(self, vb=0):
        if not hasattr(self, 'seqence'):
            self.sequence = ""
            for resid in self.resids:
                self.sequence += self.residues[resid].nt_name 
        if vb: 
            print "# sequence: %s" % self.sequence
            print "# num_nt: %3d " % len(self.sequence)
        return self.sequence
        
        
    def getUniqAtoms(self, atom_names):
        """ remove duplicate atom name in the atom_names list """
        uniq_atom = []
        for _atom_name in atom_names:
            if _atom_name not in uniq_atom:
                uniq_atom.append(_atom_name)
        return uniq_atom 


    def getOrientationSS(self, res1, res2):
        n1, n2 = self.__orientationSS(res1, res2) 
        bp_type = "S/S"
        if n1>0 and n2>0:             #   the larger suger side  
            bp_type = "s/S"
        else:
            n3, n4 = self.__orientationSS(res2, res1) 
            if n3>0 and  n4>0: bp_type = "S/s" 
        return bp_type    


    def __orientationSS(self, res1, res2):
        n1 = 0
        n2 = 0
        res1_atom_xyz = res1.getAtomVect("O2'")
        for  res2_atom in res2.atom_names:
            if res2_atom[0] not in ["N", "O"]: continue 
            res2_atom_xyz = res2.getAtomVect(res2_atom)
            dd = self.getVecLen( res1_atom_xyz - res2_atom_xyz )
            if res2_atom != "O2'":
                if dd < 3.8: n1 += 1 
            if res2_atom[0] in ["N", "O"] and res2_atom[-1] != "'":
                if dd < 3.8: n2 += 1
        return n1, n2     

    def isBaseContact(self, res1, res2):
        """ check if two bases are contacted (at least one pair of base O & N & C is within range) """
        no_contact = 1
        base1_xyz = res1.ring_xyz
        base2_xyz = res2.ring_xyz
        for i in range(res1.num_atom_in_base):
            for j in range(res2.num_atom_in_base):
                dxyz = self.getVecLen(base1_xyz[i,:] - base2_xyz[j,:])
                #print "# dxzy= %8.3f " % dxyz 
                if dxyz <= self.MAX_HB: # H-bond upper limit 
                    #print "# dxzy= %8.3f " % dxyz
                    no_contact = 0
                    break 
        return no_contact        


    def isBaseStack(self, res1,res2, dv):
        """get the distance between the ring center, get rid of the stacks """
        cc_dist = self.getVecLen(res1.getRingCenter() - res2.getRingCenter())
        if res1.nt_name[0] in ['C','U'] and res2.nt_name in ['C','U']:        
            if dv > 2.3:
                if cc_dist<5.6: return 1
            elif dv>= 1.9:
                if cc_dist<5.5: return 1
            elif dv <1.9:
                if  cc_dist<5.0: return 1
        else:
            if dv >2.3:
                if cc_dist<5.8: return 1
            elif  dv >= 1.9:
                if cc_dist<5.7: return 1
            elif dv <1.9:
                if cc_dist<5.4: return 1
        return 0 



    def getAng_Norm(self, dval):
        return 90.0 - fabs(self.dot2deg(dval)-90.0)


    def dot2deg(self, dval):
        if dval >= 1.0:
            ang_deg = 0.0
        elif dval <= -1.0:
            ang_deg = 180.0
        else:
            ang_deg = degrees(acos(dval))
        return ang_deg




    
    def __addNewResidue(self, residue):
        if residue is None: return 
        self.residues[residue.res_sid] = residue 
        self.res_sid += 1
        

    def simpleParsePDB(self):
        """ parse pdb datafile to get residues, sequence, cartesian coordinates, No duplicate resids 
        """
        self.residues = {}
        self.pdb_atom_lns = []
        residue  = None
        pre_resid = -99 
        num_ter = 0
        serial_resid = -1 # handle different chains have same residue number 
        for ln in self.pdb_lns:
            if ln[0:4] == "ATOM" or ln[0:6] == "HETATM":
                atom = self.parseAtom(ln)
                resid = atom.resid
                if pre_resid != resid: #new residue 
                    if pre_resid> 0: self.residues[pre_resid] = residue 
                    residue = Nucleotide(name=atom.resname, resid=resid, chain_id=atom.chain_id, res_sid=resid)
                    pre_resid = resid
                residue.AddAtom(atom)
            elif ln[0:3] == "TER" and num_ter:
                num_ter += 1
                break 
            elif ln[0:6] == 'ENDMDL': break 
        self.residues[resid] = residue 
        
        self.resids = self.residues.keys()
        self.resids.sort()
        self.num_res = len(self.residues)
        self.pdb_first_resid = self.resids[0]
        

    def getRNASeq(self):
        rna_seq = {}
        for resid in self.resids:
            rna_seq[resid] = self.residues[resid].name[0]
        return rna_seq


    
    def __setResids(self):
        self.resids = self.residues.keys()
        self.resids.sort()
        self.num_res = len(self.residues)
        pdb_resids = []
        for resid in self.resids:
            pdb_resids.append(self.residues[resid].resid)
        self.pdb_resids = pdb_resids 
        
    
    def parseAtom(self, ln):
        self.pdb_atom_lns.append(ln)
        ln = lstrip(ln)
        atomid = int(ln[6:11])
        atomname= strip(ln[12:16])
        resname = strip(ln[17:20])
        resid = int(ln[22:26])
        x = float(ln[30:38])
        y = float(ln[38:46])
        z = float(ln[46:54])
        chain_id = ln[21]
        #print "Chain: ", chain_id
        atom = Atom(atomname=atomname, atomid=atomid, xyz=(x,y,z), resid=resid,  resname=resname, chain_id=chain_id)
        return atom 

    
    def parsePDB(self):
        self.parsePDBw1M()
        #self.TparsePDB()
        
    def parsePDBw1M(self):
        """ parse pdb datafile to get residues, sequence, cartesian coordinates
            only parse one model if pdb file includes two or more models
        """
        self.residues = {}
        self.pdb_atom_lns = []
        residue  = None
        pre_resid = -99 
        pre_chain_id = None 
        chain_ids = []
        for ln in self.pdb_lns:
            if ln[0:4] == "ATOM" or ln[0:6] == "HETATM":
                atom = self.parseAtom(ln)
                resid = atom.resid
                if self.res_sid is None: self.res_sid = resid 
                
                if pre_resid != resid: #new residue 
                    self.__addNewResidue(residue)
                    residue = Nucleotide(name=atom.resname, resid=resid, chain_id=atom.chain_id, res_sid=self.res_sid)
                    if not residue.isRNANT(): break 
                    pre_resid = resid
                if pre_chain_id != atom.chain_id:
                    pre_chain_id = atom.chain_id
                    if pre_chain_id in chain_ids:
                        break
                    else:
                        chain_ids.append(pre_chain_id)
                residue.AddAtom(atom)
            #elif ln[0:3] == "TER":
            #    break 
            #elif ln[0:6] == 'ENDMDL': break 
        #print "# chain_ids: ", chain_ids 
        if residue.isRNANT():
            self.__addNewResidue(residue)
        self.__setResids()
        
        
    def TparsePDB(self):
        """ parse pdb datafile to get residues, sequence, cartesian coordinates
            only parse 2 models if pdb file includes two or more models
        """
        pre_resid = -99 
        num_models = 2
        model_id = 0       # 0: means no model in pdb file
        residue = None
        for ln in self.pdb_lns:
            ln = ln.lstrip()
            if not ln: continue 
            if ln[0:5]=="MODEL":
                tmp, model_id = split("\s+", ln.rstrip())
                model_id = int(model_id)
                num_models += 1
            if  ln[0:4] == "ATOM": # in ln[0:6] == or "HETATM" in ln:
                note = ln 
                atom = self.parseAtom(ln)
                resid = atom.resid
                if  pre_resid != resid:  
                    if self.res_sid > 0: self.__addNewResidue(residue)
                    residue = Nucleotide(name=atom.resname, resid=resid, chain_id=atom.chain_id, note=note, model_id=model_id, res_sid=self.res_sid)
                    #print "# pre_resid= %3d resid= %3d model_id= %3d " % (pre_resid, resid, model_id)
                    pre_resid = resid 
                residue.AddAtom(atom)
        self.__addNewResidue(residue)
        self.__setResids()
    
    
    def prnResid(self):
        """ print residue id in pdb file(*.resid) and resid in memory( *.iresid) """
        print "# number of residues = %3d " % self.num_res
        for i in range(self.num_res):
            resid = self.resids[i]
            print "# No. %3d : resid(in resids)= %3d  resid(pdb)= %3d  resid(mem)= %3d " % (i, resid, self.residues[resid].resid, self.residues[resid].res_sid)
            
            
       
    def getDihedralAngles(self,  atom_names, start_shift=0):

        dih_angles = {}
        _residues = self.residues
        if atom_names[0][0] == -1:
            start_resid = self.pdb_first_resid + 1
        else:
            start_resid = self.pdb_first_resid
        start_resid += start_shift

        if atom_names[-1][0] == 1:
            end_resid = self.pdb_last_resid
        else:
            end_resid = self.pdb_last_resid + 1 

        for resid in range(start_resid, end_resid):
        #for resid in range(2, 5):
            residue = self.residues[resid]
            atoms_xyz = []
            for resid_shift, atom_name in atom_names: 
                residue = _residues[resid+resid_shift]
                atoms_xyz.append(self.getAtomXYZ(residue, atom_name))
                #print resid, residue.name,  atom_name, atoms_xyz[-1]
            dih_angles[resid] = dihedral_vect(atoms_xyz[0],atoms_xyz[1], atoms_xyz[2],atoms_xyz[3])
        return dih_angles


    def prnDihAngles(self, angles, cmt=""):
        if cmt: print cmt 
        resids = angles.keys()
        resids.sort()
        #for resid in resids:
        #    print "%4d %4s  %6.1f " % (resid, self.residues[resid].name,  degrees(angles[resid]))


    def prnAllAngle(self):
        print "# Alpha Beta Gamma delta  Epsilon  Zeta Xhi  Mu1  Mu2 "
        angle_dts = [self.alphas, self.betas, self.gammas, self.deltas, self.epsilons, self.zetas, self.xhis, self.mu1, self.mu2]
        for resid in self.resids:
            ln = " %4d " % resid
            for ang_dt in angle_dts:
                angle = ang_dt.get(resid, None)
                if angle:
                    ln += "%6.1f " % degrees(angle)
                else:
                    ln += "%6s " % " " 
            print ln 

    def setAllDihAngles(self):
        self.prime_format = 0
        self.getAlpha()
        self.getBeta()
        self.getGamma()
        self.getDelta()
        self.getEpsilon()
        self.getZeta()
        self.getXhi()
        self.getMu1()
        self.getMu2()


    def getDistance(self, atom1, atom2, vb=0):
        r1 = getattr(self.residues[atom1[0]], atom1[1])
        r2 = getattr(self.residues[atom2[0]], atom2[1])
        dr_len = self.getAtomDist(r1, r2)
        if vb: print "# distance: %3d %2s -- %3d %2s = %8.4fA" % (atom1[0], atom1[1], atom2[0], atom2[1], dr_len)
        return dr_len



    def start(self):
        #self.prnAllAngle()
        #atom1 = (39, 'C2')  # resid, atom name 
        #atom2 = (39, 'O2')
        #self.getDistance(atom1, atom2)
        #print self.resids
        print "lib"
        pass 




    def getEpsilon(self, resid1, resid2):
        """ i C4' - C3' - O3' - (i+1)P 
        """
        _res1 = self.residues[resid1]
        _res2 = self.residues[resid2]
        #self.prnDihAngles(angles, cmt="epsilon:  C4' - C3' - O3' - P ")
        epsilon =  dihedral_vect(_res1.getAtomVect("C4'"), _res1.getAtomVect("C3'"), _res1.getAtomVect("O3'"), _res2.getAtomVect("P"))
        return degrees(epsilon) 

    def getResEpsilon(self, resid):
        return self.getEpsilon(resid, resid+1)


    def getZeta(self, resid1, resid2, vb=0):
        """ iC3' - O3' - P  
        self.prnDihAngles(angles, cmt="zeta: C3' - O3' - (i+1) P - (i+1) O5' ")
        """
        #atom_names = [(0,"C3'"), (0, "O3'"), (1,"P"), (1,"O5'")]
        _res1 = self.residues[resid1]
        _res2 = self.residues[resid2]
        zeta =  dihedral_vect(_res1.getAtomVect("C3'"), _res1.getAtomVect("O3'"), _res2.getAtomVect("P"), _res2.getAtomVect("O5'"))

        if vb: 
            print "# resid= %3d  O5': %s "  % (resid2, _res2.getAtomVect("O5'") )
            print "# resid1= %3d  resid2= %3d " % (resid1, resid2)
        return degrees(zeta)


    def getResZeta(self, resid, **rest):
        return self.getZeta(resid, resid+1, **rest)

    def getResAlpha(self, resid):
        return self.getAlpha(resid-1, resid)

    def getAlpha(self, resid1, resid2):
        """ i-1 O3' --  iP -- iO5' --C5'
        """
        # first residue has not O3' 
        _res1 = self.residues[resid1]
        _res2 = self.residues[resid2]
        alpha =  dihedral_vect(_res1.getAtomVect("O3'"), _res2.getAtomVect("P"), _res2.getAtomVect("O5'"), _res2.getAtomVect("C5'"))
        return degrees(alpha)


    def writeXpRNAPDB(self, fn_pdb, lns=None, vb=0):
        fmt_pdb = '%-6s%5d %-4s %-4s%1s%4d%1s%3s%8.3f%8.3f%8.3f%18s%-4s'
        cnt = 0
        nt_full_names = {'G':'GUA', 'U': 'URI',  'A':'ADE', 'C':"CYT"}
        #lns = []
        if lns is None: lns = []
        #print " # lns : ", lns 
        for i in self.resids:
            _amino = self.residues[i]
            res_name = nt_full_names[_amino.nt_name]
            for _atom in _amino.atom_list:
                atom_name = _atom.name 
                if atom_name == 'OP1': atom_name = "O1P"
                elif atom_name == 'OP2': atom_name = "O2P"
                if len(atom_name) <= 3: atom_name= " %s" % atom_name
                if _atom.misc is None: 
                    misc = " "
                else:
                    misc  = _atom.misc
                if _atom.segid is None:
                    segid = ""
                else:
                    segid = _atom.segid 
                ln =  fmt_pdb  % ('ATOM', _atom.atomid, atom_name, res_name," ", _atom.resid, "","", _atom.x, _atom.y, _atom.z,  misc, segid)
                lns.append(ln)
        lns = [ "# pdb was saved in Xplor format: %s " %  fn_pdb]
        self.add2log(lns, vb=vb)
        lns = self.pdb_head_lns + lns + ["END"]
        self.saveLines(lns, fn_pdb, vb=vb)

    def writePDB_wo_H(self, fn_pdb, lns=None, vb=1):
        """ remove Protons in PDB file """
        fmt_pdb = '%-6s%5d %-4s %-4s%1s%4d%1s%3s%8.3f%8.3f%8.3f%18s%-4s'
        cnt = 0
        nt_full_names = {'G':'GUA', 'U': 'URI',  'A':'ADE', 'C':"CYT"}
        #lns = []
        if lns is None: lns = []
        #print " # lns : ", lns 
        for i in self.resids:
            _amino = self.residues[i]
            res_name = nt_full_names[_amino.nt_name]
            for _atom in _amino.atom_list:
                atom_name = _atom.name 
                if atom_name[0] == 'H': continue 
                #if atom_name == 'OP1': atom_name = "O1P"
                #elif atom_name == 'OP2': atom_name = "O2P"
                if len(atom_name) <= 3: atom_name= " %s" % atom_name
                if _atom.misc is None: 
                    misc = " "
                else:
                    misc  = _atom.misc
                if _atom.segid is None:
                    segid = ""
                else:
                    segid = _atom.segid 
                ln =  fmt_pdb  % ('ATOM', _atom.atomid, atom_name, res_name," ", _atom.resid, "","", _atom.x, _atom.y, _atom.z,  misc, segid)
                lns.append(ln)
        if vb: print "# pdb was saved in < %s >: " %  fn_pdb
        lns = self.pdb_head_lns + lns + ["END"]
        self.saveLines(lns, fn_pdb, vb=vb)


    def getAtom_by_XYZ(self, xyz, nt_name, atom):
        res_atoms = nt_atoms_woh[nt_name]
        atom_id = res_atoms.index(atom)
        return xyz[atom_id, :]


    def resXYZ2Imino(self, xyz, nt_name, vb=0):
        if nt_name in ['A','C']: 
            if vb: print "# residue %s  has No imino proton" % nt_name
            return None
        
        C2 = self.getAtom_by_XYZ(xyz, nt_name, 'C2' )
        if nt_name == 'G':
            Natom = 'N1'
            vatom = "C6"
            Hatom = 'H1'
        elif nt_name == 'U':
            Natom = 'N3'
            vatom = "C4"
            Hatom = 'H3'
        N = self.getAtom_by_XYZ(xyz, nt_name, Natom) #Vect(Natom)
        C4_6 = self.getAtom_by_XYZ(xyz,nt_name,  vatom)
        u = self.normVect(C2 - N )
        v = self.normVect(C4_6 - N )
        norm = cross(u, v)
        NH_n = - self.normVect(u+v)
        if nt_name == 'G': 
            NH_n = self.new_arb_rot(NH_n, norm, -1.4)
        NH = NH_n.copy() * 1.010 
        #H_atom = NH + N
        return NH 

    def res2pdb(self, resids, atom_id=0, vb=1):
        """
        generate pdb from nucleotide's atom coordinates xyz(dict), res_seq(dict)
	"""
        lns = []
        atom_id = atom_id
        for resid in resids:
            pdb_lns, atom_id = self.residues[resid].toPdb(atom_id=atom_id)
            lns += pdb_lns
        return lns


    def resXYZ_Seq2pdb(self, xyz, res_seq, atom_id=0, vb=0):
        """
        generate pdb from nucleotide's atom coordinates xyz(dict), res_seq(dict)
	"""
        fmt_pdb = '%-6s%5d %-4s %-4s%1s%4d%1s%3s%8.3f%8.3f%8.3f%18s%-4s'
        lns = []
        atom_id = atom_id
        resids = xyz.keys()
        resids.sort()
        for resid in resids:
            nt_name = res_seq[resid] 
            if vb: print "# resid= %3d nt_name= %s " % (resid, nt_name)
            #nt_atoms_woh[nt_name]
            nt_atom_names = nt_atoms_woh[nt_name]
            res_xyz = xyz[resid]
            for i in range(len(nt_atom_names)):
                atom_name = nt_atom_names[i]
                atom_xyz = res_xyz[i]
                if len(atom_name) <= 3: atom_name= " %s" % atom_name
                atom_id += 1
                ln =  fmt_pdb  % ('ATOM', atom_id, atom_name, nt_name," ",  resid, "","", atom_xyz[0],atom_xyz[1], atom_xyz[2], "", "") #_atom.x, _atom.y, _atom.z, _atom.misc, _atom.segid)
                lns.append(ln)
                #print ln 
        return lns, atom_id  




    def getResBBLen(self, resid, vb=0):
        """ return residue backbone bond lengths """
        bb_lens = self.residues[resid].getBB_BondLengths()
        if resid+1 in self.resids:
            next_P = self.residues[resid+1].getAtomVect("P")
            O3 = self.residues[resid].getAtomVect("O3'")
            O3_P_len = self.getVecLen(next_P - O3)
        else:
            O3_P_len = -1.0
        bb_lens.append(O3_P_len)
        if vb: self.prnResBBLen(bb_lens)
        return bb_lens 

    def getResBBAngles(self, resid, vb=0):
        """ return residue backbone bond lengths """
        bb_angles = self.residues[resid].getBB_BondAngles()
        if resid+1 in self.resids:
            next_P = self.residues[resid+1].getAtomVect("P")
            next_O5 = self.residues[resid+1].getAtomVect("O5'")
            C3 = self.residues[resid].getAtomVect("C3'")
            O3 = self.residues[resid].getAtomVect("O3'")
            O3_C3 = C3 - O3 
            O3_P = next_P - O3 
            P_O5 = next_O5 - next_P
            C3_O3_P = self.getVecAngle(O3_C3, O3_P)
            O3_P_O5 = self.getVecAngle(-O3_P, P_O5)
        else:
            C3_O3_P = 0
            O3_P_O5 = 0
        bb_angles.append(C3_O3_P)
        bb_angles.append(O3_P_O5)
        if vb: self.prnResBBAngles(bb_angles)
        return bb_angles

    def prnResBBAngles(bb_angles):
        ang_strs = ["P_O5'_C5'", "O5'_C5'_C4'", "C5'-C4'-C3'", "C4'-C3'-O3'", "C3'-O3'-P", "O3'-P-O5'"]
        for i in range(6):
            ang_deg = degrees(bb_angles[i])
            print "# %s: ang= %6.1f "  % (ang_strs[i], ang_deg)



    def prnResBBLen(self, bb_lens):
        assert len(bb_lens) == 6 
        bb_atoms = ["P", "O5'", "C5'", "C4'", "C3'", "O3'","P"]
        for i in range(6):
            print "# %3s -- %3s= %8.3f " % (bb_atoms[i], bb_atoms[i+1], bb_lens[i])

    
    def chkMissingAtoms(self):
        for resid in self.resids:
            #print "# check missing atoms: res_sid= %3d " % resid
            #print "# resid= %3d " % self.residues[resid].resid
            self.residues[resid].chkMissAtom()
            
    def rmResWithMissingAtom(self, vb=0):
        rm_resids = []
        for resid in self.resids:
            rm_res = 0
            if self.residues[resid].miss_atom:
                #print "# resid= %3d  miss atom: %s " % (self.residues[resid].resid, self.residues[resid].miss_atoms)
                #for _atom in self.residues[resid].miss_atoms:
                    #if "N" in _atom:
                        #rm_res = 1
                        #break 
                    #elif "O" in _atom and "P" not in _atom:
                        #rm_res = 1
                        #break 
            #if rm_res:
                rm_resids.append(resid)
                
        for resid in rm_resids:
            if vb:
                print "# delete resid= %3d  %s \n" % (resid, self.residues[resid].nt_name)
                print "# missing atom: %s " % self.residues[resid].miss_atoms
            del self.residues[resid]
            #self.resids.remove(resid)
        resids = self.residues.keys()
        resids.sort()
        self.resids = resids
                
    def resDict2pdb(self, resid_cutoff=498, rep_resids=None, resid_shift=900):
        """ save residue xyz array into pdb file """
        fmt_pdb = '%-6s%5d %-4s %-4s%1s%4d%1s%3s%8.3f%8.3f%8.3f%18s%-4s'
        resids = self.new_pdb_xyz.keys()
        resids.sort()
        lns = []
        atom_id = 0
        for resid in resids:
            if resid >= resid_cutoff: continue
            if rep_resids and resid in rep_resids:
                nt_name = self.seq_dict[resid+resid_shift] 
                print "# replace resid= %3d nt_name= %s " % (resid,nt_name)
            else:
                #nt_name = self.seq_dict[resid] 
                nt_name = self.getRes_by_resid(resid).nt_name 
                #print "# resid= %3d nt_name= %s " % (resid,nt_name)
            nt_atom_names = nt_atoms_woh[nt_name]
            res_xyz = Npy.array(self.new_pdb_xyz[resid])
            #print res_xyz
            #res_xyz = array(res_xyz)
            for i in range(len(nt_atom_names)):
                atom_name = nt_atom_names[i]
                #print "# i= %3d nt_name=%s "  % (i, nt_name)
                if len(atom_name) <= 3: atom_name= " %s" % atom_name
                atom_xyz = res_xyz[i,:]
                atom_id += 1
                ln =  fmt_pdb  % ('ATOM', atom_id, atom_name, nt_name," ",  resid, "","", atom_xyz[0],atom_xyz[1], atom_xyz[2], "", "") #_atom.x, _atom.y, _atom.z, _atom.misc, _atom.segid)
                lns.append(ln)
                #print ln 
        return lns 
    
    def resXYZ2pdb(self,pdb_xyz=None):
        """ save residue xyz array into pdb file """
        fmt_pdb = '%-6s%5d %-4s %-4s%1s%4d%1s%3s%8.3f%8.3f%8.3f%18s%-4s'
        if pdb_xyz is None: pdb_xyz = self.new_pdb_xyz 
        resids = pdb_xyz.keys()
        resids.sort()
        lns = []
        atom_id = 0
        for resid in resids:
            nt_name = self.seq_dict[resid]
            nt_atom_names = nt_atoms_woh[nt_name]
            res_xyz = Npy.array(pdb_xyz[resid])
            for i in range(len(nt_atom_names)):
                atom_name = nt_atom_names[i]
                #print "# i= %3d nt_name=%s "  % (i, nt_name)
                if len(atom_name) <= 3: atom_name= " %s" % atom_name
                atom_xyz = res_xyz[i,:]
                atom_id += 1
                ln =  fmt_pdb  % ('ATOM', atom_id, atom_name, nt_name," ",  resid, "","", atom_xyz[0],atom_xyz[1], atom_xyz[2], "", "") #_atom.x, _atom.y, _atom.z, _atom.misc, _atom.segid)
                lns.append(ln)
        return lns 
    

    def setBB_AtomIDs(self, vb=0):
        """ return backbone atom ids in residue atom coordinates' array """
        bb_aids = {}
        for nt in ['A','U','G','C']:
            aids = []
            _atoms_woH = nt_atoms_woh[nt]
            for _atom in nt_bb_atoms:
                _aid = _atoms_woH.index(_atom)
                aids.append(_aid)
                if vb: print "# nt=%s atom= %4s  aid= %3d " % (nt, _atom, _aid)
            bb_aids[nt] = aids
        self.bb_aids = bb_aids
    
    
    def getResBB(self, res_xyz, resn, _list=1):        
        aids = self.bb_aids[resn]
        res_bb_xyz = []
        for _aid in aids:
            #print "# _aid= %3d" % _aid 
            #print res_xyz[_aid]
            res_bb_xyz.append(res_xyz[_aid].copy())
        return res_bb_xyz
    
    
    def saveDupBB2PDB(self, dup_resids, dup_xyz, seq_dt, atom_id=1, vb=0):
        """ save duplex's backbone atoms' coordinates into pdb file """
        fmt_pdb = '%-6s%5d %-4s %-4s%1s%4d%1s%3s%8.3f%8.3f%8.3f%18s%-4s'
        lns = []
        atom_id = atom_id
        num_bb_atoms = len(nt_bb_atoms)
        atom_sid = 0
        for i in range(len(dup_resids)):
            resid = dup_resids[i]
            nt_name = seq_dt[resid]
            if vb: print "# resid= %3d nt_name= %s " % (resid, nt_name)
            for j in range(num_bb_atoms):
                atom_name = nt_bb_atoms[j]
                atom_xyz = dup_xyz[atom_sid]
                atom_sid += 1
                if len(atom_name) <= 3: atom_name= " %s" % atom_name
                atom_id += 1
                ln =  fmt_pdb  % ('ATOM', atom_id, atom_name, nt_name," ",  resid, "","", atom_xyz[0],atom_xyz[1], atom_xyz[2], "", "") #_atom.x, _atom.y, _atom.z, _atom.misc, _atom.segid)
                lns.append(ln)
        return lns

    
                

