Coarse-Grained Simulation of DNA using LAMMPS
Oliver Henrich, Yair Augusto Gutierrez-Fosado, Tine Curk, Thomas E. Ouldridge
CCoarse-Grained Simulation of DNA using LAMMPS
An implementation of the oxDNA model and its applications
Oliver Henrich ∗ Department of Physics, SUPA, University of Strathclyde, Glasgow G4 0NG, Scotland, UK ∗ Yair Augusto Guti´errez Fosado
School of Physics and Astronomy, University of Edinburgh, Edinburgh EH9 3FD, Scotland, UK
Tine Curk
CAS Key Laboratory of Soft Matter Physics, Beijing National Laboratory for Condensed Matter Physics,Institute of Physics, Chinese Academy of Sciences, Beijing 100190, China &Department of Chemistry, University of Cambridge, Cambridge CB2 1EW, UK
Thomas E. Ouldridge
Department of Bioengineering & Centre of Synthetic Biology,Imperial College London, London SW7 2AZ, UK
During the last decade coarse-grained nucleotide models have emerged that allow us to studyDNA and RNA on unprecedented time and length scales. Among them is oxDNA, a coarse-grained,sequence-specific model that captures the hybridisation transition of DNA and many structuralproperties of single- and double-stranded DNA. oxDNA was previously only available as standalonesoftware, but has now been implemented into the popular LAMMPS molecular dynamics code.This article describes the new implementation and analyses its parallel performance. Practicalapplications are presented that focus on single-stranded DNA, an area of research which has been sofar under-investigated. The LAMMPS implementation of oxDNA lowers the entry barrier for usingthe oxDNA model significantly, facilitates future code development and interfacing with existingLAMMPS functionality as well as other coarse-grained and atomistic DNA models.
I. INTRODUCTION
DNA is one of the most important bio-polymers, asits sequence encodes the genetic instructions neededin the development and functioning of many living or-ganisms. While we know now the sequence of manygenomes, we still know little as to how DNA is organ-ised in 3D inside a living cell, and of how gene reg-ulation and DNA function are coupled to this struc-ture. The complexity of the DNA molecule can bebrought to mind by highlighting a few of its quantita-tive aspects. The entire DNA within a single humancell is about 2 m long, but only 2 nm wide and organ-ised at different hierarchical levels. If compressed intoa spherical ball, this ball would have a diameter ofabout 2 µ m [1].Computational modelling of DNA appears as theonly avenue to understanding its intricacies in suf-ficient detail and has been an important field inbiophysics for decades. Traditionally, most of theavailable simulation techniques have worked at theatomistic level of detail [2]. Existing atomistic forcefields can capture fast conformational fluctuations andprotein-DNA binding, but cannot deliver the neces-sary temporal and spatial resolution to describe phe-nomena that occur on larger time and length scales asthey are often limited to a few hundred base pairs and(at most) microsecond time scales. Recent years havetherefore witnessed a rapid increase of a new research ∗ Corresponding email address: [email protected] effort at a different, coarse-grained level [3]. Coarse-grained (CG) models of DNA can provide significantcomputational and conceptual advantages over atom-istic models leading often to three or more ordersof magnitude greater efficiency. The challenge con-sists in retaining the right degrees of freedom so thatthe CG model reproduces relevant emergent structuralfeatures and thermodynamic properties of DNA. CGmodelling of DNA is not only an efficient alternativeto atomistic approaches. It is indispensable for themodelling of DNA on timescales in the millisecondrange and beyond, or when long DNA strands of tensof thousands of base pairs or more have to be consid-ered, e.g. to study the dynamics of DNA supercoiling(i.e. the local over- or under-twisting of the doublehelix, which is also important for gene expression inbacteria), of genomic DNA loops and of chromatin orchromosome fragments.A small number of very promising CG DNA modelshave emerged to date. Conceptually they can be cate-gorised into top-down approaches, which use empiricalinteractions that are parameterised to match experi-mental observables, or bottom-up approaches, whicheliminate dispensable degrees of freedom systemati-cally starting from atomistic force fields. They mayalso target different applications depending on theircapabilities, such as single versus double strandedDNA (ssDNA and dsDNA), or nanotechnological ver-sus biological applications. We refer to [4] for a com-prehensive overview of the capabilities of individualmodels and recent activities in this field.From a software point of view these models are often a r X i v : . [ c ond - m a t . s o f t ] M a y based on standalone software [5–7], which has a some-what limiting effect on uptake and user communitiesgrowth. Others models use popular MD-codes as com-putational platforms, such as GROMACS [8] in caseof the SIRAH [9] and the MARTINI force field [10],or NAMD [11, 12]. Another suitable platform for CGsimulation of DNA has emerged in form of the power-ful Large-Scale Atomic/Molecular Massively ParallelSimulator (LAMMPS) for molecular dynamics [13],including the widely used 3SPN.2 model [14, 15] andothers that target even larger length scales [16, 17].This article reports the latest effort of implementingthe popular oxDNA model [18, 19] into the LAMMPScode. Until recently this model was only available asbespoke and standalone software [20]. Through theefficient parallelisation of LAMMPS it is now pos-sible to run oxDNA in parallel on multi-core CPU-architectures, extending its capabilities to unprece-dented time and length scales. The largest systemthat could be studied by oxDNA was previously lim-ited by the size of system that can be fitted onto asingle GPU.This paper is organised as follows: In SectionII we briefly introduce the details of the oxDNAand oxDNA2 models. Section III explains how theLAMMPS implementation of the oxDNA models canbe invoked and provides further information on thecode distribution and documentation. In Section IVwe describe the LAMMPS implementation of novelLangevin-type rigid body integrators which featureimproved stability and accuracy. Section V gives de-tails of the scaling performance of parallel implemen-tation. Section VI presents results on the behaviour ofsingle-stranded DNA, an area of DNA research whichso far has not been intensively investigated. One ap-plication is concerned with lambda-DNA of a bacterio-phage, whereas the other application involves a plas-mid cloning vector pUC19. In Section VII we sum-marise this work. II. THE OXDNA MODEL
The oxDNA model consists of rigid nucleotides withthree interaction sites for the effective interactions be-tween the nucleotides. These pairwise-additive forcesarise due to the excluded volume, the connectivity ofthe phosphate backbone, the stacking, cross-stackingand coaxial stacking as a consequence of the hy-drophobicity of the bases, as well as hydrogen bond-ing between complementary base pairs. Fig. 1 illus-trates these interactions schematically for the originalversion of the model, to which we refer as oxDNA[19]. In this version all three interaction sites are co-linear. The hydrogen bonding/excluded volume siteand the stacking site are separated from the back-bone/electrostatic interaction site by 0 .
74 length units(6 . . . V coaxial stack V backbone V stack V H.B. V cross stack FIG. 1. Overview of bonded and pair interactions inoxDNA: phosphate backbone connectivity and excludedvolume, hydrogen-bonding, stacking, cross-stacking andcoaxial stacking interaction. The oxDNA2 model containsan additional implicit electrostatic interaction in form of aDebye-H¨uckel potential. Reprinted from [21] with permis-sion from ACS Nano. Copyright (2013) American Chemi-cal Society. with the relative distance vectors between the interac-tion sites, the base vector and base normal vector areused to modulate the stacking, cross-stacking, coaxialstacking and hydrogen bonding interaction betweentwo consecutive nucleotides.The simplest interaction is the backbone connectiv-ity, which is modelled with FENE (finitely extensiblenonlinear elastic) springs acting between the backboneinteraction sites. The excluded volume interactionis modelled with truncated and smoothed Lennard-Jones potentials between backbone sites, base sitesand between the backbone and base sites. The hydro-gen bonding interaction consists of smoothed, trun-cated and modulated Morse potentials between thehydrogen bonding site. The stacking interaction fallsinto three individual sub-interactions: the stacking in-teraction between consecutive nucleotides on the samestrand as well as cross-stacking and coaxial stackingbetween any nucleotide in the appropriate relative po-sition. It is worth emphasising that the duplex struc-ture is not specified or imposed in any other way,but emerges naturally through this choice of inter-actions and their parameterisation. This is anotherstrength of the oxDNA model and permits an accuratedescription of both ssDNA and dsDNA. The stack-ing interactions are modelled with a combination ofsmoothed, truncated and modulated Morse, harmonicangle and harmonic distance potentials. All interac-tions have been parameterised to match key thermo-dynamic properties of ssDNA and dsDNA such as thelongitudinal and torsional persistence length or themelting temperature of the duplex [18, 22, 23].A short schematic overview of various interactionsinvolved in the definition of oxDNA model is givenin Fig. 1. More details can be found in the originalpublications [18, 19].
FIG. 2. (a): Schematic distinction between oxDNA (left)and oxDNA2 (right). In oxDNA all interaction sites areco-linear whereas in oxDNA2 the backbone interaction siteand the stacking and hydrogen-bonding interaction sitesare oriented at an angle. (b): The non-co-linear arrange-ment of the interaction sites leads to the formation of themajor and minor groove, an important structural featureof DNA. Reproduced from [24], with the permission of AIPPublishing.
The original model (oxDNA) has been further de-veloped to include sequence-specific stacking and hy-drogen bonding interaction strengths [25] (oxDNA1.5)and implicit ions, which are modelled by means of aDebye-H¨uckel potential [24] (oxDNA2). A major im-provement of the latest version is also the fact thatit shows the correct structure with major and minorgrooves (see Fig. 2 (b)). This is achieved through amodification of the relative position of the backboneand stacking/hydrogen bonding interaction sites, asschematically depicted in Fig. 2 (a).
III. THE LAMMPS IMPLEMENTATION OFOXDNAA. Code Distribution, Force Fields andCompilation
The software is open source and distributed underGNU General Public License (GPL). It is availablefor download as LAMMPS USER-package from thecentral LAMMPS repository at Sandia National Lab-oratories, USA [13]. This includes a detailed onlinedocumentation, examples and utility scripts. We re-fer also to these materials for a general introductioninto the usage of LAMMPS.To compile the code, load the LAMMPS standardpackages
MOLECULE and
ASPHERE and the
USER-CGDNA package by issuing make yes-molecule yes-asphere yes-user-cgdna in the main source code directory and compile asusual.All three versions oxDNA, oxDNA1.5 and oxDNA2are implemented in the LAMMPS code and canbe invoked through appropriate keywords in the in-put file. This allows for instance to run with-out sequence-specific interactions and without im-plicit ions (oxDNA force field and keyword seqav ≡ oxDNA), with sequence-specific interactions andwithout implicit ions (oxDNA force field and keyword seqdep ≡ oxDNA1.5) or with implicit ions and with orwithout sequence-specific interactions (oxDNA2 forcefield and keywords seqdep or seqav , respectively).The source code is also distributed via our mainrepository at CCPForge [26] under the project name Coarse-Grained DNA Simulation (cgdna) . Pleasesend a request to join the project for full access thatincludes permission to browse the repository and com-mit changes.
B. Force and Torque Calculation
Integrating the equations of motion of rigid bod-ies requires accurate information of their relative ori-entations. In simple situations this can be achievedthrough Euler angles, which describe the orientationof a rigid body and its local reference frame with re-spect to the laboratory system. Euler angles havethe disadvantage that they are not unambiguously de-fined as a singularity arises when two rotation axes fallparallel. This situation, usually referred to as gimballock, arises easily in a system that contains a largenumber of rigid bodies. Unsurprisingly, it triggers nu-merical instabilities, which is why rigid body problemsare best formulated by means of quaternions [27] in-stead of Euler angles.Computationally it is most efficient to integrate thequaternion degrees of freedom directly via a gener-alised 4-component quaternion torque (see [19] for adetailed derivation of the oxDNA forces and gener-alised 4-torques using quaternion dynamics). Unfor-tunately such an interface for generalised quaterniontorques and momenta is not provided in LAMMPS.It expects for its rigid body integrators 3-componenttorques and angular momenta as input quantities (be-sides the Newtonian force for the integration of thecoordinate degrees of freedom). To be consistent andsimplify interfacing with existing functionality, we de-cided to adhere to this convention. This, however, en-tails conversion of the unit quaternions into Cartesianunit vectors of a body frame before forces and torquescan be calculated for the integration step, thus leadingto a computational overhead (see Appendix A).Once this choice has been made, the calculation ofthe forces and torques is most conveniently formulatedfollowing Ref. [28]. If ˆ a and ˆ b are the principal axesof two rigid bodies A and B and r is the norm of therelative distance vector r = r A − r B from B to A,then the pair potential depends on a combination ofthese quantities, U = U ( r, ˆ a , ˆ b ) = U ( r, { ˆ a m · ˆ r } , { ˆ b n · ˆ r } , { ˆ a m · ˆ b n } ) (1)where ˆ r , ˆ a m and ˆ b n are the normalised relative dis-tance and orthonormal principal axes vectors. Fromthis definition the force on A due to B are straightfor-wardly written as F A = − F B = − ∂U∂ r = − ∂U∂r ˆ r − r − (cid:88) m (cid:34) ∂U∂ (ˆ a m · r ) ˆ a ⊥ m + ∂U∂ (ˆ b m · r ) ˆ b ⊥ m (cid:35) . (2)Here ˆ a ⊥ m = ˆ a m − (ˆ a m · ˆ r )ˆ r denotes the componentof ˆ a m which is perpendicular to ˆ r . The torques areslightly more involved: τ A = (cid:88) m ∂U∂ (ˆ a m · r ) (ˆ r × ˆ a m ) − (cid:88) mn ∂U∂ (ˆ a m · ˆ b n ) (ˆ a m × ˆ b n ) (3) τ B = (cid:88) n ∂U∂ (ˆ b n · r ) (ˆ r × ˆ b n )+ (cid:88) mn ∂U∂ (ˆ a m · ˆ b n ) (ˆ a m × ˆ b n ) . (4)The fact that local angular momentum conservationrequires τ A + τ B + r × f = 0 (5)can be conveniently utilised for debugging and veri-fication purposes. The implementation was verifiedagainst two independent implementations, namelyOuldridge’s own code, which is based on quaterniondynamics [19] as well as the standalone oxDNA code[20], which makes also use of the same scheme for theforce and torque calculation. To this end two bench-marks were studied, a 5-base-pair duplex and a 8-basepair nicked duplex, which are both provided as exam-ples in the USER-CGDNA package. C. Input File
In the following we discuss the structure of the inputfile and how the newly introduced oxDNA classes areinvoked.We work with Lennard-Jones reduced units, which areinvoked in LAMMPS via units lj
The system is three-dimensional. dimension 3
In LAMMPS, an oxDNA nucleotide is represented asa bonded-ellipsoidal hybrid particle with the associ-ated degrees of freedom of bonded particles in a bead-spring polymer (backbone connectivity) and aspheri-cal particles with shape (moment of inertia), quater-nion (orientation) and angular momentum. atom style hybrid bond ellipsoid
Users are required to suppress the atom sorting algo-rithm as this can lead to problems in the bond topol-ogy of the DNA. atom modify sort 0 1.0
It is important to set the skin size correctly, which con-trols the extent of the neighbour lists. Too large a skinsize and neighbour lists become unnecessarily long,leading to superfluous communication. Too short andpartners in the pair interactions will be lost. neighbor 1.0 bin
A good way to fine-tune this parameter is to run anNVE simulation with constant energy before apply-ing Langevin integrators. We recommend neighbor2.0 bin as a safe starting point. Likewise, frequentupdate of the neighbour lists can lead to an undueperformance degradation. This parameter should betuned as well so that no dangerous builds (as reportedin the standard output of LAMMPS) occur. neigh modify every 1 delay 0 check yes
The initial configuration and topology is created bymeans of an external setup tool (see Sec. III D) andread in. read data data file name
All masses are set to 3 . set atom * mass 3.1575 Note that the moment of inertia is determined throughthe shape parameter in the data file (see below Sec.III D). There are four types of nucleotides (A=1, C=2,G=3, T=4), which are grouped together into a groupnamed all for the integration. group all type 1 4
The new oxDNA classes with its parameters are in-voked as follows: bond style oxdna2/fenebond coeff * 2.0 0.25 0.7564pair style hybrid/overlay oxdna2/excv &oxdna2/stk oxdna2/hbond oxdna2/xstk &oxdna2/coaxstk oxdna2/dhpair coeff * * oxdna2/excv 2.0 0.7 0.675 2.0& 0.515 0.5 2.0 0.33 0.32pair coeff * * oxdna2/stk seqdep 0.1 6.0 0.4& 0.9 0.32 0.6 1.3 0 0.8 0.9 0 0.95 0.9 0& 0.95 2.0 0.65 2.0 0.65pair coeff * * oxdna2/hbond seqdep 0.0 8.0 &0.4 0.75 0.34 0.7 1.5 0 0.7 1.5 0 0.71.5 &
Please note that according to the LAMMPS parsingrules the ampersands (&) represent line breaks.Visit the LAMMPS online documentation and manualfor more information and for information on oxDNA2.
D. Data File and Setup Tool
The data file contains all relevant structuralparameters for the simulation, i.e. details about thenumber of atoms, the topology of the molecules,the size of the simulation box, initial velocities, etc.The LAMMPS implementation of oxDNA follows thestandard form as discussed in the LAMMPS usermanual. We outline the relevant parts below.At the beginning of the data file the total numberof particles and bonds has to be given. As we areusing hybrid particles, we need to set the same numberof ellipsoids. For a standard DNA duplex consistingof 8 complementary base pairs we need 16 atoms, 16ellipsoids and 14 bonds, 7 on each of the two singlestrands. If the strands are nicked, which we do notassume here, the number of bonds would be reduced.
16 atoms16 ellipsoids14 bonds
We use four atom types to represent the four differentnucleotides in DNA (A=1, C=2, G=3, T=4). We use only one bond type.
The dimensions of the simulation box are defined asfollows: -20.0 20.0 xlo xhi-20.0 20.0 ylo yhi-20.0 20.0 zlo zhi
Although already stated in the input file, we need toprovide again the masses of the nucleotides.
Masses1 3.15752 3.15753 3.15754 3.1575
The nucleotides are defined after the keyword
Atoms .Each row contains the atom-ID (1,2,3 in the examplebelow), the atom type (1,1,4), the position (x,y,z), themolecule ID (all 1 in this case), an ellipsoidal flag (1)and a density (1).
Atoms1 1 0.00000 0.00000 0.00000 1 1 12 1 0.13274 -0.42913 0.37506 1 1 13 4 0.48461 -0.70835 0.75012 1 1 1...
Next we set the initial velocities to the desired value,here all equal to 0. The first column contains theatom-ID (1,2,3), the following three columns thetranslational, and the last three columns the angularvelocity.
Velocities1 0.0 0.0 0.0 0.0 0.0 0.02 0.0 0.0 0.0 0.0 0.0 0.03 0.0 0.0 0.0 0.0 0.0 0.0...
Note that this is our special choice in the setup tool.The velocities can be generally initialised to any value.Large values will lead to the FENE springs becomingoverstretched and may provoke an early abortion ofthe run.The ellipsoids are defined with atom-ID, shape(1.17398 to produce the correct moment of inertia)and initial quaternion (last four columns).
Ellipsoids1 1.17398 1.17398 1.17398 1.00000 0. 0. 0.2 1.17398 1.17398 1.17398 0.95534 0. 0.0.295523 1.17398 1.17398 1.17398 0.82534 0. 0.0.56464...
Finally, we specify the bond topology. The first col-umn contains the bond-ID (1,2,3), the second one thebond type (1) and the third and fourth the IDs of thetwo bond partners.
Bonds1 1 1 2
To simplify the setup procedure we provide a simplepython tool with the example and utility files of theUSER-CGDNA package. The script allows the userto create single- and double-stranded DNA from aninput file that specifies the sequence and requires aninstallation of numpy .The syntax is very straightforward, but the systemsize has to be specified in the following way: $> python generate.py
ACGTA
If the sequence is prepended by the keyword
DOUBLE ,then a single, helical DNA duplex is created. Thebases on the second strand are complementary tothose on the first strand, which is given in the se-quence input file:
DOUBLE ACGTA
Consecutive strands are positioned and oriented ran-domly without creating any overlap in case of morethan one ssDNA or dsDNA strand. Note that theprocedure works only below a critical density as thissimple script does not feature cell lists. Besides thesesetup tools, the USER-CGDNA package contains aswell example input, data and standard output files ofshort benchmark runs of dsDNA duplexes.
E. Output and Visualisation
LAMMPS offers a multitude of possible output for-mats, including parallel HDF5 and NetCDF formats,VTK format or very basic standard trajectory data.We will summarise here how output of basic observ-ables of the oxDNA model can be invoked in the inputfile.The xyz style writes XYZ files, which is a simpletext-based coordinate format that many codes canread, which has one line per atom with the atom typeand the x-, y-, and z-coordinate of that atom. Thisstyle is invoked via dump 1 all xyz Nint trajectory.xyz where
Nint is the output frequency in timesteps. Ad-ditional output of e.g. velocity, force and torque on aper-atom basis makes some customisation necessary, dump 2 all custom Nint filename.dat id x y z& vx vy vz fx fy fz tqx tqy tqz where id is the unique atom-ID. The output ofquaternions requires a so-called compute style. Theresult of the compute style can then be retrieved inthe following way: compute quat all property/atom quatw quati &quatj quatkdump 3 all custom Nint filename.dat id &c quat[1] c quat[2] c quat[3] c quat[4] Another observable that may be of interest is theenergy, or more specifically broken down into rota-tional, kinetic and potential energy. This is also donethrough a compute style. compute erot all erotate/aspherecompute ekin all kecompute epot all pevariable erot equal c erotvariable ekin equal c ekinvariable epot equal c epotvariable etot equal c erot+c ekin+c epot
Note that the somewhat simpler thermo style command for output discards the kinetic energy ofrotation when the kinetic energy is requested.LAMMPS does not contain a direct visualisationtoolkit. There are, however, a multitude of ways howsnapshots can be visualised. ParaView [29] for in-stance, is an open source, multi-platform data analy-sis and visualisation application. The images in thiswork have been generated with the molecular visual-isation program VMD (Visual Molecular Dynamics)[30]. More information about possible visualisationpipelines can be found in the LAMMPS online man-ual [13].
IV. LANGEVIN-TYPE RIGID-BODYINTEGRATORS
Together with the USER-CGDNA package comesalso an implementation of novel Langevin-type rigid-body integrators that were developed by Davidchack,Ouldridge and Tretyakov [31]. The motivation for thiswas that previously only a limited choice of suitableLangevin integrators for rigid bodies was available inLAMMPS. Without noise all integrators A, B and C inthe above reference are identical and basically equiv-alent to the integrator presented by Miller et al. [32].Nevertheless, we refer to this case as the “DOT in-tegrator” (the other implementation of the Miller in-tegrator is only available when using the fix rigid command in LAMMPS). The DOT integrator is an al-ternative to the standard LAMMPS NVE integratorfor aspherical particles, and can be invoked by replac-ing the standard choice fix 1 all nve/asphere with fix 1 all nve/dot in the input file. This energy-conserving integrator isuseful for an analysis of the accuracy of this family ofintegrators or the integrity of the pair interactions ata given timestep size ∆ t .The C integrator in Ref. [31], to which we referas “DOT-C integrator”, is invoked by replacing thestandard NVE integrator for aspherical particles andthe fix for Langevin dynamics fix 1 all nve/aspherefix 2 all langevin 0.1 0.1 0.03 457145angmom 10 with one single fix fix 1 all nve/dotc/langevin 0.1 0.1 0.03 &457145 angmom 10 To measure the accuracy of the new integrators,we run a test case consisting of a short, nicked du-plex with 8 base pairs (16 nucleotides). Fig. 3 showsthe accuracy measured through the normalised dif-ference between the total energy E tot for this partic-ular benchmark and the total energy at the begin-ning of the run E ∗ tot . We compared the standard fixnve/asphere integrator, which is based on a Richard-son iteration in the update of the quaternion degreesof freedom, to the new DOT integrator, which uses arotation sequence to update the quaternions. Shownare results for two different timestep sizes ∆ t = 10 − and ∆ t = 10 − . Both simulations were run for thesame physical simulation time to allow direct com-parison of the deviations of a dynamical run. As thisis done in the NVE ensemble and without noise, theenergy should be exactly conserved. This correspondsto a straight, horizontal line at 0. -4-2 0 2 4 6 8 10 12 14 16 0 2000 4000 6000 8000 10000 ( E t o t / E * t o t - ) · - τ (LJ time units) nve/asphere dt=1e-3nve/asphere dt=1e-4nve/dot dt=1e-3nve/dot dt=1e-4 FIG. 3. Relative normalised accuracy ( E tot − E ∗ tot ) /E ∗ tot of the standard LAMMPS NVE integrator for asphericalparticles and the NVE DOT integrator from Ref. [31]. E ∗ tot is the total free energy at the beginning of the simu-lation runs. It is obvious that above a certain timestep size theaccuracy of the new DOT integrator is slightly infe-rior compared to the standard integrator. Up to acertain point the DOT integrator actually seems todeviate further from the correct result, whereas the standard integrator fluctuates more around the cor-rect value. This, however, is more or less a transienteffect as longer runs show there is no permanent driftaway from the correct result. For Langevin dynam-ics, it is not possible to evaluate the accuracy andstability in the same way. We opted instead for anestimate based on the average kinetic, rotational, po-tential and total energy of the benchmark. Again,we performed runs of τ = 10000 Lennard-Jones timeunits length, this time thermalised, and averaged theresults over the time interval. The number of MD-timesteps and the output frequency for each timestepsize were adapted so that the total physical simulationtime and the statistical basis of the error calculationswere consistent. The temperature in reduced LJ-unitswas set to T = 0 .
1, whereas the translational and ro-tational friction or damping coefficients were set to γ = 1 / .
03 and Γ = 1 / .
3, respectively. The resultsare summarised in Tab. I. These values were used dur-ing the verification of the LAMMPS implementationbecause they produced relatively smooth trajectoriesthat could be easily followed. For actual productionruns it may be more appropriate to use different val-ues to allow a better and more efficient sampling ofthe configuration space.Based on three translational and three rotationaldegrees of freedom per nucleotide and 8 base pairs weexpect kinetic and rotational energies E kin = E rot =2 . T = 0 .
1. This isvery well achieved for all timestep sizes and bothintegrators, the standard LAMMPS integrator fixnve/asphere & fix langevin and the DOT-C inte-grator fix nve/dotc/langevin . However, there ap-pears to be a slight decrease in the DOT-C integratorfor very large step sizes (∆ t = 2 · − ). The deviationof the total energy between all timestep sizes, admit-tedly an ad hoc criterion to quantify the stability ofthe integrators, but one that is rather hard for theintegrators to get exactly right, is in the sub-percentrange. It is actually slightly better for the DOT-C in-tegrator than for the standard LAMMPS integrator.The statistical errors, reported in Tab.I, are the stan-dard deviations of a linear least square fit and showthat the deviations are well above the uncertainty ofthe fits.Remarkably, for the DOT-C integrator the limit fora stable integration is ∆ t = 2 · − , which representsa very large timestep size. This is about 4 times largerthan the maximum timestep size for which the stan-dard LAMMPS Langevin integrator produces soundresults. Because of the more complex rotations inquaternion space and various additional transforma-tions that the DOT-C integrator requires there is asmall overhead of about 15% compared to the stan-dard LAMMPS integrator. Nevertheless, this smalloverhead of the DOT-C integrator is very well com-pensated by the computational efficiency and possi-bility to increase the timestep size by 400% (from amaximum of ∆ t = 5 · − for the standard LAMMPSintegrator to ∆ t = 2 · − for the DOT-C integrator). ∆ t E kin E rot E pot E tot standard error of E tot fitfix nve/asphere &fix langevin 10 − ± − ± · − ± − ± − ± − ± · − ± fix nve/asphere & fix langevin and the DOT-C integrator nve/dotc/langevin for different timestep sizes. V. PERFORMANCE ANALYSIS
We devised a few simple benchmarks to study theparallel performance of the LAMMPS implementa-tion. The size of each benchmark is well beyondthe current capabilities of the standalone version, soeach demonstrates as well a minimal performance re-quirement. The benchmarks consisted of arrays of
FIG. 4. The low-density benchmark consisting of a 10 × ×
40 arrayof duplexes with 960 kbp in total. The pictures show thefinal configuration the end of a performance run and wereproduced with VMD. The centre of mass of each nucleotideis represented through a sphere. double-stranded, regularly arranged DNA duplexes,each with a length of 600 base pairs. The low-density(LD) benchmark was formed by a 10 ×
10 array of du-plexes, giving a total of 60 kbp, and is shown in Fig.4. The high-density (HD) benchmark was formed bya 40 ×
40 array of duplexes with a density 16 timeslarger than the LD case and a total number of 960kbp. Whilst a regular array of double-stranded DNAstrands appears perhaps somewhat artificial, it createsa reasonably load-balanced situation and facilitatesthe performance analysis. The obtained densities ofDNA, are however very well comparable to those ofDNA gels [33] and high density states of DNA which
100 1000 10000 S p ee dup v s M P I- t a s k s MPI-tasks Ideal performance60 kbp960 kbp 0 1
100 1000 10000 P a r a ll e l E ff i c i e n c y MPI-tasks
FIG. 5. Strong scaling behaviour: Speedup of the low andhigh density benchmarks of 60 kbp and 960 kbp, respec-tively, compared to the single node performance with 24MPI-tasks. The inset shows the parallel efficiency relativeto the single node case with 24 MPI-tasks. form liquid-crystalline phases [34].Strong scaling tests were performed on ARCHERon up to 86 nodes (LD) and 683 nodes (HD), respec-tively. The benchmark cases were run for 30,000 (LD)and 10,000 (HD) MD-timesteps with a timestep sizeof ∆ t = 5 × − . We used the standard LAMMPSintegrators for Langevin dynamics, although the scal-ing behaviour was found to be virtually identical whenusing the above described rigid body integrator DOT-C. The primary reason for this was that the wallclocktime for runs with the standard integrator was still afew percent shorter, although the improved efficiencyof the DOT-C integrator would mean these runs wereshorter in physical time. The temperature in reducedLJ-units was T = 0 .
1, whereas the translational androtational friction coefficients were set to γ = 1 / . / .
3, respectively.Fig. 5 shows the parallel speedup for both bench-marks relative to the single node performance with24 MPI-tasks. The code performs well for the LDbenchmark up to about 128 MPI-tasks with a parallelefficiency around 95% (see the inset). Beyond sev-eral hundred MPI-tasks a gradual performance degra-dation is observed. At 2048 MPI-tasks the parallelefficiency has decreased to about 45% and the totalspeedup is roughly 930-fold compared to the singlecore performance (39-fold compared to the single nodeperformance).A look at the ratio of the number of local atoms,i.e. those that are inside a process boundary, to thenumber of ghost atoms, i.e. those which need to becommunicated via neighbour lists, proves that the ob-served performance degradation is due to the com-parably small size of the problem. At the largestcore counts there are on average only about 60 localatoms present on each process, whereas the numberof ghost atoms is with about 225 atoms almost fourtimes larger. LAMMPS is known to require at least afew hundred local atoms or more for a good parallelperformance [35]. The speedup is still relatively goodbecause the fraction of time that the algorithm spendsin the force calculation is still comparably large. Forthe HD benchmark, 16 times larger than the LD case,the performance degradation is more or less mirroredat core counts that are about 16 times larger. Forthe HD benchmark the total speedup at 16384 MPI-tasks is 9680-fold with respect to the single core per-formance (400-fold compared to the single node per-formance) and the parallel efficiency is still at around60%.These two examples are of course slightly idealisedin the sense that both benchmarks fulfil easily therequirement of good load-balancing, which is neces-sary to obtain a good scaling performance. LAMMPS,however, features sophisticated load-balancing algo-rithms which permit good scaling behaviour also forvery inhomogeneous systems. We are planning to ex-tend the existing implementation to benefit furtherfrom recent developments pertaining to threaded par-allelisation on shared memory architectures such asmany-core chips and general purpose graphical pro-cessing units (GPGPUs).One of the major advantages of the new LAMMPSimplementation is that it can be directly comparedwith other coarse-grained models that are also basedon the LAMMPS code. To this end, we comparedthe single core performance of oxDNA2 with that of3SPN.2 [14]. The benchmark consisted of two comple-mentary dsDNA duplexes of 8 bps with implicit ions.In order to compare both models we set the transla-tional friction coefficient γ to about (300 fs) − . Weopted for the maximum timestep size that provideda stable integration, which was ∆ t = 35 fs (3SPN.2)and ∆ t = 48 fs (oxDNA2 + DOT-C integrator), re-spectively.On a single Intel Core i7 2.8 GHz processor us-ing the latest version of LAMMPS (16 March 2018)3SPN.2 delivered a performance of about 60 µ s perday. oxDNA2 was able to surpass this by about a factor 1.6 with a performance of roughly 100 µ s perday. Note that comparing the wall times is only anapproximate way to compare the performance as thereis no guarantee that similar processes take a similarsimulation time in the two models.Apart from the enhanced stability of the rigidbody integrator, this difference in performance will becaused by the different number of degrees of freedomthat both models require: oxDNA/oxDNA2 uses only13 degrees of freedom per nucleotide (3 coordinate po-sitions, 3 translational momenta, 3 angular momentaand 4 quaternion degrees of freedom), whereas 3SPN.2uses 18 degrees of freedom per nucleotide (3 particleswith each 3 coordinate positions and translational mo-menta).Unfortunately, we could not measure the parallelperformance of 3SPN.2. But this conceptual differ-ence between the two models is very likely to entailfurther detrimental effects when running in parallel.With the larger number of degrees of freedom pernucleotide in 3SPN.2, communication overheads arelikely to build up more quickly and neighbour listsare longer and probably have to be rebuilt more fre-quently. On the other hand, the current LAMMPSimplementation of oxDNA offers further potential foroptimisation as it spends a good part its time com-puting the inverse cosine (around 12%, see AppendixA). This could be alleviated for instance through theintroduction of appropriate lookup tables for trigono-metric functions. VI. APPLICATIONS
The structural properties of DNA such as the persis-tence length, radius of gyration and torsional rigidityplay an important role in its function. Characteris-ing these properties and their dependence on differentconditions is therefore fundamental for highly complexprocesses such as DNA packaging, replication and de-naturation. Experimentally, however, making thesemeasurements is not an easy task as it requires subtlemanipulation of single molecules and direct measure-ment of their response to applied forces or displace-ments, which can then be related to the elasticity ofDNA. By using coarse-grained computational modelslike oxDNA, we can study these systems in more de-tail. These simulations can in turn provide insightsinto experimental data or the performance of othertheoretical approaches.The radius of gyration is a particularly useful de-scriptor of the structure and compactness of macro-molecules. For ssDNA the radius of gyration R g canbe defined as R g = 1 N N (cid:88) i =1 ( r i − ¯ r ) (6)where N is the number of nucleotides, r i is the po-sition of the i − th nucleotide and ¯ r = N (cid:80) Ni =1 r i is0the mean position of the ssDNA strand. For dsDNAthis definition would be modified to use the centre ofmass coordinate of a base pair (bp) and N would bereplaced with the number of base pairs.In this section we present results obtained with theoxDNA2 model for two different systems: a sequenceof ssDNA from a λ -bacteriophage that has a multi-tude of applications in microbial and molecular ge-netics and serves e.g. as cloning vector, as well ascomplete ssDNA sequence of the pUC19 plasmid, an-other model organism and cloning vector, which con-veys antibiotic resistance.We performed Langevin dynamics simulations ofthe two above mentioned ssDNA sequences at a con-stant salt concentration of 0 .
2M NaCl. For simplicitywe used linear DNA molecules, so their ends are freelyto rotate. After a sudden quench in temperature, thesystem evolved from a random initial configurationtowards a new steady state. The criterion for reach-ing this steady state was a constant radius of gyra-tion R g and number of base-pairs N c formed alongthe chain. Equilibrium values for these observableswere obtained by averaging five different configura-tions over the last 3 × τ LJ timesteps.
10 12 14 16 18 20 22 24 0 10 20 30 40 50 60 R g [ n m ] Temperature [ ° C] R g λ -ssDNApoly-A ssDNApoly-A av. stacking strengthpoly-T ssDNApoly-T av. stacking strength FIG. 6. Response of the radius of gyration R g of thessDNA λ -DNA sequence to temperature changes. Pointsshow R g computed from averages over various configura-tions of a 500 base pair (bp) long ssDNA chain using theoxDNA2 model with sequence-specific stacking strength(red full circles), poly-A (green open circles) and poly-T(magenta open squares). These results are compared tothose for poly-A and poly-T chains with sequence-averagedstacking strength, respectively (blue full squares and cyancrosses). In Fig.6 the initial 500 nucleotide long sequence ofssDNA λ -DNA is compared with different linear DNAmolecules of the same length, namely poly-A and poly-T strands. The radius of gyration as a function oftemperature is shown. For λ -DNA we observe that R g increases with temperature until a plateau is reachedat around 50 ◦ C. While the λ -DNA sequence allowshybridisation along the ssDNA (see Fig.7), the same
10 11 12 13 14 15 16 17 0 10 20 30 40 50 60 0 0.1 0.2 0.3 0.4 0.5 R g [ n m ] C on t ac t F r ac ti on f h Temperature [ ° C] Contact Fraction λ -ssDNAR g λ -ssDNA FIG. 7. Dependence of the radius of gyration R g and thefraction of direct intra-chain nucleotide contacts in the λ -ssDNA sequence on the temperature. is not true for poly-A or ploy-T sequences. This canexplain the differences in R g between the two thatwe observe at low temperatures. In contrast, poly-A shows the opposite tendency, with the largest R g at the lowest temperature setting of 0 ◦ C. The reasonfor this different behaviour is the roughly 16% largerstacking strength between consecutive A nucleotidesas compared to T nucleotides, an explanation thatis corroborated through a sequence-averaged stack-ing strength (see poly-A-avstk and poly-T-avstk). Fi-nally, for higher temperatures self-hybridisation be-comes less important and the radius of gyration ap-proaches the same plateau value for all sequences.
FIG. 8. Simulation snapshots of the λ -ssDNA sequencefor three different temperatures, at 0 ◦ C, 20 ◦ C and 60 ◦ C.The type of each nucleotide is represented by a colourscheme: A (white), T (Cyan), G (blue) and C (red).
A fraction of complementary nucleotides (A-T orG-C) on the single-stranded λ -DNA chain are close1enough to form hydrogen bonds. Due to the cooper-ativity of base pairing, long stems with many proxi-mal base pairs tend to form between regions of highcomplementarity - these are the characteristic hair-pins in Fig.8. This transition between a flexible ss-DNA and significantly more rigid hairpins of dsDNA(the persistence length of dsDNA is 50 nm, aroundthirty times larger than that of ssDNA) is mediatedby e.g. changes in the temperature, salt concentra-tion or pH value. In Fig.7 we show the radius of gy-ration R g and the contact fraction (the number ofcontacts N c normalised by half the number of nu-cleotides in the ssDNA strand, which is the maxi-mum number of possible base pairs) for the single-stranded λ -DNA versus temperature. A contact wasdefined when the hydrogen-bonding interaction sitesof any two nucleotides were less than 0.45 length unitsapart, regardless of the individual bases. In princi-ple, this criterium cannot prevent stacked, nearest-neighbour nucleotides from being counted as a con-tact. Nevertheless it proved sufficiently accurate fora perfect dsDNA duplex where the number of con-tacts N c = N/
2. Additionally, this definition willtend to include mismatched base pairs in a duplexas contacts. It will thus overestimate the number ofcorrectly-formed Watson-Crick base pairs, but for ourpurposes it is more important that 2 N c provides agood estimate of the number of bases incorporatedinto hairpin structures.At 0 ◦ C, around 44% of the nucleotides are involvedin contacts. When the temperature increases, the sys-tem destabilises and the number of contacts decreasessignificantly until it flattens out at 50 ◦ C (the sametemperature at which R g has a plateau). However,while the contact fraction changes dramatically (morethan a factor 40 from about 0.44 to 0.01) in this tem-perature range, there is only a small change in theradius of gyration (around a factor 1.45 from 11.3 nmto 16.4 nm). Related snapshots from simulations atselected temperatures are given in Fig. 8.We apply the same protocol as before to the pUC19plasmid, consisting of a ssDNA sequence of 2686 nuc-leotides. For simplicity we opted for a linear moleculewith freely rotating ends. The radius of gyration asa function of temperature is shown in Fig.9. The be-haviour is very similar to the one of λ -ssDNA, par-ticularly the minor effect that temperature changeshave on R g despite dramatic changes in the numberof contacts between nucleotides. While for λ -ssDNAthe radius of gyration at 20 ◦ C equals 4.5% of its totalcontour length, in the case of the plasmid R g repre-sents only 2.2%. Using the theoretical expression for R g in Eq. 8 below and monomer length a =0.65 nm,Kuhn segment length b =2 nm and Flory exponent ν = 0 . R g /aN = 5% ( λ -DNA) and2 .
5% (plasmid), respectively. Hence, the computa-tional values are about 10% smaller than the theoret-ical values, but generally consistent with the latter. Ataround 50 ◦ C R g reaches a plateau, which is at leastconstant within the error bars. It is interesting to seethat the λ -ssDNA exhibits the same tendency at the
30 35 40 45 50 55 0 10 20 30 40 50 60 0 0.1 0.2 0.3 0.4 0.5 R g [ n m ] C on t ac t F r ac ti on f h Temperature [ ° C] Contact Fraction pUC19 ssDNAR g pUC19 ssDNA FIG. 9. Dependence of the radius of gyration R g andthe fraction of direct intra-chain nucleotide contacts in thessDNA pUC19 plasmid sequence on the temperature.FIG. 10. Simulation snapshots of the pUC19-ssDNA se-quence for two different temperatures, at 20 ◦ C and 60 ◦ C.The type of each nucleotide is represented by a colourscheme: A (white), T (Cyan), G (blue) and C (red). same temperature. As reference we also modelled thedouble-stranded linear pUC19 plasmid, for which wemeasured values of R g in the region of 130 nm at 20 ◦ Cand 170 nm at 60 ◦ C, respectively, so about a factor 3to 4 larger than the values of R g we obtained for thessDNA sequence.In Fig.10 we can see that at 20 ◦ C several nucleotideshave hybridised, forming hairpin structures of 20-30 bp located along the plasmid. When we increasethe temperature of the system up to 60 ◦ C the hair-pins disappear as self-hybridisation is suppressed, ac-counting for the substantial reduction of intra-chain2contacts.The interpretation of these results is not en-tirely uncomplicated as several interlinked mecha-nisms are at work that all influence the radius ofgyration. When hairpins (or indeed any contact be-tween bases) form, the hydrogen-bonding between nu-cleotides short-circuits all bases that are part of thehairpin, effectively shortening the contour length ofthe biopolymer. Hence, self-hybridisation leads toa smaller radius of gyration through a reduction ofthe effective contour length. Thus the smaller radiusof gyration at lower temperatures can be partly ex-plained with basic polymer physics. On the otherhand, the contribution of hairpins to the total valueof R g is not zero, bearing in mind that even a rigidrod has a finite radius of gyration. The impact ofself-hybridisation is thus a priori not easily assessed.Moreover, regions cut out in this way are generallybulky, tending to swell the DNA strand relative to ashorter polymer with no base pairing. This consti-tutes an excluded volume effect which increases R g .The exact number of hairpins and the degree of self-hybridisation depend ultimately on sequence of thessDNA strand and are generally not quantifiable onthe sole basis of polymer physics.Nevertheless, some of the dependence of the radiusof gyration on the number of formed base pairs can berationalised using a simple and idealised physical poly-mer model. We assume that all nucleotide contactsare contained in well-defined hairpins. The singlestranded DNA can thus be modelled as a self-avoidingpolymer with attached rigid, rod-like hairpins that arecut out of the contour length of the polymer.At high temperature the base pairing can be ne-glected and the genome can be modelled as a self-avoiding walk (SAW) polymer with radius of gyration R g = b √ N ν Kuhn (7)where b is the Kuhn segment of the polymer. At saltconditions used in the oxDNA simulations, c Na =0 . b ≈ ν is the scaling exponent [39] and N Kuhn = N a/b the number of Kuhn segments in the polymer, with a = 0 .
65 nm [36, 37]. Scaling exponent of a SAWpolymer is ν = 0 .
588 which holds for poly-T ssDNAat physiological salt concentration [37]. Therefore, theradius of gyration is R g = b √ (cid:18) aNb (cid:19) ν (8)with N the number of nucleotides. For the λ -ssDNAsequence and the linear pUC19 plasmid this leadsto R g ( N = 500) = 16 . R g ( N = 2686) =43 . N c gives a good estimate of the totalnumber of bases cut out of the contour length by hy-bridisation, the effective contour length of ssDNA is reduced to N ss = N − N c . Consequently, the effec-tive radius of gyration of the ssDNA is reduced to R g,ss = b √ (cid:18) a ( N − N c ) b (cid:19) ν (9)depending on the number of contacts N c . Hairpins,however, also contribute to R g . Assuming that a hair-pin is a rigid rod with length l (justifiable for hairpinsshorter than about 100 nm) the radius of gyration ofevery hairpin is R g,h = l/ √
12. If k hairpins of equallength are formed, each hairpin will contribute R g,h = a ds N c / ( k √
12) (10)with the effective monomer length reduced due to he-licity of double stranded DNA a ds = 0.34 nm. Thisconditions applies as all hairpins combined need toadd up to the length along the contour that is in con-tact.The total radius of gyration of an object is a sumover its subparts, where each subpart contributes itsown radius of gyration plus a centre-of-mass distancesquared, weighted by the mass. The centre-of-mass ofthe total ssDNA and hairpin system is therefore c m = f h k k (cid:88) i =1 x i + l n i (11)with x i the (vector) position of the i -th hairpin base,i.e. the end where the hairpin is attached to the poly-mer. The centre-of-mass position of the i -th hairpin is x i + l ˆ n i with ˆ n i the unit vector specifying the orienta-tion of the hairpin’s major axis. Note that only hair-pins contribute because we chose the centre-of-massof the ssDNA polymer as the origin of our coordinatesystem. The weight factor f h = 2 N c /N is determinedby the fraction of total polymer mass contained inthe hairpins. The quantity f h is equal to the contactfraction shown in Figs. 7 and 9. The total radius ofgyration of the ssDNA and hairpins system becomes R g = (1 − f h )( R g,ss + c m )+ f h k k (cid:88) i =1 R g,h +( x i + l n i − c m ) (12)where the first term on the right hand side is the con-tribution of the ssDNA and the second term, the sum,is performed over all k hairpins. Note that the fractionof total mass in each hairpin is f h /k and x i + l ˆ n i − c m is the distance between the hairpin centre-of-mass andsingle-stranded polymer centre-of-mass.Assuming that the positions of hairpins are uni-formly random and uncorrelated, inserting Eq. (11)into Eq. (12) and employing some basic algebra out-lined in Appendix B, the expected value for thesquared radius of gyration is obtained3 (cid:104) R g (cid:105) = R g,ss (cid:18) − f h k (cid:19) + R g,h (cid:18) f h − f h k (cid:19) (13)with R g,ss and R g,h given by Eqs. (9) and (10), re-spectively, and the contact fraction f h = 2 N c /N .We have assumed that hairpins do not interactwith the ssDNA polymer, or with other hairpins, andthat all k hairpins are of the same length. How-ever, even this relatively simple, idealised derivationdemonstrates that the radius of gyration depends onboth the number of contacts and the hairpin length.This is shown in Fig. 11 for a sequence of N = 500nucleotides, i.e. the length of our λ -ssDNA. The de-pendence on the number of contacts is obviously non-monotonous. The values k = 1 and k = N c / (cid:104) R g (cid:105) pro-vide the upper and lower physical limit for the ex-pected value of the radius of gyration. The simula-tions, Fig. 7, result in a radius of gyration around12.6 nm and 11.3 nm at the observed contact frac-tion of around 32% and 45%, respectively, in goodagreement with the theoretical prediction shown onFig. 11. We also see that for an even larger number ofcontacts the possible range of values of (cid:104) R g (cid:105) is quitewide. This is of course a much idealised and simpli-fied reasoning, but it elucidates the non-trivial natureof these interdependencies.The theory neglects the ex-cluded volume of regions cut out of the contour lengthby hybridisation; taking this into would increase the R g in Eq. (12), while additional bases cut out by hy-bridisation but not contributing to N c would decreaseit. We speculate that the two effects cancel out, to adegree, resulting in a good agreement between theoryand simulations. � = �� = �� = � � / � ������� �������� � � � � � � � � � � �� � � �� �� ( < � � � > ) � / � [ � � ] FIG. 11. Radius of gyration (cid:112) (cid:104) R g (cid:105) as a function ofthe contact fraction f h = 2 N c /N for different number offormed hairpins k . The curves were obtained from Eq. 13using the following parameters: Kuhn segment b =2 nm,nucleotide size a =0.65 nm, double stranded nucleotide ef-fective size a ds =0.34 nm, scaling exponent ν = 0 . N = 500. VII. CONCLUSIONS
The implementation of the oxDNA model forcoarse-grained DNA modelling into a communitymolecular dynamics code such as LAMMPS reducesthe entry barrier of using the model significantly.Moreover, it allows to combine this coarse-grainedforce field with different features that are already en-abled in LAMMPS.The Langevin-type rigid-body integrators that aredistributed together with the LAMMPS USER-package, particularly the DOT-C integrator, offer ad-ditional advantages over the existing standard rigid-body integrators for Langevin dynamics. They showimproved stability at the costs of a very small over-head. This permits larger timesteps and thereforelarger physical simulation times.The parallel performance of the MPI-only imple-mentation, as demonstrated through scaling tests us-ing a simple benchmark, is excellent provided thereare at least a few dozen particles per MPI-task. Theseresults show effectively that the oxDNA model is wellsuited for large and extremely large problems in DNAand RNA modelling. It can tackle problem sizes thatwere well beyond the reach of the original standaloneimplementation of the model. It is worth mentioningthat the GPU-accelerated version of the standalonecode is also limited to speedups of typically a factor30 compared to the single core performance. Basedon the scaling analysis of the benchmarks it could besaid that this is matched by the performance of a sin-gle multi- or many core chip.The applications we opted for, a sequence of lin-ear, single-stranded λ -bacteriophage and pUC19 plas-mid DNA, are motivated primarily by currently ongo-ing projects in the under-investigated area of single-stranded DNA, rather than by an attempt to harvestthe performance of the newLAMMPS implementation. The results shows thatthe conformation of ssDNA is strongly affected bythe tendency to self-hybridise upon cooling, i.e. toform intra-chain base pairs between complementarynucleotides on the same strand that lead to hairpins,local regions of dsDNA, and less structured domainsof clustered nucleotides. The radius of gyration R g ofboth ssDNA examples is predicted to be relatively in-sensitive towards temperature changes between 0 ◦ Cand 60 ◦ C. The slight reduction of R g can be atleast partly explained with a shorter effective contourlength of the biopolymer due to hairpin formation.This explanation, however, disregards some of themore subtle intricacies of the self-hybridisation pro-cess. Hairpins contribute as well to the total value of R g . The hybridised domains of clustered nucleotidesintroduce an excluded volume effect, which increasesthe radius of gyration. Last but not least, the DNA se-quence determines whether any self-hybridisation canoccur in the first place. It should be noted that thereis a large number of possible self-hybridised bondingconfigurations. This means that the system is likelyto fall into a particular one upon quenching and to4remain there. However, by using a number of inde-pendent configurations we have presumably reachedstates that are representative, although these are notguaranteed to be the most stable ones.In the future it may be possible to focus on ring mo-lecules that contain superhelical twist and have differ-ent number of helical turns compared to their naturalform. These rings may be opened by introducing asingle strand break, which releases the superhelicaltwist, a mechanism that is known to be highly rele-vant during gene replication and expression. ACKNOWLEDGMENTS
This work was funded under the embedded CSEprogramme of the ARCHER UK National Super- computing Service (eCSE05-10). YAGF acknowl-edges support from the Mexican National Coun-cil for Science and Technology (CONACyT, PhDGrant 384582). TC acknowledges support from theHerchel Smith Scholarship and the CAS PIFI Fel-lowship. TEO acknowledges his Royal Society Uni-versity Research Fellowship. OH acknowledges sup-port from the EPSRC Early Career Fellowship Scheme(EP/N019180/2).
AUTHORS CONTRIBUTIONS
OH, YAGF, TC and TEO designed and performedresearch, analysed data, and wrote the paper. Theimplementation was undertaken in a collaboration be-tween OH and TEO. [1] C. Calladine, H. R. Drew, B. F. Luisi, and A. A.Travers,
Understanding DNA (Elsevier AcademicPress, London, UK, 2004).[2] C. A. Laughton and S. A. Harris, WIREs Comput.Mol. Sci. , 590 (2011).[3] D. A. Potoyan, A. Savelyev, and G. A. Papoian,WIREs Comput. Mol. Sci. , 69 (2013).[4] P. D. Dans, J. Walther, H. G´omez, and M. Orozco,Curr. Opin. Struct. Biol. , 29 (2016).[5] C. Maffeo, T. T. M. Ngo, T. Ha, and A. Aksimentiev,J. Chem. Theory Comput. , 2891 (2014).[6] N. Korolev, D. Luo, A. Lyubartsev, and L. Norden-ski¨old, Polymers , 1655 (2014).[7] M. Maciejczyk, A. Spasic, A. Liwo, and H. Scheraga,J. Chem. Theory Comput. , 5020 (2014).[8] S. Pronk, S. Pall, R. Schulz, P. Larsson, P. Bjelk-mar, R. Apostolov, M. R. Shirts, J. C. Smith, P. M.Kasson, D. van der Spoel, B. Hess, and E. Lindahl,Bioinformatics , 845 (2013).[9] M. R. Machado and S. Pantano, J. Chem. TheoryComput. , 1568 (2016).[10] J. J. Uusitalo, H. I. Ingolfsson, P. Akhshi, D. P. Tiele-man, and S. J. Marrink, J. Chem. Theory Comput. , 3932 (2015).[11] J. C. Phillips, R. Braun, W. Wang, J. Gumbart,E. Tajkhorshid, E. Villa, C. Chipot, R. D. Skeel,L. Kale, and K. Schulten, J. Comput. Chem. ,1781 (2005).[12] C. Markegard, I. Fu, K. Reddy, and H. Nguyen, J.Chem. Phys. B , 1823 (2015).[13] http://lammps.sandia.gov .[14] D. M. Hinckley, G. S. Freeman, J. K. Whitmer, andJ. J. de Pablo, J. Chem. Phys. , 144903 (2013).[15] D. M. Hinckley and J. J. de Pablo, J. Chem. TheoryComput. , 5436 (2015).[16] C. A. Brackley, A. N. Morozov, and D. Marenduzzo,J. Chem. Phys. , 135103 (2014).[17] Y. A. G. Fosado, D. Michieletto, and e. a. J. Allan,Soft Matter , 9458 (2016).[18] T. E. Ouldridge, A. A. Louis, and J. P. K. Doye, J.Chem. Phys. , 085101 (2011).[19] T. Ouldridge, Coarse-grained modelling of DNA andDNA self-assembly , Ph.D. thesis, University of Ox-ford, 2011 (2011). [20] https://dna.physics.ox.ac.uk .[21] T. Ouldridge, R. Hoare, A. Louis, J. Doye, J. Bath,and A. Turberfield, ACS Nano , 2479 (2013).[22] J. Holbrook, M. Capp, R. Saecker, and M. Record,Biochemistry , 8409 (1999).[23] J. SantaLucia Jr and D. Hicks, Annu. Rev. Biophys.Biomol. Struct. , 415 (2004).[24] B. Snodin, F. Randisi, M. Mosayebi, P. Sulc, J. Ro-mano, F. Romano, T. Ouldridge, R. Tsukanov, E. Nir,A. Louis, and J. Doye, J. Chem. Phys. , 234901(2015).[25] P. Sulc, F. Romano, T. Ouldridge, L. Rovigatti,J. Doye, and A. Louis, J. Chem. Phys. , 135101(2012).[26] https://ccpforge.cse.rl.ac.uk/gf/project/cgdna .[27] M. P. Allen and D. J. Tildesley, Computer Simula-tion of Liquids (Oxford University Press, Oxford, UK,1989).[28] M. Allen and G. Germano, Mol. Phys. , 3225(2006).[29] .[30] W. Humphrey, A. Dalke, and K. Schulten, J. Mol.Graphics , 33 (1996).[31] R. L. Davidchack, T. E. Ouldridge, and M. V.Tretyakov, J. Chem. Phys. , 144114 (2015).[32] T. F. Miller, M. Eleftheriou, P. Pattnaik, A. Ndi-rango, and D. Newns, J. Chem. Phys. , 8649(2002).[33] L. Rovigatti, F. Smallenburg, F. Romano, andF. Sciortino, ACS Nano , 3567 (2014).[34] C. D. Michele, L. Rovigatti, T. Bellinic, andF. Sciortino, Soft Matter , 8388 (2012).[35] “LAMMPS Documentation,” http://lammps.sandia.gov/doc/Manual.html (2016), chapter 8:Performance and Scalability.[36] N. M. Toan and C. Micheletti, Journal of Physics:Condensed Matter , S269 (2006).[37] A. Y. L. Sim, J. Lipfert, D. Herschlag, and S. Do-niach, Phys. Rev. E , 021901 (2012).[38] H. Chen, S. P. Meisburger, S. A. Pabit, J. L. Sutton,W. W. Webb, and L. Pollack, Proceedings of theNational Academy of Sciences , 799 (2012). [39] P.-G. de Gennes, Scaling concepts in polymer physics (Cornell University Press, 1979).
APPENDIXAppendix A: Profiling
Profiling allows a detailed analysis of the imple-mentation and gives an overview of how much timethe code spends in each individual subroutine. Weused the Craypat Performance Tools on the ARCHERUK National Supercomputing Service to conduct sam-pling experiments of the high and low density bench-marks. Although the experiments where actually per-formed with the oxDNA model, the results are rep-resentative as well for oxDNA2 as the only differencebetween the two is a different local geometry of theinteraction sites and an additional pair interaction inform of a Debye-H¨uckel potential.Fig. 12 shows a pie chart of the low density (LD)run. The image on the left shows the results on a sin-gle node with 24 MPI-tasks, whereas the image on theright is for 2048 MPI-tasks. Focussing first on a singlenode, calls to the MPI-library are below 5% and do notappear with an individual pie section. The total timespent in the force calculation is around 86% (accordingto the LAMMPS breakdown). Interestingly, a signifi-cant fraction of the time is spent on calculating the lo-cal body coordinate system of the nucleotide from thequaternion degrees of freedom (MathExtra::q to exy,11 . . Appendix B: Derivation of (cid:104) R g (cid:105) We assume that the position of hairpin bases, aswell as the orientation of hairpins, is uniformly ran-dom and uncorrelated along the ssDNA contour, for-mally: (cid:104) x i (cid:105) = 0, (cid:104) x i (cid:105) = R g,ss , (cid:104) x i x j (cid:105) = 0 for i (cid:54) = j ,and similarly for the orientation: (cid:104) ˆ n i (cid:105) = 0, (cid:104) ˆ n i (cid:105) = 1, (cid:104) ˆ n i ˆ n j (cid:105) = 0 for i (cid:54) = j , (cid:104) x i ˆ n j (cid:105) = 0. These propertiesresult in (cid:104) c m (cid:105) = 0and (cid:104) c m (cid:105) = f h k ( R g,ss + l / . The average R g becomes (cid:104) R g (cid:105) = (1 − f h ) R g,ss + (1 − f h ) (cid:104) c m (cid:105) + f h R g,h ++ f h k (cid:68) (cid:88) i ( x i + l n i − c m ) (cid:69) . (B1)The average of the sum is (cid:68) k (cid:88) i =1 ( x i + l n i − c m ) (cid:69) = k (cid:88) i =1 (cid:110) (cid:104) x i (cid:105) + (cid:104) c m (cid:105) + l (cid:104) ˆ n i (cid:105) +6 FIG. 12. Craypat performance analysis of a sampling experiment for the low density benchmark (60 kbp) ona single node (left, 24 MPI-tasks) and for 2048 MPI-tasks (right). Note that the assigned colour code for thefunctions is different in both cases.FIG. 13. Craypat performance analysis of a sampling experiment for the high density benchmark (960 kbp) ona single node (left, 24 MPI-tasks) and for 2048 MPI-tasks (right). Note that the assigned colour code for thefunctions is different in both cases. l (cid:104) x i ˆ n i (cid:105) − (cid:104) x i c m (cid:105) − l (cid:104) ˆ n i c m (cid:105) (cid:111) = kR g,ss + f h ( R g,ss + l /
4) + k l − f h R g,ss − f h l (cid:80) i (cid:104) x i c m (cid:105) = f h k (cid:80) ij (cid:104) x i ( x j + l ˆ n j ) (cid:105) = f h R g,ss and (cid:80) i (cid:104) ˆ n i c m (cid:105) = f h l . Furthermore, l = 12 R g,h .Using these relations the expected value for thesquared radius of gyration is obtained (cid:104) R g (cid:105) = R g,ss (cid:18) − f h k (cid:19) + R g,h (cid:18) f h − f h k (cid:19) , (B3)with R g,ss and R g,h given by Eqs. (9) and (10), re-spectively, and the contact fraction f h = 2 N c /N/N