Commit 4d0a84c7 authored by cmaffeo2's avatar cmaffeo2
Browse files

Imported old ENM code, probably breaking python script

parent 7b4b0dcf
......@@ -2209,6 +2209,245 @@ class SegmentModel(ArbdModel):
for s in self.strands:
s.update_atomic_orientations(orientation)
def._write_atomic_ENM(self, prefix, lattice_type=None):
## TODO: ensure atomic model was generated already
if lattice_type is None:
lattice_type = self.lattice_type
if lattice_type == "square":
enmTemplate = enmTemplateSQ
elif lattice_type == "honeycomb":
enmTemplate = enmTemplateHC
else:
raise Exception("Lattice type '%s' not supported" % self.latticeType)
noStackPrime = 0
noBasepair = 0
with open("%s.exb" % prefix,'w') as fh:
# natoms=0
for seg in self.segments:
for nt1 in seg:
other = []
nt2 = nt1.basepair
if nt2 is not None:
if nt2.firstAtomicIndex > nt1.firstAtomicIndex:
other.append((nt2,'pair'))
nt2 = nt2.stack3prime
if nt2 is not None and nt2.firstAtomicIndex > nt1.firstAtomicIndex:
other.append((nt2,'paircross'))
else:
noBasepair += 1
nt2 = nt1.stack3prime
if nt2 is not None:
other.append((nt2,'stack'))
nt2 = nt2.basepair
if nt2 is not None and nt2.firstAtomicIndex > nt1.firstAtomicIndex:
other.append((nt2,'cross'))
else:
noStackPrime += 1
a1 = nt1.atoms(transform=False)
for nt2,key in other:
key = ','.join((key,nt1.sequence[0],nt2.sequence[0]))
for n1, n2, d in enmTemplate[key]:
d = float(d)
a2 = nt2.atoms(transform=False)
i,j = [a.index(n)-1 for a,n in zip((a1,a2),(n1,n2))]
# try:
# i,j = [a.index(n)-1 for a,n in zip((a1,a2),(n1,n2))]
# except:
# continue
k = 0.1
if self.latticeType == 'honeycomb':
correctionKey = ','.join((key,n1,n2))
assert(correctionKey in enmCorrectionsHC)
dk,dr = enmCorrectionsHC[correctionKey]
k = float(dk)
d += float(dr)
i += nt1.firstAtomicIndex
j += nt2.firstAtomicIndex
fh.write("bond %d %d %f %.2f\n" % (i,j,k,d))
print("NO STACKS found for:", noStackPrime)
print("NO BASEPAIRS found for:", noBasepair)
props, origins = self.part.helixPropertiesAndOrigins()
origins = [np.array(o) for o in origins]
## Push bonds
pushBonds = []
fwdDict = dict() # dictionaries of nts with keys [hid][zidx]
revDict = dict()
xo = dict() # dictionary of crossovers between [hid1][hid2]
helixIds = [hid for s in self.segments for hid in s.nodes.keys()]
helixIds = sorted(list(set(helixIds)))
#- initialize dictionaries
for h1 in helixIds:
fwdDict[h1] = dict()
revDict[h1] = dict()
xo[h1] = dict()
for h2 in helixIds:
xo[h1][h2] = []
#- fill dictionaries
for seg in self.segments:
for nt in seg:
hid = nt.prop['helixId']
## Add nt to ntDict
idx = nt.prop['idx']
zidx = nt.prop['zIdx']
isFwd = nt.prop['isFwd']
if isFwd:
fwdDict[hid][idx] = nt
else:
revDict[hid][idx] = nt
## Find crossovers and ends
for nt2 in (nt.ntAt5prime, nt.ntAt3prime):
if nt2 is not None:
hid2 = nt2.prop['helixId']
if hid2 != hid:
# xo[hid][hid2].append(zidx)
xo[hid][hid2].append(idx)
xo[hid2][hid].append(idx)
props = self.part.getModelProperties()
if props.get('point_type') == PointType.ARBITRARY:
raise Exception("Not implemented")
else:
props, origins = self.part.helixPropertiesAndOrigins()
neighborHelices = atomicModel._getNeighborHelixDict(self.part)
if self.latticeType == 'honeycomb':
# minDist = 8 # TODO: test against server
minDist = 11 # Matches ENRG MD server... should it?
elif self.latticeType == 'square':
minDist = 11
for hid1 in helixIds:
for hid2 in neighborHelices[hid1]:
if hid2 not in helixIds:
# print("WARNING: writeENM: helix %d not in helixIds" % hid2) # xo[%d] dict" % (hid2,hid1))
continue
## Build chunk for helix pair
chunks = self._getDsChunksForHelixPair(hid1,hid2)
for lo,hi in chunks:
# pdb.set_trace()
nts1 = self.getNtsBetweenCadnanoIdxInclusive(hid1,lo,hi)
nts2 = self.getNtsBetweenCadnanoIdxInclusive(hid2,lo,hi)
if len(nts1) <= len(nts2):
iVals = list(range(len(nts1)))
if len(nts1) > 1:
jVals = [int(round(j*(len(nts2)-1)/(len(nts1)-1))) for j in range(len(nts1))]
else:
jVals = [0]
else:
if len(nts2) > 1:
iVals = [int(round(i*(len(nts1)-1)/(len(nts2)-1))) for i in range(len(nts2))]
else:
iVals = [0]
jVals = list(range(len(nts2)))
ntPairs = [[nts1[i],nts2[j]] for i,j in zip(iVals,jVals)]
## Skip pairs close to nearest crossover on lo side
xoIdx = [idx for idx in xo[hid1][hid2] if idx <= lo]
if len(xoIdx) > 0:
xoIdx = max(xoIdx)
assert( type(lo) is int )
xoDist = min([len(self.getNtsBetweenCadnanoIdxInclusive(hid,xoIdx,lo-1)) for hid in (hid1,hid2)])
skip=minDist-xoDist
if skip > 0:
ntPairs = ntPairs[skip:]
## Skip pairs close to nearest crossover on hi side
xoIdx = [idx for idx in xo[hid1][hid2] if idx >= hi]
if len(xoIdx) > 0:
xoIdx = min(xoIdx)
xoDist = min([len(self.getNtsBetweenCadnanoIdxInclusive(hid,hi+1,xoIdx)) for hid in (hid1,hid2)])
skip=minDist-xoDist
if skip > 0:
ntPairs = ntPairs[:-skip]
for ntPair in ntPairs:
bps = [nt.basepair for nt in ntPair]
if None in bps: continue
for nt1,nt2 in [ntPair,bps]:
i,j = [nt.firstAtomicIndex for nt in (nt1,nt2)]
if j <= i: continue
if np.linalg.norm(nt1.position-nt2.position) > 45: continue
ai,aj = [nt.atoms(transform=False) for nt in (nt1,nt2)]
try:
i += ai.index('P')-1
j += aj.index('P')-1
pushBonds.append(i)
pushBonds.append(j)
except:
pass
print("PUSH BONDS:", len(pushBonds)/2)
if not self.useTclForces:
with open("%s.exb" % prefix, 'a') as fh:
for i,j in zip(pushBonds[::2],pushBonds[1::2]):
fh.write("bond %d %d %f %.2f\n" % (i,j,1.0,31))
else:
atomList = list(set(pushBonds))
with open("%s.forces.tcl" % prefix,'w') as fh:
fh.write("set atomList {%s}\n\n" %
" ".join([str(x-1) for x in atomList]) )
fh.write("set bonds {%s}\n" %
" ".join([str(x-1) for x in pushBonds]) )
fh.write("""
foreach atom $atomList {
addatom $atom
}
proc calcforces {} {
global atomList bonds
loadcoords rv
foreach i $atomList {
set force($i) {0 0 0}
}
foreach {i j} $bonds {
set rvec [vecsub $rv($j) $rv($i)]
# lassign $rvec x y z
# set r [expr {sqrt($x*$x+$y*$y+$z*$z)}]
set r [getbond $rv($j) $rv($i)]
set f [expr {2*($r-31.0)/$r}]
vecadd $force($i) [vecscale $f $rvec]
vecadd $force($j) [vecscale [expr {-1.0*$f}] $rvec]
}
foreach i $atomList {
addforce $i $force($i)
}
}
""")
def _generate_atomic_model(self, scale=1):
self.children = self.strands
......
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