Commit a5353000 authored by cmaffeo2's avatar cmaffeo2
Browse files

Merge branch 'master' into dev

parents a8ed3cac a7f370bb
......@@ -2,3 +2,4 @@ include mrdna/RELEASE-VERSION
include mrdna/README.md
include mrdna/LICENSE
recursive-include mrdna/resources *
recursive-include mrdna/arbdmodel/resources *
......@@ -38,9 +38,9 @@ cd ../../
# Setup python dependencies (you may want to set up a virtual environement)
conda install numpy>=1.14 scipy>=1.1 appdirs >= 1.4
conda install -c conda-forge mdanalysis>=0.18
pip install cadnano >= 2.5.2.1
conda install "numpy>=1.14" "scipy>=1.1" "appdirs>=1.4"
conda install -c conda-forge "mdanalysis>=0.18"
pip install "cadnano>=2.5.2.1"
git clone https://gitlab.engr.illinois.edu/tbgl/tools/mrdna
cd mrdna
......
......@@ -14,9 +14,25 @@ parser.add_argument('-o','--output-prefix', type=str, default=None,
help="Name for your job's output")
parser.add_argument('-d','--directory', type=str, default=None,
help='Directory for simulation; does not need to exist yet')
parser.add_argument('--sequence-file', type=str, default=None,
help='Sequence of longest strand.')
parser.add_argument('--fill-sequence', choices=['A','T','C','G','random'], default='T',
help='Fill sequence where otherwise unset.')
parser.add_argument('--crossover-to-intrahelical-cutoff', type=float, default=-1,
help='Set the distance (in Angstroms) beyond which crossovers between helix ends are converted to intrahelical connections; a negative value means no crossovers will be converted')
parser.add_argument('--backbone-scale', type=float, default=1.0,
help='Factor to scale DNA backbone in atomic model; try 0.25 to avoid clashes for atomistic simulations')
parser.add_argument('--enrg-md-steps', type=float, default=1e6,
help='Number of ENRG-MD steps')
parser.add_argument('--write-pqr', action='store_true',
help='Write PQR file in addition PDB?')
parser.add_argument('input_file', type=str,
help="""Any of the following:
(1) a cadnano JSON file;
......@@ -55,6 +71,15 @@ if __name__ == '__main__':
model = read_model( str(infile) )
if args.crossover_to_intrahelical_cutoff >= 0:
model.convert_crossovers_at_ends_beyond_cutoff_to_intrahelical(
args.crossover_to_intrahelical_cutoff)
if args.sequence_file is not None:
from mrdna.model.dna_sequence import read_sequence_file
seq = read_sequence_file(args.sequence_file)
model.set_sequence( seq, args.fill_sequence )
d_orig = os.getcwd()
try:
if directory is None:
......@@ -78,7 +103,7 @@ http://dx.doi.org/10.1093/nar/gkw155
os.makedirs('output')
except FileExistsError:
...
model.atomic_simulate( output_name = prefix )
model.atomic_simulate( output_name = prefix, num_steps=args.enrg_md_steps, write_pqr = args.write_pqr, dry_run=True )
except:
raise
......
......@@ -9,6 +9,7 @@ import re
import pathlib
from mrdna import __version__ as __version__
from mrdna.simulate import multiresolution_simulation as simulate
import numpy as np
"""Easy multiresolution simulations of DNA nanotechnology objects using ARBD"""
......@@ -29,9 +30,15 @@ parser.add_argument('--debye-length', type=float, default=None,
parser.add_argument('--temperature', type=float, default=295,
help='Temperature in Kelvin.')
parser.add_argument('--dimensions', type=str, default=None,
help='Comma-separated sides of box in Angstroms.')
parser.add_argument('--sequence-file', type=str, default=None,
help='Sequence of longest strand.')
parser.add_argument('--fill-sequence', choices=['A','T','C','G','random'], default='T',
help='Fill sequence where otherwise unset.')
parser.add_argument('--output-period', type=float, default=1e4,
help='Simulation steps between DCD frames')
# parser.add_argument('--minimization-steps', type=float, default=0,
......@@ -49,18 +56,37 @@ parser.add_argument('--oxdna-steps', type=float, default=None,
help='Perform an oxDNA simulation instead of creating atomic model')
parser.add_argument('--oxdna-output-period', type=float, default=1e4,
help='Simulation steps between oxDNA configuration and energy output')
parser.add_argument('--coarse-bond-cutoff', type=float, default = 0,
help='Ignore bonds beyond this cutoff during 25%% of coarse simulation, adding a step to the simulation. Default value of zero implies bonds are not ignored')
parser.add_argument('--coarse-persistence-length', type=float, default = None,
help='Specify a custom persistence length (in nm) for the first 25%% of coarse simulation (after coarse-bond-cutoff simulation if specified), adding a step to the simulation')
parser.add_argument('--crossover-to-intrahelical-cutoff', type=float, default=-1,
help='Set the distance (in Angstroms) beyond which crossovers between helix ends are converted to intrahelical connections; a negative value means no crossovers will be converted')
parser.add_argument('--backbone-scale', type=float, default=1.0,
help='Factor to scale DNA backbone in atomic model; try 0.25 to avoid clashes for atomistic simulations')
parser.add_argument('--run-enrg-md', action='store_true',
help='Perform the ENRG-MD simulation?')
parser.add_argument('--enrg-md-steps', type=float, default=1e6,
help='Number of ENRG-MD steps')
parser.add_argument('--debug', action='store_true',
help='Run through the python debugger?')
parser.add_argument('--write-pqr', action='store_true',
help='Write PQR file in addition PDB?')
parser.add_argument('--draw-cylinders', action='store_true',
help='Whether or not to draw the cylinders')
parser.add_argument('--draw-tubes', action='store_true',
help='Whether or not to draw the tubes')
parser.add_argument('--hj-equilibrium-angle', type=float, default=0.0,
help='Specify an equilibrium dihedral angle for the Holliday junction angle potential; the default value works well for origami')
parser.add_argument('--version', action='version',
version='%(prog)s {}'.format(__version__),
help='Print the version of mrdna')
......@@ -73,6 +99,16 @@ parser.add_argument('input_file', type=str,
args = parser.parse_args()
if args.dimensions is not None:
try:
tmp = [float(f) for f in args.dimensions.split(',')]
if len(tmp) != 3:
raise ValueError
args.dimensions = tmp
except:
raise ValueError('--dimensions requires a comma separated list of floating point values, not "{}"'.format(args.dimensions))
def main():
infile = pathlib.Path(args.input_file)
try:
......@@ -87,9 +123,27 @@ def main():
from mrdna.readers import read_vhelix as read_model
elif extension == '.pdb':
from mrdna.readers import read_atomic_pdb as read_model
elif extension == '.dat' or extension == '.txt':
from mrdna.readers import read_list as read_model
else:
raise Exception("Unrecognized input file '{}'".format(infile))
model = read_model( str(infile), debye_length=args.debye_length, temperature=args.temperature )
model = read_model( str(infile), debye_length=args.debye_length, temperature=args.temperature, hj_equilibrium_angle = args.hj_equilibrium_angle )
if args.dimensions is None:
dim = [max(x,1000) for x in model.dimensions_from_structure( padding_factor=2.5 )]
model.dimensions = dim
model.origin = model.get_center()-0.5*np.array(model.dimensions)
else:
model.dimensions = args.dimensions
if args.crossover_to_intrahelical_cutoff >= 0:
model.convert_crossovers_at_ends_beyond_cutoff_to_intrahelical(
args.crossover_to_intrahelical_cutoff)
if args.sequence_file is not None:
from mrdna.model.dna_sequence import read_sequence_file
seq = read_sequence_file(args.sequence_file)
model.set_sequence( seq, args.fill_sequence )
if args.output_prefix is not None:
prefix = args.output_prefix
......@@ -105,7 +159,7 @@ def main():
model._clear_beads()
model._generate_atomic_model(scale=args.backbone_scale)
model.write_atomic_ENM( prefix )
model.atomic_simulate( output_name = prefix )
model.atomic_simulate( output_name = prefix, write_pqr=args.write_pqr, dry_run=True )
elif args.draw_tubes is True:
output_name = prefix + ".tubes.tcl"
model.vmd_tube_tcl( output_name )
......@@ -117,7 +171,7 @@ def main():
model._clear_beads()
model._generate_atomic_model(scale=args.backbone_scale)
model.write_atomic_ENM( prefix )
model.atomic_simulate( output_name = prefix )
model.atomic_simulate( output_name = prefix, write_pqr=args.write_pqr, dry_run=True )
else:
run_args = dict(
model = model,
......@@ -128,6 +182,8 @@ def main():
minimization_output_period = int(args.output_period),
coarse_local_twist = args.coarse_local_twist,
fix_linking_number = args.fix_linking_number,
bond_cutoff = args.coarse_bond_cutoff,
coarse_persistence_length = args.coarse_persistence_length,
coarse_output_period = int(args.output_period),
fine_output_period = int(args.output_period),
minimization_steps = 0, # int(args.minimization_steps),
......@@ -135,7 +191,10 @@ def main():
fine_steps = int(args.fine_steps),
backbone_scale = args.backbone_scale,
oxdna_steps = args.oxdna_steps,
oxdna_output_period = args.oxdna_output_period
oxdna_output_period = args.oxdna_output_period,
enrg_md_steps = args.enrg_md_steps,
write_pqr = args.write_pqr,
run_enrg_md = args.run_enrg_md
)
simulate( **run_args )
......
## Set up loggers
import logging
def _get_username():
import sys
try:
return sys.environ['USER']
except:
return None
logging.basicConfig(format='%(name)s: %(levelname)s: %(message)s')
logger = logging.getLogger(__name__)
logger.setLevel(logging.INFO)
_ch = logging.StreamHandler()
_ch.setFormatter(logging.Formatter('%(name)s: %(levelname)s: %(message)s'))
logger.addHandler(_ch)
devlogger = logging.getLogger(__name__+'.dev')
# devlogger.setLevel(logging.DEBUG)
if _get_username() not in ('cmaffeo2',):
devlogger.addHandler(logging.NullHandler())
## Import resources
from pathlib import Path
from .version import get_version
......@@ -8,6 +30,7 @@ _RESOURCE_DIR = Path(__file__).parent / 'resources'
def get_resource_path(relative_path):
return _RESOURCE_DIR / relative_path
## Import useful messsages
from .segmentmodel import SegmentModel, SingleStrandedSegment, DoubleStrandedSegment
from .simulate import multiresolution_simulation
from .model.dna_sequence import read_sequence_file
......@@ -15,3 +38,15 @@ from .model.dna_sequence import read_sequence_file
from .reporting import report
report("mrdna.__init__")
"""
_____ ___| \| | | _ |
| | _| | | | | | |
|_|_|_|_| |____/|_|___|__|__|
"""
print(""" _
_____ ___ _| |___ ___
| | _| . | | .'|
|_|_|_|_| |___|_|_|__,| v{}
""".format(__version__))
../LICENSE
\ No newline at end of file
../README.md
\ No newline at end of file
This diff is collapsed.
This diff is collapsed.
import numpy as np
from scipy.optimize import newton
def minimizeRmsd(coordsB, coordsA, weights=None, maxIter=100):
## Going through many iterations wasn't really needed
tol = 1
count = 0
R = np.eye(3)
comB = np.zeros([3,])
cNext = coordsB
while tol > 1e-6:
q,cB,comA = _minimizeRmsd(cNext,coordsA, weights)
R = R.dot(quaternion_to_matrix(q))
assert( np.all(np.isreal( R )) )
comB += cB
cLast = cNext
cNext = (coordsB-comB).dot(R)
tol = np.sum(((cNext-cLast)**2)[:]) / np.max(np.shape(coordsB))
if count > maxIter:
Exception("Exceeded maxIter (%d)" % maxIter)
count += 1
print("%d iterations",count)
return R, comB, comA
def minimizeRmsd(coordsB, coordsA, weights=None):
q,comA,comB = _minimizeRmsd(coordsB, coordsA, weights)
assert( np.all(np.isreal( q )) )
return quaternion_to_matrix(q),comA,comB
## http://onlinelibrary.wiley.com/doi/10.1002/jcc.21439/full
def _minimizeRmsd(coordsB, coordsA, weights=None):
A = coordsA
B = coordsB
shapeA,shapeB = [np.shape(X) for X in (A,B)]
for s in (shapeA,shapeB): assert( len(s) == 2 )
A,B = [X.T if s[1] > s[0] else X for X,s in zip([A,B],(shapeA,shapeB))] # TODO: print warning
shapeA,shapeB = [np.shape(X) for X in (A,B)]
assert( shapeA == shapeB )
for X,s in zip((A,B),(shapeA,shapeB)):
assert( s[1] == 3 and s[0] >= s[1] )
# if weights is None: weights = np.ones(len(A))
if weights is None:
comA,comB = [np.mean( X, axis=0 ) for X in (A,B)]
else:
assert( len(weights[:]) == len(B) )
W = np.diag(weights)
comA,comB = [np.sum( W.dot(X), axis=0 ) / np.sum(W) for X in (A,B)]
A = np.array( A-comA )
B = np.array( B-comB )
if weights is None:
s = A.T.dot(B)
else:
s = A.T.dot(W.dot(B))
sxx,sxy,sxz = s[0,:]
syx,syy,syz = s[1,:]
szx,szy,szz = s[2,:]
K = [[ sxx+syy+szz, syz-szy, szx-sxz, sxy-syx],
[syz-szy, sxx-syy-szz, sxy+syx, sxz+szx],
[szx-sxz, sxy+syx, -sxx+syy-szz, syz+szy],
[sxy-syx, sxz+szx, syz+szy, -sxx-syy+szz]]
K = np.array(K)
# GA = np.trace( A.T.dot(W.dot(A)) )
# GB = np.trace( B.T.dot(W.dot(B)) )
## Finding GA/GB can be done more quickly
# I = np.eye(4)
# x0 = (GA+GB)*0.5
# vals = newtoon(lambda x: np.det(K-x*I), x0 = x0)
vals, vecs = np.linalg.eig(K)
i = np.argmax(vals)
q = vecs[:,i]
# RMSD = np.sqrt( (GA+GB-2*vals[i]) / len(A) )
# print("CHECK:", K.dot(q)-vals[i]*q )
return q, comB, comA
def quaternion_to_matrix(q):
assert(len(q) == 4)
## It looks like the wikipedia article I used employed a less common convention for q (see below
## http://en.wikipedia.org/wiki/Rotation_formalisms_in_three_dimensions#Rotation_matrix_.E2.86.94_quaternion
# q1,q2,q3,q4 = q
# R = [[1-2*(q2*q2 + q3*q3), 2*(q1*q2 - q3*q4), 2*(q1*q3 + q2*q4)],
# [ 2*(q1*q2 + q3*q4), 1-2*(q1*q1 + q3*q3), 2*(q2*q3 - q1*q4)],
# [ 2*(q1*q3 - q2*q4), 2*(q1*q4 + q2*q3), 1-2*(q2*q2 + q1*q1)]]
q = q / np.linalg.norm(q)
q0,q1,q2,q3 = q
R = [[1-2*(q2*q2 + q3*q3), 2*(q1*q2 - q3*q0), 2*(q1*q3 + q2*q0)],
[ 2*(q1*q2 + q3*q0), 1-2*(q1*q1 + q3*q3), 2*(q2*q3 - q1*q0)],
[ 2*(q1*q3 - q2*q0), 2*(q1*q0 + q2*q3), 1-2*(q2*q2 + q1*q1)]]
return np.array(R)
def quaternion_from_matrix( R ):
e1 = R[0]
e2 = R[1]
e3 = R[2]
# d1 = 0.5 * np.sqrt( 1+R[0,0]+R[1,1]+R[2,2] )
# d2 = 0.5 * np.sqrt( 1+R[0,0]-R[1,1]-R[2,2] )
# d2 = 0.5 * np.sqrt( 1+R[0,0]-R[1,1]-R[2,2] )
# d2 = 0.5 * np.sqrt( 1+R[0,0]-R[1,1]-R[2,2] )
d1 = 1+R[0,0]+R[1,1]+R[2,2]
d2 = 1+R[0,0]-R[1,1]-R[2,2]
d3 = 1-R[0,0]+R[1,1]-R[2,2]
d4 = 1-R[0,0]-R[1,1]+R[2,2]
maxD = max((d1,d2,d3,d4))
d = 0.5 / np.sqrt(maxD)
if d1 == maxD:
return np.array(( 1.0/(4*d),
d * (R[2,1]-R[1,2]),
d * (R[0,2]-R[2,0]),
d * (R[1,0]-R[0,1]) ))
elif d2 == maxD:
return np.array(( d * (R[2,1]-R[1,2]),
1.0/(4*d),
d * (R[0,1]+R[1,0]),
d * (R[0,2]+R[2,0]) ))
elif d3 == maxD:
return np.array(( d * (R[0,2]-R[2,0]),
d * (R[0,1]+R[1,0]),
1.0/(4*d),
d * (R[1,2]+R[2,1]) ))
elif d4 == maxD:
return np.array(( d * (R[1,0]-R[0,1]),
d * (R[0,2]+R[2,0]),
d * (R[1,2]+R[2,1]),
1.0/(4*d) ))
def rotationAboutAxis(axis,angle, normalizeAxis=True):
if normalizeAxis: axis = axis / np.linalg.norm(axis)
angle = angle * 0.5 * np.pi/180
cos = np.cos( angle )
sin = np.sin( angle )
q = [cos] + [sin*x for x in axis]
return quaternion_to_matrix(q)
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)
def readAvgArbdCoords(psf,pdb,dcd,rmsdThreshold=3.5, first_frame=0, last_frame=-1, stride=1):
import MDAnalysis as mda
usel = mda.Universe(psf, dcd)
sel = usel.select_atoms("name D*")
# r0 = ref.xyz[0,ids,:]
ts = usel.trajectory[last_frame]
r0 = sel.positions
pos = []
for t in range(ts.frame-1,first_frame-1,-stride):
usel.trajectory[t]
R,comA,comB = minimizeRmsd(sel.positions,r0)
r = np.array( [(r-comA).dot(R)+comB for r in sel.positions] )
rmsd = np.mean( (r-r0)**2 )
r = np.array( [(r-comA).dot(R)+comB for r in usel.atoms.positions] )
pos.append( r )
if rmsd > rmsdThreshold**2:
break
t0=t+1
print( "Averaging coordinates in %s after frame %d" % (dcd, t0) )
pos = np.mean(pos, axis=0)
return pos
def unit_quat_conversions():
for axis in [[0,0,1],[1,1,1],[1,0,0],[-1,-2,0]]:
for angle in np.linspace(-180,180,10):
R = rotationAboutAxis(axis, angle)
R2 = quaternion_to_matrix( quaternion_from_matrix( R ) )
if not np.all( np.abs(R-R2) < 0.01 ):
import pdb
pdb.set_trace()
quaternion_to_matrix( quaternion_from_matrix( R ) )
if __name__ == "__main__":
unit_quat_conversions()
# -*- coding: utf-8 -*-
## Test with `python -m arbdmodel.hps_polymer_model`
import numpy as np
import sys
## Local imports
from . import ArbdModel, ParticleType, PointParticle, Group, get_resource_path
from .abstract_polymer import PolymerSection, AbstractPolymerGroup
from .interactions import NonbondedScheme, HarmonicBond, HarmonicPotential
from .coords import quaternion_to_matrix
"""Define particle types"""
type_ = ParticleType("X",
damping_coefficient = 40.9,
mass=120
)
## Bonded potentials
class FjcNonbonded(NonbondedScheme):
def __init__(self, resolution=0.1, rMin=0):
NonbondedScheme.__init__(self, typesA=None, typesB=None, resolution=resolution, rMin=rMin)
def potential(self, r, typeA, typeB):
""" Constant force excluded volume """
force = 10 # kcal_mol/AA
radius = 6
u = np.zeros(r.shape)
# s = r < 2*radius
# u[s] = (2*radius - r[s]) * force
return u
class FjcBeadsFromPolymer(Group):
def __init__(self, polymer, sequence=None,
rest_length = 3.8, spring_constant = 25,
**kwargs):
# if sequence is None:
# raise NotImplementedError
# # ... set random sequence
self.polymer = polymer
self.sequence = sequence
self.spring_constant = 25
self.rest_length = 3.8
for prop in ('segname','chain'):
if prop not in kwargs:
# import pdb
# pdb.set_trace()
try:
self.__dict__[prop] = polymer.__dict__[prop]
except:
pass
# if len(sequence) != polymer.num_monomers:
# raise ValueError("Length of sequence does not match length of polymer")
Group.__init__(self, **kwargs)
def _clear_beads(self):
...
def _generate_beads(self):
# beads = self.children
nb = self.polymer.num_monomers
for i in range(nb):
c = self.polymer.monomer_index_to_contour(i)
r = self.polymer.contour_to_position(c)
bead = PointParticle(type_, r,
resid = i+1)
self.add(bead)
## Two consecutive nts
for i in range(len(self.children)-1):
b1 = self.children[i]
b2 = self.children[i+1]
""" units "10 kJ/N_A" kcal_mol """
bond = HarmonicBond(k = self.spring_constant,
r0 = self.rest_length,
rRange = (0,500),
resolution = 0.01,
maxForce = 50)
self.add_bond( i=b1, j=b2, bond = bond, exclude=True )
class FjcModel(ArbdModel):
def __init__(self, polymers,
sequences = None,
rest_length = 3.8,
spring_constant = 25,
damping_coefficient = 40.9,
DEBUG=False,
**kwargs):
"""
[damping_coefficient]: ns
"""
print("WARNING: diffusion coefficient arbitrarily set to 100 AA**2/ns in FjcModel")
kwargs['timestep'] = 50e-6
kwargs['cutoff'] = 10
if 'decompPeriod' not in kwargs:
kwargs['decompPeriod'] = 100000
""" Assign sequences """
if sequences is None:
# raise NotImplementedError("HpsModel must be provided a sequences argument")
sequences = [None for i in range(len(polymers))]
self.polymer_group = AbstractPolymerGroup(polymers)
self.sequences = sequences
self.rest_length = rest_length
self.spring_constant = spring_constant
ArbdModel.__init__(self, [], **kwargs)
""" Update type diffusion coefficients """