Skip to content
Snippets Groups Projects
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)