# Action incudes  
#  the movement in homogeneous coordinates 
#  Rotation  of moleucles

import os
import numpy as Npy
from numpy import matrix, mat, dot, transpose, sqrt, shape, array, argsort
from numpy.linalg import det, svd  
_array = Npy.array
from math import  log, sin, cos, sqrt, pi, fabs, radians,degrees, acos, asin, atan2, pi, exp, ceil 
from libmath import Vector
from libntatoms import * 

class Rotation(Vector):
    """ using Eular angle definiations from Mathword:
             http://mathworld.wolfram.com/EulerAngles.html
             
    """
    ZAXIS = Npy.array([0.0, .0, 1.0])
    YAXIS = Npy.array([0.0, 1.0, 0.0])
    XAXIS = Npy.array([1.0, 0.0, 0.0])
    
    def getRz(self, phi):
        """ rotate around z """
        B =matrix([[ cos(phi), sin(phi), 0.], 
                    [-sin(phi), cos(phi), 0.],
                    [    0.,        0.,   1.]])
        return B 
    
    def getRMD(self, phi):
        """ rotate around Z """
        return self.getRz(phi)
    
    def getRMC(self, theta):
        """ rotate around X """ 
        return self.getRx(theta)
    
    def getRx(self, theta):
        """ rotate around x """
        B = matrix([[    1.,         0. ,    0.     ],
                    [    0.,  cos(theta), sin(theta)],
                    [    0., -sin(theta), cos(theta)]])
        return B 
    
    def getRy(self, phi):
        B = matrix([[ cos(phi),   0.,  -sin(phi)],
                    [    0.   ,   1.,      0.   ],
                    [ sin(phi),   0.,   cos(phi)]])
        return B 
    
    def getRMB(self, psi):
        """ rotate around z """
        return self.getRz(psi)
        #B = matrix([[ cos(psi), sin(psi), 0.], 
                    #[-sin(psi), cos(psi), 0.],
                    #[    0.,        0.,   1.]])
        #return B 
  
    def getBackRotM(self, phi, theta, psi):
        """ back Rotation Matrix:  AB = D C B:  z-x-z """
        B = self.getRMB(-psi)
        C = self.getRMC(-theta)
        D = self.getRMD(-phi)
        AB = C*B 
        AB = D * AB 
        return AB 
    
    def getRotM_By_Product(self, phi, theta, psi):
        """ A = B C D """
        B = self.getRMB(psi)
        C = self.getRMC(theta)
        D = self.getRMD(phi)
        AB = C * D 
        AB = B * AB 
        return AB 
        
    def getRotM(self, phi, theta, psi):
        """ A = B C D """
        cos_phi = cos(phi)
        sin_phi = sin(phi)
        cos_psi = cos(psi)
        sin_psi = sin(psi)
        cos_the = cos(theta)
        sin_the = sin(theta)
        RM = matrix([[ cos_psi*cos_phi - cos_the*sin_phi*sin_psi, 
                       cos_psi*sin_phi + cos_the*cos_phi*sin_psi,
                       sin_psi*sin_the],
                     [-sin_psi*cos_phi - cos_the*sin_phi*cos_psi,
                      -sin_psi*sin_phi + cos_the*cos_phi*cos_psi,
                       cos_psi*sin_the],
                     [ sin_the*sin_phi,
                      -sin_the*cos_phi,
                       cos_the]])
        return RM 
    
    def normVect(self, vect):
        """ return a normalized vector)"""
        return vect/sqrt(vect[0]*vect[0] + vect[1]*vect[1] + vect[2]*vect[2])    
    
    def rotVec(self, rot_M, vec):
        vec_m = matrix(vec)
        rot_vec = rot_M * vec_m.T 
        return _array(rot_vec.T)[0]
        
    
    
    def new_arb_rot(self, vect, axis, gamma):
        """ rotate vector around axis by gamma degree
        gamma in units of degree 
        """
        gamma = radians(gamma)
        cos_gamma = cos(gamma)
        sin_gamma = sin(gamma)
        #print "# gamma(rad): ", gamma
        vec_len = self.getVectNorm(vect)
        if vec_len < 0.000001: return vect 
        vect = vect/vec_len 
        axis = self.normVect(axis)
        x,y,z   = vect[0], vect[1], vect[2]
        u, v, w  = axis[0], axis[1], axis[2]
        v2 = v * v 
        w2 = w * w 
        d2 = v2 + w2  
        
        if d2 == 0.0:
            rot_vector = _array([x, y*cos_gamma - z * sin_gamma, z * cos_gamma + y * sin_gamma])
        else:
            u2 = u * u 
            ux = u * x 
            vy = v * y
            wz = w * z 
            pp = ux + vy + wz 
            #an = u2 + v2 + w2 
            #print "# an : ", an 
            nx = u *pp + ( x * d2 - u *( vy + wz )) * cos_gamma + \
                 (- w*y + v*z ) * sin_gamma
            ny = v * pp  + ( y* (u2 + w2)  - v*(ux + wz))*cos_gamma + \
                 (w*x - u*z) * sin_gamma 
            nz = w * pp  + ( z*(u2 + v2) - w * (ux + vy)) * cos_gamma + \
               (-v*x + u*y)*sin_gamma 
            rot_vector = _array([nx,ny,nz])
        return  self.normVect(rot_vector)*vec_len

    
