diff --git a/dnarbd/MANIFEST.in b/dnarbd/MANIFEST.in
new file mode 100644
index 0000000000000000000000000000000000000000..b5d192857288eb5cc0fe9ea4ba76874a14b01ac3
--- /dev/null
+++ b/dnarbd/MANIFEST.in
@@ -0,0 +1 @@
+include resources
\ No newline at end of file
diff --git a/dnarbd/__init__.py b/dnarbd/__init__.py
new file mode 100644
index 0000000000000000000000000000000000000000..493251d9e522b04a7c7e7eb21b704ed596db3fef
--- /dev/null
+++ b/dnarbd/__init__.py
@@ -0,0 +1,8 @@
+import os
+
+_RESOURCE_DIR = os.path.join(os.path.dirname(__file__), 'resources')
+def get_resource_path(relative_path):
+    return os.path.join(_RESOURCE_DIR, relative_path)
+
+from .segmentmodel import SegmentModel, SingleStrandedSegment, DoubleStrandedSegment
+from .simulate import multiresolution_simulation
diff --git a/dnarbd/coords.py b/dnarbd/coords.py
index 11833fac09e7faf771e6cd965522f908c0017928..70b6216d5e4a478da96a0bbb594d16d17278c0a8 100644
--- a/dnarbd/coords.py
+++ b/dnarbd/coords.py
@@ -1,7 +1,6 @@
 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
diff --git a/dnarbd/examples/run.py b/dnarbd/examples/run.py
index fd24ebc674b58885ca43087ae4d23536765312c0..6d5f77d6977c24be6e9b61bf5e5e97e4b8adb546 100644
--- a/dnarbd/examples/run.py
+++ b/dnarbd/examples/run.py
@@ -2,8 +2,9 @@
 import sys
 from glob import glob
 import re
-from cadnano_segments import read_json_file, read_model
-from simulate import run_simulation_protocol
+
+from dnarbd.readers import cadnano_reader()
+from dnarbd.simulate import run_simulation_protocol
 
 if __name__ == '__main__':
     if len(sys.argv) > 1:
diff --git a/dnarbd/model/CanonicalNucleotideAtoms.py b/dnarbd/model/CanonicalNucleotideAtoms.py
index f2eb3aebf9174757a54b22733d27503291d2dc23..39613912c421b93406767ae57f063bc18ce33190 100644
--- a/dnarbd/model/CanonicalNucleotideAtoms.py
+++ b/dnarbd/model/CanonicalNucleotideAtoms.py
@@ -1,8 +1,9 @@
 import json
 import numpy as np
 import copy
-from  coords import rotationAboutAxis
-from arbdmodel import Group, PointParticle, ParticleType
+from ..coords import rotationAboutAxis
+from .arbdmodel import Group, PointParticle, ParticleType
+from .. import get_resource_path
 
 seqComplement = dict(A='T',G='C')
 resnames = dict(A="ADE", G="GUA", C="CYT", T="THY")
@@ -130,7 +131,7 @@ canonicalNtFwd = dict()
 direction = 'fwd'
 for seq in seqComplement.keys():
     for keystring,suff in zip(("5%s","%s","%s3","5%s3"),("-5prime","","-3prime","-singlet")):
-        prefix = 'resources/%s-%s%s' % (seq,direction,suff)
+        prefix = get_resource_path('%s-%s%s' % (seq,direction,suff))
         key = keystring % seq
         canonicalNtFwd[key] = CanonicalNucleotideFactory( prefix, seq )
 
@@ -138,13 +139,13 @@ canonicalNtRev = dict()
 direction = 'rev'
 for seq in seqComplement.keys():
     for keystring,suff in zip(("5%s","%s","%s3","5%s3"),("-5prime","","-3prime","-singlet")):
-        prefix = 'resources/%s-%s%s' % (seq,direction,suff)
+        prefix = get_resource_path('%s-%s%s' % (seq,direction,suff))
         key = keystring % seq
         canonicalNtRev[key] = CanonicalNucleotideFactory( prefix, seq )
 
-with open("resources/enm-template-honeycomb.json") as ch:
+with open(get_resource_path("enm-template-honeycomb.json")) as ch:
     enmTemplateHC = json.load(ch)
-with open("resources/enm-template-square.json") as ch:
+with open(get_resource_path("enm-template-square.json")) as ch:
     enmTemplateSQ = json.load(ch)
-with open("resources/enm-corrections-honecomb.json") as ch:
+with open(get_resource_path("enm-corrections-honecomb.json")) as ch:
     enmCorrectionsHC = json.load(ch)
