mrdna 10.5 KB
Newer Older
1
2
#! /usr/bin/env python
# -*- coding: utf-8 -*-
3
4
5
6
import sys
if sys.version_info[0] != 3:
    raise Exception("This script must be run with python3")
    
7
8
import argparse
import re
9
import pathlib
cmaffeo2's avatar
cmaffeo2 committed
10
from mrdna import __version__ as __version__
cmaffeo2's avatar
cmaffeo2 committed
11
from mrdna.simulate import multiresolution_simulation as simulate
12
import numpy as np
13

14
15
"""Easy multiresolution simulations of DNA nanotechnology objects using ARBD"""

cmaffeo2's avatar
cmaffeo2 committed
16
parser = argparse.ArgumentParser(prog="mrdna",
17
18
19
				 description=__doc__,
                                 formatter_class=argparse.ArgumentDefaultsHelpFormatter,
)
20
21
22
23
24
25
26

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('-g','--gpu', type=int, default=0,
                    help='GPU for simulation; check nvidia-smi for availability')
27
28
29
parser.add_argument('--debye-length', type=float, default=None,
                    help='Adjust the electrostatic repulsion matching osmotic pressure data from Rau and Parsegian under 25 mM MgCl2 condition with a Debye-Hueckel correction from the default 11.1 Angstroms.')

cmaffeo2's avatar
cmaffeo2 committed
30
31
32
parser.add_argument('--temperature', type=float, default=295,
                    help='Temperature in Kelvin.')

33
34
35
parser.add_argument('--dimensions', type=str, default=None,
                    help='Comma-separated sides of box in Angstroms.')

cmaffeo2's avatar
cmaffeo2 committed
36
37
38
parser.add_argument('--sequence-file', type=str, default=None,
                    help='Sequence of longest strand.')

39
40
41
parser.add_argument('--fill-sequence', choices=['A','T','C','G','random'], default='T',
                    help='Fill sequence where otherwise unset.')

42
43
parser.add_argument('--output-period', type=float, default=1e4,
                    help='Simulation steps between DCD frames')
44
45
46
47
# parser.add_argument('--minimization-steps', type=float, default=0,
#                     help='Simulation steps between DCD frames')
parser.add_argument('--coarse-local-twist', action='store_true',
                    help='Use local twist for coarsest representation?')
48
49
parser.add_argument('--fix-linking-number', action='store_true',
                    help='Fix the linking number so that it cannot change once local twist is introduced?')
50
parser.add_argument('--coarse-steps', type=float, default=1e7,
51
                    help='Simulation steps for coarse model (200 fs/step)')
52
parser.add_argument('--fine-steps', type=float, default=1e7,
53
                    help='Simulation steps for fine model (50 fs/step)')
54
55
56
57
58

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')
59
60
parser.add_argument('--coarse-bond-cutoff', type=float, default = 0,
                    help='Ignore bonds beyond this cutoff during first step of simulation; a value of 0 implies bonds are not ignored')
61

62
63
64
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')

65
66
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')
cmaffeo2's avatar
cmaffeo2 committed
67
68
69
70

parser.add_argument('--run-enrg-md', action='store_true',
                    help='Perform the ENRG-MD simulation?')

71
72
73
parser.add_argument('--enrg-md-steps', type=int, default=1e6,
                    help='Number of ENRG-MD steps')

74
75
parser.add_argument('--debug', action='store_true',
                    help='Run through the python debugger?')
76

cmaffeo2's avatar
cmaffeo2 committed
77
78
79
parser.add_argument('--write-pqr', action='store_true',
                    help='Write PQR file in addition PDB?')

80
parser.add_argument('--draw-cylinders', action='store_true',
cmaffeo2's avatar
cmaffeo2 committed
81
                    help='Whether or not to draw the cylinders')
82
parser.add_argument('--draw-tubes', action='store_true',
cmaffeo2's avatar
cmaffeo2 committed
83
                    help='Whether or not to draw the tubes')
84

85
86
87
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')

cmaffeo2's avatar
cmaffeo2 committed
88
89
90
91
parser.add_argument('--version', action='version',
                    version='%(prog)s {}'.format(__version__),
                    help='Print the version of mrdna')

92
93
94
95
96
97
98
99
parser.add_argument('input_file', type=str,
                    help="""Any of the following:
                    (1) a cadnano JSON file;
                    (2) a vHelix Maya (.ma) file;
                    (3) an atomic PDB file""")

args = parser.parse_args()

100
101
102
103
104
105
106
107
108
109
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))


110
def main():
111
    infile = pathlib.Path(args.input_file)
112
    try:
113
        prefix = infile.stem
114
        extension = infile.suffix
115
116
117
    except:
        raise Exception("Unrecognized input file '{}'".format(infile))
        
118
    if extension == '.json':
cmaffeo2's avatar
cmaffeo2 committed
119
        from mrdna.readers import read_cadnano as read_model
120
    elif extension == '.ma':
cmaffeo2's avatar
cmaffeo2 committed
121
        from mrdna.readers import read_vhelix as read_model
122
    elif extension == '.pdb':
cmaffeo2's avatar
cmaffeo2 committed
123
        from mrdna.readers import read_atomic_pdb as read_model
