run.py 3.40 KiB
# -*- 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 = 90
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)
# return coarsestModel._writeArbdDihedralFile(outPrefix+"-1")
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=60000000, 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,
ssDnaHelixNeighborDist=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, ssDnaHelixNeighborDist=50)
fineModel.simulate( out, outputDirectory='output', numSteps=1000000,
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)
files = sorted(glob.glob('json/*.json'))
gpu=1
f = files[1]
# f = files
print(f)
out = re.match('json/(.*).json',f).group(1)
simulateDesign(f,out,gpu=gpu)