diff --git a/dnarbd/model/dna_sequence.py b/dnarbd/model/dna_sequence.py
new file mode 100644
index 0000000000000000000000000000000000000000..a189f114ed1421f13dd83fa1a9c2632f36c9f18b
--- /dev/null
+++ b/dnarbd/model/dna_sequence.py
@@ -0,0 +1,13 @@
+from .. import get_resource_path
+
+_m13_path = get_resource_path("cadnano2pdb.seq.m13mp18.dat")
+def read_sequence_file(sequenceFile=_m13_path):
+    seq = []
+    with open(sequenceFile) as ch:
+        for l in ch:
+            l = l.strip().replace(" ", "")
+            if l[0] in (";","#"): continue
+            seq.extend([c.upper() for c in l])
+    return seq
+
+m13 = read_sequence_file()
diff --git a/dnarbd/model/nbPot.py b/dnarbd/model/nbPot.py
index 20757abbcdeed481c23632e2001a0e7ed75fd17b..1b9c145b4c229a76d32eae8b9e9383b9e0db08cd 100644
--- a/dnarbd/model/nbPot.py
+++ b/dnarbd/model/nbPot.py
@@ -3,7 +3,8 @@ import scipy.optimize as opt
 from scipy.interpolate import interp1d
 from scipy.signal import savgol_filter as savgol
 
-from nonbonded import NonbondedScheme
+from .nonbonded import NonbondedScheme
+from .. import get_resource_path
 
 dxParam = -1
 
@@ -25,7 +26,7 @@ def maxForce(x,u,maxForce=40):
     return u
 
 # d = np.loadtxt("resources/jj-nar-dna-pmf-20Mg-200Na.dat")
-d = np.loadtxt("resources/jj-nar-dna-pmf-100Mg.dat")
+d = np.loadtxt(get_resource_path("jj-nar-dna-pmf-100Mg.dat"))
 
 x0 = d[:,1] + dxParam
 y0 = d[:,0]                     # kcal/mol for 10bp with 10bp, according to JY
diff --git a/dnarbd/readers/__init__.py b/dnarbd/readers/__init__.py
new file mode 100644
index 0000000000000000000000000000000000000000..b846d7228dc75cef0c6c38e2904686d12079d97c
--- /dev/null
+++ b/dnarbd/readers/__init__.py
@@ -0,0 +1,25 @@
+from ..segmentmodel import SegmentModel
+
+""" cadnano """
+from .cadnano_segments import read_json_file
+from .cadnano_segments import read_model as model_from_cadnano_json
+
+""" vHelix """
+from .polygon_mesh import parse_maya_file, convert_maya_to_segments
+
+""" Atomic PDB """
+from .segmentmodel_from_pdb import SegmentModelFromPdb
+
+## TODO: make module this package conform to a single style for input/output
+
+def read_cadnano(json_file, **model_parameters):
+    data = read_json_file
+    return model_from_cadnano_json(data)
+
+def read_vhelix(maya_file, **model_parameters):
+    data = parse_maya_file(maya_file)
+    segments, dsSegmentDict = convert_maya_to_segments( data )
+    return SegmentModel( segments,**model_parameters )
+
+def read_atomic_pdb(pdb_file, **model_parameters):
+    return SegmentModelFromPdb(pdb_file)
diff --git a/dnarbd/readers/cadnano_segments.py b/dnarbd/readers/cadnano_segments.py
index 4bb7b062892b89146e7ad1054a0008b7ca0ef5a2..f2a4836bb1df13405a9e26bcdcc855c3061b4ec2 100644
--- a/dnarbd/readers/cadnano_segments.py
+++ b/dnarbd/readers/cadnano_segments.py
@@ -5,22 +5,9 @@ import os,sys
 from glob import glob
 import re
 
-from coords import readArbdCoords, readAvgArbdCoords
-from segmentmodel import SegmentModel, SingleStrandedSegment, DoubleStrandedSegment
-
-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"
-
-def readSequenceFile(sequenceFile='resources/cadnano2pdb.seq.m13mp18.dat'):
-    seq = []
-    with open(sequenceFile) as ch:
-        for l in ch:
-            l = l.strip().replace(" ", "")
-            if l[0] in (";","#"): continue
-            seq.extend([c.upper() for c in l])
-    return seq
-
-m13seq = readSequenceFile(sequenceFile='resources/cadnano2pdb.seq.m13mp18.dat')
+from ..coords import readArbdCoords, readAvgArbdCoords
+from ..segmentmodel import SegmentModel, SingleStrandedSegment, DoubleStrandedSegment
+from ..model.dna_sequence import m13 as m13seq
 
 ## TODO: separate SegmentModel from ArbdModel so multiple parts can be combined
 ## TODO: catch circular strands in "get_5prime" cadnano calls
diff --git a/dnarbd/readers/polygon_mesh.py b/dnarbd/readers/polygon_mesh.py
index c7e5f23b4452d4bff56f39a077148eecb73e6796..15c4c88a097951e264b75dd47f1bb2608f57d878 100644
--- a/dnarbd/readers/polygon_mesh.py
+++ b/dnarbd/readers/polygon_mesh.py
@@ -1,9 +1,9 @@
 import numpy as np
 import sys
 import re
-from coords import rotationAboutAxis
+from ..coords import rotationAboutAxis
 
-from segmentmodel import SegmentModel, SingleStrandedSegment, DoubleStrandedSegment
+from ..segmentmodel import SegmentModel, SingleStrandedSegment, DoubleStrandedSegment
 
 class MayaObj():
     def __init__(self, maya_lines):
