Skip to content

Top-Level Functions

These functions are the main entry points to MDNA, available directly from import mdna.

mdna.nucleic.make(sequence=None, control_points=None, circular=False, closed=False, n_bp=None, dLk=None, bp_per_turn=10.5)

Generate a DNA structure from a given DNA sequence and control points.

Parameters:

Name Type Description Default
sequence str

DNA sequence code. If not provided, the default sequence 'CGCGAATTCGCG' will be used. (default: None)

 None
control_points ndarray

Control points of the DNA structure. Should be a numpy array of shape (n, 3) where n is the number of control points. If not provided, a straight line will be used as the default control points. (default: None)

 None
circular bool

Flag indicating if the DNA structure is circular/closed. If True, the DNA structure will be closed. If False, the DNA structure will be open. (default: False)

 False
closed bool

Flag indicating if the DNA structure is closed. If True, the DNA structure will be closed. If False, the DNA structure will be open. This argument is deprecated and will be removed in a future version. Please use the 'circular' argument instead. (default: False)

 False
n_bp int

Number of base pairs to scale the shape with. If not provided, the number of base pairs will be determined based on the length of the control points or the sequence. (default: None)

 None
dLk int

Change in twist in terms of Linking number of the DNA structure. If not provided, a neutral twist based on bp_per_turn = 10.5 will be used. (default: None)

 None
bp_per_turn int

Number of base pairs per turn of the DNA structure. (default: 10.5)

 10.5

Returns:

Name Type Description
Nucleic object

DNA structure object.
Decision matrix
  • No control_points, sequence provided: uses default shape (line/circle from circular) and sets n_bp = len(sequence).
  • No control_points, n_bp provided: uses default shape and generates a random sequence of length n_bp.
  • control_points provided, sequence and n_bp both omitted: infers n_bp from spline sampling, then generates a random sequence of that length.
  • control_points + sequence: scales spline to len(sequence).
  • control_points + n_bp: scales spline to n_bp and generates a random sequence.
  • sequence + n_bp: both are accepted only when lengths match.
Additional behavior
  • For control_points-only input (sequence=None, n_bp=None), n_bp is inferred from the spline frame count (shape-dependent), not from the default dodecamer.
Validation rules
  • n_bp must be positive when provided.
  • control_points must have shape (n, 3) with n >= 4.
  • closed is deprecated and treated as an alias of circular.
Example

Generate a DNA structure from a sequence 

dna = make(sequence='CGCGAATTCGCG', control_points=None, circular=False, closed=False, n_bp=None, dLk=None)

