Fixed problems with ENM for atomic model

atomicModel.py

@@ -378,6 +378,7 @@ class atomicModel():
                         last = s
                 strandSegments = [(i+1,j) for i,j in zip([lo-1]+split,split+[hi])]
 
+
                 ## Find zIdx for each nucleotide
                 zIdxs = []
                 for l,h in strandSegments:
@@ -391,25 +392,34 @@ class atomicModel():
 
                 assert(len(zIdxs) == numNt)
 
-                baseIds = range(numNt)
                 if isFwd:
                     baseToSeg = baseToSegFwd[hid]
                 else:
-                    baseIds = reversed(baseIds)
                     baseToSeg = baseToSegRev[hid]
                     zIdxs = reversed(zIdxs)
 
+                ## Find strandOccupancy index for each nucleotide
+                strandOccIdx = []
+                for l,h in strandSegments:
+                    for i in range(l,h+1):
+                        nts = 1
+                        for ins in strand.insertionsOnStrand(i,i):
+                            nts += ins.length()
+                        baseToSeg[i] = []
+                        for j in range(nts):
+                            strandOccIdx.append(i)
+                if not isFwd:
+                    strandOccIdx = reversed(strandOccIdx)
+
                 ## Add nucleotides
-                for i,s,zIdx in zip(baseIds,seqs,zIdxs):
+                for s,zIdx,idx in zip(seqs,zIdxs,strandOccIdx):
                     p = zIdxToPos(zIdx) 
                     a = zIdxToAngle(zIdx)
 
-                    lastNt = seg.addNucleotide(hid, i, zIdx, s, isFwd,
+                    lastNt = seg.addNucleotide(hid, idx, zIdx, s, isFwd,
                                                p, a, lastNt)
 
-                    # baseToSeg[i+1] = [lastNt]
-                    assert(zIdx not in baseToSeg)
-                    baseToSeg[zIdx] = [lastNt]
+                    baseToSeg[idx].append(lastNt)
 
         ## 2) Find basepairs and stacks
         for hid in range(numHID):
@@ -418,7 +428,6 @@ class atomicModel():
                 continue
 
             ends1,ends2 = self._helixStrandsToEnds(helixStrands)
-
             strandOccupancies = [ [x for i in range(0,len(e),2) 
                                    for x in range(e[i],e[i+1]+1)] 
                                   for e in (ends1,ends2) ]
@@ -435,7 +444,7 @@ class atomicModel():
                                    reversed(baseToSegRev[hid][i])):
                     self._pairBases(nt1,nt2)
                 
-                    if prevBasepair is not None and prevBasepair[0] == i-1:
+                    if prevBasepair is not None and i - prevBasepair[0] <= 3: # TODO: find better way of locating stacks
                         atomicModel._stackBases( prevBasepair[1], nt1)
                         atomicModel._stackBases( nt2, prevBasepair[2])
                     # else:
@@ -450,16 +459,6 @@ class atomicModel():
         below.addNodeAbove(above, sep)
         above.addNodeBelow(below, sep)
         
-    def _addCrossover(self, hid1, hid2, xo):
-        zid1, zid2, isFwd1, isFwd2 = xo
-        node1 = self.helices[hid1].nodes[zid1]
-        node2 = self.helices[hid2].nodes[zid2]
-
-        ## TODO add polarity
-        polarity = 0
-        node1.addXover(node2, isFwd1)
-        node2.addXover(node1, isFwd2)
-
     def _getNeighborHelixDict(part):
         props, origins = part.helixPropertiesAndOrigins()
         neighborHelices = dict()
@@ -860,7 +859,8 @@ class atomicModel():
             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:
@@ -875,12 +875,17 @@ class atomicModel():
                         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:
@@ -888,10 +893,11 @@ class atomicModel():
                         for n1, n2, d in enmTemplate[key]:
                             d = float(d)
                             a2 = nt2.atoms(transform=False)
-                            try:
-                                i,j = [a.index(n)-1 for a,n in zip((a1,a2),(n1,n2))]
-                            except:
-                                continue
+                            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':
@@ -905,6 +911,8 @@ class atomicModel():
                             fh.write("bond %d %d %f %.2f\n" % (i,j,k,d))
 
 
+            print("NO STACKS found for : %d" % noStackPrime)
+            print("NO BPs: %d" % noBasepair)
             props, origins = self.part.helixPropertiesAndOrigins()
             origins = [np.array(o) for o in origins]
             
@@ -956,7 +964,7 @@ class atomicModel():
                         if nt.basepair is None: continue
                         hid = nt.prop['helixId']
                         if hid not in xo:
-                            print("WARNING: writeENM: helix %d not in xo dict" % hid)
+                            # print("WARNING: writeENM: helix %d not in xo dict" % hid)
                             continue
                         zidx = nt.prop['zIdx']
                         isFwd = nt.prop['isFwd']
@@ -964,7 +972,7 @@ class atomicModel():
                         for hid2 in neighborHelices[hid]:
                             # if hid2 not in xo[hid]:
                             if hid2 not in helixIds:
-                                print("WARNING: writeENM: helix %d not in helixIds" % hid2) # xo[%d] dict" % (hid2,hid))
+                                # print("WARNING: writeENM: helix %d not in helixIds" % hid2) # xo[%d] dict" % (hid2,hid))
                                 continue
                             dist = min([np.abs(z-zidx) for z in xo[hid][hid2]])
                             key = (hid2,zidx,isFwd)

run.py

@@ -95,7 +95,7 @@ def simulateDesign(inFile, outPrefix, numCoarsestSteps=5000000, gpu=0, shiftEule
     finerModel = beadModelTwist(part, twistPersistenceLength=twistLp, 
                           maxBpsPerDNode=1, maxNtsPerSNode=1)
 
-    finestModel = atomicModel(part, ntScale=0.5)
+    finestModel = atomicModel(part, ntScale=0.25)
 
     coarsestModelNoLongBonds._removeIntrahelicalConnectionsBeyond(100)
     coarsestModelNoLongBonds._removeCrossoversBeyond(100)