124
125
    elif extension == '.dat' or extension == '.txt':
        from mrdna.readers import read_list as read_model
126
127
    else:
        raise Exception("Unrecognized input file '{}'".format(infile))
128
    model = read_model( str(infile), debye_length=args.debye_length, temperature=args.temperature, hj_equilibrium_angle = args.hj_equilibrium_angle )
129

130
    if args.dimensions is None:
131
        dim = [max(x,1000) for x in model.dimensions_from_structure( padding_factor=2.5 )]
132
133
        model.dimensions = dim
        model.origin = model.get_center()-0.5*np.array(model.dimensions)
134
135
136
    else:
        model.dimensions = args.dimensions

137
138
139
140
141
142
143
    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)
144
        model.set_sequence( seq, args.fill_sequence )
145

146
147
    if args.output_prefix is not None:
        prefix = args.output_prefix
148
    
cmaffeo2's avatar
cmaffeo2 committed
149
150
151
152
153
154
155
156
    if args.draw_cylinders is True:
        output_name = prefix + ".cylinders.tcl"
        model.vmd_cylinder_tcl( output_name )
        if prefix == infile.stem and extension == ".pdb":
            print("Are you sure you want to overwrite '{}'? (you may provide '-o' option) [y/N]".format(infile))
            if input() not in ("y","Y","yes","Y","Yes","YES"):
                import sys
                sys.exit(1)
157
        model._clear_beads()
cmaffeo2's avatar
cmaffeo2 committed
158
159
        model._generate_atomic_model(scale=args.backbone_scale)            
        model.write_atomic_ENM( prefix )
cmaffeo2's avatar
cmaffeo2 committed
160
        model.atomic_simulate( output_name = prefix, write_pqr=args.write_pqr, dry_run=True )
cmaffeo2's avatar
cmaffeo2 committed
161
162
163
164
165
166
167
168
169
170
    elif args.draw_tubes is True:
        output_name = prefix + ".tubes.tcl"
        model.vmd_tube_tcl( output_name )
        if prefix == infile.stem and extension == ".pdb":
            print("Are you sure you want to overwrite '{}'? (you may provide '-o' option) [y/N]".format(infile))
            if input() not in ("y","Y","yes","Y","Yes","YES"):
                import sys
                sys.exit(1)
        model._clear_beads()
        model._generate_atomic_model(scale=args.backbone_scale)            
171
        model.write_atomic_ENM( prefix )
cmaffeo2's avatar
cmaffeo2 committed
172
        model.atomic_simulate( output_name = prefix, write_pqr=args.write_pqr, dry_run=True )
173
174
175
176
177
178
179
    else:
        run_args = dict(
            model = model,
            output_name = prefix,
            job_id = "job-" + prefix,
            directory = args.directory,
            gpu = args.gpu,
180
181
            minimization_output_period = int(args.output_period),
            coarse_local_twist = args.coarse_local_twist,
182
            fix_linking_number = args.fix_linking_number,
183
            bond_cutoff = args.coarse_bond_cutoff,
184
185
            coarse_output_period = int(args.output_period),
            fine_output_period = int(args.output_period),
186
            minimization_steps = 0, # int(args.minimization_steps),
187
188
            coarse_steps = int(args.coarse_steps),
            fine_steps = int(args.fine_steps),
189
190
            backbone_scale = args.backbone_scale,
            oxdna_steps = args.oxdna_steps,
cmaffeo2's avatar
cmaffeo2 committed
191
            oxdna_output_period = args.oxdna_output_period,
192
            enrg_md_steps = args.enrg_md_steps,
cmaffeo2's avatar
cmaffeo2 committed
193
            write_pqr = args.write_pqr,
cmaffeo2's avatar
cmaffeo2 committed
194
            run_enrg_md = args.run_enrg_md
195
        )
196

197
        simulate( **run_args )
198
199

if __name__ == '__main__':
200

201
202
203
204
205
206
207
208
209
210
211
212
213
214
    if args.debug:

        class Restart(Exception):
            """Causes a debugger to be restarted for the debugged python program."""
            pass

        import traceback
        import sys
        from pdb import Pdb
        pdb = Pdb()

        while True:
            try:
                # pdb._runscript(mainpyfile)
215
                pdb._user_requested_quit = False
216
217
218
219
220
221
222
223
224
                pdb.runcall(main)
                if pdb._user_requested_quit:
                    break
                print("The program finished and will be restarted")
            except Restart:
                print("Restarting...")
                print("\t" + " ".join(args))
            except SystemExit:
                # In most cases SystemExit does not warrant a post-mortem session.
225
226
                # print("The program exited via sys.exit(). Exit status:", end=' ')
                print("The program exited via sys.exit(). Exit status:")
227
228
229
230
231
232
233
234
235
236
237
238
239
240
                print(sys.exc_info()[1])
            except SyntaxError:
                traceback.print_exc()
                sys.exit(1)
            except:
                traceback.print_exc()
                print("Uncaught exception. Entering post mortem debugging")
                print("Running 'cont' or 'step' will restart the program")
                t = sys.exc_info()[2]
                pdb.interaction(None, t)
                print("Post mortem debugger finished. mrdna will be restarted")

    else:
        main()