Source code in mdna/nucleic.py 
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
def make(sequence: str = None, control_points: np.ndarray = None, circular : bool = False, closed: bool = False, n_bp : int = None, dLk : int = None, bp_per_turn : int = 10.5):
    """Generate a DNA structure from a given DNA sequence and control points.

    Args:
        sequence (str, optional): DNA sequence code. If not provided, the default sequence 'CGCGAATTCGCG' will be used. (default: None)
        control_points (ndarray, optional): Control points of the DNA structure. Should be a numpy array of shape (n, 3) where n is the number of control points. If not provided, a straight line will be used as the default control points. (default: None)
        circular (bool, optional): Flag indicating if the DNA structure is circular/closed. If True, the DNA structure will be closed. If False, the DNA structure will be open. (default: False)
        closed (bool, optional): Flag indicating if the DNA structure is closed. If True, the DNA structure will be closed. If False, the DNA structure will be open. This argument is deprecated and will be removed in a future version. Please use the 'circular' argument instead. (default: False)
        n_bp (int, optional): Number of base pairs to scale the shape with. If not provided, the number of base pairs will be determined based on the length of the control points or the sequence. (default: None)
        dLk (int, optional): Change in twist in terms of Linking number of the DNA structure. If not provided, a neutral twist based on bp_per_turn = 10.5 will be used. (default: None)
        bp_per_turn (int, optional): Number of base pairs per turn of the DNA structure. (default: 10.5)

    Returns:
        Nucleic (object): DNA structure object.

    Decision matrix:
        - No `control_points`, `sequence` provided: uses default shape (line/circle from `circular`) and sets `n_bp = len(sequence)`.
        - No `control_points`, `n_bp` provided: uses default shape and generates a random sequence of length `n_bp`.
        - `control_points` provided, `sequence` and `n_bp` both omitted: infers `n_bp` from spline sampling, then generates a random sequence of that length.
        - `control_points` + `sequence`: scales spline to `len(sequence)`.
        - `control_points` + `n_bp`: scales spline to `n_bp` and generates a random sequence.
        - `sequence` + `n_bp`: both are accepted only when lengths match.

    Additional behavior:
        - For `control_points`-only input (`sequence=None`, `n_bp=None`), `n_bp` is inferred from the spline frame count (shape-dependent), not from the default dodecamer.

    Validation rules:
        - `n_bp` must be positive when provided.
        - `control_points` must have shape `(n, 3)` with `n >= 4`.
        - `closed` is deprecated and treated as an alias of `circular`.

    Example:
        Generate a DNA structure from a sequence
        ```python
        dna = make(sequence='CGCGAATTCGCG', control_points=None, circular=False, closed=False, n_bp=None, dLk=None)
        ```
    """

    # Backward compatibility: keep supporting `closed`, but funnel semantics into `circular`.
    if closed:
        warnings.warn(
            "The 'closed' argument is deprecated and will be removed in a future release. Use 'circular' instead.",
            DeprecationWarning,
            stacklevel=2,
        )
        circular = circular or closed

    # Guard against invalid explicit base-pair counts early.
    if n_bp is not None and n_bp <= 0:
        raise ValueError('n_bp should be a positive integer')

    # If control points are provided without sequence/n_bp, infer n_bp from spline arc-length sampling.
    infer_n_bp_from_control_points = False

    if control_points is not None:
        # Normalize and validate control-point geometry input.
        control_points = np.asarray(control_points)
        if control_points.ndim != 2 or control_points.shape[1] != 3:
            raise ValueError('control_points should have shape (n, 3)')
        if len(control_points) < 4:
            raise ValueError('Control points should contain at least 4 points [x, y, z]')
        infer_n_bp_from_control_points = sequence is None and n_bp is None
    elif circular:
        # Default closed template when no custom shape is provided.
        control_points = Shapes.circle(radius=1)
    else:
        # Default open template when no custom shape is provided.
        control_points = Shapes.line(length=1)

    if infer_n_bp_from_control_points:
        # Build spline first to infer n_bp, then generate a matching random sequence.
        spline = SplineFrames(control_points=control_points, n_bp=None, closed=circular, dLk=dLk, bp_per_turn=bp_per_turn)
        sequence, n_bp = _check_input(sequence=None, n_bp=spline.n_bp)
    else:
        # Resolve sequence/n_bp first (including random/default generation), then scale spline accordingly.
        sequence, n_bp = _check_input(sequence=sequence, n_bp=n_bp)
        spline = SplineFrames(control_points=control_points, n_bp=n_bp, closed=circular, dLk=dLk, bp_per_turn=bp_per_turn)

    return Nucleic(sequence=sequence, n_bp=n_bp, frames=spline.frames, chainids=[0, 1], circular=circular)

mdna.nucleic.load(traj=None, frames=None, sequence=None, chainids=[0, 1], circular=None, filename=None, top=None, stride=None)

Load DNA representation from either base step mean reference frames/spline frames or an MDtraj trajectory.

Parameters:

Name Type Description Default
traj object

MDtraj trajectory containing the DNA structure. If provided, the frames and sequence arguments are ignored. (default: None)

 None
frames array

Base step mean reference frames of shape (n_bp, n_timesteps, 4, 3) or (n_bp, 4, 3). If provided, the traj and sequence arguments are ignored. (default: None)

 None
sequence str

