Commit 39736a40 authored by cmaffeo2's avatar cmaffeo2
Browse files

Updated run.py to use average coordinates

parent 7c1baa78
......@@ -10,6 +10,8 @@ from beadModel import beadModel
from beadModelTwist import beadModelTwist
from atomicModel import atomicModel
import mdtraj as md
import multiprocessing as mp
import sys
......@@ -19,7 +21,9 @@ twistLp = 90
# arbd="/home/cmaffeo2/development/cuda/arbd.dev/src/arbd"
arbd="/home/cmaffeo2/development/cuda/arbd.dbg/src/arbd" # reduced the mem footprint cause vmd
namd="/home/cmaffeo2/development/namd-bin/NAMD_Git-2017-07-06_Linux-x86_64-multicore-CUDA/namd2"
# namd="/home/cmaffeo2/development/namd-bin/NAMD_Git-2017-07-06_Linux-x86_64-multicore-CUDA/namd2"
# namd="/home/wilsja/namd/NAMD_CVS-2016-11-09_Linux-x86_64-multicore-CUDA/namd2"
namd="/home/cmaffeo2/development/namd-bin/NAMD_2.12_Linux-x86_64-multicore-CUDA/namd2"
def readSequenceFile(sequenceFile='resources/cadnano2pdb.seq.m13mp18.dat'):
seq = []
......@@ -37,6 +41,27 @@ def readArbdCoords(fname):
coords.append([float(x) for x in line.split()[1:]])
return np.array(coords)
def readAvgArbdCoords(psf,pdb,dcd,rmsdThreshold=20):
## align trajectory to pdb
ref = md.load(pdb, top=psf)
sel = md.load(dcd, top=psf)
# ref = sel[-1]
sel.superpose(ref)
## Find where rmsd compared to last frame increases
ref = sel[-1]
rmsds = md.rmsd(sel,ref)*10
t0 = np.where(rmsds > rmsdThreshold)[0]
if len(t0) == 0:
t0 = 0
else:
t0 = t0[-1]
pos = sel.xyz[t0:]
pos = np.mean(pos, axis=0)
return 10*pos # convert to Angstroms
def simulateDesign(inFile, outPrefix, numCoarsestSteps=5000000, gpu=0, shiftEulerZ=0,
sequenceFile='resources/cadnano2pdb.seq.m13mp18.dat'):
# ----------------- #
......@@ -87,7 +112,6 @@ def simulateDesign(inFile, outPrefix, numCoarsestSteps=5000000, gpu=0, shiftEule
maxBpsPerDNode=7, maxNtsPerSNode=4)
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,
......@@ -143,42 +167,43 @@ def simulateDesign(inFile, outPrefix, numCoarsestSteps=5000000, gpu=0, shiftEule
# Simulate and backmap each model #
# ------------------------------- #
out = outPrefix + "-0"
coarsestModelNoLongBonds.simulate(out, outputDirectory='output', numSteps=int(0.2*numCoarsestSteps), timestep=200e-6,
arbd=arbd, gpu=gpu)
coarsestModelCoords = readArbdCoords( "output/%s.0.restart" % out )
out = outPrefix + "-1"
coarsestModel.backmap( coarsestModelNoLongBonds, coarsestModelCoords,
dsDnaHelixNeighborDist=100, dsDnaAllNeighborDist=50,
ssDnaHelixNeighborDist=50)
coarsestModel.simulate(out, outputDirectory='output', numSteps=int(0.8*numCoarsestSteps), timestep=200e-6,
arbd=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=100000,
arbd=arbd, gpu=gpu)
coarseModelCoords = readArbdCoords( "output/%s.0.restart" % out )
out = outPrefix + "-3"
fineModel.backmap( coarseModel, coarseModelCoords,
dsDnaHelixNeighborDist=100, dsDnaAllNeighborDist=50,
ssDnaHelixNeighborDist=50)
fineModel.simulate( out, outputDirectory='output', numSteps=100000,
arbd=arbd, gpu=gpu)
fineModelCoords = readArbdCoords( "output/%s.0.restart" % out )
out = outPrefix + "-4"
finerModel.backmap( fineModel, fineModelCoords,
dsDnaHelixNeighborDist=100, dsDnaAllNeighborDist=50,
ssDnaHelixNeighborDist=50)
finerModel.simulate(out, outputDirectory='output', numSteps=1000000, timestep=50e-6,
arbd=arbd, gpu=gpu)
finerModelCoords = readArbdCoords( "output/%s.0.restart" % out )
stepId = [i for i in range(5)]
numSteps = [int(0.2*numCoarsestSteps), int(0.8*numCoarsestSteps)]
timestep = [200e-6, 200e-6]
lastModel = None
for i in stepId:
args = dict(arbd=arbd, gpu=gpu, numSteps=10000,
outputDirectory = 'output')
if i < len(numSteps): args['numSteps'] = numSteps[i]
if i < len(numSteps): args['timestep'] = timestep[i]
if i == 4:
args['numSteps'] = 1000000
args['timestep'] = timestep=50e-6
out = "%s-%d" % (outPrefix,i)
lastOut = "%s-%d" % (outPrefix,i-1)
model = models[i]
model.writePsf(out + '.ideal.psf')
model.writePdb(out + '.ideal.pdb')
if lastModel is not None:
if i == 1:
coords = readArbdCoords( "output/%s.0.restart" % lastOut )
else:
coords = readAvgArbdCoords("%s.psf" % out, "%s.ideal.pdb" % lastOut,
"output/%s.0.dcd" % lastOut )
model.backmap( lastModel, coords,
dsDnaHelixNeighborDist=100, dsDnaAllNeighborDist=50,
ssDnaHelixNeighborDist=50)
if i > 0: model.simulate(out, **args)
lastModel = model
coords = readAvgArbdCoords(
"%s.psf" % out, "%s.pdb" % out,
"output/%s.0.dcd" % out )
out = outPrefix + "-5"
try:
......@@ -188,7 +213,7 @@ def simulateDesign(inFile, outPrefix, numCoarsestSteps=5000000, gpu=0, shiftEule
finestModel.randomizeUnassignedNtSeqs()
# finestModel.setUnassignedNtSeqs("T")
finestModel.writePdbPsf( out + ".ideal" )
finestModel.backmap( finerModel, finerModelCoords )
finestModel.backmap( lastModel, coords )
finestModel.simulate( out, numSteps=480000, gpus=gpu, numprocs=4, namd=namd )
if len(sys.argv) > 1:
......
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