diff --git a/CanonicalNucleotideAtoms.py b/CanonicalNucleotideAtoms.py index e8403c05c78dcd2523bbee1d90e17d8d3467e47f..e793888d6b0bac26316167aea4de5a53ef39f99d 100644 --- a/CanonicalNucleotideAtoms.py +++ b/CanonicalNucleotideAtoms.py @@ -8,9 +8,6 @@ seqComplement = dict(A='T',G='C') resnames = dict(A="ADE", G="GUA", C="CYY", T="THY") for x,y in list(seqComplement.items()): seqComplement[y] = x - -atomic_types = dict() - def stringToIntTuples(string, tupleLen, offset): l = string.split() # for i in range(0,len(l),tupleLen): @@ -43,10 +40,10 @@ class CanonicalNucleotideFactory(Group): for i in range(len(tns)): tn = tns[i] if tn not in atomic_types: - atomic_types[tn] = ParticleType(tn.decode('UTF-8'), - charge = d['f3'][i], - mass = d['f4'][i]) - p = PointParticle( type_ = atomic_types[tn], + type_ = ParticleType(tn.decode('UTF-8'), + charge = d['f3'][i], + mass = d['f4'][i]) + p = PointParticle( type_ = type_, position = [xs[i],ys[i],zs[i]], name = names[i].decode('UTF-8'), beta = 1) @@ -83,107 +80,6 @@ class CanonicalNucleotideFactory(Group): nt.add_improper( atoms[i], atoms[j], atoms[k], atoms[l], None ) return nt -class CanonicalNucleotide(Group): - DefaultOrientation = rotationAboutAxis([0,0,1], 80) - def __init__(self, prefix, seq): - self.sequence = seq # TODO: used? - self.resname = resnames[seq] - Group.__init__(self, orientation = CanonicalNucleotide.DefaultOrientation) - - ## idx, type, - d = np.loadtxt("%s.dat" % prefix, dtype='i4,S4,S4,f4,f4,f4,f4,f4') - - ids = d['f0'] - assert( np.all(ids == sorted(ids)) ) - - names = d['f1'] - tns = d['f2'] - xs = d['f5'] - ys = d['f6'] - zs = d['f7'] - - type_list = [] - for i in range(len(tns)): - tn = tns[i] - if tn not in atomic_types: - atomic_types[tn] = ParticleType(tn.decode('UTF-8'), - charge = d['f3'][i], - mass = d['f4'][i]) - p = PointParticle( type_ = atomic_types[tn], - position = [xs[i],ys[i],zs[i]], - name = names[i].decode('UTF-8') ) - self.add( p ) - - - atoms = self.children - self.atoms_by_name = {a.name:a for a in atoms} - - - # minIdx = np.min(self.indices) - 1 - b = np.loadtxt("%s.bonds.txt" % prefix, dtype='i4') # - minIdx - for i,j in zip(b[::2],b[1::2]): - self.add_bond( atoms[i], atoms[j], None ) - - a = np.loadtxt("%s.angles.txt" % prefix, dtype='i4') # - minIdx - for i,j,k in zip(a[::3],a[1::3],a[2::3]): - self.add_angle( atoms[i], atoms[j], atoms[k], None ) - - d = np.loadtxt("%s.dihedrals.txt" % prefix, dtype='i4')# - minIdx - for i,j,k,l in zip(d[::4],d[1::4],d[2::4],d[3::4]): - self.add_dihedral( atoms[i], atoms[j], atoms[k], atoms[l], None ) - - i = np.loadtxt("%s.impropers.txt" % prefix, dtype='i4')# - minIdx - for i,j,k,l in zip(d[::4],d[1::4],d[2::4],d[3::4]): - self.add_improper( atoms[i], atoms[j], atoms[k], atoms[l], None ) - - ## TODO: remove this! - def transformed(self, orientation, translation, scale=0.5, singleStranded=False): - # rotationAboutZ -= 68 - c = copy.copy(self) - c.props = copy.deepcopy(self.props) - - ## Find transformation - R = CanonicalNucleotide.DefaultOrientation.dot( orientation ) - - ## Find offset of nt (shift for ssDNA) - center = np.array((0,0,0)) - if singleStranded: - i0 = self.index("C3'")-1 - center = np.array([c.props[k][i0] for k in ('x','y','z')]) - - - ## Apply transformation - for i in range(len(c)): - rvec = np.array([c.props[k][i] for k in ('x','y','z')]) - center - rvec = rvec.dot( R ) - for k,r,o in zip(('x','y','z'), rvec, translation): - c.props[k][i] = r + o - - if scale != 1: - if singleStranded: - i0 = self.index("C3'")-1 - r0 = np.array([c.props[k][i0] for k in ('x','y','z')]) - for i in range(len(c)): - rvec = np.array([c.props[k][i] for k in ('x','y','z')]) - rvec = scale*(rvec-r0) + r0 - for k,r in zip(('x','y','z'), rvec): - c.props[k][i] = r - c.props["beta"][i] = 0 - else: - if self.sequence in ("A","G"): - i0 = self.index("N9")-1 - else: - i0 = self.index("N1")-1 - r0 = np.array([c.props[k][i0] for k in ('x','y','z')]) - for i in range(len(c)): - if c.props['name'][i][-1] in ("'","P","T"): - rvec = np.array([c.props[k][i] for k in ('x','y','z')]) - rvec = scale*(rvec-r0) + r0 - for k,r in zip(('x','y','z'), rvec): - c.props[k][i] = r - c.props["beta"][i] = 0 - - return c canonicalNtFwd = dict() direction = 'fwd'