DNA sequence. If provided, the traj and frames arguments are ignored. (default: None)

 None
chainids list

Chain IDs of the DNA structure. (default: [0,1])

 [0, 1]
circular bool

Flag indicating if the DNA structure is circular/closed. If not provided, it will be determined based on the input data. (default: None)

 None
filename str

The filename or filenames of the trajectory. If provided, the traj and frames arguments are ignored. (default: None)

 None
top str

The topology file of the trajectory. (default: None)

 None
stride int

The stride of the trajectory. (default: None)

 None

Returns:

Name Type Description
Nucleic object

DNA structure object.
Notes
  • filename is resolved first to an MDtraj trajectory (optionally using top and stride).
  • The resulting traj then takes precedence over frames and sequence when constructing Nucleic.
Example

Load a DNA structure from a trajectory 

traj = md.load('dna.pdb')
dna = mdna.load(traj=traj, chainids=[0, 1])

Source code in mdna/nucleic.py 
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
def load(traj=None, frames=None, sequence=None, chainids=[0,1], circular=None, filename=None, top=None, stride=None):
    """Load DNA representation from either base step mean reference frames/spline frames or an MDtraj trajectory.

    Args:
        traj (object, optional): MDtraj trajectory containing the DNA structure. If provided, the frames and sequence arguments are ignored. (default: None)
        frames (np.array, optional): Base step mean reference frames of shape (n_bp, n_timesteps, 4, 3) or (n_bp, 4, 3). If provided, the traj and sequence arguments are ignored. (default: None)
        sequence (str, optional): DNA sequence. If provided, the traj and frames arguments are ignored. (default: None)
        chainids (list, optional): Chain IDs of the DNA structure. (default: [0,1])
        circular (bool, optional): Flag indicating if the DNA structure is circular/closed. If not provided, it will be determined based on the input data. (default: None)
        filename (str, optional): The filename or filenames of the trajectory. If provided, the traj and frames arguments are ignored. (default: None)
        top (str, optional): The topology file of the trajectory. (default: None)
        stride (int, optional): The stride of the trajectory. (default: None)

    Returns:
        Nucleic (object): DNA structure object.

    Notes:
        - `filename` is resolved first to an MDtraj trajectory (optionally using `top` and `stride`).
        - The resulting `traj` then takes precedence over `frames` and `sequence` when constructing `Nucleic`.

    Example:
        Load a DNA structure from a trajectory
        ```python
        traj = md.load('dna.pdb')
        dna = mdna.load(traj=traj, chainids=[0, 1])
        ```
    """
    # Load the trajectory directly using MDtraj from a file
    if filename is not None and top is None:
        traj = md.load(filename_or_filenames=filename, stride=stride)
    elif filename is not None and top is not None:
        traj = md.load(filename_or_filenames=filename, top=top, stride=stride)

    return Nucleic(sequence=sequence, n_bp=None, traj=traj, frames=frames, chainids=chainids, circular=None)

mdna.nucleic.connect(Nucleic0, Nucleic1, sequence=None, n_bp=None, leader=0, frame=-1, margin=1, minimize=True, exvol_rad=0.0, temperature=300, control_points=None, index=0)

Connect two DNA structures by creating a new DNA structure with a connecting DNA strand.

The 3' end of the first DNA structure is connected to the 5' end of the second DNA structure. To connect the two strands, a straight line is interpolated between the two ends, and the optimal number of base pairs is distributed to achieve a neutral twist.

Parameters:

Name Type Description Default
Nucleic0 Nucleic

First DNA structure to connect.

 required
Nucleic1 Nucleic

Second DNA structure to connect.

 required
sequence str or List

DNA sequence of the connecting DNA strand. Default is None.

 None
n_bp int

Number of base pairs of the connecting DNA strand. Default is None.

 None
leader int

The leader of the DNA structure to connect. Default is 0.

 0
frame int

The time frame to connect. Default is -1.

 -1
margin int

Number of base pairs to fix at the end. Default is 1.

 1
minimize bool