class HomoMove:
    def getHM(self, conformation1, conformation2, is_center=True, vb = False):     
        """ return moving matrix:  best fit conformation2 to conformation1, move cf2 ---> cf1 """
        vecs1 = conformation1.copy()      
        vecs2 = conformation2.copy()      
        if vb: 
            print "# conformation1.shape: ", vecs1.shape 
            print "# conformation2.shape: ", vecs2.shape
        n_vec, vec_size = shape(vecs1)      
        center1 = sum(vecs1,0) / float(n_vec)      
        center2 = sum(vecs2,0) / float(n_vec)      
        if is_center:      
            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      
        # reflect along smallest principal axis      
        if is_reflection: 
            S[-1] = -S[-1]      
            V[-1,:] = V[-1,:] * (-1.0)      
        optimal_rotation = dot(V, W_trans)   
        #conformation2[:,:] = dot(vecs2, optimal_rotation)  # + center1 
        #self.testrmsd(vecs1, conformation2)
        #vb=1
        #print "# moved: ", conformation2
        if vb: 
            rmsd_sq = (E0 - 2. * sum(S)) / float(n_vec)      
            rmsd_sq = max([rmsd_sq, 0.])      
            rmsd = sqrt(rmsd_sq)
            print "#  rmsd= %8.4f " %  rmsd 
        return optimal_rotation, center1, center2 
    
    
    def tgetHM(self, conformation1, conformation2, is_center=True, vb = False):     
        """ return moving matrix:  best fit conformation2 to conformation1, move cf2 ---> cf1 """
        vecs1 = conformation1.copy()      
        vecs2 = conformation2.copy()      
        if vb: 
            print "# conformation1.shape: ", vecs1.shape 
            print "# conformation2.shape: ", vecs2.shape
        n_vec, vec_size = shape(vecs1)      
        center1 = sum(vecs1,0) / float(n_vec)      
        center2 = sum(vecs2,0) / float(n_vec)      
        if is_center:      
            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      
        # reflect along smallest principal axis      
        if is_reflection: 
            self.add2log(["# reflection" ])
            #S[-1] = -S[-1]      
            #V[-1,:] = V[-1,:] * (-1.0)      
        optimal_rotation = dot(V, W_trans)   
        #print "# origcf2: \n",
        #print conformation2
        #print "# ct1: ", center1
        #print "# cf1: \n", vecs1[:4,:] #conformation1
        conformation2[:,:] = dot(vecs2, optimal_rotation)  # + center1 
        #print "# ct2: ", center2
        #print "# cf2: \n ", conformation2[:4,:]
        #vb = 1
        #self.testrmsd(vecs1, conformation2)
        #vb=1
        #print "# moved: ", conformation2
        if vb: 
            rmsd_sq = (E0 - 2. * sum(S)) / float(n_vec)      
            rmsd_sq = max([rmsd_sq, 0.])      
            rmsd = sqrt(rmsd_sq)
            print "#  rmsd= %8.4f " %  rmsd 
        return optimal_rotation, center1, center2 

            #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) + center1   
            #rmsd = sqrt(rmsd_sq)
            #if vb: print "#  rmsd= %8.4f " %  rmsd 
            #return rmsd 
        
    
    def testrmsd(self, cf1, cf2):
        sum = 0.0
        #print cf1.shape
        #print cf2.shape 
        num_vec = len(cf1)
        for i in range(num_vec):
            for j in range(3):
                dd = cf1[i,j] - cf2[i,j]
                sum += dd *dd 
        rmsd = sqrt(sum/num_vec)
        print "# test rmsd= %8.3f " % rmsd
    
    def homoMove(self, conformation, hm):
        """ homogeneous movement: move from center2 to center 1  """
        ROTM, center1, center2  = hm 
        next_cf = conformation.copy()
        next_cf -= center2
        next_cf_TR = dot(next_cf, ROTM) + center1 
        return next_cf_TR
    
    
    def moveBase(self, conf,  hm, old_anchor, new_anchor):
        """ move new base to old base position(mutation), new_anchor --> old_anchor """
        ROTM, center1, center2  = hm 
        next_cf = conf.copy() -  new_anchor
        next_cf_TR = dot(next_cf, ROTM) + old_anchor
        return next_cf_TR
        
    
    def hm_cf2_cf1(self, cf1, cf2):
        """ homomove cf2 to cf1 """
        hm = self.getHM(cf1, cf2)
        new_cf2 = self.homoMove(cf2, hm)
        return new_cf2 
 
 
    def homoMoveRes(self, res_xyz,  hm):
        ROTM, center1, center2 = hm 
        new_res_xyz = {}
        for _atom in res_xyz.keys():
            _atom_xyz = res_xyz[_atom]
            #   next_cf = conformation.copy()
            #print "\n# atom_xyz: \n", _atom_xyz
            _atom_xyz -= center2
            new_atom_xyz =   dot(_atom_xyz, ROTM) + center1 
            new_res_xyz[_atom] = new_atom_xyz
            #print "# atom %s  xyz: %s  mxyz: %s " % (_atom, _atom_xyz, new_atom_xyz)
        return new_res_xyz
    

    def getCF_RMSD(self, conformation1, conformation2, is_center=True, vb = False):     
        """ return moving matrix:  best fit conformation2 to conformation1, move cf2 ---> cf1 """
        vecs1 = conformation1.copy()      
        vecs2 = conformation2.copy()      
        if vb: 
            print vecs1.shape 
            print vecs2.shape
        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      
        else:
            center1 = array([0,0,0])
            center2 = center1.copy()

        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) + center1# center2    
        #print "# moved: ", conformation2
        rmsd = sqrt(rmsd_sq)
        if vb: print "#  rmsd= %8.4f " %  rmsd 
        return rmsd 
    


    
    def alignTwo(self, cf1, cf2, **rest):
        self.calculate_rmsd(cf1, cf2, is_align=True, **rest)
    
    def alignTwo_byRef(self, ref1, ref2, cf2):
        """# align two structure by reference frame """
        pass 
    
        
    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       
	"""      
        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) + center1   
            rmsd = sqrt(rmsd_sq)
            if vb: print "#  rmsd= %8.4f " %  rmsd 
            return rmsd 


class Sequence:
    def getMatchSeqScore(self, target_seq, tpt_seq):
        """ """
        assert len(target_seq) == len(tpt_seq)
        _score = 0
        for i in range(len(target_seq)):
            if target_seq[i] == tpt_seq[i]:
                match_score = 20
            elif target_seq[i] in ['A','G'] and tpt_seq[i] in ['A', 'G']:
                match_score = 1
            elif target_seq[i] in ['C','U'] and tpt_seq[i] in ['C', 'U']:
                match_score = 1
            else:
                match_score = 0
            if i < 2:
                _score += match_score * 10
            else:
                _score += match_score 
        return _score
    
    def bestBySeqScore(self, scores, max_tpt=3):
        scores = array(scores)
        sort_ids = argsort(scores)
        num_best= min(max_tpt,len(sort_ids))
        best_ids = sort_ids[-num_best:]
        best_ids = best_ids.tolist()
        best_ids.reverse()
        return best_ids
    
    def matchSequence(self, _tloop_dp, rna_loop):
        tp_len, tp_seq, tp_bp, tp_resids = rna_loop
        print "# target resids: ", tp_resids
        _tpt_seq = _tloop_dp[2]
        _tpt_resids = _tloop_dp[3]
        tloop_xyz = {}
        _tpt_xyz = _tloop_dp[4]
        for i in range(tp_len):
            _tpt_nt = _tpt_seq[i]
            _tpt_resid = _tpt_resids[i] 
            rna_nt = tp_seq[i]
            rna_resid = tp_resids[i]
            print "#  tpt: resid= %3d nt= %s    --> rna_nt: resid= %3d  nt= %s  " % (_tpt_resid, _tpt_nt, rna_resid, rna_nt)
            _tpt_res_xyz = _tpt_xyz[_tpt_resid]
            tloop_xyz[rna_resid] = self.matchBase(_tpt_nt,_tpt_res_xyz, rna_nt)
            #break 
        #print tloop_xyz
        _tloop_anchor_bp = _tpt_xyz[_tloop_dp[1][0]], _tpt_xyz[_tloop_dp[1][1]]
        return (_tloop_dp[0], _tloop_anchor_bp, tloop_xyz)
  
    def resids2seq(self, resids): # pass
        tp_seq = ""
        for resid in resids:
            tp_seq += self.seq_dict[resid]
        return tp_seq
    
    
class Mutation(HomoMove, Vector):
    ref_bases_dt = { "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)}}
        
    
    comp_seq = {"U":"C", "A":"G", "G":"A", "C":"U"}
    def prnResHeavyAtoms(self, res_nt, res_xyz):
        _nt_atoms = nt_atoms_woh[res_nt]
        for i in range(len(_nt_atoms)):
            print "# id= %3d  atom= %s  (%7.3f %7.3f %7.3f)" % (i, _nt_atoms[i], res_xyz[i][0], res_xyz[i][1], res_xyz[i][2])

    
    def matchBase(self, _tpt_nt, _tpt_res_xyz, rna_nt, vb=0): # deprecated
        """ mutate base to template to target nt:  _tpt_nt --> rna_nt  """
        if _tpt_nt == rna_nt: return _tpt_res_xyz # base match
        
        
    def mutate(self, old_nt, new_nt, nt_xyz, vb=0):
        if old_nt == new_nt: return nt_xyz
        if vb: print "# mutation: %s -->  %s " % (old_nt, new_nt)
        new_tpt_ref, new_anchor = self.getRefBaseRefXYZ_UV(new_nt)  # standard ring atoms 
        old_tpt_ref, old_anchor = self.getResBaseRefXYZ_UV(old_nt, nt_xyz) 
        hm = self.getHM(old_tpt_ref, new_tpt_ref, is_center=False, vb=0)
        new_base_atoms = base_anames[new_nt]
        
        old_nt_atoms = nt_atoms_woh[old_nt]
        new_nt_atoms = nt_atoms_woh[new_nt]
        # mutate base atoms' coordinates 
        new_base = self.ref_bases_dt[new_nt]
        #new_hm = (hm[0], new_anchor, old_anchor)
        #print "# new ah: ", new_anchor
        #print "# old ah: ", old_anchor 
        new_res = {}
        for _atom in new_base_atoms:
            _atom_xyz = array(new_base[_atom]) 
            #print "# bemut: ", _atom_xyz
            #print "hm: ", hm 
            _atom_xyz = self.moveBase(_atom_xyz, hm, old_anchor, new_anchor)  
            #print "# aftmut: ", _atom_xyz
            new_res[_atom] = _atom_xyz 
        
        # copy backbone atoms' coordinates        
        for _atom in nt_bb_atoms_wS:
            atom_id = old_nt_atoms.index(_atom)
            new_res[_atom] = nt_xyz[atom_id].copy()
        return self.resDict2Arr(new_nt_atoms, new_res)
    

    
    def resDict2Arr(self, nt_atoms, res_dt):
        """ transfer nt in dictionary format to array format """
        res_xyz = []
        for _atom in nt_atoms:
            res_xyz.append(res_dt[_atom])
        return array(res_xyz)
     
      
    
    def getResRefFrame(self, res_xyz): 
        """ return residue reference frame: A,G: P, C4',C1',N9,  C,U:, P, C4',C1',N1 """
        ref_xyz = []
        for atom_id in [0,5,7,8]:
            ref_xyz.append(res_xyz[atom_id,:].copy())
        return array(ref_xyz)

    
    def chkResNumAtom(self, seq, resids, seg_xyz):
        """ check residue's heavy atoms and make sure no missing atoms """
        if seq[0] == '-': seq = seq[1:]
        for i in range(len(seq)):
            res_name = seq[i]
            resid = resids[i]
            res_xyz = seg_xyz[resid]
            num_atom ,dim = res_xyz.shape 
            num_res_heavy_atom = len(nt_atoms_woh[res_name])
            if num_res_heavy_atom != num_atom: return 0 
        return 1 
   
    
    def getRefBaseRefXYZ_UV(self, res_nt_name):
        """ return  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 """
        res_xyz = self.ref_bases_dt[res_nt_name]
        if res_nt_name in ['A','G']:
            N9 = array(res_xyz['N9'])
            C4 = array(res_xyz['C4'])
            C8 = array(res_xyz['C8'])
            u =  C4 - N9 
            v =  C8 - N9 
            anchor  = N9 
        else:
            N1 = array(res_xyz['N1'])
            C2 = array(res_xyz['C2'])
            C6 = array(res_xyz['C6'])
            u = C2 - N1 
            v = C6 - N1 
            anchor = N1 
        u  = self.normVect(u)
        v = self.normVect(v)
        y = self.normVect( u + v )
        z = self.normVect( self.cross(u,v))
        x = self.normVect( self.cross(y,z))
        return array([x,y,z]), anchor


    def getResBaseRefXYZ_UV(self, res_nt_name, res_xyz):
        """
	return base referenace frame x,y,z 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

	# nt_name: A 
	#    7 --> C1' 
	#    9 --> C4  
	#   17 --> C8  
	#    8 --> N9  
	# nt_name: G 
	#    7 --> C1' 
	#    9 --> C4  
	#   18 --> C8  
	#    8 --> N9  
	# nt_name: C 
	#    7 --> C1' 
	#   10 --> C2  
	#    9 --> C6  
	#    8 --> N1  
	# nt_name: U 
	#    7 --> C1' 
	#   10 --> C2  
	#    9 --> C6  
	#    8 --> N1  

	A: N9,C4,C8, G: N9, C4, C8,  C:N1,C2,C6, U:N1,C2,C6
	"""
        ref_atom_ids = {'A':(8,9,17), 'G':(8,9,18), 'C':(8,10,9),'U':(8,10,9)} 
        res_ref_atom_ids = ref_atom_ids[res_nt_name]
        u = res_xyz[res_ref_atom_ids[1],:]-res_xyz[res_ref_atom_ids[0]]
        v = res_xyz[res_ref_atom_ids[2],:]-res_xyz[res_ref_atom_ids[0]]
        center = res_xyz[8,:].copy()
        u  = self.normVect(u)
        v = self.normVect(v)
        y = self.normVect( u + v )
        z = self.normVect( self.cross(u,v))
        x = self.normVect( self.cross(y,z))
        return array([x,y,z]),center 



if (__name__ == "__main__"):
    pass 

