Commit f4406138 authored by cmaffeo2's avatar cmaffeo2
Browse files

Fixed bug in atomic structure builder

parent 0340c25f
......@@ -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'
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment