OxDNA - User contributions [en]

SsDNA Kuhn length

2012-04-21T17:15:32Z

Rovigatti: /* 10-base ssDNA Kuhn length */

[[Category:Examples]]
In this example we will compute the Kuhn length strand of ssDNA with two
different lengths: 10 and 30 nucleotides. Given the simplicity of the system,
the computation required is of the order on one hour to get very accurate
results.

Since ssDNA has stacked and unstacked regions, in constant evolution, the usual
definition of persistence length as a function of the correlation function
between segments is not directly applicable. Instead, we resort (see
[[Publications|Ref. 2]]) to a more general property, which is the Kuhn length defined as:

k = <'''L''' * '''L'''> / ''L''_max

where '''L''' is the end to end vector of the strand and ''L''_max is the
maximum length of the polymer. Brackets indicate thermal average. Since our ssDNA is not a rigid polymer, we use
as a proxy of the maximum length the average distance between nucleotides times the number of covalent
bonds along the backbone.

For this simulations, we will use Brownian Dynamics with the usual values for
the [[Thermostat|thermostat]. We'll run the simulations at room temperature,
and we will use poly-A sequence with the average parametrisation of the model.
We use nucleotides that are all the same to avoid the formation of secondary
structure that could affect the end-to-end distance. With the
sequence-independent parametrisation, a poly-X sequence is for all intents and
purposes the same as any sequence that does not form secondary structure.

The directories <tt>30B</tt> and <tt>10B</tt> inside <tt>EXAMPLES/SSDNA/</tt>
are ready to run the simulations of the 10-base and 30-base ssDNA,
respectively.

The files relative to this example are in the directory <tt>${oxDNA}/EXAMPLES/SSDNAKUHNLENGTH/</tt>.

==10-base ssDNA Kuhn length==
We start by creating a sequence file to
generate the initial configuration. We just need to put 30 A's in a file which
we'll call <tt>seq.txt</tt>

<pre>
AAAAAAAAAAAAAAAAAAAAAAAAAAAAAA
</pre>

and then generate the initial configuration with

<pre>
${OXDNA}/UTILS/generate-sa.py 20. seq.txt
</pre>

which generates <tt>generated.top</tt> and <tt>generated.dat</tt>. We can use
any standart input file for dynamics with the temperature set to 23C, which we
will call <tt>input</tt>. We store configurations every 50000 steps, which
should produce almost independent values for the end-to-end distance with the
standart thermostat settings. Of course this is cannot be known ''a priori'',
but we know this since we have ran this simulation before.

The end-to-end distance of ssDNA, which is the relevant parameter to monitor
for measuring the Kuhn length, has very large fluctuations in our model. A
large number of values is thus needed to have a reliable average; we store 2000
configurations, thus we use a total number of steps of 50 000 * 2 000 = 100
000 000.

The subdirectory <tt>10B</tt> contains a working sequence and input file ready to be run.

The next step is to generate the trajectory by running the program on the input
file.

<pre>
${OXDNA}/Release/oxDNA input
</pre>
Let's now discuss the analysis script <tt>sspl.py</tt> that will extract
the kuhn length from the trajectory. The core of the program is as
follows:

<pre>
line = conffile.readline ()
line = conffile.readline ()
line = conffile.readline ()
line = conffile.readline ()
while line:
# end 2 end
r1 = [float(x) for x in line.split()[0:3]]
a1 = [float(x) for x in line.split()[3:6]]
line = conffile.readline ()

r2 = [float(x) for x in line.split()[0:3]]
a2 = [float(x) for x in line.split()[3:6]]

b1 = [r1[i] - 0.4 * a1[i] for i in xrange(3)]
b2 = [r2[i] - 0.4 * a2[i] for i in xrange(3)]

for i in xrange (nnucl - 2):
line = conffile.readline()

rN = [float(x) for x in line.split()[0:3]]
aN = [float(x) for x in line.split()[3:6]]
bN = [rN[i] - 0.4 * aN[i] for i in xrange(3)]

for i in xrange (4):
line = conffile.readline()

r1N = [rN[i] - r1[i] for i in xrange(3)]
r12 = [r2[i] - r1[i] for i in xrange(3)]
b1N = [bN[i] - b1[i] for i in xrange(3)]
b12 = [b2[i] - b1[i] for i in xrange(3)]

niter += 1

Ll0 /= float (niter)
L2 /= float (niter)
l0 /= float (niter)

Kl = L2 / ((nnucl -1) * l0)

print "Kuhn length: ", Kl
</pre>

We read the position of the first, second and last nucleotide in the system,
skipping all other lines, and compute the thermal average of '''L''' * '''L'''
and ''l_0''. All distances are computed relative to backbone sites, which are
0.4 simulation units away from the centres along the principal axis of the
molecule as described in the [[Documentation#Geometry_of_the_Model|geometry]] of the model. The script produces a file called <tt>end2end.dat</tt> with the time series of the end-to-end distance, which can be used to check for good statistics.

At the end, we use those the averages to compute the Kuhn
length of the strand.

The result coming out should be around 1.9 simulation length units, which
translates to 1.6nm (see [[Model_introduction#Simulation_units|units]]).

==30-Base ssDNA==
All the steps above can be repeated specifying a longer sequence in the
<tt>seq.txt</tt> file.
The subdirectory <tt>30B</tt> contains a working sequence and input file ready to be run. Note that we expect a threefold increase in simulation time due to the system being three times as large. We can use the same python script to analyse the trajectory; the Kuhn length should come out as roughly 2.9 simulation units.

Hairpin formation

2012-04-21T17:14:39Z

Rovigatti: /* Results */

[[Category:Examples]]
== Introduction ==
In this example, you will simulate a single strand of length 18 and sequence <tt>GCGTTGCTTCTCCAACGC</tt> at 334 K (~61 °C) in three different ways:
* with a molecular dynamics (MD) simulation of the sequence-averaged (SA) model. The input file is <tt>inputMD</tt>.
* with an MD simulation of the sequence-dependent (SD) model. The input file is <code>inputMD_seq_dep</code>.
* with a Monte Carlo (MC) simulation of the SA model in which two base pairs are connected by mutual traps (i.e. additional attractive interactions between two nucleotides). The input file is <tt>inputTRAP</tt>.
The traps act between the pairs depicted in blue and red in the sequence G<span style="color:red">C</span><span style="color:blue">G</span>TTGCTTCTCCAA<span style="color:blue">C</span><span style="color:red">G</span>C. The details of the interaction associated to the traps can be changed in the file <tt>hairpin_forces.dat</tt>.

This strand, if ''T'' is sufficiently low, tends to form an hairpin with a 6-base long stem and a 6-base long loop.
The temperature has been chosen to be close to the melting temperature of such a hairpin in the SA version of the model

This document explains how to prepare the ''hairpin'' example (see [[#Preparation|Preparation]]) and how to run it ([[#Running|Running]]). Section [[#Results|Results]] contains results and plots extracted from the simulation output. In the following, <tt>$EXEC</tt> refers to the oxDNA executable.

== Preparation ==
The script <tt>run.sh</tt> generates the input files and runs all the three simulations, one after the other. With the default input files, each simulation, lasting <math>10^8</math> steps by default, takes approximately one hour on a modern CPU. The default <tt>run.sh</tt> expects <tt>$EXEC</tt> to be in the <tt>../..</tt> directory. If this is not the case, open <tt>run.sh</tt> and change the variable <tt>CODEDIR</tt> accordingly.

If you only want to generate the initial configuration, you can issue <tt>./run.sh --generate-only</tt>. Then you can run the simulations by yourself. The generated initial configuration files are <tt>initial.top</tt> (which contains the topology) and <tt>initial.conf</tt> (which contains positions and orientations of the nucleotides).

== Running ==
In order to run the whole example, it is sufficient to issue the command <tt>./run.sh</tt> (or <tt>bash run.sh</tt>). As described in [[#Preparation|Preparation]], you can generate the initial configuration and then run the simulations by hand. The three simulations, described in [[#Introduction|Introduction]], can be performed issuing <tt>$EXEC input</tt>, where <tt>input</tt> is a text file that specifies the simulation configuration. For this example, three files have been prepared: <tt>inputMD</tt>, <tt>inputMD_seq_dep</tt> and <tt>inputTRAP</tt>. The table below reports all the files associated to each simulations.

{|
|-
! Type
! Input
! energy file
! trajectory file
! last configuration file
! log file
|-
|SA model
|inputMD
|energy.dat
|trajectory.dat
|last_conf.dat
|log.dat
|-
|SD model
|inputMD_seq_dep
|energy_seq_dep.dat
|trajectory_seq_dep.dat
|last_conf_seq_dep.dat
|log_seq_dep.dat
|-
|SA model with traps
|inputTRAP
|energy_trap.dat
|trajectory_trap.dat
|last_conf_trap.dat
|log_trap.dat
|-
|}

== Results ==
As mentioned in the [[#Introduction|introduction]], temperature has been chosen so that the hairpin is near its melting temperature when simulated via SA model. Figure 1 shows the probability distribution histogram <math>P(U_{HB})</math> for the hydrogen-bonding (HB) energy <math>U_{HB}</math> of the system (in simulation units). <math>U_{HB}</math> is stored in the fifth column of the energy file as a function of MD steps. Both the minimum in <math>P(U_{HB})</math> and a direct inspection of the time series (inset) show that the hairpin is indeed near the melting temperature, since <math>U_{HB}</math> oscillates between 0 (no HBs) and ~ -3.5, which corresponds to 4-5 hydrogen-bonded nucleotides. At this temperature the fraying effect is relevant, and therefore the last base-pair is not always formed. This is clearly visible in the snapshot at the beginning of this document, in which 4 base-pairs are formed.

<div style="text-align:center">
http://kratos.phys.uniroma1.it/histo_MD.png

<h5>Figure 1: Probability of having a hydrogen-bonding interaction energy <math>U_{HB}</math> for the SA model. Inset: time series of the hydrogen-bonding energy.</h5>
</div>

Figure 2 shows <math>P(U_{HB})</math> and <math>U_{HB}</math> for the SD hairpin. In this case, one expects a different melting temperature. Indeed, at this temperature the hairpin is nearly always in its folded conformation and the average value of <math>U_{HB}</math> is lower than in the SA case. Although stacking interactions play a role, the major contribution to this effect is due to the presence of four <tt>GC</tt> only two <tt>AT</tt> base-pairs in the stem.

<div style="text-align:center">
http://kratos.phys.uniroma1.it/histo_MD_seq_dep.png

<h5>Figure 2: Same as Figure 1, but for the SD model.</h5>
</div>

If mutual traps between stem base-pairs are introduced, then the equilibrium properties of the hairpin are changed and, even if the SA model is employed, the hairpin is always (after the initial equilibration) in its folded conformation. The use of mutual traps can highly decrease the simulation time required by the folding of strands into target structures (like DNA origami or DNA constructs).

<div style="text-align:center">
http://kratos.phys.uniroma1.it/histo_TRAP.png

<h5>Figure 3: Same as Figure 1, but for the system with mutual traps. Time is expressed in Monte Carlo sweeps.</h5>
</div>

Hairpin formation

2012-04-21T17:14:15Z

Rovigatti: /* Results */

[[Category:Examples]]
== Introduction ==
In this example, you will simulate a single strand of length 18 and sequence <tt>GCGTTGCTTCTCCAACGC</tt> at 334 K (~61 °C) in three different ways:
* with a molecular dynamics (MD) simulation of the sequence-averaged (SA) model. The input file is <tt>inputMD</tt>.
* with an MD simulation of the sequence-dependent (SD) model. The input file is <code>inputMD_seq_dep</code>.
* with a Monte Carlo (MC) simulation of the SA model in which two base pairs are connected by mutual traps (i.e. additional attractive interactions between two nucleotides). The input file is <tt>inputTRAP</tt>.
The traps act between the pairs depicted in blue and red in the sequence G<span style="color:red">C</span><span style="color:blue">G</span>TTGCTTCTCCAA<span style="color:blue">C</span><span style="color:red">G</span>C. The details of the interaction associated to the traps can be changed in the file <tt>hairpin_forces.dat</tt>.

This strand, if ''T'' is sufficiently low, tends to form an hairpin with a 6-base long stem and a 6-base long loop.
The temperature has been chosen to be close to the melting temperature of such a hairpin in the SA version of the model

This document explains how to prepare the ''hairpin'' example (see [[#Preparation|Preparation]]) and how to run it ([[#Running|Running]]). Section [[#Results|Results]] contains results and plots extracted from the simulation output. In the following, <tt>$EXEC</tt> refers to the oxDNA executable.

== Preparation ==
The script <tt>run.sh</tt> generates the input files and runs all the three simulations, one after the other. With the default input files, each simulation, lasting <math>10^8</math> steps by default, takes approximately one hour on a modern CPU. The default <tt>run.sh</tt> expects <tt>$EXEC</tt> to be in the <tt>../..</tt> directory. If this is not the case, open <tt>run.sh</tt> and change the variable <tt>CODEDIR</tt> accordingly.

If you only want to generate the initial configuration, you can issue <tt>./run.sh --generate-only</tt>. Then you can run the simulations by yourself. The generated initial configuration files are <tt>initial.top</tt> (which contains the topology) and <tt>initial.conf</tt> (which contains positions and orientations of the nucleotides).

== Running ==
In order to run the whole example, it is sufficient to issue the command <tt>./run.sh</tt> (or <tt>bash run.sh</tt>). As described in [[#Preparation|Preparation]], you can generate the initial configuration and then run the simulations by hand. The three simulations, described in [[#Introduction|Introduction]], can be performed issuing <tt>$EXEC input</tt>, where <tt>input</tt> is a text file that specifies the simulation configuration. For this example, three files have been prepared: <tt>inputMD</tt>, <tt>inputMD_seq_dep</tt> and <tt>inputTRAP</tt>. The table below reports all the files associated to each simulations.

{|
|-
! Type
! Input
! energy file
! trajectory file
! last configuration file
! log file
|-
|SA model
|inputMD
|energy.dat
|trajectory.dat
|last_conf.dat
|log.dat
|-
|SD model
|inputMD_seq_dep
|energy_seq_dep.dat
|trajectory_seq_dep.dat
|last_conf_seq_dep.dat
|log_seq_dep.dat
|-
|SA model with traps
|inputTRAP
|energy_trap.dat
|trajectory_trap.dat
|last_conf_trap.dat
|log_trap.dat
|-
|}

== Results ==
As mentioned in the [[#Introduction|Introduction]], temperature has been chosen so that the hairpin is near its melting temperature when simulated via SA model. Figure 1 shows the probability distribution histogram <math>P(U_{HB})</math> for the hydrogen-bonding (HB) energy <math>U_{HB}</math> of the system (in simulation units). <math>U_{HB}</math> is stored in the fifth column of the energy file as a function of MD steps. Both the minimum in <math>P(U_{HB})</math> and a direct inspection of the time series (inset) show that the hairpin is indeed near the melting temperature, since <math>U_{HB}</math> oscillates between 0 (no HBs) and ~ -3.5, which corresponds to 4-5 hydrogen-bonded nucleotides. At this temperature the fraying effect is relevant, and therefore the last base-pair is not always formed. This is clearly visible in the snapshot at the beginning of this document, in which 4 base-pairs are formed.

<div style="text-align:center">
http://kratos.phys.uniroma1.it/histo_MD.png

<h5>Figure 1: Probability of having a hydrogen-bonding interaction energy <math>U_{HB}</math> for the SA model. Inset: time series of the hydrogen-bonding energy.</h5>
</div>

Figure 2 shows <math>P(U_{HB})</math> and <math>U_{HB}</math> for the SD hairpin. In this case, one expects a different melting temperature. Indeed, at this temperature the hairpin is nearly always in its folded conformation and the average value of <math>U_{HB}</math> is lower than in the SA case. Although stacking interactions play a role, the major contribution to this effect is due to the presence of four <tt>GC</tt> only two <tt>AT</tt> base-pairs in the stem.

<div style="text-align:center">
http://kratos.phys.uniroma1.it/histo_MD_seq_dep.png

<h5>Figure 2: Same as Figure 1, but for the SD model.</h5>
</div>

If mutual traps between stem base-pairs are introduced, then the equilibrium properties of the hairpin are changed and, even if the SA model is employed, the hairpin is always (after the initial equilibration) in its folded conformation. The use of mutual traps can highly decrease the simulation time required by the folding of strands into target structures (like DNA origami or DNA constructs).

<div style="text-align:center">
http://kratos.phys.uniroma1.it/histo_TRAP.png

<h5>Figure 3: Same as Figure 1, but for the system with mutual traps. Time is expressed in Monte Carlo sweeps.</h5>
</div>

Hairpin formation

2012-04-21T17:13:43Z

Rovigatti:

[[Category:Examples]]
== Introduction ==
In this example, you will simulate a single strand of length 18 and sequence <tt>GCGTTGCTTCTCCAACGC</tt> at 334 K (~61 °C) in three different ways:
* with a molecular dynamics (MD) simulation of the sequence-averaged (SA) model. The input file is <tt>inputMD</tt>.
* with an MD simulation of the sequence-dependent (SD) model. The input file is <code>inputMD_seq_dep</code>.
* with a Monte Carlo (MC) simulation of the SA model in which two base pairs are connected by mutual traps (i.e. additional attractive interactions between two nucleotides). The input file is <tt>inputTRAP</tt>.
The traps act between the pairs depicted in blue and red in the sequence G<span style="color:red">C</span><span style="color:blue">G</span>TTGCTTCTCCAA<span style="color:blue">C</span><span style="color:red">G</span>C. The details of the interaction associated to the traps can be changed in the file <tt>hairpin_forces.dat</tt>.

This strand, if ''T'' is sufficiently low, tends to form an hairpin with a 6-base long stem and a 6-base long loop.
The temperature has been chosen to be close to the melting temperature of such a hairpin in the SA version of the model

This document explains how to prepare the ''hairpin'' example (see [[#Preparation|Preparation]]) and how to run it ([[#Running|Running]]). Section [[#Results|Results]] contains results and plots extracted from the simulation output. In the following, <tt>$EXEC</tt> refers to the oxDNA executable.

== Preparation ==
The script <tt>run.sh</tt> generates the input files and runs all the three simulations, one after the other. With the default input files, each simulation, lasting <math>10^8</math> steps by default, takes approximately one hour on a modern CPU. The default <tt>run.sh</tt> expects <tt>$EXEC</tt> to be in the <tt>../..</tt> directory. If this is not the case, open <tt>run.sh</tt> and change the variable <tt>CODEDIR</tt> accordingly.

If you only want to generate the initial configuration, you can issue <tt>./run.sh --generate-only</tt>. Then you can run the simulations by yourself. The generated initial configuration files are <tt>initial.top</tt> (which contains the topology) and <tt>initial.conf</tt> (which contains positions and orientations of the nucleotides).

== Running ==
In order to run the whole example, it is sufficient to issue the command <tt>./run.sh</tt> (or <tt>bash run.sh</tt>). As described in [[#Preparation|Preparation]], you can generate the initial configuration and then run the simulations by hand. The three simulations, described in [[#Introduction|Introduction]], can be performed issuing <tt>$EXEC input</tt>, where <tt>input</tt> is a text file that specifies the simulation configuration. For this example, three files have been prepared: <tt>inputMD</tt>, <tt>inputMD_seq_dep</tt> and <tt>inputTRAP</tt>. The table below reports all the files associated to each simulations.

{|
|-
! Type
! Input
! energy file
! trajectory file
! last configuration file
! log file
|-
|SA model
|inputMD
|energy.dat
|trajectory.dat
|last_conf.dat
|log.dat
|-
|SD model
|inputMD_seq_dep
|energy_seq_dep.dat
|trajectory_seq_dep.dat
|last_conf_seq_dep.dat
|log_seq_dep.dat
|-
|SA model with traps
|inputTRAP
|energy_trap.dat
|trajectory_trap.dat
|last_conf_trap.dat
|log_trap.dat
|-
|}

== Results ==
As mentioned in the [[#Introduction|Introduction]], temperature has been chosen so that the hairpin is near its melting temperature when simulated via SA model. Figure 1 shows the probability distribution histogram <math>P(U_{HB})</math> for the hydrogen-bonding (HB) energy <math>U_{HB}</math> of the system (in simulation units. Both the minimum in <math>P(U_{HB})</math> and a direct inspection of the time series (inset) show that the hairpin is indeed near the melting temperature, since <math>U_{HB}</math> oscillates between 0 (no HBs) and ~ -3.5, which corresponds to 4-5 hydrogen-bonded nucleotides. At this temperature the fraying effect is relevant, and therefore the last base-pair is not always formed. This is clearly visible in the snapshot at the beginning of this document, in which 4 base-pairs are formed.

<div style="text-align:center">
http://kratos.phys.uniroma1.it/histo_MD.png

<h5>Figure 1: Probability of having a hydrogen-bonding interaction energy <math>U_{HB}</math> for the SA model. Inset: time series of the hydrogen-bonding energy.</h5>
</div>

Figure 2 shows <math>P(U_{HB})</math> and <math>U_{HB}</math> for the SD hairpin. In this case, one expects a different melting temperature. Indeed, at this temperature the hairpin is nearly always in its folded conformation and the average value of <math>U_{HB}</math> is lower than in the SA case. Although stacking interactions play a role, the major contribution to this effect is due to the presence of four <tt>GC</tt> only two <tt>AT</tt> base-pairs in the stem.

<div style="text-align:center">
http://kratos.phys.uniroma1.it/histo_MD_seq_dep.png

<h5>Figure 2: Same as Figure 1, but for the SD model.</h5>
</div>

If mutual traps between stem base-pairs are introduced, then the equilibrium properties of the hairpin are changed and, even if the SA model is employed, the hairpin is always (after the initial equilibration) in its folded conformation. The use of mutual traps can highly decrease the simulation time required by the folding of strands into target structures (like DNA origami or DNA constructs).

<div style="text-align:center">
http://kratos.phys.uniroma1.it/histo_TRAP.png

<h5>Figure 3: Same as Figure 1, but for the system with mutual traps. Time is expressed in Monte Carlo sweeps.</h5>
</div>

DsDNA persistence length

2012-04-21T17:08:46Z

Rovigatti: /* Persistence length of a double-stranded DNA */

[[Category:Examples]]
== Persistence length of a double-stranded DNA ==

The example shows how to calculate a persistence length of a double stranded DNA molecule.
dsDNA persistence length. The persistence length in this example is calculated using the following formula (see [http://jcp.aip.org/resource/1/jcpsa6/v134/i8/p085101_s1?bypassSSO=1] for details):

http://www-thphys.physics.ox.ac.uk/people/PetrSulc/images/eqn.png

In the <tt> EXAMPLES/PERSISTENCE_LENGTH</tt> directory, you will find a setup for calculating the persistence length of a 202 base pairs long dsDNA.
Note that for calculating a persistence length of a dsDNA, one needs a large number of decorrelated states. To obtain the states (which will be saved into a trajectory file), run the simulation program using the prepared input_persistence file:

<pre>oxDNA input_persistence</pre>

The program will run a molecular dynamics simulation at 23 °C and record the individual configurations. By default, they are saved in the <tt>trajectory.dat</tt> file. To analyze the data, use the python script <tt>dspl.py</tt>:

<pre>dspl.py trajectory.dat init.top 10 50</pre>

This program will produce a table of correlations between helical vectors, http://www-thphys.physics.ox.ac.uk/people/PetrSulc/images/eqn2.png.

The program <tt>dspl.py</tt> requires <tt>base.py</tt> in the UTILS directory to be present and to have UTILS directory set in your PYTHONPATH environment variable. The program calculates the local helical axis vector (n_k) as a unit vector pointing from the midpoint of hydrogen bonding sites of k-th base pair to the midpoint between (k+1)-th base pair.
The init.top file contains topology of the 202 base pairs long strand (included in the <tt> EXAMPLES/PERSISTENCE_LENGTH</tt>). In the example above, the program starts at the 10-th base pair and calculates correlations of n_10 with n_11, n_12, etc. up to n_60. It then prints out the correlations (one per line). Using an exponential fit to these data, one can find the persistence length, as illustrated in the following picture:

http://www-thphys.physics.ox.ac.uk/people/PetrSulc/images/ds.png

The exponential fit shows, in this particular example, a persistence length of 124.8 base pairs.

DsDNA persistence length

2012-04-21T17:07:47Z

Rovigatti: /* Persistence length of a double-stranded DNA */

[[Category:Examples]]
== Persistence length of a double-stranded DNA ==

The example shows how to calculate a persistence length of a double stranded DNA molecule.
dsDNA persistence length. The persistence length in this example is calculated using the following formula (see [http://jcp.aip.org/resource/1/jcpsa6/v134/i8/p085101_s1?bypassSSO=1] for details):

http://www-thphys.physics.ox.ac.uk/people/PetrSulc/images/eqn.png

In the <tt> EXAMPLES/PERSISTENCE_LENGTH</tt> directory, you will find a setup for calculating the persistence length of a 202 base pairs long dsDNA.
Note that for calculating a persistence length of a dsDNA, one needs a large number of decorrelated states. To obtain the states (which will be saved into a trajectory file), run the simulation program using the prepared input_persistence file:

<pre>oxDNA input_persistence</pre>

The program will run a molecular dynamics simulation at 23 °C and record the individual configurations. They are saved in trajectory.dat file. To analyze the data, use the python script dspl.py:

<pre>dspl.py trajectory.dat init.top 10 50</pre>

This program will produce a table of correlations between helical vectors, http://www-thphys.physics.ox.ac.uk/people/PetrSulc/images/eqn2.png.

The program <tt>dspl.py</tt> requires <tt>base.py</tt> in the UTILS directory to be present and to have UTILS directory set in your PYTHONPATH environment variable. The program calculates the local helical axis vector (n_k) as a unit vector pointing from the midpoint of hydrogen bonding sites of k-th base pair to the midpoint between (k+1)-th base pair.
The init.top file contains topology of the 202 base pairs long strand (included in the <tt> EXAMPLES/PERSISTENCE_LENGTH</tt>). In the example above, the program starts at the 10-th base pair and calculates correlations of n_10 with n_11, n_12, etc. up to n_60. It then prints out the correlations (one per line). Using an exponential fit to these data, one can find the persistence length, as illustrated in the following picture:

http://www-thphys.physics.ox.ac.uk/people/PetrSulc/images/ds.png

The exponential fit shows, in this particular example, a persistence length of 124.8 base pairs.

DsDNA persistence length

2012-04-21T17:06:59Z

Rovigatti: /* Persistence length of a double-stranded DNA */

[[Category:Examples]]
== Persistence length of a double-stranded DNA ==

The example shows how to calculate a persistence length of a double stranded DNA molecule.
dsDNA persistence length. The persistence length in this example is calculated using the following formula (see [http://jcp.aip.org/resource/1/jcpsa6/v134/i8/p085101_s1?bypassSSO=1] for details):

http://www-thphys.physics.ox.ac.uk/people/PetrSulc/images/eqn.png

In the <tt> EXAMPLES/PERSISTENCE_LENGTH </tt> directory, you will find a setup for calculating the persistence length of a 202 base pairs long dsDNA.
Note that for calculating a persistence length of a dsDNA, one needs a large number of decorrelated states. To obtain the states (which will be saved into a trajectory file), run the simulation program using the prepared input_persistence file:

<pre>oxDNA input_persistence</pre>

The program will run a molecular dynamics simulation at 23 °C and record the individual configurations. They are saved in trajectory.dat file. To analyze the data, use the python script dspl.py:

<pre>dspl.py trajectory.dat init.top 10 50</pre>

This program will produce a table of correlations between helical vectors, http://www-thphys.physics.ox.ac.uk/people/PetrSulc/images/eqn2.png.

The program dspl.py requires base.py in the UTILS directory to be present and to have UTILS directory set in your PYTHONPATH environment variable. The program calculates the local helical axis vector (n_k) as a unit vector pointing from the midpoint of hydrogen bonding sites of k-th base pair to the midpoint between (k+1)-th base pair.
The init.top file contains topology of the 202 base pairs long strand (included in the <tt> EXAMPLES/PERSISTENCE_LENGTH </tt>). In the example above, the program starts at the 10-th base pair and calculates correlations of n_10 with n_11, n_12, etc. up to n_60. It then prints out the correlations (one per line). Using an exponential fit to these data, one can find the persistence length, as illustrated in the following picture:

http://www-thphys.physics.ox.ac.uk/people/PetrSulc/images/ds.png

The exponential fit shows, in this particular example, a persistence length of 124.8 base pairs.

Thermostat

2012-04-20T09:43:45Z

Rovigatti:

The best thermostat implemented in the simulation code (''john'' thermostat) is a simple thermostat
that emulates Brownian dynamics. The system is evolved integrating Newton's
equations of motion ('NVE' ensemble) for a given (small) number of steps.
Then the velocity and momentum of each particle are refreshed, with a given
fixed probability. The new velocities and momenta are chosen according to
the Maxwell distribution of the temperature at which the simulation is ran.
This approximates a Brownian dynamics on time scales much longer than the
refresh interval.

* '''dt''' time steps of the integration
* '''newtonian steps''' number of steps between refresh attempts
* '''pt''' the probability with which each particle gets its velocity and momentum refreshed at each attempt.
* '''diff_coeff''' the overall monomer diffusion coefficient resulting from the thermostat. The code internally sets <tt>pt</tt> to get this value. Specifying <tt>pt</tt> will override this, regardless which comes first in the input file.

The algorithm works as follows: the system is evolved for a number of steps equal to <tt>newtonian_steps</tt> according to Newton's equations of motion. Than for each particle a random number is extracted; if it is larger than the value for <tt>pt</tt> (either set explicitly or derived
from <tt>diff_coeff</tt>) the particle is left untouched. If the random number extracted is lower than the value of <tt>pt</tt>, each of the components of the the velocity and angular momentum of the particle are refreshed according to
the Maxwell distribution dictated by the value of the temperature.

A completely Browninan dynamics (on the time scale set by <tt>dt</tt>) can be obtained setting <tt>pt = 1</tt> and <tt>newtonian_steps = 0</tt>. Of course, this makes little sense.

Increasing the value of <tt>pt</tt> makes the dynamics more dumped, so that overall diffusion is slower but local motion is somewhat better explored. We found that a good thermostat setting to study diffusion-limited events is to set <tt>diff_coeff = 2.5</tt> and <tt>newtonian_steps = 103</tt>. If one is not limited by diffusion, internal relaxation can be speeded up by lowering the value of <tt>diff_coeff</tt> by a factor 2 or 4.

Thermostat

2012-04-20T09:38:42Z

Rovigatti:

The thermostat implemented in the simulation code is a simple thermostat
that emulates Brownian dynamics. The system is evolved integrating Newton's
equations of motion ('NVE' ensemble) for a given (small) number of steps.
Then the velocity and momentum of each particle are refreshed, with a given
fixed probability. The new velocities and momenta are chosen according to
the Maxwell distribution of the temperature at which the simulation is ran.
This approximates a Brownian dynamics on time scales much longer than the
refresh interval.

* '''dt''' time steps of the integration
* '''newtonian steps''' number of steps between refresh attempts
* '''pt''' the probability with which each particle gets its velocity and momentum refreshed at each attempt.
* '''diff_coeff''' the overall monomer diffusion coefficient resulting from the thermostat. The code internally sets <tt>pt</tt> to get this value. Specifying <tt>pt</tt> will override this, regardless which comes first in the input file.

The algorithm works as follows: the system is evolved for a number of steps equal to <tt>newtonian_steps</tt> according to Newton's equations of motion. Than for each particle a random number is extracted; if it is larger than the value for <tt>pt</tt> (either set explicitly or derived
from <tt>diff_coeff</tt>) the particle is left untouched. If the random number extracted is lower than the value of <tt>pt</tt>, each of the components of the the velocity and angular momentum of the particle are refreshed according to
the Maxwell distribution dictated by the value of the temperature.

A completely Browninan dynamics (on the time scale set by <tt>dt</tt>) can be obtained setting <tt>pt = 1</tt> and <tt>newtonian_steps = 0</tt>. Of course, this makes little sense.

Increasing the value of <tt>pt</tt> makes the dynamics more dumped, so that overall diffusion is slower but local motion is somewhat better explored. We found that a good thermostat setting to study diffusion-limited events is to set <tt>diff_coeff = 2.5</tt> and <tt>newtonian_steps = 103</tt>. If one is not limited by diffusion, internal relaxation can be speeded up by lowering the value of <tt>diff_coeff</tt> by a factor 2 or 4.

Documentation

2012-04-20T09:35:49Z

Rovigatti: /* Geometry of the Model */

==Compile options==

Compiling oxDNA requires that you change the first rows in the makefile to match your machine configuration. The following parameters can be passed to make:

*'''dbg=1''' oxDNA will be compiled with debug flags (both for nvcc and gcc). The resulting executable will be put in the Debug directory.
*'''g=1''' oxDNA will be compiled with both debug and optimization flags. The resulting executable will be put in the Release directory.
*'''intel=1''' oxDNA will be compiled using the Intel icpc compiler. The resulting executable will be named oxDNA_intel.

==Usage==

<pre>oxDNA input_file</pre>

The input file contains all the relevant information for the program to run, such as what initial configuration to use, the topology of the system, how often to print the energies to a file, etc. Please make sure you read the [[Thermostat|thermostat]] page if you use molecular dynamics.

==Input file==

As always in UNIX environments, everything is case sensitive.
The options are in the form key = value. There can be arbitrary spaces before and after both key and value. Line with a leading # will be treated as comments.
In this part | (pipe) is the separator between the different values that can be used to specify a value for the key.
Keys between [ and ] are optional, the value after the equal sign is the default value.

===Generic options===
The options listed here define the generic behavior of the entire program.
;[sim_type=MD]: MD|MC
:MD = Molecular Dynamics, MC = Monte Carlo
;backend: CPU
;backend_precision: float|double
;[debug=0]: 0|1
: 1 if you want verbose logs, 0 otherwise.

===Simulation options===
The options listed here specify the behaviour of the simulation.

;steps: number of steps to be performed.

;[restart_step_counter=0]: 0|1
:0 means that the step counter will start from the value read in the configuration file; if 1, the step counter will be reset to 0. The total duration of the simulation is unchanged.

;[seed=time(NULL)]: seed for the random number generator. On Unix systems, it will use by default a number from /dev/urandom + time(NULL)

;T: temperature of the simulation. It can be expressed in simulation units or kelvin (append a k or K after the value) or celsius (append a c or C after the value).
:Examples:
{|
|-
! Value
! Simulation Units
|-
| 0.1
| 0.1
|-
| 300 K
| 0.1
|-
| 300k
| 0.1
|-
| 26.85c
| 0.1
|-
| 26.85 C
| 0.1
|-
|}

;verlet_skin: if a particle moves more than verlet_skin then the lists will be updated. Its name is somewhat misleading: the actual verlet skin is 2*verlet_skin.

;[use_average_seq=1]: 0|1
: specifies whether to use the default hard-coded average parameters for base-pairing and stacking interaction strengths or not. If sequence dependence is to be used, set this to 0 and specify seq_dep_file.

;[seq_dep_file]: specifies the file from which the sequence dependent parameters should be read. Mandatory if use_average_seq=no, ignored otherwise. A sample file is provided (sequence_dependent_parameters.txt).

;[external_forces=0]: 0|1
: specifies whether there are external forces acting on the nucleotides or not. If it is set to 1, then a file which specifies the external forces' configuration has to be provided (see external_forces_file).

;[external_forces_file]: specifies the file containing all the external forces' configurations. Currently there are six supported force types (see EXAMPLES/TRAPS for some examples):
:*string
:*twist
:*trap
:*repulsion_plane
:*repulsion_plane_moving
:*mutual_trap

====Molecular dynamics simulations options====

;dt: time step of the integration.

;thermostat: no|refresh|john
:no means no thermostat will be used. refresh will refresh all the particle's velocities from a maxwellian every newtonian_steps steps. john is an Anderson-like thermostat (see pt). Make sure you read [[Thermostat|thermostat]].

;newtonian_steps: required if thermostat != no
:number of steps after which a procedure of thermalization will be performed.

;pt: used if thermostat == john. It's the probability that a particle's velocity will be refreshed during a thermalization procedure.

;diff_coeff: required if pt is not specified
:used internally to automatically compute the pt that would be needed if we wanted such a self diffusion coefficient. Not used if pt is set.

====Monte Carlo simulations options====

;[check_energy_every=10]: this number times print_energy_every gives the number of steps after which the energy will be computed from scratch and checked against the actual value computed adding energy differences.

;[check_energy_threshold=1e-4]: if abs((old_energy - new_energy)/old_energy) > check_energy_threshold then the program will die and warn the user.

;ensemble: NVT
:ensemble of the simulation. More ensembles could be added in future versions.

;delta_translation: maximum displacement (per dimension) for translational moves in simulation units.

;delta_translation: maximum displacement for rotational moves in simulation units.

===Input/output===
The options listed here are used to manage the I/O (read and write configurations, energies and so on)

;conf_file: initial configuration file.

;topology: file containing the system's topology.

;trajectory_file: the main output of the program. All the configurations will be appended to this file as they are printed.

;[confs_to_skip=0]: valid only if conf_file is a trajectory. Skip the first confs_to_skip configurations and then load in memory the (confs_to_skip+1)th.

;[lastconf_file=last_conf.dat]: this is the file where the last configuration is saved (when the program finishes or is killed). Set to last_conf.dat by default

;[refresh_vel=0]: 0|1
:if 1 the initial velocities will be refreshed from a maxwellian.

;energy_file: energy output file.

;[print_energy_every=1000]: this will make the program print the energies every print_energy_every steps.

;[no_stdout_energy=0]: 0|1
:if 1 the energy will be printed just to the energy_file.

;[time_scale=linear]: linear|log_lin
:using linear configurations will be saved every print_conf_interval.
:using log_lin configurations will be saved logarithmically for print_conf_ppc times. After that the logarithmic sequence will restart.

;print_conf_interval: linear interval if time_scale == linear. First step of the logarithmic scale if time_scale == log_lin.

;print_conf_ppc: used if time_scale == log_lin
:points per logarithmic cycle.

;[print_reduced_conf_every=0]: every print_reduced_conf_every steps the program will print out the reduced configurations (i.e. confs containing only the centers of mass of strands).

;reduced_conf_output_dir: used if print_reduced_conf_every > 0
:output directory for reduced_conf files.

;[log_file=stderr]: file where generic and debug informations will be logged. If not specified then stderr will be used.

;[print_timings=0]: 0|1
:if 1 the MD step timing have be printed to a file.

;timings_filename: used if print_timings == 1
:output file where the MD step timing will be appended to.

==Output files==
*The log file contains all relevant informations about the simulation (specified options, activated external forces, warnings about misconfiguratios, critical errors, etc.). If the log file is omitted, all these informations will be displayed on the standard output.

*The energy file layout for MD simulations is

:{|
| time
| potential energy
| kinetic energy
| total energy
| hydrogen bonding energy
|}

:while for MC simulations is

:{|
| time
| potential energy
| hydrogen bonding energy
| acceptance ratio for translational moves
| acceptance ratio for rotational moves
|}

:Mind that potential, kinetic and total energies are divided by the number of particles whereas the hydrogen bonding energy is not.

*Configurations are saved in the trajectory file.

==Configuration and topology files==
The current state of a system, as by oxDNA, is described by two files: a configuration file and a topology file. The configuration file contains all the general informations (timestep, energy and box size) and orientations and positions of each nucleotide. The topology file, on the other hand, keeps track of the backbone-backbone bonds between nucleotides in the same strand. Working configuration and topology files can be found in the <tt>[[Examples|EXAMPLES]]</tt> directory.

===Configuration file===
The first three rows of a configuration file contain the timestep <tt>T</tt> at which the configuration has been printed, the length of the box sides <tt>Lx</tt>, <tt>Ly</tt> and <tt>Lz</tt> and the total, potential and kinetic energies, <tt>Etot</tt>, <tt>U</tt> and <tt>K</tt>, respectively:

<pre>
t = T
b = Lz Ly Lz
E = Etot U K
</pre>

after this header, each row contains position of the centre of mass, orientation, velocity and angular velocity of a single nucleotide in the following order:

<math>\overbrace{r_x r_y r_z}^{\rm Position} \overbrace{b_x b_y b_z}^{\rm Backbone-base versor} \overbrace{n_x n_y n_z}^{\rm Normal versor} \overbrace{v_x v_y v_z}^{\rm Velocity} \overbrace{L_x L_y L_z}^{\rm Angular velocity}</math>

===Topology file===
The topology file stores the intra-strand, fixed bonding topology (i.e. which nucleotides share backbone links). The first row contains the total number of nucleotides <tt>N</tt> and the number of strands <tt>Ns</tt>:

<pre>
N Ns
</pre>

After this header, the <tt>i</tt>-th row specifies strand, base and 3' and 5' neighbors of the <tt>i</tt>-th nucleotide in this way:

<pre>
S B 3' 5'
</pre>

where S is the number of the strand (starting from 1) which the nucleotide belongs to, B is the base and 3' and 5' specify the index of the nucleotides with which the <tt>i</tt>-th nucleotide is bonded in the 3' and 5' direction, respectively. A <tt>-1</tt> signals that the nucleotide terminates the strand in either 3' or 5' direction. The topology file of a strand of sequence <tt>GCGTTG</tt> would be:

<pre>
6 1
1 G -1 1
1 C 0 2
1 G 1 3
1 T 2 4
1 T 3 5
1 G 4 -1
</pre>

Specifying the topology in this way can simplify the process of simulating, for example, circular DNA.

===Generation of initial configurations===
In order to generate initial configuration and topology files, we provide the <tt>${oxDNA}/UTILS/generate-sa.py</tt> script. The usage of the script is

<pre>generate-sa.py <box side> <file with sequence></pre>

where <tt><box side></tt> specifies the length of the box side in simulation units and <tt><file with sequence></tt> contains the sequence of the strands to be generated, one row per strand. If double strands are needed, each sequence must be preceded by <tt>DOUBLE</tt>. For example, a file containing

<pre>
DOUBLE AGGGCT
CCTGTA
</pre>

would generate a double strand with a sequence <tt>AGGGCT</tt> and a single strand with a sequence <tt>CCTGTA</tt>.

Positions and orientations of the strands are all chosen at random in such a way that the resulting initial configuration does not contain significant excluded volume interactions between nucleotides belonging to different strands. Generated single- and double-strands have helical conformations (i.e. they are in the minimum of the intra-strand interaction energy).

The output configuration and topology are stored in <tt>generated.dat</tt> and <tt>generated.top</tt>, respectively.
Since this script will initialize nucleotides' velocities and angular velocities to 0, when performing molecular (or Brownian) dynamics simulation remember to put <tt>refresh_vel = 1</tt> in the [[Documentation#Input_file|input]] file.

==Analysis of configurations==
The configurations produced by oxDNA can be analysed with the <tt>output_bonds</tt> program in <tt>${oxDNA}/UTILS/process_data/</tt> directory. This program takes an input the input file (to recover the temperature and topology file), a configuration/trajectory file and an optional number. Since <tt>output_bonds</tt> reads analyses a single configuration, the optional number selects the configuration which it needs to analyse in the trajectory. Analysing a whole trajectory can be done by looping over a counter.

Please note that <tt>output_bonds</tt> is not compiled automatically. If you never compiled it, do so as described in the [[Download_and_Installation#Installation|installation instructions]].

<tt>output_bonds</tt> can be used as follows:
<pre>
${oxDNA}/UTILS/process_data/output_bonds <input_file> <trajectory_file> [counter]
</pre>

The program outputs some debugging information to the standard error and information regarding the interaction energies to the standard output. The contributions arising from each of the terms in the potential (see the appendix of [[Publications|Ref. 2]]) are reported for each pair of nucleotides that have non-zero total interactions.

This output can be easily parsed to analyse the configurations.

For each pair of nucleotides that do interact in the configuration, the program prints out a line containing:
* The id of the two particles (starting from 0)
* The total interaction energy
* The hydrogen bonding (base pairing) energy
* The stacking energy
* The cross stacking energy
* The excluded volume energy
* The FENE interaction energy
* A letter indicating a status code. This will be <tt>N</tt> for pairs that interact through bonded interactions (i.e. they are neighbors along a strand) and it will be <tt>H</tt> when a base pair is present. Our definition of base pair is when two nucleotides have a hydrogen bonding energy less than 0.1 in simulation units (see [[Publications|Ref. 2]]).

===Geometry of the Model===
In the configuration/trajectory files only the positions and orientations of the nucleotides are stored. If one wants to recover the positions of the individual interaction sites in the model, some maths need to be done.

The position of the base, stacking and backbone sites can be recovered as follows:

base site: (center) + 0.40 * (axis vector)

stacking site: (center) + 0.34 * (axis vector)

backbone site: (center) - 0.40 * (axis_vector)

The picture in the [[Model_introduction|introduction]] might help understanding where the sites are.

==External Forces==
The code implements several types of external forces that can be imposed on the system that can be used either to simulate tension exerted on DNA or simply to accelerate the formation of secondary or tertiary structure. External forces can be tricky to treat, especially in a dynamics simulation, since they are an external source of work. Care should be taken in adjusting the time step, thermostat parameters and such.

To enable external forces, one needs to specify <tt>external_forces = 1</tt> in the input file and also supply an external force file to read from with the key <tt>external_forces_file = <file></tt>.

The syntax of the external forces file is quite simple. Examples of such files can be found in the [[Hairpin_formation|hairpin formation]] and [[Pseudoknot|Pseudoknot formation]] examples. Each force is specified within a block contained in curly brackets. Empty lines and lines beginning with an hash symbol (<tt>#</tt>) are ignored. Different forces require different keys to be present. If the file has the wrong syntax, oxDNA should spit out a sensible error message while parsing the file.

The different types of forces implemented at the moment are:
* harmonic trap
* string
* repulsion plane
* mutual trap

All forces act on the centre of the particle.

Forces of different kinds can be combined in the same simulation. There is a maximum number of 10 external forces per particle for memory reasons. This can be manually overridden recompiling the code with a different value of the macro <tt>MAX_EXT_FORCES</tt> (currently 10) in <tt>defs.h</tt>.

===String===
A string is implemented as a force that does not depend on the particle position. Its value can be constant or can change linearly with time. It is useful as it does not fluctuate with time.

A force of this kind is specified with <tt>type = string</tt>. The relevant keys are:
* '''particle''' (int) the particle on which to exert the force
* '''F0''' (float) the value of the force at time = 0 in simulation units (please note that the value of the time may or may not be reset when starting a simulation, depending on the input file)
* '''rate''' (float) growing rate of the force (simulation units/time steps). Typical values are very small (< 10^(-5))
* '''dir''' (3 floats separated by commas) direction of the force (automatically normalised by the code)

The following bit of code will create an external force on the first nucleotide in the system starting at 1 simulation units (48.15 pN) and growing linearly with time at the rate of 48.15pN every million time steps. The force will pull the nucleotide along the <tt>z</tt> direction.

<pre>
{
type = string
particle = 0
F0 = 1.
rate = 1e-6
dir = 0., 0., 1.
}
</pre>

===Harmonic trap===
This type of force implements an harmonic trap, of arbitrary stiffness, that can move linearly with time. It can be useful to fix the position of the nucleotides to simulate attachment to something or to implement (quasi) constant extension simulations.

A force of this kind is specified with <tt>type = trap</tt>. The relevant keys are:
* '''particle''' (int) the particle on which to exert the force
* '''pos0''' (3 floats separated by commas) rest position of the trap
* '''stiff''' (float) stiffness of the trap (the force is stiff * dx)
* '''rate''' (float) speed of the trap (length simulation units/time steps)
* '''dir''' (3 floats separated by commas) direction of movement of the trap

Here is an example input for a harmonic trap acting on the third nucleotide constraining it to stay close to the origin. In this example the trap does not move (<tt>rate=0</tt>), but one could have it move at a constant speed along the direction specified by <tt>dir</tt>, in this case the <tt>x</tt> direction.

<pre>
{
type = trap
particle = 2
pos0 = 0., 0., 0.
stiff = 1.0
rate = 0.
dir = 1.,0.,0.
}
</pre>

Please note that the trap does not comply with periodic boundary conditions. This is most likely what you want.

===Repulsion plane===
This kind of external force implements a repulsion plane that constrains a particle (or the whole system) to stay on one side of it. It is implemented as a harmonic repulsion, but the stiffness can be made arbitrarily high to mimic a hard repulsion.

A force of this kind is specified with <tt>type = repulsion_plane</tt>. The relevant keys are:
* '''particle''' (int) the particle on which to exert the force. If set to the special value -1, the force will be exerted on all particles.
* '''stiff''' (float) stiffness of the trap (the force is stiff * D, where D is distance from the plane. The force is exerted only if the nucleotide is below the plane)
* '''dir''' (3 floats separated by commas) a direction normal to the plane
* '''position''' (1 float number) specifies the position of the plane

If direction is <tt> direction = u,v,w </tt> , then the plane contains all the points (x,y,z) that satisfy the equation: u*x + v*y + w*z + position = 0.
Only nucleotides with coordinates (x,y,z) that satisfy u*x + v*y + w*z + position < 0 will feel the force.
The force exerted on a nucleotide is equal to stiff * D, where D is the distance of the nucleotide from the plane, where <tt> D = | u*x + v*y + w*z + position | / \sqrt(v^2 + u^2 + z^2 ).</tt>
For nucleotides for which u*x + v*y + w*z + position >= 0, no force will be exerted.

Here is an example. This plane acts on the whole system and will not exert any force on nucleotides with a positive <tt>x</tt> coordinate. A force proportional to 96.3pN * (<tt>x</tt> coordinate) will be exerted on all particles .

<pre>
{
type = repulsion_plane
#whole system
particle = -1
stiff = 1. #96.3 pN in simulation units
dir = 1, 0, 0
position = 0
}
</pre>

If in the above example you would specify position = 3, then the force would be exerted on all nucleotides with coordinate x > -3.

===Mutual trap===
This force is useful to form initial configurations. It is a harmonic force that at every moment pulls a particle towards a reference particle. It is possible to specify the separation at which the force will be 0.

A force of this kind is specified with <tt>type = mutual_trap</tt>. The relevant keys are:
* '''particle''' (int) the particle on which to exert the force.
* '''ref_particle''' (int) particle to pull towards. Please note that this particle will not feel any force (the name mutual trap is thus misleading).
* '''stiff''' (float) stiffness of the trap
* '''r0''' (float) equilibrium distance of the trap.
* '''position''' (3 floats separated by commas) one point belonging to the plane.

Here is an example, extracted from the [[Pseudoknot|pseudoknot formation example]]. This will pull particle 14 towards particle 39, favouring an equilibrium distance of 1.4 (which corresponds roughly to the minimum of the hydrogen bonding potential, not a coincidence). The same force with opposite sign is exerted on particle 39 through a separate force. It is not necessary to have both particles feel the force, but it is usually works much better.

<pre>
{
type = mutual_trap
particle = 14
ref_particle = 39
stiff = 1.
r0 = 1.2
}

{
type = mutual_trap
particle = 39
ref_particle = 14
stiff = 1.
r0 = 1.2
}
</pre>

==Visualisation of structures==
oxDNA produces a trajectory file where all the relevant information is
stored. A converter is provided (<tt>traj2vis.py</tt>) in the
<tt>UTILS</tt> directory that is able to produce files in the <tt>xyz</tt>
and <tt>pdb</tt> formats. The same program can be used on a configuration
file and it will produce a snapshot.

Since the model is coearse-grained, we have to "trick" the visualisers into
thinking that the interaction sites in the model are actually atoms.
Advanced nucleic acids representations such as ribbons will not work on the
outputs.

All the images in the [[Screenshots]] page were produced with the pdb representation using UCSF chimera (see later on).

===xyz format===

just run

<pre>$oxDNA/UTILS/traj2vis.py xyz <trajectory> <topology> </pre>

(where <tt>$oxDNA</tt> is the oxDNA source directory) to get the xyz representation in a file called the same as the trajectory
file with <tt>.xyz</tt> appended. Please note that boundary conditions are
implemented strand-wise, so strands that are bound might appear at two
different sizes of the box. Also, the center of mass of the system (where
each strand is weighted the same regardless of the length) is set to 0 at
each frame. Carbons represent the backbone sites and oxygens the base sites.

The resulting file can be read with a variety of programs. Here we will
explain how to visualise it sensibly in [http://www.ks.uiuc.edu/Research/vmd/ VMD].

* Run VMD and load the xyz file.
* In the graphics menu, go to Representations.
* In the Selected Atoms line, input <tt>name C</tt>. Also select Drawing method CPK, sphere scale 0.8 and Bond Radius 0.
* In the Selected Atoms line, input <tt>name O</tt>. Also select Drawing method CPK, sphere scale 0.6 and Bond Radius 0.

This should produce a ball representation of our model DNA. Bonds
automatically produced by VMD are NOT meaningful in our context.

===pdb format===
Run

<pre>$oxDNA/UTILS/traj2vis.py xyz <trajectory> <topology> </pre>

to produce a trajectory/configuration in the pdb format. A further file
called <tt>chimera.com</tt> will be produced (more on this later). All
comments above about periodic boundaries and centre of mass apply here as
well.

The pdb file can be visualised in VMD just like the xyz format, but a nicer
output can be produced with [http://www.cgl.ucsf.edu/chimera/ UCSF Chimera] (although only for snapshots at
the moment) as follows:

Run chimera and load the pdb file. An ugly output will be displayed.

Bring up the command line under the <tt>Tools → General Controls</tt> menu.
Input <tt>read chimera.com</tt> in the command line and press enter. You
should get a nicer visualisation with different bases in different colors,
all the covalent bonds in the right place, etc.

On large configurations, the production of ellipsoids will be extremely
slow. You can remove it by removing the line

<code>aniso scale 0.75 smoothing 4</code>

from the commands file. Loading the resulting file should be much faster.

UCSF chimera can in turn export the scene in a variety of formats.

Documentation

2012-04-20T09:34:49Z

Rovigatti: /* Topology file */

Model introduction

2012-04-20T09:32:47Z

Rovigatti: /* Simulation units */

The model treats DNA as a string of rigid nucleotides, which interact through potentials which depend on the position and orientation of the nucleotides. The interactions are:
#Sugar-phosphate backbone connectivity,
#Excluded volume,
#Hydrogen bonding,
#Nearest-neighbour stacking,
#Cross-stacking between base-pair steps in a duplex,
#Coaxial stacking.

This interactions are illustrated below. Orientational modulations of the stacking potential encourage the bases to form coplanar stacks, the twist arising from the different length scales of the backbone separation and the optimum stacking separation. The possibility of unstacking allows single strands to be very flexible. Hydrogen bonding can occur between complementary bases when they are anti-aligned, leading to the formation of double helical structures.

http://www-thphys.physics.ox.ac.uk/people/ThomasOuldridge/Site/The_model.png

(a) Model interaction sites with their interaction ranges (the typical range of an interaction is twice the radius of the sphere shown).

(b) Representation of these interaction site in a visualisation that makes the planarity of the base clear.

(c) A duplex in this representation.

http://www-thphys.physics.ox.ac.uk/people/ThomasOuldridge/Site/interactions.png

Indication of the interactions which hold together a typical duplex. V(b.b.) indicates the phosphate-sugar backbone connectivity.

In the original model, all complementary base pairs and stacking partners interact with the same strength (there is no attractive interaction between non-complementary bases). A preliminary sequence-dependent parameterisation of the hydrogen-bonding and stacking interactions is included as an option in the code release: a paper discussing this parameterisation and its effects is currently in preparation.

The model does not incorporate the differentiation between the major and minor grooves of DNA double helices, and incorporates electrostatics only through the short-ranged excluded volume. For this reason, it is only appropriate for the study of systems at high salt concentration, when electrostatic interactions are strongly screened.

===Simulation units===
The code uses dimensionless energy, mass, length and timescales for convenience. The relationship between simulation units (SU) and SI units is given below.

{|
|-
! Simulation unit
! Physical unit
|-
| 1 unit of length
| 8.518x10<math>^{-10}</math> m
|-
| 1 unit of energy
| 4.142<math>^{-20}</math> J
|-
| 1 unit of Temperature
| 3000 K
|-
| 1 unit of Force
| 4.863<math>^{-11}</math> N
|-
| 1 unit of Mass
| 1.66<math>^{-25}</math> kg
|-
| 1 unit of Time
| 1.71<math>^{-12}</math> s
|-
|}

The model and its performance is discussed in detail in the following references (the thesis provides the most complete analysis):

T. E. Ouldridge, D.Phil. Thesis, University of Oxford, 2011.
[http://ora.ox.ac.uk/objects/uuid:b2415bb2-7975-4f59-b5e2-8c022b4a3719 Coarse-grained modelling of DNA and DNA self-assembly]

T. E. Ouldridge, A. A. Louis and J. P. K. Doye, J. Chem. Phys, 134, 085101 (2011)
[http://link.aip.org/link/?JCP/134/085101 Structural, mechanical and thermodynamic properties of a coarse-grained DNA model] ([http://arxiv.org/abs/arXiv:1009.4480 arXiv])

Hairpin formation

2012-04-19T15:23:07Z

Rovigatti:

[[Category:Examples]]
== Introduction ==
In this example, you will simulate a single strand of length 18 and sequence <tt>GCGTTGCTTCTCCAACGC</tt> at 334 K (~61 °C) in three different ways:
* with a molecular dynamics (MD) simulation of the sequence-averaged (SA) model. The input file is <tt>inputMD</tt>.
* with an MD simulation of the sequence-dependent (SD) model. The input file is <code>inputMD_seq_dep</code>.
* with a Monte Carlo (MC) simulation of the SA model in which two base pairs are connected by mutual traps (i.e. additional attractive interactions between two nucleotides). The input file is <tt>inputTRAP</tt>.
The traps act between the pairs depicted in blue and red in the sequence G<span style="color:red">C</span><span style="color:blue">G</span>TTGCTTCTCCAA<span style="color:blue">C</span><span style="color:red">G</span>C. The details of the interaction associated to the traps can be changed in the file <tt>hairpin_forces.dat</tt>.

This strand, if ''T'' is sufficiently low, tends to form an hairpin with a 6-base long stem and a 6-base long loop.
The temperature has been chosen to be close to the melting temperature of such a hairpin in the SA version of the model

This document explains how to prepare the ''hairpin'' example (see [[#Preparation|Preparation]]) and how to run it ([[#Running|Running]]). Section [[#Results|Results]] contains results and plots extracted from the simulation output. In the following, <tt>$EXEC</tt> refers to the oxDNA executable.

== Preparation ==
The script <tt>run.sh</tt> generates the input files and runs all the three simulations, one after the other. With the default input files, each simulation, lasting <math>10^8</math> steps by default, takes approximately one hour on a modern CPU. The default <tt>run.sh</tt> expects <tt>$EXEC</tt> to be in the <tt>../..</tt> directory. If this is not the case, open <tt>run.sh</tt> and change the variable <tt>CODEDIR</tt> accordingly.

If you only want to generate the initial configuration, you can issue <tt>./run.sh --generate-only</tt>. Then you can run the simulations by yourself. The generated initial configuration files are <tt>initial.top</tt> (which contains the topology) and <tt>initial.conf</tt> (which contains positions and orientations of the nucleotides).

== Running ==
In order to run the whole example, it is sufficient to issue the command <tt>./run.sh</tt> (or <tt>bash run.sh</tt>). As described in [[#Preparation|Preparation]], you can generate the initial configuration and then run the simulations by hand. The three simulations, described in [[#Introduction|Introduction]], can be performed issuing <tt>$EXEC input</tt>, where <tt>input</tt> is a text file that specifies the simulation configuration. For this example, three files have been prepared: <tt>inputMD</tt>, <tt>inputMD_seq_dep</tt> and <tt>inputTRAP</tt>. The table below reports all the files associated to each simulations.

{|
|-
! Type
! Input
! energy file
! trajectory file
! last configuration file
! log file
|-
|SA model
|inputMD
|energy.dat
|trajectory.dat
|last_conf.dat
|log.dat
|-
|SD model
|inputMD_seq_dep
|energy_seq_dep.dat
|trajectory_seq_dep.dat
|last_conf_seq_dep.dat
|log_seq_dep.dat
|-
|SA model with traps
|inputTRAP
|energy_trap.dat
|trajectory_trap.dat
|last_conf_trap.dat
|log_trap.dat
|-
|}

== Results ==
As mentioned in the [[#Introduction|Introduction]], temperature has been chosen so that the hairpin is near its melting temperature when simulated via SA model. Figure 1 shows the probability distribution histogram <math>P(U_{HB})</math> for the hydrogen-bonding (HB) energy <math>U_{HB}</math> of the system (in simulation units. Both the minimum in <math>P(U_{HB})</math> and a direct inspection of the time series (inset) show that the hairpin is indeed near the melting temperature, since <math>U_{HB}</math> oscillates between 0 (no HBs) and ~ -3.5, which corresponds to 4-5 hydrogen-bonded nucleotides. At this temperature the fraying effect is relevant, and therefore the last base-pair is not always formed. This is clearly visible in the snapshot at the beginning of this document, in which 4 base-pairs are formed.

<div style="text-align:center">
http://kratos.phys.uniroma1.it/histo_MD.png

<h5>Figure 1: Probability of having a hydrogen-bonding interaction energy <math>U_{HB}</math> for the SA model. Inset: time series of the hydrogen-bonding energy (stored in the fifth column of the energy file as a function of MD steps.</h5>
</div>

Figure 2 shows <math>P(U_{HB})</math> and <math>U_{HB}</math> for the SD hairpin. In this case, one expects a different melting temperature. Indeed, at this temperature the hairpin is nearly always in its folded conformation and the average value of <math>U_{HB}</math> is lower than in the SA case. Although stacking interactions play a role, the major contribution to this effect is due to the presence of four <tt>GC</tt> only two <tt>AT</tt> base-pairs in the stem.

<div style="text-align:center">
http://kratos.phys.uniroma1.it/histo_MD_seq_dep.png

<h5>Figure 2: Same as Figure 1, but for the SD model.</h5>
</div>

If mutual traps between stem base-pairs are introduced, then the equilibrium properties of the hairpin are changed and, even if the SA model is employed, the hairpin is always (after the initial equilibration) in its folded conformation. The use of mutual traps can highly decrease the simulation time required by the folding of strands into target structures (like DNA origami or DNA constructs).

<div style="text-align:center">
http://kratos.phys.uniroma1.it/histo_TRAP.png

<h5>Figure 3: Same as Figure 1, but for the system with mutual traps. Time is expressed in Monte Carlo sweeps.</h5>
</div>

Hairpin formation

2012-04-19T15:21:49Z

Rovigatti: /* Results */

[[Category:Examples]]
== Introduction ==
In this example, you will simulate a single strand of length 18 and sequence <tt>GCGTTGCTTCTCCAACGC</tt> at 334 K (~61 °C) in three different ways:
* with a molecular dynamics (MD) simulation of the sequence-averaged (SA) model. The input file is <tt>inputMD</tt>.
* with an MD simulation of the sequence-dependent (SD) model. The input file is <code>inputMD_seq_dep</code>.
* with a Monte Carlo (MC) simulation of the SA model in which two base pairs are connected by mutual traps (i.e. additional attractive interactions between two nucleotides). The input file is <tt>inputTRAP</tt>.
The traps act between the pairs depicted in blue and red in the sequence G<span style="color:red">C</span><span style="color:blue">G</span>TTGCTTCTCCAA<span style="color:blue">C</span><span style="color:red">G</span>C. The details of the interaction associated to the traps can be changed in the file <tt>hairpin_forces.dat</tt>.

This strand, if ''T'' is sufficiently low, tends to form an hairpin with a 6-base long stem and a 6-base long loop.
The temperature has been chosen to be close to the melting temperature of such a hairpin in the SA version of the model

This document explains how to prepare the ''hairpin'' example (see [[#Preparation|Preparation]]) and how to run it ([[#Running|Running]]). Section [[#Results|Results]] contains results and plots extracted from the simulation output. In the following, <tt>$EXEC</tt> refers to the oxDNA executable.

== Preparation ==
The script <tt>run.sh</tt> generates the input files and runs all the three simulations, one after the other. With the default input files, each simulation, lasting <math>10^8</math> steps by default, takes approximately one hour on a modern CPU. The default <tt>run.sh</tt> expects <tt>$EXEC</tt> to be in the <tt>../..</tt> directory. If this is not the case, open <tt>run.sh</tt> and change the variable <tt>CODEDIR</tt> accordingly.

If you only want to generate the initial configuration, you can issue <tt>./run.sh --generate-only</tt>. Then you can run the simulations by yourself. The generated initial configuration files are <tt>initial.top</tt> (which contains the topology) and <tt>initial.conf</tt> (which contains positions and orientations of the nucleotides).

== Running ==
In order to run the whole example, it is sufficient to issue the command <tt>./run.sh</tt> (or <tt>bash run.sh</tt>). As described in [[#Preparation|Preparation]], you can generate the initial configuration and then run the simulations by hand. The three simulations, described in [[#Introduction|Introduction]], can be performed issuing <tt>$EXEC input</tt>, where <tt>input</tt> is a text file that specifies the simulation configuration. For this example, three files have been prepared: <tt>inputMD</tt>, <tt>inputMD_seq_dep</tt> and <tt>inputTRAP</tt>. The table below reports all the files associated to each simulations.

{|
|-
! Type
! Input
! energy file
! trajectory file
! last configuration file
! log file
|-
|SA model
|inputMD
|energy.dat
|trajectory.dat
|last_conf.dat
|log.dat
|-
|SD model
|inputMD_seq_dep
|energy_seq_dep.dat
|trajectory_seq_dep.dat
|last_conf_seq_dep.dat
|log_seq_dep.dat
|-
|SA model with traps
|inputTRAP
|energy_trap.dat
|trajectory_trap.dat
|last_conf_trap.dat
|log_trap.dat
|-
|}

== Results ==
As mentioned in the [[#Introduction|Introduction]], temperature has been chosen so that the hairpin is near its melting temperature when simulated via SA model. Figure 1 shows the probability distribution histogram <math>P(U_{HB})</math> for the hydrogen-bonding (HB) energy <math>U_{HB}</math> of the system (in simulation units. Both the minimum in <math>P(U_{HB})</math> and a direct inspection of the time series (inset) show that the hairpin is indeed near the melting temperature, since <math>U_{HB}</math> oscillates between 0 (no HBs) and ~ -3.5, which corresponds to 4-5 hydrogen-bonded nucleotides. At this temperature the fraying effect is relevant, and therefore the last base-pair is not always formed. This is clearly visible in the snapshot at the beginning of this document, in which 4 base-pairs are formed.

<div style="text-align:center">
http://kratos.phys.uniroma1.it/histo_MD.png

<h5>Figure 1: Probability of having a hydrogen-bonding interaction energy <math>U_{HB}</math> for the SA model. Inset: time series of the hydrogen-bonding energy (stored in the fifth column of the energy file as a function of MD steps.</h5>
</div>

Figure 2 shows <math>P(U_{HB})</math> and <math>U_{HB}</math> for the SD hairpin. In this case, one expects a different melting temperature. Indeed, at this temperature the hairpin is nearly always in its folded conformation and the average value of <math>U_{HB}</math> is lower than in the SA case. Although stacking interactions play a role, the major contribution to this effect is due to the presence of four <tt>GC</tt> only two <tt>AT</tt> base-pairs in the stem.

<div style="text-align:center">
http://kratos.phys.uniroma1.it/histo_MD_seq_dep.png

<h5>Figure 2: Same as Figure 1, but for the SD model.</h5>
</div>

If mutual traps between stem base-pairs are introduced, then the equilibrium properties of the hairpin are changed and, even if the SA model is employed, the hairpin is always (after the initial equilibration) in its folded conformation. The use of mutual traps can highly decrease the simulation time required by the folding of strands into target structures (like DNA origami or DNA constructs).

<div style="text-align:center">
http://kratos.phys.uniroma1.it/histo_TRAP.png

<h5>Figure 3: Same as Figure 1, but for the system with mutual traps. Time is expressed in Monte Carlo sweeps.</h5>
</div>

Hairpin formation

2012-04-19T15:21:39Z

Rovigatti: /* Running */

[[Category:Examples]]
== Introduction ==
In this example, you will simulate a single strand of length 18 and sequence <tt>GCGTTGCTTCTCCAACGC</tt> at 334 K (~61 °C) in three different ways:
* with a molecular dynamics (MD) simulation of the sequence-averaged (SA) model. The input file is <tt>inputMD</tt>.
* with an MD simulation of the sequence-dependent (SD) model. The input file is <code>inputMD_seq_dep</code>.
* with a Monte Carlo (MC) simulation of the SA model in which two base pairs are connected by mutual traps (i.e. additional attractive interactions between two nucleotides). The input file is <tt>inputTRAP</tt>.
The traps act between the pairs depicted in blue and red in the sequence G<span style="color:red">C</span><span style="color:blue">G</span>TTGCTTCTCCAA<span style="color:blue">C</span><span style="color:red">G</span>C. The details of the interaction associated to the traps can be changed in the file <tt>hairpin_forces.dat</tt>.

This strand, if ''T'' is sufficiently low, tends to form an hairpin with a 6-base long stem and a 6-base long loop.
The temperature has been chosen to be close to the melting temperature of such a hairpin in the SA version of the model

This document explains how to prepare the ''hairpin'' example (see [[#Preparation|Preparation]]) and how to run it ([[#Running|Running]]). Section [[#Results|Results]] contains results and plots extracted from the simulation output. In the following, <tt>$EXEC</tt> refers to the oxDNA executable.

== Preparation ==
The script <tt>run.sh</tt> generates the input files and runs all the three simulations, one after the other. With the default input files, each simulation, lasting <math>10^8</math> steps by default, takes approximately one hour on a modern CPU. The default <tt>run.sh</tt> expects <tt>$EXEC</tt> to be in the <tt>../..</tt> directory. If this is not the case, open <tt>run.sh</tt> and change the variable <tt>CODEDIR</tt> accordingly.

If you only want to generate the initial configuration, you can issue <tt>./run.sh --generate-only</tt>. Then you can run the simulations by yourself. The generated initial configuration files are <tt>initial.top</tt> (which contains the topology) and <tt>initial.conf</tt> (which contains positions and orientations of the nucleotides).

== Running ==
In order to run the whole example, it is sufficient to issue the command <tt>./run.sh</tt> (or <tt>bash run.sh</tt>). As described in [[#Preparation|Preparation]], you can generate the initial configuration and then run the simulations by hand. The three simulations, described in [[#Introduction|Introduction]], can be performed issuing <tt>$EXEC input</tt>, where <tt>input</tt> is a text file that specifies the simulation configuration. For this example, three files have been prepared: <tt>inputMD</tt>, <tt>inputMD_seq_dep</tt> and <tt>inputTRAP</tt>. The table below reports all the files associated to each simulations.

{|
|-
! Type
! Input
! energy file
! trajectory file
! last configuration file
! log file
|-
|SA model
|inputMD
|energy.dat
|trajectory.dat
|last_conf.dat
|log.dat
|-
|SD model
|inputMD_seq_dep
|energy_seq_dep.dat
|trajectory_seq_dep.dat
|last_conf_seq_dep.dat
|log_seq_dep.dat
|-
|SA model with traps
|inputTRAP
|energy_trap.dat
|trajectory_trap.dat
|last_conf_trap.dat
|log_trap.dat
|-
|}

== Results ==
As mentioned in the [[#Introduction|Introduction]], temperature has been chosen so that the hairpin is near its melting temperature when simulated via SA model. Figure 1 shows the probability distribution histogram <math>P(U_{HB})</math> for the hydrogen-bonding (HB) energy <math>U_{HB}</math> of the system (in simulation units. Both the minimum in <math>P(U_{HB})</math> and a direct inspection of the time series (inset) show that the hairpin is indeed near the melting temperature, since <math>U_{HB}</math> oscillates between 0 (no HBs) and ~ -3.5, which corresponds to 4-5 hydrogen-bonded nucleotides. At this temperature the fraying effect is relevant, and therefore the last base-pair is not always formed~\cite{ouldridge_jcp}. This is clearly visible in the snapshot at the beginning of this document, in which 4 base-pairs are formed.

<div style="text-align:center">
http://kratos.phys.uniroma1.it/histo_MD.png

<h5>Figure 1: Probability of having a hydrogen-bonding interaction energy <math>U_{HB}</math> for the SA model. Inset: time series of the hydrogen-bonding energy (stored in the fifth column of the energy file as a function of MD steps.</h5>
</div>

Figure 2 shows <math>P(U_{HB})</math> and <math>U_{HB}</math> for the SD hairpin. In this case, one expects a different melting temperature. Indeed, at this temperature the hairpin is nearly always in its folded conformation and the average value of <math>U_{HB}</math> is lower than in the SA case. Although stacking interactions play a role, the major contribution to this effect is due to the presence of four <tt>GC</tt> only two <tt>AT</tt> base-pairs in the stem.

<div style="text-align:center">
http://kratos.phys.uniroma1.it/histo_MD_seq_dep.png

<h5>Figure 2: Same as Figure 1, but for the SD model.</h5>
</div>

If mutual traps between stem base-pairs are introduced, then the equilibrium properties of the hairpin are changed and, even if the SA model is employed, the hairpin is always (after the initial equilibration) in its folded conformation. The use of mutual traps can highly decrease the simulation time required by the folding of strands into target structures (like DNA origami or DNA constructs).

<div style="text-align:center">
http://kratos.phys.uniroma1.it/histo_TRAP.png

<h5>Figure 3: Same as Figure 1, but for the system with mutual traps. Time is expressed in Monte Carlo sweeps.</h5>
</div>

Hairpin formation

2012-04-19T15:21:10Z

Rovigatti:

[[Category:Examples]]
== Introduction ==
In this example, you will simulate a single strand of length 18 and sequence <tt>GCGTTGCTTCTCCAACGC</tt> at 334 K (~61 °C) in three different ways:
* with a molecular dynamics (MD) simulation of the sequence-averaged (SA) model. The input file is <tt>inputMD</tt>.
* with an MD simulation of the sequence-dependent (SD) model. The input file is <code>inputMD_seq_dep</code>.
* with a Monte Carlo (MC) simulation of the SA model in which two base pairs are connected by mutual traps (i.e. additional attractive interactions between two nucleotides). The input file is <tt>inputTRAP</tt>.
The traps act between the pairs depicted in blue and red in the sequence G<span style="color:red">C</span><span style="color:blue">G</span>TTGCTTCTCCAA<span style="color:blue">C</span><span style="color:red">G</span>C. The details of the interaction associated to the traps can be changed in the file <tt>hairpin_forces.dat</tt>.

This strand, if ''T'' is sufficiently low, tends to form an hairpin with a 6-base long stem and a 6-base long loop.
The temperature has been chosen to be close to the melting temperature of such a hairpin in the SA version of the model

This document explains how to prepare the ''hairpin'' example (see [[#Preparation|Preparation]]) and how to run it ([[#Running|Running]]). Section [[#Results|Results]] contains results and plots extracted from the simulation output. In the following, <tt>$EXEC</tt> refers to the oxDNA executable.

== Preparation ==
The script <tt>run.sh</tt> generates the input files and runs all the three simulations, one after the other. With the default input files, each simulation, lasting <math>10^8</math> steps by default, takes approximately one hour on a modern CPU. The default <tt>run.sh</tt> expects <tt>$EXEC</tt> to be in the <tt>../..</tt> directory. If this is not the case, open <tt>run.sh</tt> and change the variable <tt>CODEDIR</tt> accordingly.

If you only want to generate the initial configuration, you can issue <tt>./run.sh --generate-only</tt>. Then you can run the simulations by yourself. The generated initial configuration files are <tt>initial.top</tt> (which contains the topology) and <tt>initial.conf</tt> (which contains positions and orientations of the nucleotides).

== Running ==
In order to run the whole example, it is sufficient to issue the command <tt>./run.sh</tt> (or <tt>bash run.sh</tt>). As described in [[#Preparation|Preparation]], you can generate the initial configuration and then run the simulations by hand. The three simulations, described in [[#Introduction|Introduction]], can be performed issuing <tt>$EXEC input</tt>, where <tt>input</tt> is a text file that specifies the simulation configuration. For this example, three files have been prepared: <tt>inputMD</tt>, <tt>inputMD_seq_dep</tt> and <tt>inputTRAP</tt>. Table~\ref{tbl:sim} report all the files associated to each simulations.

{|
|-
! Type
! Input
! energy file
! trajectory file
! last configuration file
! log file
|-
|SA model
|inputMD
|energy.dat
|trajectory.dat
|last_conf.dat
|log.dat
|-
|SD model
|inputMD_seq_dep
|energy_seq_dep.dat
|trajectory_seq_dep.dat
|last_conf_seq_dep.dat
|log_seq_dep.dat
|-
|SA model with traps
|inputTRAP
|energy_trap.dat
|trajectory_trap.dat
|last_conf_trap.dat
|log_trap.dat
|-
|}

== Results ==
As mentioned in the [[#Introduction|Introduction]], temperature has been chosen so that the hairpin is near its melting temperature when simulated via SA model. Figure 1 shows the probability distribution histogram <math>P(U_{HB})</math> for the hydrogen-bonding (HB) energy <math>U_{HB}</math> of the system (in simulation units. Both the minimum in <math>P(U_{HB})</math> and a direct inspection of the time series (inset) show that the hairpin is indeed near the melting temperature, since <math>U_{HB}</math> oscillates between 0 (no HBs) and ~ -3.5, which corresponds to 4-5 hydrogen-bonded nucleotides. At this temperature the fraying effect is relevant, and therefore the last base-pair is not always formed~\cite{ouldridge_jcp}. This is clearly visible in the snapshot at the beginning of this document, in which 4 base-pairs are formed.

<div style="text-align:center">
http://kratos.phys.uniroma1.it/histo_MD.png

<h5>Figure 1: Probability of having a hydrogen-bonding interaction energy <math>U_{HB}</math> for the SA model. Inset: time series of the hydrogen-bonding energy (stored in the fifth column of the energy file as a function of MD steps.</h5>
</div>

Figure 2 shows <math>P(U_{HB})</math> and <math>U_{HB}</math> for the SD hairpin. In this case, one expects a different melting temperature. Indeed, at this temperature the hairpin is nearly always in its folded conformation and the average value of <math>U_{HB}</math> is lower than in the SA case. Although stacking interactions play a role, the major contribution to this effect is due to the presence of four <tt>GC</tt> only two <tt>AT</tt> base-pairs in the stem.

<div style="text-align:center">
http://kratos.phys.uniroma1.it/histo_MD_seq_dep.png

<h5>Figure 2: Same as Figure 1, but for the SD model.</h5>
</div>

If mutual traps between stem base-pairs are introduced, then the equilibrium properties of the hairpin are changed and, even if the SA model is employed, the hairpin is always (after the initial equilibration) in its folded conformation. The use of mutual traps can highly decrease the simulation time required by the folding of strands into target structures (like DNA origami or DNA constructs).

<div style="text-align:center">
http://kratos.phys.uniroma1.it/histo_TRAP.png

<h5>Figure 3: Same as Figure 1, but for the system with mutual traps. Time is expressed in Monte Carlo sweeps.</h5>
</div>

Hairpin formation

2012-04-19T15:19:51Z

Rovigatti: /* Results */

[[Category:Examples]]
== Introduction ==
In this example, you will simulate a single strand of length 18 and sequence <tt>GCGTTGCTTCTCCAACGC</tt> at 334 K (~61 °C) in three different ways:
* with a molecular dynamics (MD) simulation of the sequence-averaged (SA) model. The input file is <tt>inputMD</tt>.
* with an MD simulation of the sequence-dependent (SD) model. The input file is <code>inputMD_seq_dep</code>.
* with a Monte Carlo (MC) simulation of the SA model in which two base pairs are connected by mutual traps (i.e. additional attractive interactions between two nucleotides). The input file is <tt>inputTRAP</tt>.
The traps act between the pairs depicted in blue and red in the sequence G<span style="color:red">C</span><span style="color:blue">G</span>TTGCTTCTCCAA<span style="color:blue">C</span><span style="color:red">G</span>C. The details of the interaction associated to the traps can be changed in the file <tt>hairpin_forces.dat</tt>.

This strand, if ''T'' is sufficiently low, tends to form an hairpin with a 6-base long stem and a 6-base long loop.
The temperature has been chosen to be close to the melting temperature of such a hairpin in the SA version of the model

This document explains how to prepare the ''hairpin'' example (see [[#Preparation|Preparation]]) and how to run it ([[#Running|Running]]). Section [[#Results|Results]] contains results and plots extracted from the simulation output. In the following, <tt>$EXEC</tt> refers to the oxDNA executable.

== Preparation ==
The script <tt>run.sh</tt> generates the input files and runs all the three simulations, one after the other. With the default input files, each simulation, lasting <math>10^8</math> steps by default, takes approximately one hour on a modern CPU. The default <tt>run.sh</tt> expects <tt>$EXEC</tt> to be in the <tt>../..</tt> directory. If this is not the case, open <tt>run.sh</tt> and change the variable <tt>CODEDIR</tt> accordingly.

If you only want to generate the initial configuration, you can issue <tt>./run.sh --generate-only</tt>. Then you can run the simulations by yourself. The generated initial configuration files are <tt>initial.top</tt> (which contains the topology) and <tt>initial.conf</tt> (which contains positions and orientations of the nucleotides).

== Running ==
In order to run the whole example, it is sufficient to issue the command <tt>./run.sh</tt> (or <tt>bash run.sh</tt>). As described in [[#Preparation|Preparation]], you can generate the initial configuration and then run the simulations by hand. The three simulations, described in [[#Introduction|Introduction]], can be performed issuing <tt>$EXEC input</tt>, where <tt>input</tt> is a text file that specifies the simulation configuration. For this example, three files have been prepared: <tt>inputMD</tt>, <tt>inputMD_seq_dep</tt> and <tt>inputTRAP</tt>. Table~\ref{tbl:sim} report all the files associated to each simulations.

{|
|-
! Type
! Input
! energy file
! trajectory file
! last configuration file
! log file
|-
|SA model
|inputMD
|energy.dat
|trajectory.dat
|last_conf.dat
|log.dat
|-
|SD model
|inputMD_seq_dep
|energy_seq_dep.dat
|trajectory_seq_dep.dat
|last_conf_seq_dep.dat
|log_seq_dep.dat
|-
|SA model with traps
|inputTRAP
|energy_trap.dat
|trajectory_trap.dat
|last_conf_trap.dat
|log_trap.dat
|-
|}

== Results ==
As mentioned in the [[#Introduction|Introduction]], temperature has been chosen so that the hairpin is near its melting temperature when simulated via SA model. Figure 1 shows the probability distribution histogram $P(U_{HB})$ for the hydrogen-bonding (HB) energy $U_{HB}$ of the system (in simulation units~\cite{ouldridge_jcp}). Both the minimum in $P(U_{HB})$ and a direct inspection of the time series (inset) show that the hairpin is indeed near the melting temperature, since $U_{HB}$ oscillates between $0$ (no HBs) and ~ -3.5, which corresponds to 4-5 hydrogen-bonded nucleotides. At this temperature the fraying effect is relevant, and therefore the last base-pair is not always formed~\cite{ouldridge_jcp}. This is clearly visible in the snapshot at the beginning of this document, in which $4$ base-pairs are formed.

<div style="text-align:center">
http://kratos.phys.uniroma1.it/histo_MD.png

<h5>Figure 1: Probability of having a hydrogen-bonding interaction energy <math>U_{HB}</math> for the SA model. Inset: time series of the hydrogen-bonding energy (stored in the fifth column of the energy file as a function of MD steps.</h5>
</div>

Figure 2 shows $P(U_{HB})$ and $U_{HB}$ for the SD hairpin. In this case, one expects a different melting temperature. Indeed, at this temperature the hairpin is nearly always in its folded conformation and the average value of $U_{HB}$ is lower than in the SA case. Although stacking interactions play a role, the major contribution to this effect is due to the presence of four <tt>GC</tt> only two <tt>AT</tt> base-pairs in the stem.

<div style="text-align:center">
http://kratos.phys.uniroma1.it/histo_MD_seq_dep.png

<h5>Figure 2: Same as Figure 1, but for the SD model.</h5>
</div>

If mutual traps between stem base-pairs are introduced, then the equilibrium properties of the hairpin are changed and, even if the SA model is employed, the hairpin is always (after the initial equilibration) in its folded conformation. The use of mutual traps can highly decrease the simulation time required by the folding of strands into target structures (like DNA origami or DNA constructs).

<div style="text-align:center">
http://kratos.phys.uniroma1.it/histo_TRAP.png

<h5>Figure 3: Same as Figure 1, but for the system with mutual traps. Time is expressed in Monte Carlo sweeps.</h5>
</div>

Hairpin formation

2012-04-19T15:13:19Z

Rovigatti: /* Results */

[[Category:Examples]]
== Introduction ==
In this example, you will simulate a single strand of length 18 and sequence <tt>GCGTTGCTTCTCCAACGC</tt> at 334 K (~61 °C) in three different ways:
* with a molecular dynamics (MD) simulation of the sequence-averaged (SA) model. The input file is <tt>inputMD</tt>.
* with an MD simulation of the sequence-dependent (SD) model. The input file is <code>inputMD_seq_dep</code>.
* with a Monte Carlo (MC) simulation of the SA model in which two base pairs are connected by mutual traps (i.e. additional attractive interactions between two nucleotides). The input file is <tt>inputTRAP</tt>.
The traps act between the pairs depicted in blue and red in the sequence G<span style="color:red">C</span><span style="color:blue">G</span>TTGCTTCTCCAA<span style="color:blue">C</span><span style="color:red">G</span>C. The details of the interaction associated to the traps can be changed in the file <tt>hairpin_forces.dat</tt>.

This strand, if ''T'' is sufficiently low, tends to form an hairpin with a 6-base long stem and a 6-base long loop.
The temperature has been chosen to be close to the melting temperature of such a hairpin in the SA version of the model

This document explains how to prepare the ''hairpin'' example (see [[#Preparation|Preparation]]) and how to run it ([[#Running|Running]]). Section [[#Results|Results]] contains results and plots extracted from the simulation output. In the following, <tt>$EXEC</tt> refers to the oxDNA executable.

== Preparation ==
The script <tt>run.sh</tt> generates the input files and runs all the three simulations, one after the other. With the default input files, each simulation, lasting <math>10^8</math> steps by default, takes approximately one hour on a modern CPU. The default <tt>run.sh</tt> expects <tt>$EXEC</tt> to be in the <tt>../..</tt> directory. If this is not the case, open <tt>run.sh</tt> and change the variable <tt>CODEDIR</tt> accordingly.

If you only want to generate the initial configuration, you can issue <tt>./run.sh --generate-only</tt>. Then you can run the simulations by yourself. The generated initial configuration files are <tt>initial.top</tt> (which contains the topology) and <tt>initial.conf</tt> (which contains positions and orientations of the nucleotides).

== Running ==
In order to run the whole example, it is sufficient to issue the command <tt>./run.sh</tt> (or <tt>bash run.sh</tt>). As described in [[#Preparation|Preparation]], you can generate the initial configuration and then run the simulations by hand. The three simulations, described in [[#Introduction|Introduction]], can be performed issuing <tt>$EXEC input</tt>, where <tt>input</tt> is a text file that specifies the simulation configuration. For this example, three files have been prepared: <tt>inputMD</tt>, <tt>inputMD_seq_dep</tt> and <tt>inputTRAP</tt>. Table~\ref{tbl:sim} report all the files associated to each simulations.

{|
|-
! Type
! Input
! energy file
! trajectory file
! last configuration file
! log file
|-
|SA model
|inputMD
|energy.dat
|trajectory.dat
|last_conf.dat
|log.dat
|-
|SD model
|inputMD_seq_dep
|energy_seq_dep.dat
|trajectory_seq_dep.dat
|last_conf_seq_dep.dat
|log_seq_dep.dat
|-
|SA model with traps
|inputTRAP
|energy_trap.dat
|trajectory_trap.dat
|last_conf_trap.dat
|log_trap.dat
|-
|}

== Results ==
As mentioned in the [[#Introduction|Introduction]], temperature has been chosen so that the hairpin is near its melting temperature when simulated via SA model. The figure below shows the probability distribution histogram $P(U_{HB})$ for the hydrogen-bonding (HB) energy $U_{HB}$ of the system (in simulation units~\cite{ouldridge_jcp}). Both the minimum in $P(U_{HB})$ and a direct inspection of the time series (inset) show that the hairpin is indeed near the melting temperature, since $U_{HB}$ oscillates between $0$ (no HBs) and ~ -3.5, which corresponds to 4-5 hydrogen-bonded nucleotides. At this temperature the fraying effect is relevant, and therefore the last base-pair is not always formed~\cite{ouldridge_jcp}. This is clearly visible in the snapshot at the beginning of this document, in which $4$ base-pairs are formed.

http://kratos.phys.uniroma1.it/histo_MD.png

The figure below shows $P(U_{HB})$ and $U_{HB}$ for the SD hairpin. In this case, one expects a different melting temperature. Indeed, at this temperature the hairpin is nearly always in its folded conformation and the average value of $U_{HB}$ is lower than in the SA case. Although stacking interactions play a role, the major contribution to this effect is due to the presence of four <tt>GC</tt> only two <tt>AT</tt> base-pairs in the stem.

http://kratos.phys.uniroma1.it/histo_MD_seq_dep.png

If mutual traps between stem base-pairs are introduced, then the equilibrium properties of the hairpin are changed and, even if the SA model is employed, the hairpin is always (after the initial equilibration) in its folded conformation. The use of mutual traps can highly decrease the simulation time required by the folding of strands into target structures (like DNA origami or DNA constructs).

http://kratos.phys.uniroma1.it/histo_TRAP.png

File:Histo MD.png

2012-04-19T10:55:41Z

Rovigatti: uploaded a new version of "Image:Histo MD.png"