Whether to minimize the new DNA structure. Default is True.

 True
exvol_rad float

Radius for excluded volume interactions during minimization. Default is 0.0.

 0.0
temperature int

Temperature for minimization. Default is 300.

 300

Returns:

Name Type Description
Nucleic object

DNA structure with the two DNA structures connected.

Raises:

Type Description
ValueError

If either of the DNA structures is circular.
Notes
  • The minimization does not use excluded volume interactions by default.This is because the excluded volume interactions require the EV beads to have no overlap. However, in the initial configuration, the EV beads are likely to have overlap. If desired, the resulting Nucleic object can be further minimized with the excluded volume interactions.
Example

Connect two DNA structures 

dna = connect(Nucleic0, Nucleic1, margin=5)

Source code in mdna/nucleic.py 
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
def connect(Nucleic0, Nucleic1, sequence: Union[str|List] = None, n_bp: int = None, leader: int = 0, frame: int = -1, margin: int = 1, minimize: bool = True, exvol_rad: float = 0.0, temperature: int = 300, control_points : np.ndarray = None, index : int = 0):  
    """Connect two DNA structures by creating a new DNA structure with a connecting DNA strand.

    The 3' end of the first DNA structure is connected to the 5' end of the second DNA structure.
    To connect the two strands, a straight line is interpolated between the two ends,
    and the optimal number of base pairs is distributed to achieve a neutral twist.

    Args:
        Nucleic0 (Nucleic): First DNA structure to connect.
        Nucleic1 (Nucleic): Second DNA structure to connect.
        sequence (str or List, optional): DNA sequence of the connecting DNA strand. Default is None.
        n_bp (int, optional): Number of base pairs of the connecting DNA strand. Default is None.
        leader (int, optional): The leader of the DNA structure to connect. Default is 0.
        frame (int, optional): The time frame to connect. Default is -1.
        margin (int, optional): Number of base pairs to fix at the end. Default is 1.
        minimize (bool, optional): Whether to minimize the new DNA structure. Default is True.
        exvol_rad (float, optional): Radius for excluded volume interactions during minimization. Default is 0.0.
        temperature (int, optional): Temperature for minimization. Default is 300.

    Returns:
        Nucleic (object): DNA structure with the two DNA structures connected.

    Raises:
        ValueError: If either of the DNA structures is circular.

    Notes:
        - The minimization does not use excluded volume interactions by default.This is because the excluded volume interactions require the EV beads to have no overlap. However, in the initial configuration, the EV beads are likely to have overlap. If desired, the resulting Nucleic object can be further minimized with the excluded volume interactions.

    Example:
        Connect two DNA structures
        ```python
        dna = connect(Nucleic0, Nucleic1, margin=5)
        ```
    """
    if Nucleic0.circular or Nucleic1.circular:
        raise ValueError('Cannot connect circular DNA structures')

    if (sequence is not None and n_bp is None) or (sequence is None and n_bp is not None) or (sequence is not None and n_bp is not None):
        sequence, n_bp = _check_input(sequence=sequence, n_bp=n_bp)

    # Connect the two DNA structures
    connector = Connector(Nucleic0, Nucleic1, sequence=sequence, n_bp=n_bp, leader=leader, frame=frame, margin=margin, control_points=control_points, index=index)
    if minimize:
        connector.connected_nuc.minimize(exvol_rad=exvol_rad, temperature=temperature, fixed=connector.fixed)
    return connector.connected_nuc

mdna.nucleic.compute_rigid_parameters(traj, chainids=[0, 1], fit_reference=False)

Compute the rigid base parameters of the DNA structure.

Parameters:

Name Type Description Default
traj object

MDtraj trajectory containing the DNA structure.

 required
chainids list

List of chain IDs of the DNA structure. Default is [0, 1].

 [0, 1]
fit_reference bool

Fit each base to canonical reference bases before frame extraction. Default is False.

 False

Returns:

Name Type Description
NucleicFrames object

Object representing the rigid base parameters of the DNA structure.

