Commit 3fe4bfc3 authored by cmaffeo2's avatar cmaffeo2
Browse files

Added support for tclForces

parent 1e16cc4d
......@@ -6,7 +6,7 @@ from math import pi,sqrt,exp,floor
import numpy as np
import random
import os, sys, subprocess
import pdb
from pdb import set_trace
from coords import minimizeRmsd, rotationAboutAxis
from CanonicalNucleotideAtoms import canonicalNtFwd, canonicalNtRev, seqComplement
......@@ -202,7 +202,7 @@ class Segment():
class atomicModel():
def __init__(self, part, ntScale=1.0, randomSeed=1):
self.chunks = dict()
self.dsChunks = dict()
self.numNodes = 0
self.segments = []
......@@ -227,9 +227,11 @@ class atomicModel():
# self._dihedralParamFiles = set()
random.seed(randomSeed)
self.useTclForces = False
self.latticeType = atomicModel._determineLatticeType(part)
self._buildModel(part)
def _strandOccupancies(self, helixId):
try:
......@@ -323,7 +325,7 @@ class atomicModel():
for hid in range(numHID):
self.helixNtsFwd[hid] = dict()
self.helixNtsRev[hid] = dict()
self.chunks[hid] = []
self.dsChunks[hid] = []
## Loop over oligos/segments
segments = []
......@@ -375,12 +377,18 @@ class atomicModel():
## We place nts as pairs; spacing is not always uniform along strand.
## Here we find the locations between which spacing will be uniform.
compIdxs = [comp.idxs() for comp in strand.getComplementStrands()]
chunks = self.combineCompactRegionLists(
chunks = self.combineRegionLists(
[[lo,hi]],
compIdxs,
intersect=True
)
self.chunks[hid].extend( chunks ) # Add running list of chunks for push bonds
self.dsChunks[hid].extend( chunks ) # Add running list of chunks for push bonds
chunks = self.combineRegionLists(
[[lo,hi]],
chunks,
intersect=False
)
## Find zIdx for each nucleotide
zIdxs = []
......@@ -390,10 +398,14 @@ class atomicModel():
for ins in strand.insertionsOnStrand(l,h):
nts += ins.length()
assert( nts > 0 )
for i in range(nts):
zIdxs.append( l + (h-l) * float(i)/(nts-1) )
if nts == 1:
zIdxs.append( l )
else:
zIdxs.append( l + (h-l) * float(i)/(nts-1) )
assert(len(zIdxs) == numNt)
assert( len(zIdxs) == numNt )
if isFwd:
helixNts = self.helixNtsFwd[hid]
......@@ -645,7 +657,7 @@ class atomicModel():
nt.firstAtomicIndex = natoms
natoms += len(nt.atoms(transform=False))
def combineRegionLists(self,loHi1,loHi2):
def combineRegionLists(self,loHi1,loHi2,intersect=False):
"""Combines two lists of (lo,hi) pairs specifying integer
regions a single list of regions. """
......@@ -657,6 +669,17 @@ class atomicModel():
assert(len(pair) == 2)
assert(pair[0] <= pair[1])
if len(loHi1) == 0:
if intersect:
return []
else:
return loHi2
if len(loHi2) == 0:
if intersect:
return []
else:
return loHi1
## Break input into lists of compact regions
compactRegions1,compactRegions2 = [[],[]]
for compactRegions,loHi in zip(
......@@ -712,6 +735,11 @@ class atomicModel():
# print("loHi2:",loHi2)
# print(result,"\n")
if intersect:
lo = max( [loHi1[0][0], loHi2[0][0]] )
hi = min( [loHi1[-1][1], loHi2[-1][1]] )
result = [r for r in result if r[0] >= lo and r[1] <= hi]
return result
def combineCompactRegionLists(self,loHi1,loHi2,intersect=False):
......@@ -739,6 +767,16 @@ class atomicModel():
for pair1,pair2 in zip(l[::2],l[1::2]):
assert(pair1[1]+1 == pair2[0])
if len(loHi1) == 0:
if intersect:
return []
else:
return loHi2
if len(loHi2) == 0:
if intersect:
return []
else:
return loHi1
## Find the ends of the region
lo = min( [loHi1[0][0], loHi2[0][0]] )
......@@ -780,9 +818,9 @@ class atomicModel():
# print(returnList,"\n")
return returnList
def _getChunksForHelixPair(self,hid1,hid2):
def _getDsChunksForHelixPair(self,hid1,hid2):
chunks = []
chunks1,chunks2 = [sorted(list(set(self.chunks[hid]))) for hid in [hid1,hid2]]
chunks1,chunks2 = [sorted(list(set(self.dsChunks[hid]))) for hid in [hid1,hid2]]
return self.combineRegionLists(chunks1,chunks2)
i,j = [0,0]
......@@ -1043,7 +1081,6 @@ class atomicModel():
def writeENM(self, prefix):
pushCount = 0
if self.latticeType == "square":
enmTemplate = enmTemplateSQ
elif self.latticeType == "honeycomb":
......@@ -1103,128 +1140,172 @@ class atomicModel():
print("NO STACKS found for:", noStackPrime)
print("NO BPs:", noBasepair)
props, origins = self.part.helixPropertiesAndOrigins()
origins = [np.array(o) for o in origins]
## Push bonds
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._getChunksForHelixPair(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)))
print("NO BASEPAIRS found for:", noBasepair)
ntPairs = [[nts1[i],nts2[j]] for i,j in zip(iVals,jVals)]
props, origins = self.part.helixPropertiesAndOrigins()
origins = [np.array(o) for o in origins]
## 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)])
## Push bonds
pushBonds = []
skip=minDist-xoDist
if skip > 0:
ntPairs = ntPairs[skip:]
fwdDict = dict() # dictionaries of nts with keys [hid][zidx]
revDict = dict()
xo = dict() # dictionary of crossovers between [hid1][hid2]
## 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)])
helixIds = [hid for s in self.segments for hid in s.nodes.keys()]
helixIds = sorted(list(set(helixIds)))
skip=minDist-xoDist
if skip > 0:
ntPairs = ntPairs[:-skip]
#- 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']
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
## 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
ai,aj = [nt.atoms(transform=False) for nt in (nt1,nt2)]
i += ai.index('P')-1
j += aj.index('P')-1
pushBonds.append(i)
pushBonds.append(j)
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}
}
ai,aj = [nt.atoms(transform=False) for nt in (nt1,nt2)]
i += ai.index('P')-1
j += aj.index('P')-1
fh.write("bond %d %d %f %.2f\n" % (i,j,1,31))
pushCount += 1
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)
}
}
""")
print("PUSH BONDS:", pushCount)
def writeNamdFile(self,prefix,numSteps=4800000):
with open("%s.namd" % prefix,'w') as fh:
......@@ -1266,8 +1347,10 @@ switching on
switchdist 8
cutoff 10
pairlistdist 12
margin 30
""")
if not self.useTclForces:
fh.write("margin 30\n")
fh.write("""
# Integrator Parameters
timestep 2.0 ;# 2fs/step
rigidBonds all ;# needed for 2fs steps
......@@ -1300,7 +1383,13 @@ restartfreq 48000
# ENM and intrahelical extrabonds
extraBonds on
extraBondsFile $prefix.exb
""")
if self.useTclForces:
fh.write("""
tclForces on
tclForcesScript $prefix.forces.tcl
""")
fh.write("""
#############################################################
## RUN ##
#############################################################
......
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