diff --git a/dnarbd/readers/segmentmodel_from_pdb.py b/dnarbd/readers/segmentmodel_from_pdb.py
index cb3c9e4ff489a3e26c8b912316007dc0ad664e01..cc9708d85926126c0cc75159d518425bce803d7d 100644
--- a/dnarbd/readers/segmentmodel_from_pdb.py
+++ b/dnarbd/readers/segmentmodel_from_pdb.py
@@ -10,21 +10,8 @@ import MDAnalysis as mda
 import MDAnalysis.analysis.align as align
 from MDAnalysis.analysis.distances import distance_array
 
-from segmentmodel import SegmentModel, SingleStrandedSegment, DoubleStrandedSegment
-
-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"
-
-def readSequenceFile(sequenceFile='resources/cadnano2pdb.seq.m13mp18.dat'):
-    seq = []
-    with open(sequenceFile) as ch:
-        for l in ch:
-            l = l.strip().replace(" ", "")
-            if l[0] in (";","#"): continue
-            seq.extend([c.upper() for c in l])
-    return seq
-
-m13seq = readSequenceFile(sequenceFile='resources/cadnano2pdb.seq.m13mp18.dat')
+from ..segmentmodel import SegmentModel, SingleStrandedSegment, DoubleStrandedSegment
+from .. import get_resource_path
 
 resnames = dict(A='ADE DA DA5 DA3', T='THY DT DT5 DT3',
                 C='CYT DC DC5 DC3', G='GUA DG DG5 DG3')
@@ -34,7 +21,7 @@ bases = resnames.keys()
 resname_to_key = {n:k for k,v in resnames.items() for n in v.split()}
 complement = dict(A='T', C='G', T='A', G='C')
 
-refUnis = {b:mda.Universe("resources/generate-nts/1-w3dna/{}{}.pdb".format(b,complement[b]))
+refUnis = {b:mda.Universe(get_resource_path("generate-nts/1-w3dna/{}{}.pdb".format(b,complement[b])))
            for b in bases}
 
 stack_text = "not name O5' P O1P O2P OP1 OP2 O3' H*"
diff --git a/dnarbd/segmentmodel.py b/dnarbd/segmentmodel.py
index 88f152cc660f1bd95332ed310eacf04d71f2ebb2..657ab574da405a60aa895f191af151545b83b8b2 100644
--- a/dnarbd/segmentmodel.py
+++ b/dnarbd/segmentmodel.py
@@ -1,18 +1,18 @@
 import pdb
 import numpy as np
 import random
-from arbdmodel import PointParticle, ParticleType, Group, ArbdModel
-from coords import rotationAboutAxis, quaternion_from_matrix, quaternion_to_matrix
-from nonbonded import *
+from .model.arbdmodel import PointParticle, ParticleType, Group, ArbdModel
+from .coords import rotationAboutAxis, quaternion_from_matrix, quaternion_to_matrix
+from .model.nonbonded import *
 from copy import copy, deepcopy
-from nbPot import nbDnaScheme
+from .model.nbPot import nbDnaScheme
 
 from scipy.special import erf
 import scipy.optimize as opt
 from scipy import interpolate
 
-from CanonicalNucleotideAtoms import canonicalNtFwd, canonicalNtRev, seqComplement
-from CanonicalNucleotideAtoms import enmTemplateHC, enmTemplateSQ, enmCorrectionsHC
+from .model.CanonicalNucleotideAtoms import canonicalNtFwd, canonicalNtRev, seqComplement
+from .model.CanonicalNucleotideAtoms import enmTemplateHC, enmTemplateSQ, enmCorrectionsHC
 
 # import pdb
 """
diff --git a/dnarbd/simulate.py b/dnarbd/simulate.py
index 933db19e0cd13cdc92e25c5841df6862a9e01a3a..fe8e0b0d42ced786b76dcd48b11cd45678eb3781 100644
--- a/dnarbd/simulate.py
+++ b/dnarbd/simulate.py
@@ -1,22 +1,22 @@
 import os
-from coords import readArbdCoords, readAvgArbdCoords
+import tempfile
+from .coords import readArbdCoords, readAvgArbdCoords
 
 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"
 
-def run_simulation_protocol( output_name, job_id, model,
-                             remove_long_bonds=False,
-                             gpu = 0,
-                             directory=None,
-                             coarse_steps = 5e7+1,
-                             fine_steps = 5e7+1
-                         ):
+def multiresolution_simulation( output_name, job_id, model,
+                                remove_long_bonds=False,
+                                gpu = 0,
+                                directory=None,
+                                coarse_steps = 5e7+1,
+                                fine_steps = 5e7+1
+                            ):
 
     ret = None
     d_orig = os.getcwd()
     try:
         if directory is None:
-            import tempfile
             directory = tempfile.mkdtemp(prefix='origami-%s-' % job_id, dir='/dev/shm/')
         elif not os.path.exists(directory):
             os.makedirs(directory)