Raises:

Type Description
ValueError

If the traj argument is not provided.
Notes
  • The returned NucleicFrames object contains information about the rigid base parameters of the DNA structure, such as the positions and orientations of the base steps.
Example

Compute the rigid base parameters of a DNA structure ```python traj = md.load('dna.pdb') rigid_params = mdna.compute_rigid_parameters(traj, chainids=[0, 1]) ````

Source code in mdna/nucleic.py 
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
def compute_rigid_parameters(traj, chainids=[0,1], fit_reference=False):
    """Compute the rigid base parameters of the DNA structure.

    Args:
        traj (object): MDtraj trajectory containing the DNA structure.
        chainids (list, optional): List of chain IDs of the DNA structure. Default is [0, 1].
        fit_reference (bool, optional): Fit each base to canonical reference bases before
            frame extraction. Default is False.

    Returns:
        NucleicFrames (object): Object representing the rigid base parameters of the DNA structure.

    Raises:
        ValueError: If the traj argument is not provided.

    Notes:
        - The returned NucleicFrames object contains information about the rigid base parameters of the DNA structure, such as the positions and orientations of the base steps.

    Example:
        Compute the rigid base parameters of a DNA structure
        ```python
        traj = md.load('dna.pdb')
        rigid_params = mdna.compute_rigid_parameters(traj, chainids=[0, 1])
        ````
    """
    if traj is None:
        raise ValueError("The traj argument must be provided.")
    return NucleicFrames(traj, chainids, fit_reference=fit_reference)

mdna.nucleic.sequence_to_pdb(sequence='CGCGAATTCGCG', filename='my_dna', save=True, output='GROMACS', shape=None, n_bp=None, circular=False, dLk=None, save_location='./')

Generate a DNA structure from a DNA sequence code.

Parameters:

Name Type Description Default
sequence str

The DNA sequence code. Default is 'CGCGAATTCGCG'.

 'CGCGAATTCGCG'
filename str

The filename for the pdb output. Default is 'my_dna'.

 'my_dna'
save bool

Whether to save the pdb file. Default is True.

 True
output str

The type of pdb DNA format. Default is 'GROMACS'.

 'GROMACS'
shape ndarray

Control points of shape (n,3) with n > 3 that is used for spline interpolation to determine DNA shape. Default is None, which is a straight line.

 None
n_bp int

Number of base pairs to scale shape with. Default is None, then the sequence is used to determine n_bp.

 None
circular bool

Indicates if the structure is circular/closed. Default is False.

 False
dLk int

Change in twist in terms of Linking number of DNA structure to output. Default is None, which corresponds to a neutral twist based on bp_per_turn = 10.5.

 None
save_location str

Location to save the trajectory. Default is './'.

 './'

Returns:

Type Description
Trajectory

md.Trajectory: An MDtraj trajectory object of the DNA structure (containing only a single frame).

Raises:

Type Description
ValueError

If the sequence is not provided.
Notes
  • The pdb file is saved in the current directory with the specified filename.
Example

Generate a DNA structure from a sequence 

traj = mdna.sequence_to_pdb(sequence='CGCGAATTCGCG', filename='my_dna', save=True, output='GROMACS', shape=None, n_bp=None, circular=False, dLk=None)

Source code in mdna/nucleic.py 
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
def sequence_to_pdb(sequence: str = 'CGCGAATTCGCG', filename: str = 'my_dna', save: bool = True, output: str = 'GROMACS', shape: np.ndarray = None, n_bp: int = None, circular: bool = False, dLk: int = None, save_location : str = './') -> md.Trajectory:
    """Generate a DNA structure from a DNA sequence code.

    Args:
        sequence (str, optional): The DNA sequence code. Default is 'CGCGAATTCGCG'.
        filename (str, optional): The filename for the pdb output. Default is 'my_dna'.
        save (bool, optional): Whether to save the pdb file. Default is True.
        output (str, optional): The type of pdb DNA format. Default is 'GROMACS'.
        shape (np.ndarray, optional): Control points of shape (n,3) with n > 3 that is used for spline interpolation to determine DNA shape. Default is None, which is a straight line.
        n_bp (int, optional): Number of base pairs to scale shape with. Default is None, then the sequence is used to determine n_bp.
        circular (bool, optional): Indicates if the structure is circular/closed. Default is False.
        dLk (int, optional): Change in twist in terms of Linking number of DNA structure to output. Default is None, which corresponds to a neutral twist based on bp_per_turn = 10.5.
        save_location (str, optional): Location to save the trajectory. Default is './'.

    Returns:
        md.Trajectory: An MDtraj trajectory object of the DNA structure (containing only a single frame).

    Raises:
        ValueError: If the sequence is not provided.

    Notes:
        - The pdb file is saved in the current directory with the specified filename.

    Example:
        Generate a DNA structure from a sequence
        ```python
        traj = mdna.sequence_to_pdb(sequence='CGCGAATTCGCG', filename='my_dna', save=True, output='GROMACS', shape=None, n_bp=None, circular=False, dLk=None)
        ```
    """

    # Check if the sequence is provided
    if sequence is None:
        raise ValueError("The sequence argument must be provided.")

    # TODO: Update with make function
    sequence, n_bp = _check_input(sequence=sequence, n_bp=n_bp)

    # Linear strand of control points
    if shape is None:
        shape = Shapes.line(length=1)

    # Convert the control points to a spline
    spline = SplineFrames(control_points=shape, n_bp=n_bp, closed=circular, dLk=dLk)

    # Generate the DNA structure
    generator = StructureGenerator(sequence=sequence, spline=spline, circular=circular)

    # Edit the DNA structure to make it compatible with the AMBER force field
    traj = generator.get_traj(remove_terminal_phosphates=(output == 'GROMACS'))

    # Save the DNA structure as a pdb file
    if save:
        traj.save(f'{save_location}{filename}.pdb')

    return traj

mdna.nucleic.sequence_to_md(sequence=None, time=10, time_unit='picoseconds', temperature=310, solvated=False, filename='my_dna', save=True, output='GROMACS', shape=None, n_bp=None, circular=False, dLk=None, save_location='./')

Simulate DNA sequence using OpenMM.

Parameters:

Name Type Description Default
sequence str

DNA sequence code.

 None
time int

Simulation time.

 10
time_unit str

Time unit (picoseconds or nanoseconds).

 'picoseconds'
temperature int

Temperature in Kelvin.

 310
solvated bool

Solvate DNA with water and ions.

 False
filename str

Filename for pdb output.

 'my_dna'
save bool

Save the trajectory.

 True
output str

Output format for the trajectory (GROMACS or HDF5).

 'GROMACS'
shape str

Shape of the DNA structure (linear or circular).

 None
n_bp int

Number of base pairs in the DNA structure.

 None
circular bool

Flag indicating if the DNA structure is circular.

 False
dLk int

Change in linking number of the DNA structure.

 None
save_location str

Location to save the trajectory.

 './'

Returns:

Name Type Description
MDTraj object

MDtraj trajectory object of DNA structure.
Notes
  • This function uses the OpenMM library to simulate the behavior of a DNA sequence.
  • The simulation can be performed for a specified time period at a given temperature.
  • The DNA structure can be solvated with water and ions.
  • The trajectory of the simulation can be saved in either GROMACS or HDF5 format.
Example

Simulate a linear DNA structure for 100 picoseconds at 300 K 

trajectory = mdna.sequence_to_md(sequence='ATCGATA', time=100, time_unit='picoseconds', temperature=300, shape='linear')

Source code in mdna/nucleic.py 
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
def sequence_to_md(sequence=None, time=10, time_unit='picoseconds',temperature=310, solvated=False,  filename='my_dna', save=True, output='GROMACS',shape=None,n_bp=None,circular=False,dLk=None,save_location='./'):
    """Simulate DNA sequence using OpenMM.

        Args:
            sequence (str): DNA sequence code.
            time (int): Simulation time.
            time_unit (str): Time unit (picoseconds or nanoseconds).
            temperature (int): Temperature in Kelvin.
            solvated (bool): Solvate DNA with water and ions.
            filename (str): Filename for pdb output.
            save (bool): Save the trajectory.
            output (str): Output format for the trajectory (GROMACS or HDF5).
            shape (str): Shape of the DNA structure (linear or circular).
            n_bp (int): Number of base pairs in the DNA structure.
            circular (bool): Flag indicating if the DNA structure is circular.
            dLk (int): Change in linking number of the DNA structure.
            save_location (str): Location to save the trajectory.

        Returns:
            MDTraj (object): MDtraj trajectory object of DNA structure.

        Notes:
            - This function uses the OpenMM library to simulate the behavior of a DNA sequence.
            - The simulation can be performed for a specified time period at a given temperature.
            - The DNA structure can be solvated with water and ions.
            - The trajectory of the simulation can be saved in either GROMACS or HDF5 format.

        Example:
            Simulate a linear DNA structure for 100 picoseconds at 300 K
            ```python
            trajectory = mdna.sequence_to_md(sequence='ATCGATA', time=100, time_unit='picoseconds', temperature=300, shape='linear')
            ```
        """

    # TODO update with make function
    try:
        import openmm as mm
        import openmm.app as app
        import openmm.unit as unit
        from mdtraj.reporters import HDF5Reporter
        import mdtraj as md
        openmm_available = True
    except ImportError:
        openmm_available = False
        print("Openmm is not installed. You shall not pass.")
    if filename is None:
        filename = sequence

    pdb = sequence_to_pdb(sequence=sequence, filename=f'{save_location}{filename}', save=True, output='GROMACS',shape=shape,n_bp=n_bp,circular=circular,dLk=dLk)

    if not openmm_available:
        print('But here is your DNA structure')
        return pdb
    else:
        if time_unit == 'picoseconds':
            simulation_time = time * unit.picoseconds
        elif time_unit == 'nanoseconds':
            simulation_time = time * unit.nanoseconds
        else:
            raise ValueError("time_unit should be 'picoseconds' or 'nanoseconds'")

        time_step = 2 * unit.femtoseconds
        temperature = temperature * unit.kelvin
        steps = int(simulation_time / time_step)

        print(f'Initialize DNA openMM simulation at {temperature._value} K for', simulation_time, 'time units')
        topology = pdb.topology.to_openmm()
        modeller = app.Modeller(topology, pdb.xyz[0])

        forcefield = app.ForceField('amber14-all.xml', 'amber14/tip3pfb.xml')
        modeller.addHydrogens(forcefield)
        if solvated:
            print('Solvate DNA with padding of 1.0 nm and 0.1 M KCl')
            modeller.addSolvent(forcefield, padding=1.0*unit.nanometers, ionicStrength=0.1*unit.molar, positiveIon='K+')

        system = forcefield.createSystem(modeller.topology, nonbondedMethod=app.CutoffNonPeriodic)
        integrator = mm.LangevinIntegrator(temperature, 1.0/unit.picoseconds, time_step)

        simulation = app.Simulation(modeller.topology, system, integrator)
        simulation.context.setPositions(modeller.positions)
        simulation.reporters.append(HDF5Reporter(f'{save_location}{filename}'+'.h5', 100))
        simulation.reporters.append(app.StateDataReporter(f'{save_location}/output_{filename}.csv', 100, step=True, potentialEnergy=True, temperature=True,speed=True))

        print('Minimize energy')
        simulation.minimizeEnergy()

        print('Run simulation for', steps, 'steps')
        simulation.step(steps)
        simulation.reporters[0].close()
        print('Simulation completed')
        print('Saved trajectory as:', f'{save_location}{filename}'+'.h5')
        traj = md.load_hdf5(f'{save_location}{filename}'+'.h5')
        return traj