diff --git a/beadModel.py b/beadModel.py
index 9d0edbb4252509b0ac22abfe0810961628a81422..47bca5b497bc5899331fa72c6241347458395c35 100644
--- a/beadModel.py
+++ b/beadModel.py
@@ -418,12 +418,9 @@ class beadModel():
## Find crossovers for this helix
xoZids = [x[0] for h0 in range(hid) if hid in xoDicts[h0] for x in xoDicts[h0][hid]]
xoZids.extend([x[0] for hid2,xos in xoDicts[hid].items() for x in xos])
- xoZids.extend([x[0] for h0 in range(hid) if hid in ssXoDicts[h0] for x in xoDicts[h0][hid]])
+ xoZids.extend([x[0] for h0 in range(hid) if hid in ssXoDicts[h0] for x in ssXoDicts[h0][hid]])
xoZids.extend([x[0] for hid2,xos in ssXoDicts[hid].items() for x in xos])
-
reqNodeZids = sorted(list(set( ends1 + ends2 + xoZids ) ) )
- # [x for x in xoDict.keys()] +
- # [x for x in xoDictRev.keys()] )))
## Build lists of which nt sites are occupied in the helix
strandOccupancies = [ [x for i in range(0,len(e),2)
@@ -988,10 +985,15 @@ component "data" value 3
## units "3e-19 erg cm/ 295 k K" "nm" =~ 73
Lp = self.twistPersistenceLength/0.34 # set semi-arbitrarily as there is a large spread in literature
- # k = exp(-float(sep)/Lp) * kT * 0.00030461742
+ # fitFun = lambda x: np.real(erf( (4*np.pi*x + 1j)/(2*np.sqrt(x)) )) * np.exp(-1/(4*x)) / erf(2*np.sqrt(x)*np.pi) - exp(-sep/Lp)
+ # k = opt.leastsq( fitFun, x0=exp(-sep/Lp) )
+ # k = k[0][0] * 2*kT*0.00030461742
+
+ intrinsicDegrees=20
+ fitFun = lambda x: (1.0/(2*x) - 2*np.sqrt(np.pi)*np.exp(-4*np.pi**2*x) / (np.sqrt(x)*erf(2*np.pi*np.sqrt(x))) ) - \
+ ( (intrinsicDegrees*np.pi/180)**2 + 2*(1-exp(-sep/Lp)) )
- fitFun = lambda x: np.real(erf( (4*np.pi*x + 1j)/(2*np.sqrt(x)) )) * np.exp(-1/(4*x)) / erf(2*np.sqrt(x)*np.pi) - exp(-sep/Lp)
- k = opt.leastsq( fitFun, x0=exp(-sep/Lp) )
+ k = opt.leastsq( fitFun, x0=1/(1-exp(-sep/Lp)) )
k = k[0][0] * 2*kT*0.00030461742
t0 = sep*(360.0/10.5)
diff --git a/encode-1bp.py b/encode-1bp.py
deleted file mode 100644
index b9f30a5e6c99fa8ac5cc61a9bdf6375ee2dbd4e7..0000000000000000000000000000000000000000
--- a/encode-1bp.py
+++ /dev/null
@@ -1,85 +0,0 @@
-# -*- coding: utf-8 -*-
-import json
-import cadnano
-from cadnano.document import Document
-import numpy as np
-
-from beadModel import beadModel
-from atomicModel import atomicModel
-
-
-twistLp = 40
-outPrefix0 = 'twist-%d' % twistLp
-inFile = 'input.json'
-
-# ----------------- #
-# Read cadnano file #
-# ----------------- #
-
-with open(inFile) as ch:
- data = json.load(ch)
-try:
- doc = Document()
- cadnano.fileio.v3decode.decode(doc, data)
-except:
- doc = Document()
- cadnano.fileio.v2decode.decode(doc, data)
-
-
-def readArbdCoords(fname):
- coords = []
- with open(fname) as fh:
- for line in fh:
- coords.append([float(x) for x in line.split()[1:]])
- return np.array(coords)
-
-
-parts = [p for p in doc.getParts()]
-if len(parts) != 1:
- raise Exception("Only documents containing a single 'part' are implemented at this time.")
-part = parts[0]
-
-## Manually fix eulerZ, which is improperly set by cadnano2.5 for some structures
-props, origins = part.helixPropertiesAndOrigins()
-for t,i in zip( props['eulerZ'], range(part.getIdNumMax()+1) ):
- part.setVirtualHelixProperties(i,'eulerZ',t+180.0)
-
-
-# ------------- #
-# Create models #
-# ------------- #
-
-coarsestModel = beadModel(part, twistPersistenceLength=twistLp,
- maxBpsPerDNode=7, maxNtsPerSNode=4)
-coarseModel = beadModel(part, twistPersistenceLength=twistLp,
- maxBpsPerDNode=4, maxNtsPerSNode=2)
-fineModel = beadModel(part, twistPersistenceLength=twistLp,
- maxBpsPerDNode=1, maxNtsPerSNode=1)
-# atomicModel = atomicModel(part, ntScale=0.5)
-
-
-# ------------------------------- #
-# Simulate and backmap each model #
-# ------------------------------- #
-
-gpu = 0
-outPrefix = outPrefix0 + "1"
-coarsestModel.simulate(outPrefix, outputDirectory='output', numSteps=5000000, timestep=300e-6,
- arbd="/home/cmaffeo2/development/cuda/arbd.dev/src/arbd", gpu=gpu)
-coarsestModelCoords = readArbdCoords( "output/%s.0.restart" % outPrefix )
-
-outPrefix = outPrefix0 + "2"
-coarseModel.backmap( coarsestModel, coarsestModelCoords,
- dsDnaHelixNeighborDist=100, dsDnaAllNeighborDist=50)
-coarseModel.simulate( outPrefix, outputDirectory='output', numSteps=1000000,
- arbd="/home/cmaffeo2/development/cuda/arbd.dev/src/arbd", gpu=gpu)
-coarseModelCoords = readArbdCoords( "output/%s.0.restart" % outPrefix )
-
-fineModel.backmap( coarseModel, coarseModelCoords )
-fineModel.simulate( 'fine', outputDirectory='output', numSteps=1000000,
- arbd="/home/cmaffeo2/development/cuda/arbd.dev/src/arbd", gpu=gpu)
-# fineModelCoords = readArbdCoords( "output/%s.0.restart" % 'fine' )
-
-# atomicModel.randomizeUnassignedNtSeqs()
-# atomicModel.backmap( fineModel, fineModelCoords )
-# atomicModel.writePdbPsf( 'atomic' )
diff --git a/run.py b/run.py
new file mode 100644
index 0000000000000000000000000000000000000000..46e2453b18dab1765d3d0cdc162984e5193558c2
--- /dev/null
+++ b/run.py
@@ -0,0 +1,128 @@
+# -*- coding: utf-8 -*-
+import json
+import cadnano
+from cadnano.document import Document
+import numpy as np
+
+import glob,re
+
+from beadModel import beadModel
+from atomicModel import atomicModel
+
+
+# twistLp = 1e-6
+twistLp = 40
+def simulateDesign(inFile, outPrefix, gpu=0, shiftEulerZ=0):
+ # ----------------- #
+ # Read cadnano file #
+ # ----------------- #
+
+
+ with open(inFile) as ch:
+ data = json.load(ch)
+ try:
+ doc = Document()
+ cadnano.fileio.v3decode.decode(doc, data)
+ except:
+ doc = Document()
+ cadnano.fileio.v2decode.decode(doc, data)
+
+
+ def readArbdCoords(fname):
+ coords = []
+ with open(fname) as fh:
+ for line in fh:
+ coords.append([float(x) for x in line.split()[1:]])
+ return np.array(coords)
+
+
+ parts = [p for p in doc.getParts()]
+ if len(parts) != 1:
+ raise Exception("Only documents containing a single 'part' are implemented at this time.")
+ part = parts[0]
+
+ if shiftEulerZ != 0:
+ ## Manually fix eulerZ, which is improperly set by cadnano2.5 for some structures
+ props, origins = part.helixPropertiesAndOrigins()
+ for t,i in zip( props['eulerZ'], range(part.getIdNumMax()+1) ):
+ part.setVirtualHelixProperties(i,'eulerZ',t+180.0)
+
+
+ # ------------- #
+ # Create models #
+ # ------------- #
+
+ coarsestModel = beadModel(part, twistPersistenceLength=twistLp,
+ maxBpsPerDNode=7, maxNtsPerSNode=4)
+ coarseModel = beadModel(part, twistPersistenceLength=twistLp,
+ maxBpsPerDNode=4, maxNtsPerSNode=2)
+ fineModel = beadModel(part, twistPersistenceLength=twistLp,
+ maxBpsPerDNode=1, maxNtsPerSNode=1)
+ finestModel = atomicModel(part, ntScale=0.5)
+
+
+ # ------------------------------- #
+ # Simulate and backmap each model #
+ # ------------------------------- #
+
+ out = outPrefix + "-1"
+ coarsestModel.simulate(out, outputDirectory='output', numSteps=20000000, timestep=200e-6,
+ arbd="/home/cmaffeo2/development/cuda/arbd.dev/src/arbd", gpu=gpu)
+ coarsestModelCoords = readArbdCoords( "output/%s.0.restart" % out )
+
+ out = outPrefix + "-2"
+ coarseModel.backmap( coarsestModel, coarsestModelCoords,
+ dsDnaHelixNeighborDist=100, dsDnaAllNeighborDist=50)
+ coarseModel.simulate( out, outputDirectory='output', numSteps=1000000,
+ arbd="/home/cmaffeo2/development/cuda/arbd.dev/src/arbd", gpu=gpu)
+ coarseModelCoords = readArbdCoords( "output/%s.0.restart" % out )
+
+ out = outPrefix + "-3"
+ fineModel.backmap( coarseModel, coarseModelCoords )
+ fineModel.simulate( out, outputDirectory='output', numSteps=2000000,
+ arbd="/home/cmaffeo2/development/cuda/arbd.dev/src/arbd", gpu=gpu)
+ fineModelCoords = readArbdCoords( "output/%s.0.restart" % out )
+
+ out = outPrefix + "-4"
+ finestModel.randomizeUnassignedNtSeqs()
+ finestModel.backmap( fineModel, fineModelCoords )
+ finestModel.writePdbPsf( out )
+
+
+# def simulateDesign(f,out,gpu=0,shiftEulerZ=180):
+# print(out)
+
+gpu=0
+for f in list(reversed(glob.glob('json/bend-*')))[1::2]:
+ out = re.match('json/(.*).json',f).group(1)
+ # simulateDesign(f,out,gpu=gpu,shiftEulerZ=180)
+gpu+=1
+
+for f in list(reversed(glob.glob('json/bend-*')))[::4]:
+ out = re.match('json/(.*).json',f).group(1)
+ # simulateDesign(f,out,gpu=gpu,shiftEulerZ=180)
+gpu+=1
+
+gpu=0
+for f in list(reversed(glob.glob('json/bend-*')))[2::4]:
+ out = re.match('json/(.*).json',f).group(1)
+ # simulateDesign(f,out,gpu=gpu,shiftEulerZ=180)
+gpu+=1
+
+for f in glob.glob('json/*twist*.json'):
+ out = re.match('json/(.*).json',f).group(1)
+ # simulateDesign(f,out,gpu=gpu)
+gpu+=1
+
+gpu=0
+for f in ('json/right_twist_mod.json',):
+ out = re.match('json/(.*).json',f).group(1)
+ simulateDesign(f,out,gpu=gpu)
+gpu+=1
+
+
+# for f in ('json/Lowerring_V2_sticky_2_opt_allT.json',:
+# out = re.match('json/(.*).json',f).group(1)
+# simulateDesign(f,out,gpu=gpu,shiftEulerZ=180)
+# gpu+=1
+