Migration von DNA auf strukturierten

Oberflachen in einem außeren Feld

Diplomarbeit

zur Erlangung des Grades eines Diplom-Physikers

vorgelegt von

Martin Andreas StreekTheorie der kondensierten Materie

Fakultat fur PhysikUniversitat Bielefeld

begutachtet durch

Prof. Dr. Friederike SchmidProf. Dr. Dario Anselmetti

vorgelegt am

5. April 2002

2

Contents

Contents 1

1 Introduction 31.1 Properties of the DNA molecule . . . . . . . . . . . . . . . . . . . 31.2 Separation techniques . . . . . . . . . . . . . . . . . . . . . . . . 5

1.2.1 Gel electrophoresis . . . . . . . . . . . . . . . . . . . . . . 51.2.2 Entropic traps . . . . . . . . . . . . . . . . . . . . . . . . 61.2.3 Ratchets . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71.2.4 Other techniques . . . . . . . . . . . . . . . . . . . . . . . 8

2 Simulation Technique 92.1 The Verlet algorithm . . . . . . . . . . . . . . . . . . . . . . . . . 102.2 Langevin dynamics . . . . . . . . . . . . . . . . . . . . . . . . . . 10

2.2.1 Selection of the time step . . . . . . . . . . . . . . . . . . 12

3 The model 143.1 The Gaussian chain model . . . . . . . . . . . . . . . . . . . . . . 14

3.1.1 End-to-End vector . . . . . . . . . . . . . . . . . . . . . . 153.1.2 Radius of gyration . . . . . . . . . . . . . . . . . . . . . . 16

3.2 Modifications of the model . . . . . . . . . . . . . . . . . . . . . . 163.2.1 The repulsive Lennard-Jones potential . . . . . . . . . . . 163.2.2 Bond length potential . . . . . . . . . . . . . . . . . . . . 173.2.3 Bond-angle potential . . . . . . . . . . . . . . . . . . . . . 193.2.4 Electric field . . . . . . . . . . . . . . . . . . . . . . . . . 193.2.5 Further aspects . . . . . . . . . . . . . . . . . . . . . . . . 20

4 Simulation program 214.1 Programming details . . . . . . . . . . . . . . . . . . . . . . . . . 21

4.1.1 Boxing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 214.1.2 Prefactors . . . . . . . . . . . . . . . . . . . . . . . . . . . 244.1.3 Random number generator . . . . . . . . . . . . . . . . . 25

4.2 Testing the program . . . . . . . . . . . . . . . . . . . . . . . . . 254.2.1 End-to-end vector . . . . . . . . . . . . . . . . . . . . . . 254.2.2 Rouse-mode analysis of a Gaussian chain . . . . . . . . . 26

1

5 The free Chain 305.1 Parameter set for the chain in simulation . . . . . . . . . . . . . 305.2 Radius of gyration . . . . . . . . . . . . . . . . . . . . . . . . . . 315.3 Mobility in free solution . . . . . . . . . . . . . . . . . . . . . . . 325.4 Distribution of the bond length . . . . . . . . . . . . . . . . . . . 345.5 Persistence length . . . . . . . . . . . . . . . . . . . . . . . . . . 36

6 The entropic trap 396.1 Simulation of the wall . . . . . . . . . . . . . . . . . . . . . . . . 396.2 Geometric setup of the trap . . . . . . . . . . . . . . . . . . . . . 406.3 First runs in the trap . . . . . . . . . . . . . . . . . . . . . . . . . 426.4 Tilted force . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 436.5 The trap with tilted walls . . . . . . . . . . . . . . . . . . . . . . 45

7 Results 477.1 Tilt dependence of the mobility . . . . . . . . . . . . . . . . . . . 477.2 Force dependence of the mobility . . . . . . . . . . . . . . . . . . 487.3 Separation by chain length in the trap . . . . . . . . . . . . . . . 49

8 Conclusions 55

A Rouse-mode analysis 56A.1 Rouse-modes of the bead-spring model . . . . . . . . . . . . . . . 56A.2 Time evolution of the different modes . . . . . . . . . . . . . . . 58

A.2.1 Movement of the center of mass . . . . . . . . . . . . . . . 59A.2.2 Correlation of the higher modes . . . . . . . . . . . . . . . 59

B All mobilities 61

C Source code 69

List of Figures 87

Bibliography 89

Acknowledgements 91

2

Chapter 1

Introduction

DNA (DeoxyriboNucleic Acid) is probably the most famous polyelectrolyte innature because it contains all information necessary to reproduce almost everyknown living organism. The information stored in the DNA is encrypted in asequence of four bases. In a cell, DNA is stored in a pair of strands that areentwined in a double helix and the bases are arranged such that two comple-mentary bases always face each other:

• (A)denine ↔ (T)hymine

• (G)uanine ↔ (C)ytosine

This provides a rather easy method of replication: when a cell divides, thetwo DNA strands separate, and the complementary partner strand of each isrebuilt from surrounding base acids. This is shown in fig. 1.1.

Figure 1.1: The recombination of double stranded DNA [19]

1.1 Properties of the DNA molecule

DNA is easy to model for a number of reasons:

3

• DNA is a polyelectrolyte with two negative charges per base pair. Thus,an electric force acting on the polymer acts uniformly on each monomer.

• Usually, one would expect that DNA segments to repel each other due toelectrostatic interactions. However, in most cases, DNA is examined in asalty solvent. The solvent provides enough ions to screen the electrostaticinteractions for distances greater than the Debye screening length. In aphysiological environment, the Debye screening length is only a few nm,which is of the same order of magnitude as the diameter of the DNA doublehelix, which in turn is found to be 2nm1 and much smaller than the per-sistence length. Thus, when simulating long DNA molecules in a solvent,the electrostatic interactions can be neglected in good approximation.

• The folding of the DNA molecule is independent of the sequence of basesstored in it. This is much easier to model than, for example, proteinfolding, where changing one single amino acid can change the folding ofthe protein.

• The radius of gyration is typically up to a few µm, depending on thenumber of monomers.

• DNA is a rather stiff polymer. The persistence length of DNA depends onthe ion concentration in the buffer and is usually around 50 nm [17, 18].

• The number of base pairs extends to 3 · 109 for human DNA, resulting ina contour length of up to 1 m, but it is arranged in 24 chromosomes.

As a result, it becomes clear that the simulation of DNA in a solvent can beeasily executed, in good approximation, using stiff bead spring chains. This isdescribed in chapter 3 and [19].

Thus, it is obvious that DNA is interesting for biologists and physicists. Forbiologists, the capability of handling and understanding DNA is key to success-fully doing genetic engineering. For physicists, DNA is the classical example ofa polyelectrolyte, because the folding is independent of the sequence of bases.Furthermore, a dragging force can be be applied in a simple manner, but elec-trostatic interactions are screened.

In this thesis I will examine a new method of separating DNA strands bylength. The most famous application of separating DNA segments by length wasthe sequencing of the human genome [15]. Although the sequencing was finishedyears ahead of the original estimates, it is still necessary to find procedureswhich are easy to use and are able to separate chains effectively and quicklyby length. Once this task has been accomplished, it will become possible tosequence additional genomes and learn to understand the genetic code. As forapplications, a thorough understanding of the genetic code is set to spark hugeadvances in agricultural science and other fields.

1Modeling the DNA helix, to be found athttp://wug.physics.uiuc.edu/courses/phys450/Lectures2001/L15/ModelingDNAhelix.pdf(03-26-2002)

4

1.2 Separation techniques

DNA is a polyelectrolyte. This suggests that a separation by chain length cansimply be performed by applying an electric field to the DNA in solvent, e.g.water. Unfortunately, the diffusion coefficient of the DNA in water behaves likeO(N−1), where N denotes the number of monomers. When applying a forceonto each monomer, this leads to a force that behaves like O(N). Thus, themobility of DNA in solvent is constant for all chain lengths, and no separationis possible.

1.2.1 Gel electrophoresis

The common method currently used for DNA separation is electrophoresis in apolymer gel. The gel introduces a sieving matrix for the DNA. A sketch can befound in fig. 1.2.

Figure 1.2: Electrophoresis in a gel [19]

Here three regimes of the constant-field electrophoresis are shown:

(a) Ogston Sieving: in this case, the polymer is rather short. Thus, the polymerhas to pass through the barriers as a whole.

(b) Reptation without orientation: in this regime, the polymer length is onthe characteristic length of the sieving medium (∼pore size). When thepolymer passes through the obstacles, it has to unfold, and one end hasto find the way through the sieving matrix. The term reptation reflectsthe movement of the polymer, which is very similar to the movement ofa snake or reptile. That’s where the term is derived from. On the otherhand, the polymer is too short to identify an orientation.

(c) Reptation with orientation: in this regime, the polymer is sufficiently longso an orientation can be identified.

5

In [19], the mobility depending on the polymer length is computed in thebiased reptation model. The result is:(

µ

µ0

)BRM

≈ 13

(1N

+ε2

3

)(1.1)

where µ0 is the mobility in free solution, N is the number of monomers and εdenotes the reduced electric field. For large values of N , the mobility obviouslydoes not change as a function of the number of monomers. Thus only shortpolymers can be separated in gel electrophoresis.

Another problem arises when examining the pore size of the gel, which isnot uniform. Thus, all experimental results are softened, which complicates theseparation by chain length. Typical times for electrophoresis in a gel are around12− 24h [8].

1.2.2 Entropic traps

A rather new method used in electrophoresis is entropic trapping. These artifi-cial sieving matrices consist of alternating wide and thin volumes. A schematicplot of an entropic trap is shown in fig. 1.3. The applied electric field is per-pendicular to the transition lines between the thin and wide regions. The ideabehind entropic traps is that the polymer has to form a thread when passingthe thin region and can form a coil again upon reaching the wide region.

F

Figure 1.3: The entropic trap

The mobility of long polymers in entropic traps is expected to be lowerthan the mobility of short polymers, because long polymers lose more entropicfree energy than short polymers when unfolding. Surprisingly, it is observedexperimentally that long polymers travel at higher speeds than short polymers.This is described and explained in [8]. The explanation states that the importantprocess during trapping is when the first loop enters the thin area. Obviously,longer chains have a greater chance of forming the initial loop. Typical timesfor successful separation in entropic traps are around 15 minutes, which is veryshort [8] compared to electrophoresis in a gel.

In this thesis, I will focus on entropic traps.

6

1.2.3 Ratchets

Another promising idea is presented by ratchets. Ratchets use Brownian motionto accomplish a movement with applied no net force. A schematic drawing of aratchet is given in fig. 1.4.

Figure 1.4: Schematic drawing of a ratchet

In a ratchet, there is an asymmetric potential in space and time which is aperiodical sequence of short, steep gradients on the one hand and long, ratherflat gradients on the other hand. The potential is periodically turned on andoff. A particle inside the ratchet will move as follows (see fig. 1.5):

• First, the potential is turned on, the particle moves to a minimum of thepotential.

• Then the potential is turned off. The particle diffuses freely and theprobability distribution is Gaussian.

• When the potential is switched on again, the particle moves accordingto the gradient of the potential. The probability to travel the distancenecessary to pass the short region is greater than diffusing through thelong region. Using a proper setup, no particles will have the time todiffuse to the maximum on the right, whereas it is possible to diffuse tothe maximum on the left. Thus some particles have moved to the left, butnone to the right.

Figure 1.5: A ratchet device [13]

7

Apparently, a net transport to the left has occurred. The diffusion constantof a polymer behaves like N−1. It is expected that a separation by lengthis possible using a ratchet device because long chains diffuse slower. Furtherinformation concerning ratchets can be found in [1, 13].

1.2.4 Other techniques

The separation methods presented above are exemplarily and do not representa complete list. Other approaches (like capillary electrophoresis or pulsed fieldgel electrophoresis) have been made, as well. The interested reader is referredto [16, 19].

8

Chapter 2

Simulation Technique

A common method for computer simulations, the Monte Carlo (MC) algorithmproposes moves, i.e. changes, for the current configuration which are eitheraccepted or rejected according to a “Metropolis” criterion that depends on thechange in energy and on the temperature. Although this method is fast, it isdesigned for the study of equilibria, not that of dynamics, for it lacks a physicaltime parameter. It can thus not be used to study mobility.

Another common method is the molecular-dynamics (MD) method. Here,the equations of motion are integrated numerically:

mi~ri = −∂V∂~ri

(2.1)

The crucial point with MD simulations is to keep the computational effortas low as possible. This means in particular that it shall suffice to compute theforces only once per time step because their computation may involve compli-cated functions. Moreover, it is important to avoid too many errors during thesimulation as they would require smaller time steps and, in the long run, caneven destroy the results:

• Short-time stabilityIt is essential for the algorithm to be stable and very accurate in any singletime step.

• Long-time stabilityJust as important as short-time stability is the requirement that the con-servation laws and symmetries be fulfilled. In the micro-canonical ensem-ble, it is especially important that the energy of the system be conserved:∣∣∣E(t)−E(0)

E(0)

∣∣∣� 1.

9

2.1 The Verlet algorithm

The most commonly used MD algorithm is the Verlet algorithm. It is based onthe following time step:

~ri(t+ ∆t) = 2~ri(t)− ~ri(t−∆t) +h2

mi

~fi(t) +O(∆4t ) (2.2)

~vi(t) =1

2∆t[~ri(t+ ∆t)− ~ri(t−∆t)] +O(∆2

t ) (2.3)

which can be derived from the Taylor series of the real trajectory. ~fi is theforce acting on each particle. In this form, the algorithm has two major dis-advantages: there is a numerical rounding error in eqn. 2.2 when computing2~ri(t)−~ri(t− δt), and velocities can only be computed one time step later (eqn.2.3). These problems can be avoided by using the Velocity-Verlet algorithm,which is algebraically equivalent to the Verlet algorithm:

~ri(t+ ∆t) = ~ri(t) + h~vi(t) +∆2t

2mi

~fi(t) +O(∆3t ) (2.4)

~vi(t+ ∆t) = ~vi(t) +∆t

2mi

(~fi(t) + ~fi(t+ ∆t)

)+O(∆3

t ) (2.5)

This algorithm has a lot of advantages:

• Local error of order ∆3t

• Global error of order ∆2t

• The two equations are symmetric in ∆t, hence the algorithm is time-reversible

• Exact phase-space conservation (”symplectic”)

Benefiting from time-reversibility and symplecticity, extremely good long-time stability is achieved. This also guarantees very good conservation of energy,so this algorithm is very convenient. Today, it is standard in MD simulations.For a detailed discussion of the Verlet-algorithm, see [5].

2.2 Langevin dynamics

When simulating a DNA molecule in solvent with MD, it is inevitable to takethe surrounding solvent molecules into consideration. In order to get a realisticbehavior of the DNA in dilute solution, much more solvent particles than DNAmonomers have to be simulated. When doing so, most computing-time is spenton the solvent particles. However, the motion of the solvent molecules is ofminor interest.

10

The idea behind Langevin dynamics is to replace the solvent particles withBrownian motion acting on each DNA monomer. To do so, a random force~ηi is introduced. To prevent the system from heating up, a friction ζ is alsointroduced. Hence the equations of motion become:

~ri = ~vi (2.6)

~vi = ~fi − ζ~vi + ~ηi. (2.7)

The random force has to fulfill [2]:

〈~ηi〉 = 0 (2.8)〈ηiα(t)ηjβ(t′)〉 = δijδαβδ(t− t′)2kBTζ (2.9)

Substituting the Langevin equations into the Verlet algorithm, the followingtime step is obtained:

~v(t+ ∆t

2 ) = ~v(t) + ∆t

2m~f(~r(t))− ζ

m∆t~v(t) +√

2kbT ζm∆tη

~r(t+ ∆t) = ~r(t) + ∆t~v(t+ ∆t

2 )

~v(t+ ∆t) = ~v(t+ ∆t

2 ) + ∆t

2~f(~r(t+ ∆t))

(2.10)

where η is Gaussian distributed and fulfills

〈~η〉 = 0

〈η2i 〉 = 1

(2.11)

Generating a Gaussian distributed random number requires a lot of com-putating time. In [4], it was shown that Gaussian random numbers are notnecessary for a successful simulation. Uniformly distributed random numbersare just as good, as long as they fulfill 〈~η2〉 = 3, which is equivalent to eqn.(2.11). When using uniformly distributed random numbers in three dimensions,it is necessary to use uniformly distributed random numbers on the unit sphere.I found that the fastest algorithm was to simply pick random numbers out ofthe unit cube, then check if the numbers picked were within the unit sphere. Ifthe numbers were outside the unit sphere, they were discarded and a new setwas generated which of course had to be checked again. Nevertheless, this kindof random number generation was roughly a factor of 2.8 faster than one whichuses a Gaussian distribution.

11

2.2.1 Selection of the time step

In computing eqn. 2.10, it is mandatory to take a small time step ∆t in order tokeep local errors low. It turns out that a crucial parameter for simulation is ζ∆t.This parameter influences the long-time averages, as well. For simplicity, themass is set to unity in the following calculations. Computing a single particle,one gets:

Diffusion

At each time step, an additional velocity is picked up by the particle:

~v0 =√

2kBTζ∆tη.

After n time steps the velocity has decreased to

~vn = (1− ζ∆t)n~v0,

due to friction. Using equations 2.11 and 2.8, it is possible to compute thediffusion in a single time step as follows:

⟨(~x(∆t)− ~x(0))2

⟩=⟨(~v(t)∆t)2

⟩=⟨(∑n ~vn∆t)2

⟩= (∆t)2

⟨(∑n ~v0(1− ζ∆t)n)2

⟩= (∆t)2〈~v2

0〉 (∑n(1− ζ∆t)n)2

= (∆t)22kBTζ∆t〈η2〉 1(ζ∆t)2

= 6kBTζ ∆t

(2.12)

Comparing this result to eqn. 3.6, it is easy to see that the equation ofdiffusion is fulfilled exactly for a chain of one monomer and at every time step.

Expectation value of the energy

Statistical mechanics tells us that

〈E〉 =32〈v2〉,

and hence

〈v2〉 = 3kBT.

which should naturally be fulfilled by a meaningful algorithm. However, wesee that these equations are not fulfilled for finite time steps. To remedy this

12

situation, we introduce a corrected temperature, T ′, at which the simulation isperformed:

〈v2(t+ ∆t)〉 = 〈v2〉(1− ζ∆t)2 + 2kBT ′ζ∆t〈η2〉

= 3kBT (1− 2ζ∆t + (ζ∆t)2) + 6kBT ′ζ∆t

= 3kBT + 6kB(T ′ − T )ζ∆t + 3kBT (ζ∆t)2

(2.13)

Demanding that 〈v2(t+ ∆t)〉 = 3kBT leads to a corrected temperature definedby

T ′ := T

(1− ζ∆t

2

). (2.14)

Now we face the problem that when using an uncorrected temperature, thediffusion is correct and when correcting the temperature the equation of dif-fusion is not fulfilled. The only way out is to reduce the critical parameterζ∆t. For example: when using an uncorrected temperature and simulating atζ∆t = 0.01, the energy will be 0.5% higher than the correct expectation value.Furthermore, it becomes clear that this algorithm is only useful for low frictionbecause otherwise the necessary time steps would become very small.

13

Chapter 3

Modeling a DNA Molecule

As described in [2], the properties of a sufficiently long chain do not dependon the interaction between two neighboring monomers. Thus, the interactionbetween two neighboring monomers becomes negligible for long chains, andthe properties depend merely on the number of monomers. Consequently, anysimple model reproduces the properties of a long chain.

3.1 The Gaussian chain model

The simplest model for a DNA molecule is the so-called bead-spring model.In this model, the molecule is represented as a chain of beads with mass m.For simplicity, m is set to unity. The monomers are connected by springs withpotential VSP = k

2 |~ri − ~ri+1|2 (see fig. 3.1).

Figure 3.1: The bead-spring model

The beads do not interact with each other, and the springs always try tocontract the bonds to zero length. Thermal fluctuations will keep the elongationof the bonds on average non-zero, however. This leads to the following Langevinequation for N monomers:

~rn = −k (2~rn − ~rn+1 − ~rn−1)− ζ~rn(t) + ~ηn , n = 1 .. N (3.1)

with the artificial end beads:

~r0 ≡ ~r1 (3.2)~rN+1 ≡ ~rN (3.3)

14

where ~ηi is the random force, which is Gaussian distributed and obeys thefollowing equations:

〈~ηi(t)〉 = 0 (3.4)〈~ηi,α(t)~ηj,β(t′)〉 = 2ζkBT δij δαβ δ(t− t′) (3.5)

Equation 3.4 simply states that the random force is not biased in any direc-tion. Equation 3.5 demands that different components of the random force aswell as the random force at different moments be uncorrelated. Additionally,this equation ensures that the diffusion of the center of mass ~RG of a chain withN monomers behaves like⟨

(~RG(t)− ~RG(0))2⟩

= 6kBT

Nζt =: Dt. (3.6)

D is called the diffusion constant.This model also benefits us for being solvable analytically, which is described

in app. A. Hence it provides a good test case in developing the simulationprogram.

3.1.1 End-to-End vector

In this model, it is straight forward to compute 〈(~R)2〉 ≡ 〈(~r1 − ~rN )2〉. Themonomers do not interact with each other, so each bond only depends on thetwo monomers at its ends. Equation A.26 describes the mean excitation of eachbond. One finds:

⟨~R2⟩

=⟨(~r1 − ~rN )2

⟩=∑N−1n=1

⟨(~rn − ~rn+1)2

⟩= (N − 1)

⟨(~r1 − ~r2)2

⟩= (N − 1) kBTζ

(3.7)

Put another way, one has: ⟨~R2⟩∝ (N − 1)2ν , (3.8)

whereN−1 is the number of bonds. The exponent ν describes how the entangledchain grows when N is increased. ν is called the Flory exponent.

ν =12

(3.9)

This exponent is characteristic for the Gaussian chain, because the Gaussianchain belongs to the universality class of random walks.

15

3.1.2 Radius of gyration

In experiments, it is next to impossible to track a single monomer. Instead,a meaningful observable is the radius of gyration, which describes the meandiameter of the entangled polymer and can be determined by light scattering.It is defined as:

R2g =

12N2

N∑n=1

N∑m=1

⟨(~rn − ~rm)2

⟩. (3.10)

In [3] it is shown that

R2g =

1N

∑n=1

N⟨

(~rn − ~RG)2⟩, (3.11)

which is easier to compute. Furthermore, it is derived that

R2g =

16

(N − 1)kBT

ζ=

16⟨(~r1 − ~rN )2

⟩, (3.12)

which is only valid for sufficiently long chains. Take a chain consisting of onlytwo monomers as an example of a short chain:

R2g = 1

2·22

∑2n=1

∑2m=1

⟨(~rn − ~rm)2

⟩= 1

4

⟨(~r2 − ~r1)2

⟩,

(3.13)

so the ratio of⟨

(~R)2⟩

and R2g is constant only for large chains. Thus for long

chains the scaling law (eqn. 3.8) holds, as well.

3.2 Modifications of the model

Whereas the above model has the advantage of being solvable analytically, thereare some aspects in which it is unrealistic. As a consequence, some modifica-tions have to be made in order to obtain a more realistic description. Thesemodifications will be explained now.

3.2.1 The repulsive Lennard-Jones potential

The most obvious error of the model is that the monomers do not interact witheach other. Thus, they can overlap in any way. In nature, the DNA monomerscannot overlap with each other due to the Pauli interaction. The simplest ideawould be to demand that |~ri − ~rj | ≥ rmin for all i 6= j. However, this is notpossible when simulating a dynamics that includes velocities because one hasto calculate forces acting on each particle, so the potentials have to be soft,because derivatives have to be computed. The most widely used potential is theLennard-Jones potential, which we shall use to implement interactions between

16

any two monomers of the chain [12]. To reduce the computational expense, thepotential is cut off at the equilibrium point and shifted so that the transitionbetween both is smooth:

VLJ(r) =

{εLJ

[(σLJr

)12 −(σLJr

)6]− Vcut : r ≤ Rcut0 : r > Rcut

(3.14)

where

Rcut = 216σLJ (3.15)

Vcut = −εLJ4. (3.16)

The interaction between monomers affects the folding of the whole chainbecause monomers which are removed from each other in terms of the chaincontour will still interact if they come too close in real space (see fig. 3.2).

rnrm

Figure 3.2: The excluded volume

This effect is called the excluded-volume effect. In [3] a brief estimate showsν = 3

5 , a better calculation gives [7, 10]:

ν = 0.588± 0.001. (3.17)

Comparing this exponent to the one of the non-interacting chain (ν = 12 ), it is

evident that the folding of very long chains cannot be neglected.

3.2.2 Bond length potential

The next unrealistic behavior of the model lies with the springs: the potentialdoes not forbid thermal motion to elongate the bond to a great extent, becausethe Boltzmann distribution does not prevent the energy from rising above acertain value. The probability only becomes small. Thus, the Lennard-Jonespotential cannot keep the polymer strings from crossing each other (see fig. 3.3).If you imagine DNA as a thin string, it is obviously forbidden for the chains tocross each other since the Pauli interaction forbids such a behavior. This effectis called the entanglement interaction.

The first idea to solve this problem would be to introduce a repulsive inter-action between all bonds as it was done for the monomers, but this is almost

17

Figure 3.3: The entanglement interaction

impossible, because we would have to calculate the distance of two bonds comingclose to one another.

A far better idea is to limit the length of the bonds and to use the repulsiveinteraction between the monomers to avoid that the chains cross each other.Since the Boltzmann distribution does not forbid high energies, it is necessaryto replace the spring potential by a potential featuring a singularity at themaximum elongation of each bond. The potential used in the simulation wasthe FENE potential which was also used in [12]:

VBL(r) = −εBL2d2BL · ln

(1−

(d− d0

dBL

)2)

(3.18)

This potential limits the bonds to a maximum of dmax = d0 + dBL. Fur-thermore, it also offers the option of introducing an equilibrium length d0 forthe bonds. However, when using a repulsive Lennard-Jones potential betweenall monomers, the repulsive interaction between two neighboring monomers willestablish an equilibrium length, as well, because the bond length potential triesto contract the bond while the Lennard-Jones potential tries to extract it. Adisadvantage of this potential is the singularity in combination with a finite timestep: each monomer moves a finite distance at each time step, so for large timesteps, it is possible that two neighboring monomers move too far away fromeach other. Thus, the bond is broken because the potential is repulsive for bondlengths greater than the limit. In the simulation program, there is a check forbroken bonds which stops the program immediately when necessary.

Another possibility is to introduce very strong springs, so that the minimumenergy for the chains to cross each other is much greater than the thermal energy.This was done in the simulation to speed up calculations. In the simulations,the minimum energy for the chains to cross illegally was roughly 80 times thethermal energy. This means that the Boltzmann factor of the probability forthe chains to pass each other is roughly 10−35, which is very improbable. Thisis described in chapter 5.4.

The bond length potential (eqn. 3.18) has been implemented and tested butnot yet simulated.

18

3.2.3 Bond-angle potential

The last modification concerns the fact that DNA is an elastic polymer. Thismeans that neighboring bonds always try to align in a way that the monomersform a straight line. Thus the polymer becomes stretched on a comparativelysmall length scale, the so-called persistence length. In fig. 3.4, it is obvious thatthe persistence length of the chain is around one bond length on the left-handside whereas on the right it is several bond lengths.

Figure 3.4: The persistence length

The potential used in the program was

VBA = εBA (1− cos θab) (3.19)

where cos θab can be calculated from

cos θab =~a ·~b|~a||~b|

. (3.20)

To compute this potential, it is inevitable to calculate the square root ofboth bonds ~a and ~b. When the bond length potential is used, as well, it ispossible to reduce the computational effort because each bond length has to becalculated only once. In [11] other potentials were discussed, as well, but therecommended solution by Kolb is to use equation 3.19.

Experimental data shows that the persistence length of DNA is not constantfor all situations. It depends on the ion concentration in the buffer[20]. Whenchanging εBA, it is possible to adjust the persistence length to the experimentalsituation. For a detailed discussion of the persistence length, see [17, 18].

The bond-angle potential (eqn. 3.19) has been implemented and tested butnot yet simulated.

3.2.4 Electric field

Because DNA is a charged polyelectrolyte, it is easy to apply a dragging force byusing an electric field. The charge of the DNA is equally distributed along thechain. This leads to a dragging force that scales with the number of monomers.Considering eqn. 3.6, one sees:

µ ∝ D · f ∝ 1N·N = const. (3.21)

19

This means that in free solution the mobility is independent of the chainlength. This effect is reproduced by the model and is checked in chapter 5.3.

3.2.5 Further aspects

When examining DNA in water, the DNA strand is always surrounded by watermolecules. Therefore, when the DNA is moving, it will always produce somedrag on the nearby water molecules, and, if present, DNA monomers. Thisis called the hydrodynamic interaction. It is described in [2]. Another effectconcerning the mobility is caused by the counterions surrounding the DNApolymer. When applying an electric field, the counterions are always acceleratedin the opposite direction of the DNA monomer. Thus, they hinder the movementof the DNA. Neither the hydrodynamic nor the counterion interaction has beenimplemented yet.

Nevertheless, experimental data shows that the mobility of DNA in freesolution does not depend on the chain length, so these two effects seem tocompensate each other. This is described in [19].

The Debye screening length of DNA in water is small compared to the radiusof gyration, which is described in chapter 1.1. Thus, electrostatic interactionsof the monomers along the chain can be neglected. However, when applying anelectric field, the simulation is very simple because the electric field simply actsas another force interacting with each monomer.

20

Chapter 4

Development of theProgram

4.1 Programming details

The main objective in developing the program was to obtain a high degree ofefficiency.

4.1.1 Boxing

The main problem during the development of the program was the efficientcalculation of the Lennard-Jones potential described in chapter 3.2.1. This po-tential acts between all monomers of the chain, so a simple algorithm comparingthe locations of all monomers would take O(N2) comparisons. Doing so, thesimulation of a chain that is twice as long would take four times much time,which is very bad.

The next idea was to divide the location space into cubes of equal size andlinear length equal to the cutoff radius of the interaction. For each box thereexists a list containing all monomers within the box. From the position of eachmonomer it is possible to compute the box coordinates and then check the boxitself and all direct neighbors for nearby monomers. Because the density islimited due to the repulsive potential, this takes only constant time for onemonomer. It is convenient not to have boundary conditions when studying themovement of a particle. Then it is easy to see that it is no longer possible todivide the whole space into cubes and store a list for each cube, because memoryis limited. Hence, it was necessary to store only occupied boxes.

Key generation

When introducing boxing, it is necessary to create a key from the coordinatesof the monomer in order to find the appropriate box. In this case, it is possible

21

to compute the a key consisting of three integers in the following way:

ix = bx/rcutc

iy = by/rcutc

iz = bz/rcutc

(4.1)

Here, rcut is the cutoff radius of the Lennard-Jones potential (eqn. 3.14).When computing the Lennard-Jones potential, it is necessary to take all 26direct neighboring boxes into account, but no additional boxes.

Binary search tree

It is necessary to access the boxes from each monomer’s position. The simplestidea would be to use a sorted list of all boxes, but then searching within the listwould be O(N) for a single monomer, leading to a time scaling equal to whenno boxes are used at all. Therefore, in the first implementations of the program,a binary search tree was used.

00

00

00

00

0000

00

00

Root of tree

44

22

11

33

66

77

88

Figure 4.1: The binary search tree

Each node of a binary search tree contains the search key, links to the furthernodes, and the data (see fig. 4.1). Furthermore, it is necessary to have anordering relation for the keys. In the program the ordering used was simplyto sort the integer coordinates first by the x coordinate, then by the y andat last by the z coordinate; this is an ordering relation for three dimensionalcoordinates. In a binary search tree data, is arranged such that from each node,the upward link leads only to keys with smaller values and the downward link togreater values of the key. Searching in a binary tree always starts from the rootof the tree. By comparing the keys within each node and deciding whether tostop the search (keys are equal) or following the links up or down, it is possibleto tell whether a node is stored in the tree. When a link does not point to anode, it has to be marked as void and the search has to be stopped.

It is easy to see that searching in an optimal binary tree takes time likeO(log2N), which is far better than a simple list. Still, the binary search treehas two major disadvantages:

22

• The search time goes like O(lnN), which is clearly better than O(N).Nevertheless, for large numbers of stored boxes, most time is spent onsearching nodes in the tree, which is unsatisfying.

• The time scaling like O(lnN) is only valid for a perfectly balanced tree.In the worst case, the tree behaves like a single linked list and the searchtime rises to O(N). In this case it becomes necessary to rebuild the tree.This takes time like O(N), but only needs to be done when the tree hasbecome unbalanced.

For a detailed description of the binary search tree, see [21].

Hash map

In [21], a better method of organizing links to the boxes is introduced:Assume a large key space with only very few keys that are actually occupied.

Assume that there exists a function (the hash function) which can be easilycomputed, gives different values for almost all keys and has values in a small,limited range of integers. Then it is possible to store the occupied keys in acomparatively small array with the following algorithm:

1. For the key, compute the value of the hash-function.

2. Use the hash value of the key as the starting point in the array.

3. If the spot is occupied, proceed to the next until an empty spot is found.

4. Store the key and the data.

This way of organizing data is called a hash map. The data for a specific keycan be obtained in a similar way. If a key is not in the hash map, then an emptyplace will be found before finding the key.

Assuming that searching in step 3 is constant for a given hash map, it iseasy to see that an access via the hash map shows a constant time behaviorindependent of the number of keys stored. The crucial point in finding theobject is step 3: when the hash map is very full, it is likely that many otherkeys have to be scanned to find the proper key. It is important to keep an eyeon the filling level of the hash map: if it gets too full, it has to be rebuilt witha larger array.

The standard template library includes an implementation of a hash map,which was used in the simulation program. All you need is an equal key anda hash function. The hash function simply has to give an integer value foreach key, independent of the size of the hash map. The integer value should beequally distributed on the range of all valid integers. Hence the hash functioncan be computed very easily. The hash function used in the program was:

fhash(ix, iy, iz) = ix % 220 +(iy % 220

)· 220 +

(iz % 220

)· 240, (4.2)

23

where % denotes the modulo division for integers. Considering that the monomersof the polymer are quite close to each other, it is easy to see that this functiondoes not produce many collisions for different keys because the polymer willform a coil which is small compared to the repetition length of the hash func-tion. Hence this hash function only takes those bits of the key into considerationwhich change along the chain.

Further optimization of the boxing

With a binary tree, most time is spent on searching keys in the tree, because allother steps have constant time behavior per monomer and time step. Looking atthe computation of the Lennard-Jones potential, it is apparent that only neigh-boring boxes are regarded. So the most obvious optimization is to introducelinks to the nearest occupied neighbors in each box, reducing the time to searchkeys in the tree. This optimization is kept in the hash map, too. Because it issometimes necessary to rebuild the hash map, only links to the box objects arestored in the hash map. This made sure no link to a neighbor is broken duringrebuilding.

Even though searching in the tree could be reduced, it takes a lot of time,even for the hash map. When searching a cube, at first the key has to begenerated from the coordinates of the monomer, then the hash function has tobe computed and at last the box has to be located in the hash map. Duringthe simulation, the monomers rarely change their box. Thus, much computationtime is spent on searching links to boxes which have been computed many timesbefore. One sees that a list containing the links to the boxes for each monomersaves computation time.

4.1.2 Prefactors

In the first versions of the program, the changes of the velocity of each monomerwere computed by simply multiplying the time step and the force acting on themonomer. As an example, let us assume the potential is of the form V =εV · Vi(~r), so that εV is the prefactor of the potential which can be adjustedfrom the outside:

∆~v =[−εV ~∇Vi(~r)

]·∆t (4.3)

Computing this potential the following effort is required:

1. ~∇Vi(~r) is computed. In this general form, this cannot be simplified orestimated.

2. −εV ~∇Vi(~r) is computed. This takes three floating point operations.

3.[−εV ~∇Vi(~r)

]·∆t is computed.

24

Changing the order of computation, equation 4.3 becomes

∆~v = [−εV ∆t] · ~∇Vi(~r). (4.4)

It is easy to see that−εV ∆t is constant during the simulation. Thus, it is efficientto compute this prefactor during initialization of the program. Calculating theeffort, one sees that three floating point operations have been eliminated. Usingpotentials which only depend on the distance of two particles, the computingtime reduction becomes even bigger.

4.1.3 Random number generator

In this program, the random numbers were obtained by using the GNU ScientificLibrary (gsl) of version 0.7. The random number generator used was the ”taus”algorithm. For more detailed information on the gsl in general and the randomnumber generator in particular, see [6].

4.2 Testing the program

During the development of the program, the code was tested regularly in orderto find errors as early as possible. All potentials implemented had to fulfillconservation of energy and momentum when friction and Brownian motion weredisabled.

In testing the program, the Gaussian chain model proved very valuable,because it is solvable analytically. Thus, almost all tests concerning expectationvalues were carried out using the Gaussian chain model.

When applying thermal fluctuations and friction, the energy had to be pro-portional to the temperature. For a Gaussian chain of length n the thermalenergy was found to be

〈En〉 =(n− 1

2

)· 3kBT, (4.5)

which can be easily found by using the equipartition theorem and the fact thatthere are n monomers and (n − 1) bonds and each of which representing onedegree of freedom. For great values of the time step ∆t the error in the energydescribed in chapter 2.2.1 became visible; using the corrected temperature theenergy was right but the equation of the diffusion (eqn. 3.6) was no longerfulfilled.

4.2.1 End-to-end vector

An important property of the Gaussian chain is the scaling law presented in eqn.3.8. This was checked during development of the program. The expectationvalue of the end-to-end vector plotted against the number of monomers is givenin fig. 4.2 for k = 1, T = 1 and kB = 1. The line shown was obtained by linearregression of the data.

25

1 10 100Number of bonds

1

10

100E

nd-t

o-en

d ve

ctor

SimulationRegression ν = 0.498 ± 0.002

Figure 4.2: End-to-end vector plotted against chain length

Apparently, the values form a straight line in a double logarithmic plot.Furthermore, the regression shows

ν = 0.498± 0.02, (4.6)

which is in perfect agreement with the prediction ν = 12 for a Gaussian chain.

Looking at the length of the polymer consisting of only one bond, one findsthat the length is given by

lbond = 1.73± 0.02. (4.7)

This corresponds to an energy of the bond equal to

Ebond =12kl2bond = 1.50± 0.02, (4.8)

which is in perfect agreement with the prediction given above.

4.2.2 Rouse-mode analysis of a Gaussian chain

Another important point that had to be checked were the long time correlationsof the Gaussian chain computed in appendix A. The following set of parameterswas used for the simulations:

26

T = 1 (4.9)kB = 1 (4.10)k = 2 (4.11)ζ = 4 (4.12)

∆t = 0.01 (4.13)

The correlations were computed using a block analysis over ten blocks. Theresults obtained from simulation as well as the expected correlation function aregiven in figures 4.3, 4.4 and 4.5 for modes 1, 4 and 9, respectively. Note thatall correlations were normalized such that 〈X2

p,i(0)〉 = 1. Using the equationsderived in appendix A, ones sees:

• Mode 1 (fig. 4.3) is damped so strongly that no oscillation is possible.The correlation decays only weakly as a function of the time.

• Mode 4 (fig. 4.4) is damped just strongly enough to represent the aperiodiclimit. The correlation decays very fast.

• Mode 9 (fig. 4.5) is only weakly damped because the effective springconstant rises as a function of the mode. Unfortunately, the damping isquite strong, so the correlation function does not show much more thanone oscillation.

All observed correlations are in good agreement with the expected ones. Thus,the program is capable of reproducing the dynamics of a Gaussian chain.

27

0 2 4 6 8 10Time

0.4

0.5

0.6

0.7

0.8

0.9

1

1.1C

orre

latio

n of

mod

e 1

SimulationExpected

Figure 4.3: Correlation of mode 1 plotted against time

0 2 4 6 8 10Time

0

0.2

0.4

0.6

0.8

1

Cor

rela

tion

of m

ode

4

SimulationExpected

Figure 4.4: Correlation of mode 4 plotted against time

28

0 2 4 6 8 10Time

0

0.5

1

Cor

rela

tion

of m

ode

9

SimulationExpected

Figure 4.5: Correlation of mode 9 plotted against time

29

Chapter 5

Properties of the Chain infree Solution

5.1 Parameter set for the chain in simulation

During the simulation, the parameters of the chain had to be kept constant inorder to get results which can be compared to each other. The parameters usedduring the simulation had the following values:

T = 1 (5.1)∆t = 0.01 (5.2)kB = 1 (5.3)ζ = 1 (5.4)

εLJ = 1 (5.5)σLJ = 1 (5.6)k = 100 (5.7)

Some of these parameters were initially set to unity during the developmentof the program. If those parameters turned out to give sensible results in thesimulation, they were kept in order to reduce the effort for searching new sets.This set of parameters showed the following properties:

• With this set one finds ζ∆t = 0.01, which is quite small (see chapter 2.2.1for details). Thus, a correction of the temperature can be neglected, andthe equation of diffusion (3.6) is fulfilled.

• The tight springs introduce entanglement interactions. This is describedin chapter 3.2.2 and is checked in chapter 5.4.

30

• The mobility of a chain can be computed via

µ =1ζ

= 1 (5.8)

and has to be constant for all chains. This is checked in chapter 5.3.

• By introducing excluded volume and entanglement interactions a smallpersistence length is introduced as well. This is described in chapter 5.5.

Figure 5.1: Snapshot of the polymer in free solution

A snapshot of 500 monomers in free solution is given if fig. 5.1.

5.2 Radius of gyration

Comparing experimental and simulation data, one has to identify observableswhich can be adjusted. In our case, the radius of gyration is very important,as it allows for a comparison of the length scales not only of the polymers, butalso of the geometry of the trap (chapter 6.2). The Radius of Gyration was alsoused for finding a proper geometry of the trap.

The radius of gyration was computed from chains in free diffusion with timesteps ranging from 50 · 106 up to 109. During the simulation, a snapshot wassaved regularly and appended to a file. This file was analysed later and theradius of gyration was obtained by a block analysis over ten blocks. The resultsare shown in fig. 5.2.

31

1 10 100 1000Number of Bonds

0

1

10

100R

adiu

s of

Gyr

atio

n

MeasurementExpected ν = 0.588

Figure 5.2: Radius of gyration for different chain lengths

For long chains, R2g should behave like N2ν , where ν = 0.588. By looking

at fig. 5.2, one notices that the radius of gyration does not give a straight linefor the small values of N in a double logarithmic plot. This can also be seenin fig. 5.3. The ratio changes for short chains and too many errors occur forlong chains. Thus the ratio cannot be assumed constant and a regression forestimating ν was not attempted.

Nevertheless, the radius of gyration could be determined even for long chainswith sufficiently low errors. This is used in chapter 6.2. Furthermore, the linefor the expected values, which was derived from the radius of the longest chainand ν = 0.588 fits very well to the values obtained from simulation.

5.3 Mobility in free solution

One important aspect which has to be reproduced by the model is describedin chapter 3.2.4. The mobility of the chain has to be independent of the chainlength and has to fulfill eqn. 5.8. The mobility is defined as:

µ =v(E)E

, (5.9)

where E denotes the electrical field. In fig. 5.4 the movements of chains withdifferent length in free solution with electrical field ~E = 0.01~ex are plotted.

32

1 10 100 1000Number of Bonds

3

4

5

6

7

8

Rat

io o

f End

−to−

End

vs.

Gyr

atio

n

Figure 5.3: Ratio of radius of gyration vs. end-to-end vector

Longer chains diffuse more steadily, which is easy to see. In fig. 5.5 themobility of the chains is computed using block analysis with five blocks. Thedotted line is the expected mobility (µ = 1). The diffusion of the shorter chainsleads to greater errors, as expected from fig. 5.4. Nevertheless, all chains showa chain length independent mobility, as demanded by eqn. 3.21. The mobilityof the longest chain shows a deviation greater than the error from the expectedvalue, which is due to the algorithm. This can be seen as follows:

1. The friction is computed via

~vi = (1− ζ∆t)~vi

2. The forces and Brownian motion is added:

~vi = ~vi +∆t

2~f(~r) +

√2kbT

ζ

m∆tη

3. The positions of the particles are updated:

~ri = ~ri + ∆t~vi

4. The forces are computed and a ∆t

2 time step is applied:

~vi = ~vi +∆t

2~f(~r)

33

0 20 40 60 80 100Number of timesteps (in million)

0

5

10

Pos

ition

of t

he c

ente

r of m

ass

(a. u

.) N = 10N = 100N = 1000

Figure 5.4: Movement in free solution

In step 1, the friction for the whole time step is computed and in step 2 onlyhalf a time step is computed for the force. Thus, the velocities of the particlesare slightly smaller than expected when updating the location of the particles.All chains are affected in the same way. This effect does not affect the results,because only differences between the mobilities are of interest.

5.4 Distribution of the bond length

In chapter 3.2.2 a different possibility than using eqn. 3.18 was presented.By introducing very strong springs (see chapter 5.1) the computation of thesquare-root in eqn. 3.18 can be avoided. It is now necessary to check if this newparameter set allows chains to cross or not. In order to do so, I first looked atthe potential:

Vbond(r) =k

2r2 +

{εLJ

[(σLJr

)12 −(σLJr

)6]− Vcut : r ≤ Rcut0 : r > Rcut

(5.10)

This potential is plotted in fig. 5.6 for the parameter set described above.A brief numerical calculation shows that the equilibrium of the length for twosingle bonds is given by

req ≈ 0.847 (5.11)Eeq ≈ 2 · 40.75 = 81.5 (5.12)

In order to find the minimum energy for two bonds crossing each other, thebonds have to stretch (see fig. 5.7).

34

10 100 1000Number of monomers

0.8

0.85

0.9

0.95

1

1.05

1.1

1.15

1.2

Mob

ility

ζ∆t

Figure 5.5: Mobility µ in free solution

The energy of the configuration presented in fig. 5.7 can be calculated from

E(r) = 4VLJ

(r√2

)+ 2Vbond(r), (5.13)

where r is the distance between the monomers connected by the springs. Theminimum of this potential was found to be:

rcross ≈ 1.198 (5.14)Ecross ≈ 163 (5.15)

Simulating at a thermal energy of kBT ≡ 1, the energy needed for crossingis roughly 80 times the thermal energy. This leads to a Boltzmann factor of10−35, which is very small. To make sure that no entanglement occurred insimulation, the lengths of the simulation of the chain with 1000 monomers infree solution were analyzed, giving 30 · 106 sample lengths. The distribution ofthe lengths and the correctly scaled Boltzmann factor are given in fig. 5.8.

The distribution reflects well the properties of the potential and the pre-diction of the Boltzmann factor. The minimum in the potential correspondsdirectly to the maximum in the distribution. Furthermore, the sharp raise ofthe potential for small values of r leads to a strong decay in the distribution.Values of r with r > 1 were also found in the samples. Within all samples,only 937 values of r were found to be larger than 1.0, but none larger than 1.1.This means that the elongation enabling minimum crossing energy was neverreached in a simulation, which is in good agreement with the prediction of theBoltzmann factor.

35

0.7 0.75 0.8 0.85 0.9 0.95 1Bond length

35

40

45

50

55

Pot

entia

l

Figure 5.6: Bond length potential with springs

Figure 5.7: Two bonds crossing each other

5.5 Persistence length

The chains are unable to cross each other, which implies a small persistencelength because the two neighbors of a single monomer cannot overlap. Thusboth bonds are correlated. This effect is shown in fig. 5.9.

In [9] it is shown that the projection HN of the end-to-end vector ~RE ontothe first bond ~R1 can be computed by

HN = ~e1 · ~RE =~R1

R1·N∑j=1

~Rj . (5.16)

The persistence length is then defined as

lp := limN→∞

〈HN 〉. (5.17)

Using the configuration data from chapter 5.2, it is possible to computeHN for various chain lengths and then estimate the limit for N → ∞. Theconfigurations were again analysed using block analysis over ten blocks. Thevalues found for various chain lengths are presented in fig. 5.10. It is easy tosee that HN reaches a limit for large values of N . This limit can roughly beidentified as

36

0.7 0.8 0.9 1Bond Length

0

5

10

15

Pro

pabi

lity

dens

ity fr

om s

imul

atio

n

SimulationBoltzmann

Figure 5.8: Distribution of the bond lengths

Figure 5.9: Folding of two neighboring bonds

lp ≈ 1.6 ≈ 1.9 · lbond, (5.18)

which is little less than two bond lengths but greater than the persistence lengthof an ideal chain which is equal to the bond length. For long chains the errorsincrease, because the first bond is barely correlated to the end-to-end vector.

Thus by enabling excluded volume and entanglement interactions a littlepersistence length was introduced as well.

37

10 100Number of monomers

1

1.5

2

2.5

Proj

ectio

n of

end

-to-

end

vect

or o

n fi

rst b

ond

Figure 5.10: Projection of the end-to-end vector on the first bond

38

Chapter 6

Setting up the entropicTrap

6.1 Simulation of the wall

In order to implement an entropic trap, the interaction potential of the monomerswith the walls has to be set up. The simplest case is an infinitely long, straightwall. The potential chosen was simply the repulsive Lennard-Jones potential:

Vwall(r) =

{εwall

[(σwallr

)12 −(σwallr

)6]− Vcut : r ≤ Rcut0 : r > Rcut

, (6.1)

where Vcut is adjusted to Rcut so that the change of potential at the cut offpoint is zero. The parameter set used was as follows:

εwall = 1; (6.2)σwall = 1; (6.3)Rcut = 21/6 · σwall (6.4)

Thus the walls were purely repulsive. A difficulty arises from the fact thatthe walls in entropic traps are not simply straight, but also have corners. Thisis shown in fig. 6.1.

In fig. 6.1(a) the walls are easy to implement: the effects of both walls simplysum up, assuring a smooth transition.

In the opposite case (fig. 6.1(b)), it is not useful to sum up the effects ofboth walls because they both have ended. The best one can do is use the pointat the corner as a interaction point for the monomer and then compute thepotential according to eqn. (6.1). This provides for a smooth transition at bothends.

39

++

++

(a) (b)

Figure 6.1: The corners of the walls

6.2 Geometric setup of the trap

Electric Field

LLBox

HHBox

LLCh

HHChXX

ZZYY

Figure 6.2: The entropic trap

Using chapter 6.1, it is possible to implement an idealized entropic trap. Tobe precise, the idealizations are:

• The walls are perfectly smooth and no extra friction occurs at the walls.

• All walls are perfectly aligned. This means in particular that the walls onthe left and right of each box are perfectly upright.

• The electric force applied is perfectly perpendicular to the walls on eachside of the box. Thus, any particle at those walls diffuses freely in z-direction, see fig. 6.2.

• There are no restrictions in the y direction. Thus the trap is infinitelystretched. In experiments, all traps have finite size.

• The trap is infinitely and periodically stretched in x-direction. This is nota great problem because the applied force as well as the simulation timeare finite, limiting the movement in the trap.

40

In fig. 6.2 a schematic drawing of the trap as well as the coordinate systemused are presented. With the idealizations given only four parameters are neededfor a full description of the trap. The first set of parameters used was:

HBox = 10 (6.5)LBox = 10 (6.6)HCh = 4 (6.7)LCh = 4 (6.8)

Figure 6.3: 1000 monomers in the initial trap - side view

0 10 20 30−100

−50

0

50

100

150

Figure 6.4: 1000 monomers in the initial trap - view from above

41

A simulation snapshot is presented in figures 6.3 and 6.4. In the side view,the boxes look very crowded and the polymer fills two boxes. This is furtherconfirmed in the view from above: looking at the scaling of the axes, it becomesclear that the boxes are too small for the polymer to form a coil. This setup isnot usable in an experiment because it is impossible to inject the polymer intothe trap. The next trial setup was:

HBox = 80 (6.9)LBox = 80 (6.10)HCh = 10 (6.11)LCh = 20 (6.12)

In order to find a proper geometry for the trap, it was inevitable to adaptthe sizes of the trap to the size of the polymer. In chapter 5.2 the radius ofgyration for various chain lengths has been computed. For a chain consisting ofN = 1000 monomers it was found that Rg = 28 ± 2, giving a diameter of thepolymer of roughly 60. Thus, a box with dimension 80 × 80 was expected todeliver good results. The channel has to have a height smaller than the radiusof gyration and was set to 10. Large chains are forced to unfold in order topass through. For simplicity, the length of the channel was set to 20, givingan overall length of 100 of the whole trap. A snapshot of a chain consisting ofN = 1000 monomers in the trap is shown in fig. 6.5: The chain is just about topass the channel. Note that only a part of the box is shown.

Figure 6.5: 1000 monomers in the revised trap

6.3 First runs in the trap

Now the time had come to see whether the geometry of the revised trap is usefulor not. In fig. 6.6, the movement of the center of mass for 108 time steps with

42

an electric field ~Eel = 0.01~ex is plotted. In free solution, all chains should havepassed 100 boxes.

0 20 40 60 80 100Number of timesteps (in million)

0

20

40

60

80

Pos

ition

of t

he c

ente

r of m

ass

(box

leng

ths)

N = 10N = 100N = 1000

Figure 6.6: Movement of the center of mass in the trap

The first issue to be noted is the fact that the long chain (N = 1000) diffusesquicker in the trap than all other chains. It is strange, because it was expectedthat long chains take longer time to stretch in order to pass the channel. Suchbehavior is observed experimentally as well (see [8]). The explanation given isthat longer chains have a greater chance of forming a loop that reaches into thechannel. This can be identified as the beginning of getting in line of the polymerfor passing the channel. Such a loop can also be seen in fig. 6.4.

Looking at the diffusion of the chain with N = 100, one notices that thereare passages in which the chain is hindered in motion completely and in otherpassages the chain moves over many boxes at almost maximum speed. Doing ablock analysis of such a curve in order to get a mobility will lead to great errors.

6.4 Tilted force

In order to get rid of chains which are trapped for very long times, it is necessaryto understand why the chains are trapped for such long times. This can beunderstood quite easily when looking at the idealizations: a chain located atthe wall does not feel the electric force any longer, because the electric field iscompensated by the wall. The only way to find the channel is to diffuse upwards.This is a purely random process which can happen in any direction. For chainsof medium length, one sees that the diffusion while passing the box is sufficientfor the chain to diffuse into the box. When the chain has reached the wall onthe right side, it takes quite a long time to find the channel by pure diffusion.For short chains, this effect occurs as well, but the diffusion is greater. Thus,

43

short chains diffuse quicker upwards to the channel and leave the box earlier.For long chains the diffusion constant breaks down like N−1. This induces theydo not have the time to diffuse into the box and travel very quickly in the trap.

To reduce trapping times, it is necessary to introduce some drag to themonomers in order to find the channel. The easiest way to accomplish such adrag is to tilt the electric field:

~Eel = 0.01~ex + 0.001~ez. (6.13)

Thus, the electric field in z direction will drag the polymer upwards to thechannel. On the other hand, the field is quite small and does not stick thepolymer to the upper ceiling of the trap. The drag certainly increases themobility of the chains. In order to compensate this effect the channel heightwas reduced. The following geometry was found to be a useful setting:

HBox = 80 (6.14)LBox = 80 (6.15)HCh = 7 (6.16)LCh = 20 (6.17)

0 20 40 60 80 100Number of timesteps (in million)

0

20

40

60

80

Pos

ition

of t

he c

ente

r of m

ass

(box

leng

ths) N = 10

N = 100N = 1000

Figure 6.7: Movement of the center of mass in the trap with tilted force

The movement of the center of mass driven by a tilted force is plotted in fig.6.7. One sees that the long trapping times of chains with medium length havevanished. Furthermore, it seems that chains of medium length now diffuse atmedium speed.

44

Unfortunately, I was informed by my experimental co-workers, AlexandraRos1 and Than-Tu Duong2, that it is almost impossible to benefit from a tiltedforce because one usually applies simply a voltage, thereby creating a non-tiltedforce on the polymer. Thus this idea was promising but unrealistic and had tobe discarded.

6.5 The trap with tilted walls

The production of an entropic trap is usually done by etching the structuresinto a surface. It is easy to see that perfectly upright structures are almostimpossible to manufacture. In order to add some reality to the model, the wallsof the trap are tilted slightly. This is shown in fig. 6.8.

αα

Figure 6.8: Trap with tilted walls

The black line shows the original setup of the trap. For simplicity, in simu-lations the tilt was represented by the inverse slope of the wall, twall. One seesthat

tanα = twall, (6.18)

so for example a tilt of twall = 0.1 corresponds to an angle of α ≈ 5.71◦. Inthe first simulation of the trap with tilted walls, the tilt twall was set to 0.1,which was expected to reproduce the results from chapter 6.4. In fig. 6.9 themovement of the center of mass for an electric force ~Eel = 0.01~ex is shown.

It is apparent that the trap equipped with tilted walls shows properties whichare almost similar to the properties of the trap described in chapter 6.4. Thus,this set of parameters is promising and can be reproduced in experiments.

1Dr. Alexandra Ros, Experimentelle Biophysik and angewandte Nanowissenschaften,Fakultat fur Physik, Universitat Bielefeld, Germany

2Than-Tu Duong, Experimentelle Biophysik and angewandte Nanowissenschaften,Fakultat fur Physik, Universitat Bielefeld, Germany

45

0 20 40 60 80 100Number of timesteps (in million)

0

20

40

60

80

Pos

ition

of t

he c

ente

r of m

ass

(box

leng

ths) N = 10

N = 100N = 1000

Figure 6.9: Movement of the center of mass in the trap with tilted walls

46

Chapter 7

Mobilities in the Trap

In this chapter we shall now discuss the diffusion of the polymer in the trap withtilted walls for various chain lengths N , wall tilt parameters twall and electricforces ~Eel. To be precise, the settings simulated covered all combinations of

• Chain lengths of N = 10, 20, 50, 100, 200, 500 and 1000. Thus the chainlengths are uniformly distributed on a logarithmic scale.

• Electric forces of Ex = 0.0025, 0.005, 0.01, 0.02 and 0.04. Thus the electricforces are distributed uniformly on a logarithmic scale as well.

• Tilt parameters of twall = 0, 0.01, 0.02, 0.05, 0.1, 0.15 and 0.2.

The mobilities were then obtained by block analysis over the movement ofthe center of mass and the time (for example, please have a look at fig. 6.6).For simplicity the number of blocks was set to five for each run. Each runwas simulated over at least 2 · 108 time steps, the maximum number of timesteps simulated was 109. For a better understanding, only parts of the resultsare given in this chapter, all results obtained from simulation can be found inappendix B.

7.1 Tilt dependence of the mobility

The first thing that already became obvious in chapter 6.5 is the change ofmobility for chains of medium length when tilting the walls because the longtrapping times vanished for tilted walls. In figures 7.1 and 7.2 the mobilityfor various chain lengths is plotted against the tilt of the wall at electric forcesEx = 0.01 and Ex = 0.04, respectively. The mobility depending on the tilt foran electric force of Ex = 0.0025 was not plotted because the force is too weakto pull any long or medium chain into the channel.

Looking at both figures, one sees:

47

• The mobility of the long chain depends only weakly on the tilt of the wall.Only the transition from twall = 0 to twall = 0.01 at Ex = 0.01 shows arather small change in mobility.

• For chains of length N = 100, the change in the mobility when varyingthe tilt of the wall is very big. When the electric force Ex is set to 0.01and the tilt is set to twall = 0, the mobility is smaller than the mobility ofthe shortest chain simulated. When increasing the tilt to twall = 0.2, themobility is higher than the mobility of the short chain. For this tilt, themobility of the chain of length N = 100 changes by a factor of roughly 7,only by changing the tilt of the walls.

This can also be seen when looking at the behavior of the mobility de-pending on the tilt at a force of Ex = 0.04: For zero tilt the chain is onlylittle faster than the shortest chain, but for high tilts the chain travels atspeeds almost equal to those of the long chain.

• Chains of length equal to N = 10 do not show such a strong dependenceon the tilt of the walls because the diffusion in z direction is greater. ForEx = 0.01 almost no dependence is visible and for Ex = 0.04 the mobilityrises by a factor of roughly 3, which is small compared to the raise of themobility of the chains with N = 100.

As a result, one can say that by tilting the walls the mobility of the chainsof medium length is most affected. The mobility of short chains shows only arather weak and that of long chains almost no dependence on the tilt of thewalls.

7.2 Force dependence of the mobility

Another parameter varied in simulation was the electric force applied. As thetime needed for passing a box decreases when the force applied increases, it wasexpected that the mobility rises as well, because the time for diffusion into thebox breaks down. In figures 7.3, 7.4 and 7.5 the mobility is plotted as a functionof the electric force for various tilts and chain lengths. Summarizing the plotsfor different chain lengths, one sees:

• The mobility of the short chains depends only weakly on the electric fieldapplied. For untilted walls the mobility does not vary much more thanthe errors observed. At the maximum tilt, the mobility shows only a weakdependence on the electric field.

• Chains of medium length have no mobility for the weakest force appliedand show a medium dependence on the force applied. For untilted wallsthe mobility is rather similar to the one shown by the short chains. Atgreat tilts, the mobility behaves more like the one of the long chains. Thiswas already observed in chapter 7.1.

48

• Long chains show a great dependence on the electric field applied. Forlow forces, the chains do not move at all. When the force is increased,the mobility changes quickly and for the greatest forces applied the chainstravel at almost maximum speed. Thus almost no part of the chain hasthe time to diffuse into the box.

As a conclusion one can say that the effect of the strength of the electricfield increases with the length of the chain. Furthermore, in order to move longchains, it is necessary to apply a sufficient force, because otherwise the entropicforce produced by the temperature is stronger and prevents the polymer fromdiffusing into the channel.

7.3 Separation by chain length in the trap

Using the effects we learned in the previous chapters, it is now possible to analyzethe mobility of chains inside the trap as a function of their length. A few plotsare shown in figures 7.6, 7.7, 7.8 and 7.9, all plots are given in appendix B.When looking at the plots, it becomes evident that all traps have a maximumchain length for effective separation. The setting twall = 0, Ex = 0.01 shows aminimum of the mobility for a chain length equal to N = 100. This effect isexplained in chapter 6.4. The effective ranges of the different trap setups aregiven in table 7.1. Note that the great error in the increase for the first setsimply arises from the fact that the mobility of the short chain is almost zero.Within the ranges given the mobility rises as a function of the chain length andthe errors are comparably small.

Examining the effective sets of parameters, one sees that the low forces donot yield useful values. Introducing tilted walls shifted the effective ranges toshorter chains, because the mobility is increased. Furthermore, the increase inmobility decreases as a function of the tilt of the wall. On the other hand, tiltedwalls lead to smoother diffusion inside the trap and therefore quicker results.By applying higher electric fields, the effective range is shifted to shorter chainlengths and the increase in mobility decreases as well.

Nevertheless, all effective settings allow to separate chains by length over atleast one order of magnitude in the the chain length. Thus, these settings areuseful for separating chains by length.

49

twall Ex min. length max. length increase of mobility0 0.01 100 1000 (500 ± 320)%0 0.02 100 1000 (290 ± 110)%0 0.04 50 500 (700 ± 300)%0.01 0.02 100 1000 (220 ± 50)%0.01 0.04 50 500 (290 ± 70)%0.02 0.02 20 500 (290 ± 80)%0.02 0.04 20 500 (250 ± 80)%0.05 0.01 50 500 (150 ± 40)%0.05 0.02 20 500 (190 ± 40)%0.05 0.04 10 200 (180 ± 40)%0.1 0.01 50 500 (110 ± 30)%0.1 0.02 10 500 (160 ± 40)%0.1 0.04 10 200 (111 ± 9)%0.15 0.01 10 500 (160 ± 50)%0.15 0.02 10 200 (110 ± 20)%0.15 0.04 10 100 (60 ± 4)%0.2 0.01 10 200 (110 ± 30)%0.2 0.02 10 200 (70 ± 20)%0.2 0.04 10 100 (48 ± 5)%

Table 7.1: Effective ranges of the traps

0 0.05 0.1 0.15 0.2Tilt of the wall

0

0.1

0.2

0.3

0.4

0.5

0.6

Mob

ility

N = 10N = 100N = 1000

Figure 7.1: Mobility plotted against the tilt of the wall for Ex = 0.01

50

0 0.05 0.1 0.15 0.2Tilt of the wall

0

0.2

0.4

0.6

0.8

Mob

ility

N = 10N = 100N = 1000

Figure 7.2: Mobility plotted against the tilt of the wall for Ex = 0.04

0.001 0.01 0.1Electric force

0

0.2

0.4

0.6

0.8

Mob

ility

N = 10N = 100N = 1000

Figure 7.3: Mobility plotted against the force for twall = 0

51

0.001 0.01 0.1Electric force

0

0.2

0.4

0.6

0.8

Mob

ility

N = 10N = 100N = 1000

Figure 7.4: Mobility plotted against the force for twall = 0.1

0.001 0.01 0.1Electric force

0

0.2

0.4

0.6

0.8

Mob

ility

N = 10N = 100N = 1000

Figure 7.5: Mobility plotted against the force for twall = 0.2

52

10 100 1000Number of monomers

0

0.2

0.4

0.6

0.8

Mob

ility

Ex = 0.04Ex = 0.01Ex = 0.0025

Figure 7.6: Mobility plotted against the chain length for twall = 0.

10 100 1000Number of monomers

0

0.2

0.4

0.6

0.8

Mob

ility

Ex = 0.04Ex = 0.01Ex = 0.0025

Figure 7.7: Mobility plotted against the chain length for twall = 0.02

53

10 100 1000Number of monomers

0

0.2

0.4

0.6

0.8

Mob

ility

Ex = 0.04Ex = 0.01Ex = 0.0025

Figure 7.8: Mobility plotted against the chain length for twall = 0.1

10 100 1000Number of monomers

0

0.2

0.4

0.6

0.8

1

Mob

ility

Ex = 0.04Ex = 0.01Ex = 0.0025

Figure 7.9: Mobility plotted against the chain length for twall = 0.2

54

Chapter 8

Conclusions and Outlook

In this thesis, it is shown that chains with the properties described in chapter 5can efficiently be separated by length in the trap introduced in chapter 6. Theranges of effective separation cover at least one order of magnitude and reachup to almost two orders of magnitude in the length of the chain.

Furthermore, it is shown that by varying the electric field the movement ofthe long chains is most affected, whereas varying the tilt of the walls affectsprimarily the movement of the chains of medium length. Thus, it is possible toadjust the movement by varying the electric force and the tilt of the walls. Byreducing the tilt of the walls to zero, it is possible to find chains of intermediatelength which travel at lower speed than all other chains. So it is possible to sortout chains of intermediate length.

The setting of the trap can be adjusted over a wide range of parameters.When looking for an effective way to separate long chains, it is now necessaryto find settings which are useful for longer chains.

DNA is a rather stiff polymer with persistence lengths of typically lp = 45nmin double-stranded form[18], which is much more than the two bond lengthsshown by the settings used in this simulation. Thus, the bond angle potentialhas to be enabled in order to adjust the persistence length to experimentalvalues in future simulations.

55

Appendix A

Rouse-mode analysis of thebead spring model

A.1 Rouse-modes of the bead-spring model

In order to solve equation 3.1 analytically, it is common to transform normalcoordinates so that each mode is independent of each other’s motion. Obviously,the original equation (3.1) does not fulfill this requirement since each monomerfeels attracted to both neighbors. Thus, the problem is to find the eigenstatesof the equation set. The common way to handle this problem is to translate thediscrete polymer chain into a continuous one. This is described in [2], p. 131f.for the following differential equation:

ζ~rn = −k (2~rn − ~rn+1 − ~rn−1) + ~ηn (A.1)

with the boundary conditions equal to 3.2 and 3.3 where the random forcessatisfy

〈~ηi(t)〉 = 0 (A.2)〈~ηi,α(t)~ηj,β(t′)〉 = 2ζkBT δij δαβ δ(t− t′) (A.3)

Transforming into continuous variables gives

ζ~rn = εSP∂2~rn∂n2

+ ~ηn (A.4)

with the boundary conditions ∂~rn∂n

∣∣n=0

= 0 and ∂~rn∂n

∣∣n=N

= 0. Finally, theRouse-modes are defined as:

~Xp ≡1N

∫ N

0

dn cos(pπnN

)~rn(t) , p = 0, 1, 2 . . . (A.5)

Then eqn. A.4 can be rewritten as (see [2], p. 94):

56

ζp∂

∂t~Xp = −kp ~Xp + ~ηp (A.6)

where

ζp = (2− δp0)Nζ (A.7)kp = 2π2kp2/N (A.8)

Unfortunately, the differential equation simulated is of second order and thetransition to continuous variables is not possible. Nevertheless, equation A.5hints to the correct Rouse-modes of a discrete chain:

~Xp =1N

∑n

cos(p(n+ 1

2 )πN

)~rn (A.9)

It follows:

~Xp = 1N

∑n cos

(p(n+ 1

2 )π

N

)~rn

= 1N

∑n cos

(p(n+ 1

2 )π

N

){−ζ~rn − k (2~rn − ~rn+1 − ~rn−1) + ~ηn

}= − ζ

N

∑n cos

(p(n+ 1

2 )π

N

)~rn + 1

N

∑n cos

(p(n+ 1

2 )π

N

)~ηn

− kN

∑n cos

(p(n+ 1

2 )π

N

)(2~rn − ~rn+1 − ~rn−1)

= −ζ ~Xp + 1N

∑n cos

(p(n+ 1

2 )π

N

)~ηn

− kN

∑n ~rn

{2 cos

(p(n+ 1

2 )π

N

)− cos

(p(n+ 3

2 )π

N

)− cos

(p(n− 1

2 )π

N

)}= −ζ ~Xp − 4k

N sin2(pπ2N

)∑n cos

(p(n+ 1

2 )π

N

)~rn + 1

N

∑n cos

(p(n+ 1

2 )π

N

)~ηn

= −ζ ~Xp − 4kN sin2

(pπ2N

)~Xp + 1

N

∑n cos

(p(n+ 1

2 )π

N

)~ηn

(A.10)Since it is not to be expected that the effective mass of a Rouse-mode is equalto unity, it is necessary to have:

~Xp = − ζpmp

~Xp −kpmp

~Xp +1mp

~ηp (A.11)

From equations A.10 and A.11, it is easy to see:

1mp

~ηp =1N

∑n

cos(p(n+ 1

2 )πN

)~ηn (A.12)

This equation results to:

57

〈ηpα(t)ηqβ(0)〉 = m2p

N2

∑n

∑m cos

(p(n+ 1

2 )π

N

)cos(q(m+ 1

2 )π

N

)〈ηnα(t), ηmβ(0)〉

= m2p

N2 2ζkBTδαβδ(t)∑n cos

(q(n+ 1

2 )π

N

)cos(q(n+ 1

2 )π

N

)= m2

p

N2 2ζkBTδαβδ(t)∑n

12

{cos(

(p+q)(n+ 12 )π

N

)+ cos

((p−q)(n+ 1

2 )π

N

)}= m2

p

N2 2ζkBTδαβδ(t)N 12δpq (1 + δp0)

(A.13)where 3.5 was substituted. Furthermore, it is necessary to demand:

〈~ηp〉 = 0 (A.14)〈ηp,α(t)ηq,β(0)〉 = 2ζpkbTδpqδαβδ(t) (A.15)

Equation A.14 is obviously satisfied due to equation 3.4. By comparing equa-tions A.10, A.11, A.13 and A.15, it is possible to see:

ζp =m2p

Nζ

1 + δp02

(A.16)

kpmp

= 4k sin2( pπ

2N

)(A.17)

ζpmp

= ζ (A.18)

These formulae can easily be solved:

mp = (2− δp0)N (A.19)ζp = (2− δp0)Nζ (A.20)

kp = 8kN sin2( pπ

2N

)(A.21)

mp is called the effective mass, ζp is called the effective friction and kp is calledthe effective spring constant for each Rouse-mode.

A.2 Time evolution of the different modes

Using equations A.11, A.14 and A.15, the computation of the time dependenceof the different Rouse-modes is possible.

58

A.2.1 Movement of the center of mass

First, mode 0 is investigated. ~Xp now becomes

~X0 =1N

∑n

cos(

0(n+ 12 )π

N

)~rn =

1N

∑n

~rn, (A.22)

which is simply the center of mass ~RG. Since kp = 0, there is no effective springconstant for the center of mass. Thus all random motion is summed up andleads to diffusion:

〈(X0α(t)−X0α(0)) ((X0β(t)−X0β(0))〉 = δαβ2kBTζ0

t = δαβ2kBTNζ

t (A.23)

With this equation it is simple to calculate the diffusion of the center of mass:

⟨(~RG(t)− ~RG(0)

)2⟩

=∑

α=x,y,z

⟨(X0α(t)−X0α(0))2

⟩= 6

kBT

Nζt (A.24)

which is in agreement with eqn. 3.6.

A.2.2 Correlation of the higher modes

The differential equation for each Rouse-mode has the following form:

mp~Xp = −ζp ~Xp − kp ~Xp + ~ηp (A.25)

This is a damped harmonic oscillator, driven by a random force. This is solvedin [2], p. 62f. For a potential U = 1

2kx2, the mean excitation from ground state

is given by ⟨x2⟩

=kBT

k. (A.26)

One now finds

〈Epot〉 =12k⟨x2⟩

=12kBT, (A.27)

which agrees with the Boltzmann statistics for a particle with one degree offreedom. It is now possible to use eqn. A.26 to determine

⟨X2p

⟩for A.25:⟨

X2pα

⟩=kBT

kp(A.28)

The modes are independent of each other, so the equation becomes

〈Xpα(0)Xqβ(0)〉 = δpqδαβkBT

kp(A.29)

59

Because the spring constant changes with the mode (eqn. A.21), there areoscillating, damped and aperiodically damped modes to be expected. A solutionof a harmonic damped oscillator can be found in [14]. With the assumptionthat all random-forces of the future are uncorrelated with the current state(eqn. A.15) and ∂

⟨~X2p

⟩/∂t = 0 because the mean is constant with time, it is

now possible to find the evolution of 〈Xpα(t)Xqβ(0)〉. Three cases have to beexamined:

a) Weak damping: oscillation

This is the case for ζ2p < 4kpmp. Defining

ωp =

√kpmp−

ζ2p

4m2p

, (A.30)

the result becomes:

〈Xpα(t)Xqβ(0)〉 = δpqδαβkBT

kpe− ζp

2mpt(

cosωpt+ζp

2mpωpsinωpt

)(A.31)

b) Aperiodic limit

In this case equation ζ2p = 4kpmp becomes true. Thus ωp becomes zero.

The result now becomes:

〈Xpα(t)Xqβ(0)〉 = δpqδαβkBT

kpe− ζp

2mpt(

1 +ζp

2mpt

)(A.32)

c) Strong damping

The last case becomes true for ζ2p > 4kpmp. Here ωp becomes imaginary

and no oscillation is possible. Defining:

γp =

√ζ2p

4m2p

− kpmp

(A.33)

a1 =12

(1 +

ζp2mpγp

)(A.34)

a2 =12

(1− ζp

2mpγp

)(A.35)

the result can be written as:

〈Xpα(t)Xqβ(0)〉 = δpqδαβkBT

kpe− ζp

2mpt (a1e

γpt + a2e−γpt

)(A.36)

60

Appendix B

All mobilities obtained fromsimulation

In this appendix, all mobilities obtained from simulation are plotted as a func-tion of the number of monomers in the chain for various tilts and forces.

61

1010

010

00N

umbe

r of m

onom

ers

0

0.2

0.4

0.6

0.8

Mobility

Ex =

0.0

4E

x = 0

.02

Ex =

0.0

1E

x = 0

.005

Ex =

0.0

025M

obili

ty fo

r tw

all =

0Sc

an o

ver c

hain

leng

th

1010

010

00N

umbe

r of m

onom

ers

0

0.2

0.4

0.6

0.8

Mobility

Ex =

0.0

4E

x = 0

.02

Ex =

0.0

1E

x = 0

.005

Ex =

0.0

025

Mob

ility

for t

wal

l = 0

.01

Scan

ove

r cha

in le

ngth

1010

010

00N

umbe

r of m

onom

ers

0

0.2

0.4

0.6

0.8

Mobility

Ex =

0.0

4E

x = 0

.02

Ex =

0.0

1E

x = 0

.005

Ex =

0.0

025

Mob

ility

for t

wal

l = 0

.02

Scan

ove

r cha

in le

ngth

1010

010

00N

umbe

r of m

onom

ers

0

0.2

0.4

0.6

0.8

Mobility

Ex =

0.0

4E

x = 0

.02

Ex =

0.0

1E

x = 0

.005

Ex =

0.0

025

Mob

ility

for t

wal

l = 0

.05

Scan

ove

r cha

in le

ngth

1010

010

00N

umbe

r of m

onom

ers

0

0.2

0.4

0.6

0.8

Mobility

Ex =

0.0

4E

x = 0

.02

Ex =

0.0

1E

x = 0

.005

Ex =

0.0

025M

obili

ty fo

r tw

all =

0.1

Scan

ove

r cha

in le

ngth

1010

010

00N

umbe

r of m

onom

ers

0

0.2

0.4

0.6

0.81

Mobility

Ex =

0.0

4E

x = 0

.02

Ex =

0.0

1E

x = 0

.005

Ex =

0.0

025M

obili

ty fo

r tw

all =

0.1

5Sc

an o

ver c

hain

leng

th

1010

010

00N

umbe

r of m

onom

ers

0

0.2

0.4

0.6

0.81

Mobility

Ex =

0.0

4E

x = 0

.02

Ex =

0.0

1E

x = 0

.005

Ex =

0.0

025M

obili

ty fo

r tw

all =

0.2

Scan

ove

r cha

in le

ngth

Appendix C

Source codes of thesimulation program

On the following pages, the source code of the simulation program and all ana-lysis programs used for this thesis can be found. A sample constants file and asample start configuration are also included.

69

Te

mp

era

ture

1

Tim

e−

ste

p 1

e−

2N

um

be

r_o

f_st

ep

s 1

0N

um

be

r_o

f_o

utp

uts

1

0e

psi

lon

_L

J 1

sig

ma

_L

J 1

ep

silo

n_

BL

1

0d

_B

L 1

d_

0_

BL

1

ep

silo

n_

BA

5

zeta

_F

R 1

ep

silo

n_

BM

1

wh

ich

_B

M 1

take

_n

ew

_se

ed

0

ep

silo

n_

SP

1

00

thro

w_

aw

ay

1

00

00

0e

psi

lon

_L

J_tr

ap

1

sig

ma

_L

J_tr

ap

1

l_cu

t_L

J_tr

ap

_*_

sig

ma

_L

J_tr

ap

1

.12

24

62

04

83

h_

top

_tr

ap

0

h_

bo

tto

m_

tra

p 8

0l_

bo

tto

m_

tra

p 8

0l_

top

_tr

ap

2

0w

all_

tilt_

tra

p 0

.1e

psi

lon

_e

l_x

−

0.0

4e

psi

lon

_e

l_y

0

ep

silo

n_

el_

z 0

en

ab

le_

LJ

0

en

ab

le_

BL

0

en

ab

le_

BA

0

en

ab

le_

FR

1

en

ab

le_

BM

1

en

ab

le_

SP

0

en

ab

le_

tra

p 1

en

ab

le_

el

1

dis

pla

y_co

nst

an

ts 1

con

stan

ts#include

<io

stre

am

>#include

<fs

tre

am

>#include

<st

rin

g>

#include

<cm

ath

>#include

<st

dio

.h>

#ifndef

CO

NS

TA

NT

S_

INC

LU

DE

D#define

CO

NS

TA

NT

S_

INC

LU

DE

D

// T

his

file

re

ad

s th

e c

on

sta

nts

fo

r th

e m

ain

pro

gra

m a

nd

a fu

nct

ion

fo

r d

isp

layi

ng

// th

em

, if

de

sire

d.

// T

he

se a

re th

e c

on

sta

nts

fo

r th

e L

en

na

rd−

Jon

es

po

ten

tial b

etw

ee

n tw

o m

on

om

ers

do

ub

le e

psi

lon

_L

J;d

ou

ble

sig

ma

_L

J;d

ou

ble

v_

cuto

ff_

LJ;

do

ub

le l_

cuto

ff_

LJ;

do

ub

le l_

cuto

ff_

LJ_

sqr;

do

ub

le s

igm

a_

6;

do

ub

le p

refa

cto

r_L

J;

// T

he

se a

re th

e c

on

sta

nts

fo

r th

e b

on

d−

len

gth

−e

ne

rgy

be

twe

en

tw

o m

on

om

ers

// o

n th

e c

ha

ind

ou

ble

ep

silo

n_

BL

;d

ou

ble

d_

BL

;d

ou

ble

d_

BL

_sq

r;d

ou

ble

d_

0_

BL

;d

ou

ble

pre

fact

or_

BL

;

// T

his

is th

e c

on

sta

nt fo

r th

e b

on

d−

an

gle

−e

ne

rgy

do

ub

le e

psi

lon

_B

A;

do

ub

le p

refa

cto

r_B

A;

// T

his

is th

e frict

ion

−co

nst

an

td

ou

ble

ze

ta_

FR

;d

ou

ble

de

cay_

v;

// T

his

is th

e te

mp

era

ture

do

ub

le T

;d

ou

ble

sq

rt_

T;

// T

his

is fo

r th

e b

row

nia

n m

otio

nd

ou

ble

ep

silo

n_

BM

;d

ou

ble

pre

fact

or_

BM

;u

nsi

gn

ed

lo

ng

ra

nd

_se

ed

;in

t w

hic

h_

BM

;

// T

his

is th

e n

um

be

r o

f st

ep

sin

t s

tep

s;// T

his

is th

e n

um

be

r o

f d

isp

laye

d c

on

figu

ratio

ns

int

nu

mb

er_

of_

ou

tpu

ts;

// T

his

is th

e n

um

be

r o

f th

wro

w−

aw

ay−

ste

ps

for

initi

alis

atio

nin

t th

row

_a

wa

y;

// T

his

is th

e tim

e−

ste

pd

ou

ble

de

lta_

t;d

ou

ble

sq

rt_

de

lta_

t;

// T

his

is th

e s

prin

g−

con

sta

nt

do

ub

le e

psi

lon

_S

P;

do

ub

le p

refa

cto

r_S

P;

//T

his

is th

e e

ntr

op

ic−

tra

p L

en

na

rd−

Jon

es−

po

ten

tial

do

ub

le e

psi

lon

_L

J_tr

ap

;d

ou

ble

sig

ma

_L

J_tr

ap

;d

ou

ble

sig

ma

_6

_tr

ap

;d

ou

ble

l_cu

toff_

LJ_

tra

p;

do

ub

le l_

cuto

ff_

LJ_

tra

p_

sqr;

do

ub

le v

_cu

toff_

LJ_

tra

p;

do

ub

le h

_to

p_

tra

p;

do

ub

le h

_b

otto

m_

tra

p;

do

ub

le l_

bo

tto

m_

tra

p;

do

ub

le l_

top

_tr

ap

;d

ou

ble

wa

ll_til

t_tr

ap

;d

ou

ble

l_su

rfa

ce_

tra

p;

do

ub

le p

refa

cto

r_L

J_tr

ap

;

// T

his

is th

e e

xte

rna

l ele

ctric

field

do

ub

le e

psi

lon

_e

l_x;

do

ub

le e

psi

lon

_e

l_y;

do

ub

le e

psi

lon

_e

l_z;

do

ub

le p

refa

cto

r_e

l_x;

do

ub

le p

refa

cto

r_e

l_y;

do

ub

le p

refa

cto

r_e

l_z;

// T

he

se a

re th

e e

n−

/dis

ab

lers

fo

r th

e p

ote

ntia

ls

con

stan

ts.c

c

1/17

cons

tant

s, c

onst

ants

.cc

int

en

ab

le_

LJ;

int

en

ab

le_

BL

;in

t e

na

ble

_B

A;

int

en

ab

le_

FR

;in

t e

na

ble

_B

M;

int

en

ab

le_

SP

;in

t e

na

ble

_tr

ap

;in

t e

na

ble

_e

l;in

t d

isp

lay_

con

sta

nts

;in

t s

ave

_n

ew

_se

ed

;in

t a

pp

en

d_

con

figs;

// −

−−

−−

−−

−−

−−

−−

−−

−−

−−

−−

−−

−−

−−

−−

−−

−−

−−

−−

−−

−−

−−

−−

−−

−−

−−

−−

−−

−−

−−

−−

−−

−−

−−

−−

−−

−−

−−

−vo

id in

it_co

nst

an

ts(

cha

r*

file

na

me

) {

F

ILE

* in

pu

t =

fo

pe

n(f

ilen

am

e, "

r");

if

(!in

pu

t) {

ce

rr <

< "

Inva

lid c

onst

ants

file

−na

me!

\n"; e

xit(

1);

}

// In

itia

lisa

tion

of te

mp

era

ture

fs

can

f(in

pu

t, "

%*s

%lf"

, &

T);

sq

rt_

T =

std

::sq

rt(T

);// In

itia

lisa

tion

of d

elta

_t

fs

can

f(in

pu

t, "

%*s

%lf"

, &

de

lta_

t);

sq

rt_

de

lta_

t =

sq

rt(d

elta

_t)

;// In

itia

lisa

tion

of st

ep

s +

ou

tpu

ts fs

can

f(in

pu

t, "

%*s

%i",

&st

ep

s);

fs

can

f(in

pu

t, "

%*s

%i",

&n

um

be

r_o

f_o

utp

uts

);

// In

itia

lisa

tion

of L

en

na

rd−

Jon

es

fs

can

f(in

pu

t, "

%*s

%lf"

, &

ep

silo

n_

LJ)

; fs

can

f(in

pu

t, "

%*s

%lf"

, &

sig

ma

_L

J);

v_

cuto

ff_

LJ

= e

psi

lon

_L

J *

0.2

5;

l_

cuto

ff_

LJ

= si

gm

a_

LJ

* 1

.12

24

62

04

83

09

3;

l_

cuto

ff_

LJ_

sqr

= l_

cuto

ff_

LJ

* l_

cuto

ff_

LJ;

si

gm

a_

6 =

po

w(s

igm

a_

LJ,

6);

p

refa

cto

r_L

J =

ep

silo

n_

LJ

* si

gm

a_

6 *

6 *

de

lta_

t / 2

;// In

itia

lisa

tion

of b

on

d−

len

gth

fs

can

f(in

pu

t, "

%*s

%lf"

, &

ep

silo

n_

BL

); fs

can

f(in

pu

t, "

%*s

%lf"

, &

d_

BL

); fs

can

f(in

pu

t, "

%*s

%lf"

, &

d_

0_

BL

); d

_B

L_

sqr

= d

_B

L *

d_

BL

; p

refa

cto

r_B

L =

ep

silo

n_

BL

* d

elta

_t / 2

;// In

itia

lisa

tion

of b

on

d−

an

gle

fs

can

f(in

pu

t, "

%*s

%lf"

, &

ep

silo

n_

BA

); p

refa

cto

r_B

A =

ep

silo

n_

BA

* d

elta

_t / 2

;// In

itia

lisa

tion

of fr

ictio

n fs

can

f(in

pu

t, "

%*s

%lf"

, &

zeta

_F

R);

d

eca

y_v

= (

1−

de

lta_

t *

zeta

_F

R);

// In

itia

lisa

tion

of b

row

nia

n m

otio

n fs

can

f(in

pu

t, "

%*s

%lf"

, &

ep

silo

n_

BM

);

// e

psi

lom

_B

M e

qu

als

th

e s

tep

ha

n−

bo

ltzm

an

n−

con

sta

nt

fs

can

f(in

pu

t, "

%*s

%i",

&w

hic

h_

BM

);

p

refa

cto

r_B

M =

std

::sq

rt(2

4 *

ep

silo

n_

BM

* T

* (

1−

de

lta_

t*ze

ta_

FR

/2)

* ze

ta_

FR

* d

elta

_t)

;

if

(w

hic

h_

BM

)

// w

hic

h_

BM

==

1 e

na

ble

s ra

nd

om

nu

mb

ers

pic

ked

ou

t p

refa

cto

r_B

M *

= s

td::sq

rt(5

.0/3

.0);

// o

f th

e u

nit−

sph

ere

fs

can

f(in

pu

t, "

%*s

%i",

&sa

ve_

ne

w_

see

d);

F

ILE

* rs

ee

d =

fo

pe

n("

/hom

e/m

artin

/.ran

dom

_see

d",

"r"

);

if

(!r

see

d)

{ ce

rr <

< "

Ran

dom

−se

ed n

ot fo

und!

\n"; e

xit(

1);

}

fs

can

f(rs

ee

d, "

%lu

", &

ran

d_

see

d);

fc

lose

(rse

ed

);// In

itia

lisa

tion

of sp

rin

g fs

can

f(in

pu

t, "

%*s

%lf"

, &

ep

silo

n_

SP

); fs

can

f(in

pu

t, "

%*s

%i",

&th

row

_a

wa

y);

p

refa

cto

r_S

P =

ep

silo

n_

SP

* d

elta

_t / 2

;// In

itia

lisa

tion

of e

ntr

op

ic tra

p L

en

na

rd−

Jon

es

fs

can

f(in

pu

t, "

%*s

%lf"

, &

ep

silo

n_

LJ_

tra

p);

fs

can

f(in

pu

t, "

%*s

%lf"

, &

sig

ma

_L

J_tr

ap

); fs

can

f(in

pu

t, "

%*s

%lf"

, &

l_cu

toff_

LJ_

tra

p);

l_

cuto

ff_

LJ_

tra

p *

= s

igm

a_

LJ_

tra

p;

fs

can

f(in

pu

t, "

%*s

%lf"

, &

h_

top

_tr

ap

); fs

can

f(in

pu

t, "

%*s

%lf"

, &

h_

bo

tto

m_

tra

p);

fs

can

f(in

pu

t, "

%*s

%lf"

, &

l_b

otto

m_

tra

p);

fs

can

f(in

pu

t, "

%*s

%lf"

, &

l_to

p_

tra

p);

fs

can

f(in

pu

t, "

%*s

%lf"

, &

wa

ll_til

t_tr

ap

);

con

stan

ts.c

c si

gm

a_

6_

tra

p =

po

w(s

igm

a_

LJ_

tra

p, 6

); v_

cuto

ff_

LJ_

tra

p =

ep

silo

n_

LJ_

tra

p *

(p

ow

(sig

ma

_L

J_tr

ap

/l_cu

toff_

LJ_

tra

p, 6

) −

po

w(s

igm

a_

LJ_

tra

p/l_

cuto

ff_

LJ_

tra

p, 1

2))

; l_

cuto

ff_

LJ_

tra

p_

sqr

= p

ow

(l_

cuto

ff_

LJ_

tra

p, 2

);

l_

surf

ace

_tr

ap

= l_

bo

tto

m_

tra

p +

l_to

p_

tra

p;

p

refa

cto

r_L

J_tr

ap

= e

psi

lon

_L

J_tr

ap

* s

igm

a_

6_

tra

p *

6 *

de

lta_

t / 2

;// In

itia

lisa

tion

of e

xte

rna

l ele

ctric

field

fs

can

f(in

pu

t, "

%*s

%lf"

, &

ep

silo

n_

el_

x);

fs

can

f(in

pu

t, "

%*s

%lf"

, &

ep

silo

n_

el_

y);

fs

can

f(in

pu

t, "

%*s

%lf"

, &

ep

silo

n_

el_

z);

p

refa

cto

r_e

l_x

= e

psi

lon

_e

l_x

* d

elta

_t / 2

; p

refa

cto

r_e

l_y

= e

psi

lon

_e

l_y

* d

elta

_t / 2

; p

refa

cto

r_e

l_z

= e

psi

lon

_e

l_z

* d

elta

_t / 2

;// In

itia

lisa

tion

of e

n−

/dis

ab

lers

fs

can

f(in

pu

t, "

%*s

%i",

&e

na

ble

_L

J);

fs

can

f(in

pu

t, "

%*s

%i",

&e

na

ble

_B

L);

fs

can

f(in

pu

t, "

%*s

%i",

&e

na

ble

_B

A);

fs

can

f(in

pu

t, "

%*s

%i",

&e

na

ble

_F

R);

fs

can

f(in

pu

t, "

%*s

%i",

&e

na

ble

_B

M);

fs

can

f(in

pu

t, "

%*s

%i",

&e

na

ble

_S

P);

fs

can

f(in

pu

t, "

%*s

%i",

&e

na

ble

_tr

ap

); fs

can

f(in

pu

t, "

%*s

%i",

&e

na

ble

_e

l);

fs

can

f(in

pu

t, "

%*s

%i",

&d

isp

lay_

con

sta

nts

); fc

lose

(in

pu

t);

} void

dis

pla

yco

nst

an

ts (

cha

r*

std

ou

t_fil

en

am

e)

{

FIL

E*

std

ou

t_fil

e =

fo

pe

n(s

tdo

ut_

file

na

me

, "

a");

fp

rin

tf(s

tdo

ut_

file

, "

# T

= %

lf\n", T

); fp

rin

tf(s

tdo

ut_

file

, "

# sq

rt_T

=

%lf\

n", s

qrt

_T

); fp

rin

tf(s

tdo

ut_

file

, "

# de

lta_t

= %

lf\n", d

elta

_t)

; fp

rin

tf(s

tdo

ut_

file

, "

# st

eps

=

%i\n",

ste

ps)

; fp

rin

tf(s

tdo

ut_

file

, "

# ou

tput

s

=

%i\n",

nu

mb

er_

of_

ou

tpu

ts);

fp

rin

tf(s

tdo

ut_

file

, "

# ep

silo

n_LJ

=

%lf\

n", e

psi

lon

_L

J);

fp

rin

tf(s

tdo

ut_

file

, "

# si

gma_

LJ

=

%lf\

n", s

igm

a_

LJ)

; fp

rin

tf(s

tdo

ut_

file

, "

# v_

cuto

ff_LJ

=

%lf\

n", v

_cu

toff_

LJ)

; fp

rin

tf(s

tdo

ut_

file

, "

# l_

cuto

ff_LJ

=

%lf\

n", l_

cuto

ff_

LJ)

; fp

rin

tf(s

tdo

ut_

file

, "

# l_

cuto

ff_LJ

_sqr

= %

lf\n", l_

cuto

ff_

LJ_

sqr)

; fp

rin

tf(s

tdo

ut_

file

, "

# si

gma_

6

=

%lf\

n", s

igm

a_

6);

fp

rin

tf(s

tdo

ut_

file

, "

# pr

efac

tor_

LJ

= %

lf\n", p

refa

cto

r_L

J);

fp

rin

tf(s

tdo

ut_

file

, "

# ep

silo

n_B

L

= %

lf\n", e

psi

lon

_B

L);

fp

rin

tf(s

tdo

ut_

file

, "

# d_

BL

=

%lf\

n", d

_B

L);

fp

rin

tf(s

tdo

ut_

file

, "

# d_

0_B

L

=

%lf\

n", d

_0

_B

L);

fp

rin

tf(s

tdo

ut_

file

, "

# d_

BL_

sqr

= %

lf\n", d

_B

L_

sqr)

; fp

rin

tf(s

tdo

ut_

file

, "

# pr

efac

tor_

BL

= %

lf\n", p

refa

cto

r_B

L);

fp

rin

tf(s

tdo

ut_

file

, "

# ep

silo

n_B

A

= %

lf\n", e

psi

lon

_B

A);

fp

rin

tf(s

tdo

ut_

file

, "

# pr

efac

tor_

BA

=

%lf\

n", p

refa

cto

r_B

A);

fp

rin

tf(s

tdo

ut_

file

, "

# ze

ta_F

R

=

%lf\

n", z

eta

_F

R);

fp

rin

tf(s

tdo

ut_

file

, "

# de

cay_

v

=

%lf\

n", d

eca

y_v)

; fp

rin

tf(s

tdo

ut_

file

, "

# ep

silo

n_S

P

= %

lf\n", e

psi

lon

_S

P);

fp

rin

tf(s

tdo

ut_

file

, "

# pr

efac

tor_

SP

=

%lf\

n", p

refa

cto

r_S

P);

fp

rin

tf(s

tdo

ut_

file

, "

# th

row

_aw

ay

= %

i\n", th

row

_a

wa

y);

fp

rin

tf(s

tdo

ut_

file

, "

# ep

silo

n_B

M

= %

lf\n", e

psi

lon

_B

M);

fp

rin

tf(s

tdo

ut_

file

, "

# pr

efac

tor_

BM

=

%lf\

n", p

refa

cto

r_B

M);

fp

rin

tf(s

tdo

ut_

file

, "

# w

hich

_BM

= %

i\n", w

hic

h_

BM

); fp

rin

tf(s

tdo

ut_

file

, "

# ra

nd−

seed

= %

u\n

", r

an

d_

see

d);

fp

rin

tf(s

tdo

ut_

file

, "

# ep

silo

n_LJ

_tra

p

= %

lf\n", e

psi

lon

_L

J_tr

ap

); fp

rin

tf(s

tdo

ut_

file

, "

# si

gma_

LJ_t

rap

= %

lf\n", s

igm

a_

LJ_

tra

p);

fp

rin

tf(s

tdo

ut_

file

, "

# si

gma_

6_tr

ap

= %

lf\n", s

igm

a_

6_

tra

p);

fp

rin

tf(s

tdo

ut_

file

, "

# l_

cuto

ff_LJ

_tra

p

= %

lf\n", l_

cuto

ff_

LJ_

tra

p);

fp

rin

tf(s

tdo

ut_

file

, "

# l_

cuto

ff_LJ

_tra

p_sq

r =

%lf\

n", l_

cuto

ff_

LJ_

tra

p_

sqr)

; fp

rin

tf(s

tdo

ut_

file

, "

# v_

cuto

ff_LJ

_tra

p

= %

lf\n", v

_cu

toff_

LJ_

tra

p);

fp

rin

tf(s

tdo

ut_

file

, "

# h_

top_

trap

=

%lf\

n", h

_to

p_

tra

p);

fp

rin

tf(s

tdo

ut_

file

, "

# h_

botto

m_t

rap

= %

lf\n", h

_b

otto

m_

tra

p);

fp

rin

tf(s

tdo

ut_

file

, "

# l_

botto

m_t

rap

= %

lf\n", l_

bo

tto

m_

tra

p);

fp

rin

tf(s

tdo

ut_

file

, "

# l_

top_

trap

=

%lf\

n", l_

top

_tr

ap

); fp

rin

tf(s

tdo

ut_

file

, "

# w

all_

tilt_

trap

= %

lf\n",

wa

ll_til

t_tr

ap

); fp

rin

tf(s

tdo

ut_

file

, "

# l_

surf

ace_

trap

= %

lf\n", l_

surf

ace

_tr

ap

); fp

rin

tf(s

tdo

ut_

file

, "

# pr

efac

tor_

LJ_t

rap

=

%lf\

n", p

refa

cto

r_L

J_tr

ap

); fp

rin

tf(s

tdo

ut_

file

, "

# ep

silo

n_el

_x

= %

lf\n", e

psi

lon

_e

l_x)

; fp

rin

tf(s

tdo

ut_

file

, "

# ep

silo

n_el

_y

= %

lf\n", e

psi

lon

_e

l_y)

; fp

rin

tf(s

tdo

ut_

file

, "

# ep

silo

n_el

_z

= %

lf\n", e

psi

lon

_e

l_z)

; fp

rin

tf(s

tdo

ut_

file

, "

# pr

efac

tor_

el_x

= %

lf\n", p

refa

cto

r_e

l_x)

; fp

rin

tf(s

tdo

ut_

file

, "

# pr

efac

tor_

el_y

= %

lf\n", p

refa

cto

r_e

l_y)

; fp

rin

tf(s

tdo

ut_

file

, "

# pr

efac

tor_

el_z

= %

lf\n", p

refa

cto

r_e

l_z)

; fp

rin

tf(s

tdo

ut_

file

, "

# en

able

_LJ

=

%i\n",

en

ab

le_

LJ)

; fp

rin

tf(s

tdo

ut_

file

, "

# en

able

_BL

=

%i\n",

en

ab

le_

BL

); fp

rin

tf(s

tdo

ut_

file

, "

# en

able

_BA

= %

i\n", e

na

ble

_B

A);

fp

rin

tf(s

tdo

ut_

file

, "

# en

able

_FR

= %

i\n", e

na

ble

_F

R);

fp

rin

tf(s

tdo

ut_

file

, "

# en

able

_BM

= %

i\n", e

na

ble

_B

M);

fp

rin

tf(s

tdo

ut_

file

, "

# en

able

_SP

= %

i\n", e

na

ble

_S

P);

fp

rin

tf(s

tdo

ut_

file

, "

# en

able

_tra

p

= %

i\n", e

na

ble

_tr

ap

); fp

rin

tf(s

tdo

ut_

file

, "

# en

able

_el

=

%i\n",

en

ab

le_

el);

fp

rin

tf(s

tdo

ut_

file

, "

# ap

pend

_con

figs

= %

i\n", a

pp

en

d_

con

figs)

;

con

stan

ts.c

c

2/17

cons

tant

s.cc

fclose(stdout_file);

} #endif

con

stan

ts.c

c#include<iostream>

#include<fstream>

#include<string>

#include<cmath>

#include<stdio.h>

#ifndef CONSTANTS_INCLUDED

#define CONSTANTS_INCLUDED

// T

his

file

re

ad

s th

e c

on

sta

nts

fo

r th

e m

ain

pro

gra

m a

nd

a fu

nct

ion

fo

r d

isp

layi

ng

// th

em

, if

de

sire

d.

// T

he

se a

re th

e c

on

sta

nts

fo

r th

e L

en

na

rd−

Jon

es

po

ten

tial b

etw

ee

n tw

o m

on

om

ers

extern

do

ub

le epsilon_LJ;

extern

do

ub

le sigma_LJ;

extern

do

ub

le v_cutoff_LJ;

extern

do

ub

le l_cutoff_LJ;

extern

do

ub

le l_cutoff_LJ_sqr;

extern

do

ub

le sigma_6;

extern

do

ub

le prefactor_LJ;

// T

he

se a

re th

e c

on

sta

nts

fo

r th

e b

on

d−

len

gth

−e

ne

rgy

be

twe

en

tw

o m

on

om

ers

// o

n th

e c

ha

inextern

do

ub

le epsilon_BL;

extern

do

ub

le d_BL;

extern

do

ub

le d_BL_sqr;

extern

do

ub

le d_0_BL;

extern

do

ub

le prefactor_BL;

// T

his

is th

e c

on

sta

nt fo

r th

e b

on

d−

an

gle

−e

ne

rgy

extern

do

ub

le epsilon_BA;

extern

do

ub

le prefactor_BA;

// T

his

is th

e frict

ion

−co

nst

an

textern

do

ub

le zeta_FR;

extern

do

ub

le decay_v;

// T

his

is th

e te

mp

era

ture

extern

do

ub

le T;

extern

do

ub

le sqrt_T;

// T

his

is fo

r th

e b

row

nia

n m

otio

nextern

do

ub

le epsilon_BM;

extern

do

ub

le prefactor_BM;

extern

un

sig

ne

d

lon

g rand_seed;

extern

int which_BM;

// T

his

is th

e n

um

be

r o

f st

ep

sextern

int steps;

// T

his

is th

e n

um

be

r o

f d

isp

laye

d c

on

figu

ratio

ns

extern

int number_of_outputs;

// T

his

is th

e n

um

be

r o

f th

wro

w−

aw

ay−

ste

ps

for

initi

alis

atio

nextern

int throw_away;

// T

his

is th

e tim

e−

ste

pextern

do

ub

le delta_t;

extern

do

ub

le sqrt_delta_t;

// T

his

is th

e s

prin

g−

con

sta

nt

extern

do

ub

le epsilon_SP;

extern

do

ub

le prefactor_SP;

//T

his

is th

e e

ntr

op

ic−

tra

p L

en

na

rd−

Jon

es−

po

ten

tial

extern

do

ub

le epsilon_LJ_trap;

extern

do

ub

le sigma_LJ_trap;

extern

do

ub

le sigma_6_trap;

extern

do

ub

le l_cutoff_LJ_trap;

extern

do

ub

le l_cutoff_LJ_trap_sqr;

extern

do

ub

le v_cutoff_LJ_trap;

extern

do

ub

le h_top_trap;

extern

do

ub

le h_bottom_trap;

extern

do

ub

le l_bottom_trap;

extern

do

ub

le l_top_trap;

extern

do

ub

le wall_tilt_trap;

extern

do

ub

le l_surface_trap;

extern

do

ub

le prefactor_LJ_trap;

// T

his

is th

e e

xte

rna

l ele

ctric

field

extern

do

ub

le epsilon_el_x;

extern

do

ub

le epsilon_el_y;

extern

do

ub

le epsilon_el_z;

extern

do

ub

le prefactor_el_x;

extern

do

ub

le prefactor_el_y;

extern

do

ub

le prefactor_el_z;

// T

he

se a

re th

e e

n−

/dis

ab

lers

fo

r th

e p

ote

ntia

ls

con

stan

ts.h

h

3/17

cons

tant

s.cc

, con

stan

ts.h

h

extern

in

t e

na

ble

_L

J;extern

in

t e

na

ble

_B

L;

extern

in

t e

na

ble

_B

A;

extern

in

t e

na

ble

_F

R;

extern

in

t e

na

ble

_B

M;

extern

in

t e

na

ble

_S

P;

extern

in

t e

na

ble

_tr

ap

;extern

in

t e

na

ble

_e

l;extern

in

t d

isp

lay_

con

sta

nts

;extern

in

t s

ave

_n

ew

_se

ed

;extern

in

t a

pp

en

d_

con

figs;

// −

−−

−−

−−

−−

−−

−−

−−

−−

−−

−−

−−

−−

−−

−−

−−

−−

−−

−−

−−

−−

−−

−−

−−

−−

−−

−−

−−

−−

−−

−−

−−

−−

−−

−−

−−

−−

−−

−vo

id in

it_co

nst

an

ts(

cha

r*

file

na

me

);

void

dis

pla

yco

nst

an

ts (

cha

r*

std

ou

t_fil

en

am

e);

#endif

con

stan

ts.h

h#include

<cm

ath

>#include

<io

stre

am

>#include

"po

int.h

h"

#ifndef

IP

OIN

T_

INC

LU

DE

D#define

IP

OIN

T_

INC

LU

DE

D

// T

his

is a

cla

ss o

f 3

−d

ime

nsi

on

al i

nte

ge

rs. It in

clu

de

s a

co

nve

rsio

n−

op

era

tor

// ip

oin

t(d

ou

ble

, d

ou

ble

, d

ou

ble

) a

nd

als

o ip

oin

t(p

oin

t) fo

r e

asy

co

nve

rsio

n o

f// 3

−d

ime

nsi

on

al p

oin

ts. It a

lso

incl

ud

es

ord

er−

op

era

tors

(<

, >

, <

=, >

=)

wh

ich

first

// s

ort

th

e ip

oin

ts b

y th

e first

co

mp

on

en

t, th

en

by

the

se

con

d c

om

po

ne

nt a

nd

// a

t la

st b

y th

e th

ird

en

try.

Th

is w

as

use

d fo

r so

rtin

g th

e ip

oin

ts in

th

e// in

de

x−tr

ee

. class

ipo

int {

in

t x

; in

t y

; in

t z

;public

: ip

oin

t()

{}

ip

oin

t(in

t x

x,

int

yy,

in

t z

z){

x

= x

x; y

= y

y; z

= z

z; }

ip

oin

t(d

ou

ble

xx,

d

ou

ble

yy,

d

ou

ble

zz)

{ x

=

int

(std

::flo

or(

xx))

; y

=

int

(std

::flo

or(

yy))

; z

=

int

(std

::flo

or(

zz))

; }

ip

oin

t(p

oin

t p

){ x

=

int

(std

::flo

or(

p.c

om

po

ne

nt(

1))

);

y

=

int

(std

::flo

or(

p.c

om

po

ne

nt(

2))

); z

=

int

(std

::flo

or(

p.c

om

po

ne

nt(

3))

); }

friend

ost

rea

m&

operator

<<

(ost

rea

m&

os,

ipo

int p

){

return

os

<<

"("

<<

p.x

<<

","

<<

p.y

<<

","

<<

p.z

<<

")"

; }

friend

b

oo

l operator

<(ip

oin

t p

1, ip

oin

t p

2){

if

((p

1.x

==

p2

.x)

&&

(p

1.y

==

p2

.y))

return

(p

1.z

< p

2.z

);

if

(p

1.x

==

p2

.x)

return

(p

1.y

< p

2.y

);

return

(p

1.x

< p

2.x

); }

friend

b

oo

l operator

>(ip

oin

t p

1, ip

oin

t p

2){

return

(p

2 <

p1

); }

friend

b

oo

l operator

<=

(ip

oin

t p

1, ip

oin

t p

2){

return

(p

1 =

= p

2)

|| (p

1 <

p2

); }

friend

b

oo

l operator

>=

(ip

oin

t p

1, ip

oin

t p

2){

return

(p

1 =

= p

2)

|| (p

1 >

p2

); }

friend

b

oo

l operator

==

(ip

oin

t p

1, ip

oin

t p

2){

return

((p

1.x

==

p2

.x)

&&

(p

1.y

==

p2

.y)

&&

(p

1.z

==

p2

.z))

; }

friend

b

oo

l operator

!=(ip

oin

t p

1, ip

oin

t p

2){

return

!(p

1=

=p

2);

}

friend

ipo

int

operator

+(ip

oin

t p

1, ip

oin

t p

2){

return

ipo

int(

p1

.x +

p2

.x, p

1.y

+ p

2.y

, p

1.z

+ p

2.z

); }

ip

oin

t operator

+=

(ip

oin

t p

){ x

+=

p.x

; y

+=

p.y

; z

+=

p.z

;

return

ipo

int(

x,y,

z);

}

friend

ipo

int

operator

−(ip

oin

t p

){

return

ipo

int(

−p

.x, −

p.y

, −

p.z

); }

friend

ipo

int

operator

−(ip

oin

t p

1, ip

oin

t p

2){

return

ipo

int(

p1

.x−

p2

.x, p

1.y

−p

2.y

, p

1.z

−p

2.z

); }

ip

oin

t operator

−=

(ip

oin

t p

){ x

−=

p.x

; y

−=

p.y

; z

−=

p.z

;

return

ipo

int(

x,y,

z);

}

in

t c

om

po

ne

nt(

int

ch

oic

e){

return

(ch

oic

e=

=1

? x

: (

cho

ice

==

2 ?

y : z

));

}

in

t x

_(

void

){

return

x;

}

in

t y

_(

void

){

return

y;

}

in

t z

_(

void

){

return

z;

}

}; #endif

ipo

int.

hh

4/17

cons

tant

s.hh

, ipo

int.h

h

#include

<io

stre

am

>#include

<h

ash

_m

ap

>#include

<st

dio

.h>

#include

"ip

oint

.hh"

#include

"po

int.h

h"#include

"po

tent

ials

.hh"

#include

"co

nsta

nts.

hh"

#ifndef

HA

SH

_IP

OIN

T_

INC

LU

DE

D#define

HA

SH

_IP

OIN

T_

INC

LU

DE

D

// T

his

is th

e im

ple

me

nta

tion

of th

e h

ash

tab

le c

on

tain

ing

th

e in

dic

es

// o

f e

ach

cu

be

. It c

on

tain

s fu

nct

ion

s fo

r in

sert

ing

/ re

mo

vin

g in

dic

es

// (

the

ap

pro

pria

te c

ub

es

are

inse

rte

d/ d

ele

ted

au

tom

atic

ally

to

sa

ve

// m

em

ory

), d

isp

layi

ng

a c

ert

ain

cu

be

plu

s its

ind

ice

s a

nd

ne

igh

bo

urs

,// c

ou

ntin

g th

e in

dic

es

in a

ce

rta

in c

ub

e a

nd

fu

nct

ion

s to

cre

ate

an

// a

rra

y co

nta

inin

g a

ll in

dic

es

with

in a

cu

be

an

d o

f a

ce

rta

in a

mo

un

t// o

f n

eig

hb

ou

rs. F

or

the

ne

are

st n

eig

hb

ou

rs a

nd

an

incl

ud

ed

ce

nte

r−// c

ub

e, a

sp

eci

aliz

ed

fu

nct

ion

fo

r re

ad

ou

t e

xist

s. It g

ain

s sp

ee

d b

y// a

cce

ssin

g th

e n

eig

hb

ou

rs d

ire

ctly

th

rou

gh

th

e r

efe

ren

ces.

// // In

ma

ny

fun

ctio

ns,

a p

oin

ter

’mo

use

’ is

use

d. It is

use

d to

po

int

// a

t a

ce

rta

in c

ub

e, lik

e a

co

mp

ute

r−m

ou

se is

use

d to

po

int o

n th

e s

cre

en

.// // S

tru

ctu

res

use

d:

// −

no

de

: c

on

tain

s th

e c

urr

en

t cu

be

, th

e n

um

be

r o

f in

dic

es

in th

e c

urr

en

t // cu

be

, a

n a

rra

y o

f in

dic

es,

th

e s

ize

of th

e a

rra

y, th

e r

efe

ren

ce// −

field

to

th

e n

ea

rest

ne

igh

bo

urs

an

d a

sta

nd

ard

−co

nst

ruct

or

//

// −

ha

sh<

ipo

int>

: th

e (

first

pa

rt)

of th

e h

ash

−fu

nct

ion

// −

eq

key

: th

e e

qu

alit

y−re

latio

n n

ece

ssa

ry fo

r h

ash

ing

// // C

on

sta

nts

use

d:

// −

ma

x_n

eig

hb

ou

rs : th

e m

axi

mu

m n

um

be

r o

f n

ea

rest

ne

igh

bo

urs

of a

cu

be

// (

eq

ua

ls 3

*3*3

, a

cu

be

with

sid

e−

len

gth

3)

//

// F

un

ctio

ns

use

d (

fun

ctio

ns

in b

rack

ets

are

no

t in

ten

de

d to

be

ca

lled

dire

ctly

// b

y th

e u

ser

an

d a

re n

ot in

th

e h

ea

de

r−fil

e):

// −

(cr

ea

te_

ne

w)

: T

his

fu

nct

ion

cre

ate

s a

ne

w n

od

e. S

tore

s th

e c

ub

e in

th

e

// n

od

e, cr

ea

tes

an

arr

ay

for

two

ind

ice

s, s

ets

th

e n

um

be

r// o

f va

lid e

ntr

ies

to z

ero

. C

rea

tes

the

re

fere

nce

−a

rra

y// a

nd

fill

s it

with

th

e v

alid

re

fere

nce

s. U

pd

ate

s th

e// re

fre

nce

s o

f th

e n

eig

hb

ou

rin

g c

ub

es

as

we

ll.// −

ha

sh_

inse

rt_

ind

ex

: In

sert

s a

n in

de

x in

th

e h

ash

−ta

ble

.// If n

ece

ssa

ry, it

cre

ate

s a

ne

w n

od

e fo

r th

e c

ub

e.

// C

he

cks,

if th

e in

de

x−a

rra

y is

big

en

ou

gh

to

sto

re// a

no

the

r in

de

x. If n

ot, th

e s

ize

of th

e in

de

x is

// d

ou

ble

d. S

tore

s th

e n

ew

ind

ex

at th

e e

nd

of th

e a

rra

y// a

nd

incr

ea

ses

the

nu

mb

er

of va

lid e

ntr

ies

by

on

e.

// −

ha

sh_

ou

tpu

t : D

isp

lays

th

e d

esi

red

cu

be

, a

nd

, if

valid

, th

e in

dic

es

with

in// th

e c

ub

e a

nd

th

e v

alid

ne

are

st−

ne

igh

bo

ur

cub

es.

// −

(d

est

roy)

: T

his

fu

nct

ion

de

stro

ys a

no

de

. D

ele

tes

the

re

fere

nce

s a

nd

th

e// b

ack

po

intin

g o

ne

s (f

rom

th

e n

eig

hb

ou

rs)

as

we

ll. D

ele

tes

the

// in

de

x−a

rra

y a

nd

at la

st th

e n

od

e it

self.

// −

ha

sh_

rem

ove

_in

de

x : R

em

ove

s a

n in

de

x fr

om

th

e h

ash

_ta

ble

. C

he

cks,

if th

e c

ub

e// is

va

lid a

nd

co

nta

ins

the

co

rre

ct in

de

x. S

wa

ps

the

ind

ex

// to

be

de

lete

d to

th

e e

nd

of th

e fie

ld. In

cre

ase

s th

e n

um

be

r// o

f va

lid e

ntr

ies

by

on

e. If n

o v

alid

en

trie

s a

re le

ft,

// th

e n

od

e is

de

stro

yed

.// −

ind

ex_

cou

nt : C

ou

nts

th

e in

dic

es

with

in a

ce

rta

in c

ub

e. R

etu

rns

zero

fo

r// a

no

n−

exi

stin

g c

ub

e.

// −

ind

ice

s_in

_cu

be

: R

etu

rns

an

arr

ay

con

tain

ing

all

valid

ind

ice

s w

ithin

cu

be

.// C

rea

tes

a n

ew

inte

ge

r−a

rra

y, s

o d

on

’t fo

rge

t to

fre

e th

e// m

em

ory

fo

r it

afte

rwa

rds.

Als

o r

etu

rns

the

nu

mb

er

of

// v

alid

ind

ice

s.// −

ne

igh

bo

urs

_o

f_cu

be

: R

ea

ds

ou

t th

e in

dic

es

in th

e c

urr

en

t cu

be

an

d it

s n

ea

rest

// n

eig

hb

ou

rs, w

hic

h a

re a

cce

sse

d d

ire

ctly

via

th

e r

efe

ren

ce−

// lis

t. C

rea

tes

an

arr

ay

con

tain

ing

th

e in

dic

es,

so

do

n’t

// fo

rge

t to

fre

e th

e m

em

ory

afte

rwa

rds.

// −

fa

st_

LJ_

up

da

te : S

can

s th

rou

gh

th

e h

ash

fo

r a

ll m

on

om

ers

, a

nd

eva

lua

tes

the

// L

en

na

rd−

Jon

es−

inte

ract

ion

po

ten

tial.

Alm

ost

sim

ilar

to

// n

eig

hb

ou

rs_

of_

cub

e a

nd

LJ_

up

da

te, b

ut a

void

es

cop

yin

g o

f // in

dic

es.

st

ruct

no

de

{ ip

oin

t cu

be

;

// s

tore

s th

e (

lattic

e−

)cu

be

in

t m

ax_

ind

ice

s_in

_cu

be

;

// th

e s

ize

of th

e in

de

x−a

rra

y

int

ind

ice

s_in

_cu

be

;

// th

e n

um

be

r o

f in

dic

es

in th

e c

ub

e

int

* in

de

x;

// lo

catio

n o

f th

e in

dic

es

with

in th

e c

ub

e (

arr

ay)

n

od

e *

* re

fere

nce

;

// s

tore

s th

e r

efe

ren

ces

to th

e n

ea

rest

ne

igh

bo

urs

in

t n

um

be

r_o

f_re

fere

nce

s;

// th

e n

um

be

r o

f o

ccu

pie

d r

efe

ren

ces

n

od

e()

{}

// D

efa

ult−

con

stru

cto

r}; st

ruct

ha

sh<

ipo

int>

{ si

ze_

t operator

()(ip

oin

t p

) const

{

ipo

int_

has

h.c

c

return

(p

.x_

()%

10

48

57

6)

+ (

p.y

_()

%1

04

85

76

)*1

04

85

76

+

(p

.z_

()%

10

48

57

6)*

10

48

57

6*1

04

85

76

; }

}; stru

ct e

qke

y {

b

oo

l operator

()(ip

oin

t p

1, ip

oin

t p

2)

const

{

return

(p

1 =

= p

2);

}}; // M

axi

mu

m n

um

be

r o

f n

ea

rest

ne

igh

bo

urs

of a

cu

be

const

in

t m

ax_

ne

igh

bo

urs

= 2

7;

typedef

ha

sh_

ma

p<

ipo

int, n

od

e*,

ha

sh<

ipo

int>

, e

qke

y> h

ash

_ip

oin

t;

no

de

* cr

ea

te_

ne

w(h

ash

_ip

oin

t &

idx_

ha

sh, ip

oin

t cu

be

) {

n

od

e*

mo

use

=

new

no

de

;

// c

rea

te th

e n

ew

no

de

m

ou

se−

>cu

be

= c

ub

e;

m

ou

se−

>m

ax_

ind

ice

s_in

_cu

be

= 2

;

// c

rea

te a

n a

rra

y fo

r 2

ind

ice

s m

ou

se−

>in

dic

es_

in_

cub

e =

0;

// a

nd

se

t th

e n

um

be

r o

f va

lid e

ntr

ies

m

ou

se−

>in

de

x =

new

in

t[m

ou

se−

>m

ax_

ind

ice

s_in

_cu

be

];

// to

ze

ro m

ou

se−

>re

fere

nce

=

new

(n

od

e*)

[ma

x_n

eig

hb

ou

rs];

m

ou

se−

>n

um

be

r_o

f_re

fere

nce

s =

0;

for

(in

t i

= 0

; i <

ma

x_n

eig

hb

ou

rs; i+

+)

m

ou

se−

>re

fere

nce

[i] =

0;

// s

can

th

e tre

e fo

r a

ll n

ea

rest

for

(in

t d

x =

−1

; d

x <

= 1

; d

x++

)

// n

eig

hb

ou

rs

for

(in

t d

y =

−1

; d

y <

= 1

; d

y++

)

for

(in

t d

z =

−1

; d

z <

= 1

; d

z++

)if

(ip

oin

t(d

x, d

y, d

z) !=

ipo

int(

0,0

,0))

{ h

ash

_ip

oin

t::it

era

tor

sea

rch

= id

x_h

ash

.fin

d(ip

oin

t(d

x, d

y, d

z) +

cu

be

);

if

(se

arc

h !=

idx_

ha

sh.e

nd

())

{ n

od

e*

sea

rch

_m

ou

se =

idx_

ha

sh[ip

oin

t(d

x, d

y, d

z) +

cu

be

]; m

ou

se−

>re

fere

nce

[mo

use

−>

nu

mb

er_

of_

refe

ren

ces]

= s

ea

rch

_m

ou

se;

m

ou

se−

>n

um

be

r_o

f_re

fere

nce

s++

;

se

arc

h_

mo

use

−>

refe

ren

ce[s

ea

rch

_m

ou

se−

>n

um

be

r_o

f_re

fere

nce

s] =

mo

use

; se

arc

h_

mo

use

−>

nu

mb

er_

of_

refe

ren

ces+

+;

}

}

return

mo

use

;} vo

id h

ash

_in

sert

_in

de

x(h

ash

_ip

oin

t& id

x_h

ash

, ip

oin

t cu

be

, n

od

e*&

bo

x_re

f,

int

ind

ex)

{ n

od

e*

mo

use

= id

x_h

ash

[cu

be

];

if

(m

ou

se =

= 0

) {

// c

ub

e is

no

t in

th

e h

ash

m

ou

se =

cre

ate

_n

ew

(id

x_h

ash

, cu

be

);

// a

nd

ha

s to

be

cre

ate

d id

x_h

ash

[cu

be

] =

mo

use

; }

if

(m

ou

se−

>in

dic

es_

in_

cub

e =

= m

ou

se−

>m

ax_

ind

ice

s_in

_cu

be

) {

m

ou

se−

>m

ax_

ind

ice

s_in

_cu

be

*=

2;

int

* n

ew

_in

de

x =

new

in

t[m

ou

se−

>m

ax_

ind

ice

s_in

_cu

be

];

// c

rea

te a

ne

w a

rra

y ca

pa

ble

for

(in

t i

= 0

; i <

mo

use

−>

ind

ice

s_in

_cu

be

; i+

+)

// o

f st

orin

g a

ll o

ld in

dic

es

n

ew

_in

de

x[i]

= m

ou

se−

>in

de

x[i];

// a

nd

th

e n

ew

on

e

{

int

* tm

p =

mo

use

−>

ind

ex;

// s

wa

p in

de

x−fie

lds

m

ou

se−

>in

de

x =

ne

w_

ind

ex;

n

ew

_in

de

x =

tm

p;

}

delete

[] n

ew

_in

de

x; }

m

ou

se−

>in

de

x[m

ou

se−

>in

dic

es_

in_

cub

e] =

ind

ex;

// e

nte

r th

e n

ew

ind

ex

into

th

e a

rra

y m

ou

se−

>in

dic

es_

in_

cub

e+

+;

b

ox_

ref =

mo

use

;} vo

id h

ash

_o

utp

ut(

ha

sh_

ipo

int&

idx_

ha

sh, ip

oin

t cu

be

, ch

ar

* st

do

ut_

file

na

me

) {

F

ILE

* st

do

ut_

file

= fo

pe

n(s

tdo

ut_

file

na

me

, "

a");

n

od

e*

mo

use

= id

x_h

ash

[cu

be

];

if

(m

ou

se =

= 0

) {

fp

rin

tf(s

tdo

ut_

file

, "

>>

(%

i,%i,%

i) <

<\n

", c

ub

e.x

_()

, cu

be

.y_

(), cu

be

.z_

());

id

x_h

ash

.era

se(c

ub

e);

}

else

{ fp

rin

tf(s

tdo

ut_

file

, "

(%i,%

i,%i)

{",

cu

be

.x_

(), cu

be

.y_

(), cu

be

.z_

());

for

(in

t i

= 0

; i <

mo

use

−>

ind

ice

s_in

_cu

be

; i+

+)

fp

rin

tf(s

tdo

ut_

file

, "

%i,"

, m

ou

se−

>in

de

x[i])

; fp

rin

tf(s

tdo

ut_

file

, "

%c}

, ", 8

); fp

rin

tf(s

tdo

ut_

file

, "

%i N

eigh

bour

s: >

> "

, m

ou

se−

>n

um

be

r_o

f_re

fere

nce

s);

for

(in

t i

= 0

; i <

mo

use

−>

nu

mb

er_

of_

refe

ren

ces;

i++

)

if

(m

ou

se−

>re

fere

nce

[i] !=

0)

fprin

tf(s

tdo

ut_

file

, "

(%i,%

i,%i)

", m

ou

se−

>re

fere

nce

[i]−

>cu

be

.x_

(),

mo

use

−>

refe

ren

ce[i]

−>

cub

e.y

_()

, m

ou

se−

>re

fere

nce

[i]−

>cu

be

.z_

());

fp

rin

tf(s

tdo

ut_

file

, "

<<

\n")

; }

fc

lose

(std

ou

t_fil

e);

} void

de

stro

y(n

od

e*&

mo

use

) {

ipo

int_

has

h.c

c

5/17

ipoi

nt_h

ash.

cc

// T

his

fu

nct

ion

act

ua

lly d

ele

tes

the

mo

use

, a

lso

de

lete

s th

e r

efe

ren

ces

// in

bo

th d

ire

ctio

ns

for

(in

t i

= 0

; i <

mo

use

−>

nu

mb

er_

of_

refe

ren

ces;

i++

) {

// s

can

all

valid

re

fere

nce

s

for

(in

t j

= 0

; j <

mo

use

−>

refe

ren

ce[i]

−>

nu

mb

er_

of_

refe

ren

ces;

j++

) // lo

cate

th

e

if

(m

ou

se−

>re

fere

nce

[i]−

>re

fere

nce

[j] =

= m

ou

se)

{

//b

ack

−p

oin

ting

re

fere

nce

mo

use

−>

refe

ren

ce[i]

−>

refe

ren

ce[j]

=

m

ou

se−

>re

fere

nce

[i]−

>re

fere

nce

[(m

ou

se−

>re

fere

nce

[i]−

>n

um

be

r_o

f_re

fere

nce

s) −

1

];

// a

nd

de

lete

it(m

ou

se−

>re

fere

nce

[i]−

>n

um

be

r_o

f_re

fere

nce

s)−

−;

}

}

delete

[] m

ou

se−

>re

fere

nce

;

// d

ele

te a

ll p

oin

ters

an

d a

rra

ys

delete

[] m

ou

se−

>in

de

x;

delete

mo

use

;} vo

id h

ash

_re

mo

ve_

ind

ex(

ha

sh_

ipo

int&

idx_

ha

sh, ip

oin

t cu

be

, n

od

e*&

bo

x_re

f,

int

ind

ex)

{ n

od

e*

mo

use

= id

x_h

ash

[cu

be

];

if

(m

ou

se =

= 0

) {

// c

ub

e h

as

to b

e w

ithin

th

e h

ash

id

x_h

ash

.era

se(c

ub

e);

ce

rr <

< "

Inva

lid D

elet

ion

− c

ube

is n

ot in

tree

!\n"; }

else

{

bo

ol

fo

un

d =

fa

lse

;

for

(in

t i

= 0

; i<

mo

use

−>

ind

ice

s_in

_cu

be

; i+

+)

// s

ea

rch

th

e r

igh

t in

de

x

if

(m

ou

se−

>in

de

x[i]

==

ind

ex)

{ fo

un

d =

tr

ue

;

int

tm

p =

mo

use

−>

ind

ex[

i];

// s

wa

p th

e in

de

x w

hic

h is

to

be

m

ou

se−

>in

de

x[i]

= m

ou

se−

>in

de

x[m

ou

se−

>in

dic

es_

in_

cub

e−

1];

// d

ele

ted

to

th

e m

ou

se−

>in

de

x[m

ou

se−

>in

dic

es_

in_

cub

e−

1] =

tm

p;

// e

nd

of th

e fie

ld i =

mo

use

−>

ind

ice

s_in

_cu

be

; }

if

(!fo

un

d)

{ ce

rr <

< "

Inva

lid D

elet

ion

− in

dex

is n

ot in

cub

e!\n";

}

else

{ m

ou

se−

>in

dic

es_

in_

cub

e−

−;

// a

nd

re

du

ce th

e n

um

be

r o

f va

lid in

dic

es

by

on

e }

if

(m

ou

se−

>in

dic

es_

in_

cub

e =

= 0

) {

// r

em

ove

th

e c

ub

e, if

ne

cess

ary

d

est

roy(

mo

use

); id

x_h

ash

.era

se(c

ub

e);

}

}

b

ox_

ref =

0;

} int

ind

ex_

cou

nt(

no

de

*& b

ox_

ref)

// C

ou

nts

th

e n

um

be

r o

f in

dic

es

with

in th

e n

od

e, re

turn

s 0

wh

en

cu

be

is n

ot in

// th

e h

ash

.{

if

(b

ox_

ref =

= 0

) {

return

0;

}

else

{

return

bo

x_re

f−>

ind

ice

s_in

_cu

be

; }

} int

* in

dic

es_

in_

cub

e(n

od

e*&

bo

x_re

f,

int

& n

um

be

r)// R

etu

rns

an

arr

ay

of th

e in

dic

es

loca

ted

in c

ub

e. D

oe

s n

ot re

turn

em

pty

// in

dic

es.

Th

e s

ize

of th

e a

rra

y is

als

o r

etu

rne

d in

nu

mb

er.

Th

is fu

nct

ion

// a

lloca

tes

the

me

mo

ry fo

r th

e a

rra

y, s

o d

on

’t fo

rge

t to

fre

e th

e m

em

ory

// a

fte

r u

sag

e.

{ if

(b

ox_

ref =

= 0

) {

n

um

be

r =

0;

return

0;

}

else

{ n

um

be

r =

bo

x_re

f−>

ind

ice

s_in

_cu

be

;

int

* id

x_lis

t =

new

in

t[n

um

be

r];

for

(in

t i

= 0

; i <

nu

mb

er;

i++

) id

x_lis

t[i]

= b

ox_

ref−

>in

de

x[i];

return

idx_

list;

}

} int

* n

eig

hb

ou

rs_

of_

cub

e(n

od

e*&

bo

x_re

f,

int

& n

um

be

r) {

// R

ea

ds

ou

t th

e in

dic

es

in th

e c

urr

en

t cu

be

an

d it

s n

ea

rest

ne

igh

bo

urs

, w

hic

h a

re

// a

cce

sse

d d

ire

ctly

via

th

e r

efe

ren

ce−

list. C

rea

tes

an

arr

ay

con

tain

ing

th

e in

dic

es,

// s

o d

on

’t fo

rge

t to

fre

e th

e m

em

ory

afte

rwa

rds.

if

(b

ox_

ref =

= 0

) {

ce

rr <

< "

Err

or: B

ox−

refe

renc

es a

re fa

ulty

!\n"; e

xit(

1);

}

n

um

be

r =

bo

x_re

f−>

ind

ice

s_in

_cu

be

;

// c

ou

nt th

e n

um

be

r o

f

for

(in

t i

= 0

; i <

bo

x_re

f−>

nu

mb

er_

of_

refe

ren

ces;

i++

) // in

dic

es

to b

e r

etu

rne

d n

um

be

r +

= b

ox_

ref−

>re

fere

nce

[i]−

>in

dic

es_

in_

cub

e;

in

t*

idx_

list =

new

in

t[n

um

be

r];

n

um

be

r =

0;

ipo

int_

has

h.c

c

for

(in

t i

= 0

; i <

bo

x_re

f−>

ind

ice

s_in

_cu

be

; i+

+)

{ id

x_lis

t[n

um

be

r] =

bo

x_re

f−>

ind

ex[

i]; n

um

be

r++

; }

for

(in

t i

= 0

; i <

bo

x_re

f−>

nu

mb

er_

of_

refe

ren

ces;

i++

)

for

(in

t j

= 0

; j <

bo

x_re

f−>

refe

ren

ce[i]

−>

ind

ice

s_in

_cu

be

; j+

+)

{ id

x_lis

t[n

um

be

r] =

bo

x_re

f−>

refe

ren

ce[i]

−>

ind

ex[

j]; n

um

be

r++

; }

return

idx_

list;

} void

fa

st_

LJ_

up

da

te(p

oin

t* r

, p

oin

t* d

v, n

od

e**

& b

ox_

ref,

int

N)

{// S

can

s th

rou

gh

th

e tre

e fo

r a

ll in

dic

es,

eva

lua

tes

the

Le

nn

ard

−Jo

ne

s−in

tera

ctio

n.

// A

lmo

st s

imila

r to

ne

igh

bo

urs

_o

f_cu

be

an

d L

J_u

pd

ate

, b

ut a

void

es

cop

yin

g o

f in

dic

es.

for

(in

t i

= 0

; i <

N; i+

+)

{ n

od

e*

mo

use

= b

ox_

ref[i];

if

(m

ou

se =

= 0

) {

// if

it is

no

t fo

un

d, e

xit th

e p

rog

ram

ce

rr <

< "

Inva

lid a

cces

s in

fast

_LJ_

upda

te −

cub

e is

not

in tr

ee!\n

"; e

xit(

1);

}

for

(in

t j

= 0

; j <

mo

use

−>

ind

ice

s_in

_cu

be

; j+

+)

{ // s

can

th

e c

ub

e in

th

e c

en

ter

int

ind

ex

= m

ou

se−

>in

de

x[j];

if

(in

de

x <

i) {

p

oin

t d

_V

_L

J =

d_

V_

LJ_

dr_

time

s_d

elta

_t_

2(r

[ind

ex]

− r

[i]);

d

v[i]

+

= d

_V

_L

J; d

v[in

de

x] +

= −

d_

V_

LJ;

}

}

for

(in

t k

= 0

; k

< m

ou

se−

>n

um

be

r_o

f_re

fere

nce

s; k

++

) {

// a

nd

all

ne

igh

bo

urs

n

od

e*

sea

rch

_m

ou

se =

mo

use

−>

refe

ren

ce[k

];

//th

rou

gh

th

e r

efe

ren

ces

for

(in

t j

= 0

; j <

se

arc

h_

mo

use

−>

ind

ice

s_in

_cu

be

; j+

+)

{in

t in

de

x =

se

arc

h_

mo

use

−>

ind

ex[

j];if

(in

de

x <

i) {

p

oin

t d

_V

_L

J =

d_

V_

LJ_

dr_

time

s_d

elta

_t_

2(r

[ind

ex]

− r

[i]);

d

v[i]

+

= d

_V

_L

J; d

v[in

de

x] +

= −

d_

V_

LJ;

}

}

}

}

} #endif

ipo

int_

has

h.c

c

6/17

ipoi

nt_h

ash.

cc

#include<iostream>

#include<hash_map>

#include<stdio.h>

#include"

ipoi

nt.h

h"#include"

poin

t.hh"

#include"

pote

ntia

ls.h

h"#include"

cons

tant

s.hh"

#ifndef HASH_IPOINT_INCLUDED

#define HASH_IPOINT_INCLUDED

// T

his

is th

e im

ple

me

nta

tion

of th

e h

ash

tab

le c

on

tain

ing

th

e in

dic

es

// o

f e

ach

cu

be

. It c

on

tain

s fu

nct

ion

s fo

r in

sert

ing

/ re

mo

vin

g in

dic

es

// (

the

ap

pro

pria

te c

ub

es

are

inse

rte

d/ d

ele

ted

au

tom

atic

ally

to

sa

ve

// m

em

ory

), d

isp

layi

ng

a c

ert

ain

cu

be

plu

s its

ind

ice

s a

nd

ne

igh

bo

urs

,// c

ou

ntin

g th

e in

dic

es

in a

ce

rta

in c

ub

e a

nd

fu

nct

ion

s to

cre

ate

an

// a

rra

y co

nta

inin

g a

ll in

dic

es

with

in a

cu

be

an

d o

f a

ce

rta

in a

mo

un

t// o

f n

eig

hb

ou

rs. F

or

the

ne

are

st n

eig

hb

ou

rs a

nd

an

incl

ud

ed

ce

nte

r−// c

ub

e, a

sp

eci

aliz

ed

fu

nct

ion

fo

r re

ad

ou

t e

xist

s. It g

ain

s sp

ee

d b

y// a

cce

ssin

g th

e n

eig

hb

ou

rs d

ire

ctly

th

rou

gh

th

e r

efe

ren

ces.

// // S

tru

ctu

res

use

d:

// −

no

de

: c

on

tain

s th

e c

urr

en

t cu

be

, th

e n

um

be

r o

f in

dic

es

in th

e c

urr

en

t // cu

be

, a

n a

rra

y o

f in

dic

es,

th

e s

ize

of th

e a

rra

y, th

e r

efe

ren

ce// −

field

to

th

e n

ea

rest

ne

igh

bo

urs

an

d a

sta

nd

ard

−co

nst

ruct

or

//

// −

ha

sh<

ipo

int>

: th

e (

first

pa

rt)

of th

e h

ash

−fu

nct

ion

// −

eq

key

: th

e e

qu

alit

y−re

latio

n n

ece

ssa

ry fo

r h

ash

ing

// // C

on

sta

nts

use

d:

// −

ma

x_n

eig

hb

ou

rs : th

e m

axi

mu

m n

um

be

r o

f n

ea

rest

ne

igh

bo

urs

of a

cu

be

// (

eq

ua

ls 3

*3*3

, a

cu

be

with

sid

e−

len

gth

3)

//

// F

un

ctio

ns

use

d (

fun

ctio

ns

in b

rack

ets

are

no

t in

ten

de

d to

be

ca

lled

dire

ctly

// b

y th

e u

ser

an

d a

re n

ot in

th

e h

ea

de

r−fil

e):

// −

(cr

ea

te_

ne

w)

: T

his

fu

nct

ion

cre

ate

s a

ne

w n

od

e. S

tore

s th

e c

ub

e in

th

e

// n

od

e, cr

ea

tes

an

arr

ay

for

two

ind

ice

s, s

ets

th

e n

um

be

r// o

f va

lid e

ntr

ies

to z

ero

. C

rea

tes

the

re

fere

nce

−a

rra

y// a

nd

fill

s it

with

th

e v

alid

re

fere

nce

s. U

pd

ate

s th

e// re

fre

nce

s o

f th

e n

eig

hb

ou

rin

g c

ub

es

as

we

ll.// −

ha

sh_

inse

rt_

ind

ex

: In

sert

s a

n in

de

x in

th

e h

ash

−ta

ble

.// If n

ece

ssa

ry, it

cre

ate

s a

ne

w n

od

e fo

r th

e c

ub

e.

// C

he

cks,

if th

e in

de

x−a

rra

y is

big

en

ou

gh

to

sto

re// a

no

the

r in

de

x. If n

ot, th

e s

ize

of th

e in

de

x is

// d

ou

ble

d. S

tore

s th

e n

ew

ind

ex

at th

e e

nd

of th

e a

rra

y// a

nd

incr

ea

ses

the

nu

mb

er

of va

lid e

ntr

ies

by

on

e.

// −

ha

sh_

ou

tpu

t : D

isp

lays

th

e d

esi

red

cu

be

, a

nd

, if

valid

, th

e in

dic

es

with

in// th

e c

ub

e a

nd

th

e v

alid

ne

are

st−

ne

igh

bo

ur

cub

es.

// −

(d

est

roy)

: T

his

fu

nct

ion

de

stro

ys a

no

de

. D

ele

tes

the

re

fere

nce

s a

nd

th

e// b

ack

po

intin

g o

ne

s (f

rom

th

e n

eig

hb

ou

rs)

as

we

ll. D

ele

tes

the

// in

de

x−a

rra

y a

nd

at la

st th

e n

od

e it

self.

// −

ha

sh_

rem

ove

_in

de

x : R

em

ove

s a

n in

de

x fr

om

th

e h

ash

_ta

ble

. C

he

cks,

if th

e c

ub

e// is

va

lid a

nd

co

nta

ins

the

co

rre

ct in

de

x. S

wa

ps

the

ind

ex

// to

be

de

lete

d to

th

e e

nd

of th

e fie

ld. In

cre

ase

s th

e n

um

be

r// o

f va

lid e

ntr

ies

by

on

e. If n

o v

alid

en

trie

s a

re le

ft,

// th

e n

od

e is

de

stro

yed

.// −

ind

ex_

cou

nt : C

ou

nts

th

e in

dic

es

with

in a

ce

rta

in c

ub

e. R

etu

rns

zero

fo

r// a

no

n−

exi

stin

g c

ub

e.

// −

ind

ice

s_in

_cu

be

: R

etu

rns

an

arr

ay

con

tain

ing

all

valid

ind

ice

s w

ithin

cu

be

.// C

rea

tes

a n

ew

inte

ge

r−a

rra

y, s

o d

on

’t fo

rge

t to

fre

e th

e// m

em

ory

fo

r it

afte

rwa

rds.

Als

o r

etu

rns

the

nu

mb

er

of

// v

alid

ind

ice

s.// −

ne

igh

bo

urs

_o

f_cu

be

: C

he

cks,

if c

ub

e is

loca

ted

in th

e h

ash

an

d if

th

e n

um

be

r// o

f d

esi

red

ne

are

st n

eig

hb

ou

rs is

se

t to

on

e (

the

de

fau

lt).

// If n

ot, a

wa

rnin

g is

writte

n o

n s

tde

rr a

nd

old

_n

eig

hb

ou

rs_

of_

cub

e

// is

ca

lled

, it

is a

ble

to

ha

nd

le th

ese

diff

icu

ltie

s. If th

e// ch

eck

is O

K, it

acc

ess

es

the

cu

rre

nt cu

be

dire

ctly

an

d

// th

e n

ea

rest

ne

igh

bo

urs

by

the

re

fere

nce

−lis

t, w

hic

h is

in// e

ach

no

de

an

d th

us

savi

ng

a lo

t o

f tim

e.

// −

fa

st_

LJ_

up

da

te : S

can

s th

rou

gh

th

e h

ash

fo

r a

ll m

on

om

ers

, a

nd

eva

lua

tes

the

Le

nn

ard

−Jo

ne

s−// in

tera

ctio

n p

ote

ntia

l. A

lmo

st s

imila

r to

ne

igh

bo

urs

_o

f_cu

be

an

d

// L

J_u

pd

ate

, b

ut a

void

es

cop

yin

g o

f in

dic

es.

st

ruct

node {

ipoint cube;

// s

tore

s th

e (

lattic

e−

)cu

be

int max_indices_in_cube;

// th

e s

ize

of th

e in

de

x−a

rra

y

int indices_in_cube;

// th

e n

um

be

r o

f in

dic

es

in th

e c

ub

e

int* index;

// lo

catio

n o

f th

e in

dic

es

with

in th

e c

ub

e (

arr

ay)

node ** reference;

// s

tore

s th

e r

efe

ren

ces

to th

e n

ea

rest

ne

igh

bo

urs

int number_of_references;

// th

e n

um

be

r o

f o

ccu

pie

d r

efe

ren

ces

node() {}

// D

efa

ult−

con

stru

cto

r};

stru

ct hash<ipoint> {

size_t operator()(ipoint p) const {

return (p.x_()%1048576) + (p.y_()%1048576)*1048576 + (p.z_()%1048576)*1048576*1048576;

ipo

int_

has

h.h

h }

};

stru

ct eqkey {

bo

ol operator()(ipoint p1, ipoint p2) const

{ return (p1 == p2); }

};

typedef hash_map<ipoint, node*, hash<ipoint>, eqkey> hash_ipoint;

void

hash_insert_index(hash_ipoint& idx_hash, ipoint cube, node*& box_ref,

int index);

void

hash_output(hash_ipoint& idx_hash, ipoint cube,

cha

r* stdout_filename);

void

hash_remove_index(hash_ipoint& idx_hash, ipoint cube, node*& box_ref,

int index);

int* neighbours_of_cube(node*& box_ref,

int& number);

void

fast_LJ_update(point* r, point* dv, node**& box_ref,

int N);

#endif

ipo

int_

has

h.h

h

7/17

ipoi

nt_h

ash.

hh

#include

<io

stre

am

>#include

<io

ma

nip

>#include

<cm

ath

>#include

<fs

tre

am

>#include

<cs

tdlib

>#include

<st

dio

.h>

#include

<st

rin

g.h

>#include

"po

int.h

h"#include

"co

nsta

nts.

hh"

const

in

t n

um

be

r_o

f_b

lock

s =

10

;

po

int ro

use

_m

od

e(p

oin

t* r

, in

t N

, in

t m

od

e)

{// c

om

pu

tes

the

ro

use

−m

od

e o

f a

giv

en

co

nfig

ura

tion

p

oin

t h

elp

= p

oin

t(0

,0,0

);

for

(in

t i

= 0

; i <

N; i+

+)

h

elp

+=

co

s(m

od

e *

M_

PI *

(i+

0.5

) / N

) *

r[i];

h

elp

/=

(N

);

return

he

lp;

} int

de

lta(

int

i,

int

j) {

return

(i=

=j) ?

1 : 0

;

// s

imp

ly th

e K

ron

eck

er−

de

lta} d

ou

ble

ra

diu

s(p

oin

t* r

, in

t n

um

be

r,

int

co

mp

on

en

t) {

// c

om

pu

tes

the

ra

diu

s o

f g

yra

tion

of a

giv

en

co

nfig

ura

tion

p

oin

t ce

nte

r_m

ass

= r

ou

se_

mo

de

(r, n

um

be

r, 0

);

do

ub

le s

um

_d

ist_

sqr

= 0

;

if

(co

mp

on

en

t =

= 0

)

for

(in

t i

= 0

; i <

nu

mb

er;

i++

) su

m_

dis

t_sq

r +

= a

bs_

sqr(

r[i]

− c

en

ter_

ma

ss);

else

for

(in

t i

= 0

; i <

nu

mb

er;

i++

) su

m_

dis

t_sq

r +

= (

r[i]

− c

en

ter_

ma

ss).

com

po

ne

nt(

com

po

ne

nt)

* (

r[i]

− c

en

ter_

ma

ss).

com

po

ne

nt(

com

po

ne

nt)

;

return

su

m_

dis

t_sq

r / n

um

be

r;} vo

id r

ea

d_

con

fig(

FIL

E*&

inp

ut, p

oin

t* r

, p

oin

t* v

, in

t n

um

be

r,

int

tim

est

ep

s,

do

ub

le&

en

erg

y)// r

ea

ds

the

co

nfig

ura

tion

fro

m a

file

{ in

t n

um

be

r2;

fs

can

f(in

pu

t, "

%*s

%i %

*s %

i", &

time

ste

ps,

&n

um

be

r2);

for

(in

t i

= 0

; i <

nu

mb

er;

i++

) {

do

ub

le x

, y,

z, vx

, vy

, vz

; fs

can

f(in

pu

t,"

%lf

%lf

%lf

%lf

%lf

%lf

", &

x, &

y, &

z, &

vx, &

vy, &

vz);

r[

i] =

po

int(

x, y

, z

);

v[

i] =

po

int(

vx, vy

, vz

); }

fs

can

f(in

pu

t, "

%*s

%lf"

,&

en

erg

y);

} int

ma

in(

int

arg

c,

cha

r*

arg

v[])

{

if

((a

rgc

!= 4

) &

& (

arg

c !=

5))

{ ce

rr <

< "

Inva

lid n

umbe

r of

arg

umen

ts −

\n"; ce

rr <

< "

Cal

ling

sequ

ence

: <pr

ogra

m n

ame>

con

stan

t−fil

e ou

tput

−fil

e st

dout

_file

[com

pone

nt]\n

"; e

xit(

1);

}

in

it_co

nst

an

ts(a

rgv[

1])

;

// r

ea

d c

on

sta

nts

in

t n

um

be

r;

cha

r c

;

do

ub

le e

ne

rgy;

in

t tim

est

ep

;

int

co

mp

on

en

t =

0;

if

(a

rgc

==

5)

co

mp

on

en

t =

ato

i(a

rgv[

4])

;

// d

isca

rd in

valid

se

ttin

gs

if

((c

om

po

ne

nt <

0)

|| (c

om

po

ne

nt >

3))

{ ce

rr <

< "

Inva

lid c

ompo

nent

!\n";

e

xit(

1);

}

F

ILE

* in

pu

t =

fo

pe

n(a

rgv[

2], "

r");

if

(!in

pu

t) {

ce

rr <

< "

Inva

lid o

utpu

t−fil

e−na

me!

\n"; e

xit(

1);

}

fs

can

f(in

pu

t, "

# %

i # %

i", &

time

ste

p, &

nu

mb

er)

; p

oin

t *r

=

new

po

int[n

um

be

r];

p

oin

t *v

=

new

po

int[n

um

be

r];

fc

lose

(in

pu

t);

// R

ese

t th

e r

ea

d−

bu

ffe

r in

pu

t =

fo

pe

n(a

rgv[

2], "

r");

re

ad

_co

nfig

(in

pu

t, r

, v,

nu

mb

er,

tim

est

ep

, e

ne

rgy)

;

do

ub

le e

nd

=

0;

d

ou

ble

su

m_

en

d =

0;

d

ou

ble

su

m_

en

d_

sqr

= 0

;

do

ub

le g

yr =

0;

d

ou

ble

su

m_

gyr

=

0;

mai

n_a

nal

yse_

len

gth

.cc

d

ou

ble

su

m_

gyr

_sq

r =

0;

for

(in

t i

= 1

; i <

= n

um

be

r_o

f_o

utp

uts

; i+

+)

{

if

(i%

(nu

mb

er_

of_

ou

tpu

ts/n

um

be

r_o

f_b

lock

s) =

= 0

) {

e

nd

= e

nd

/nu

mb

er_

of_

ou

tpu

ts*n

um

be

r_o

f_b

lock

s; su

m_

en

d +

= e

nd

; su

m_

en

d_

sqr

+=

en

d*e

nd

; e

nd

= 0

; g

yr =

gyr

/nu

mb

er_

of_

ou

tpu

ts*n

um

be

r_o

f_b

lock

s; su

m_

gyr

+=

gyr

; su

m_

gyr

_sq

r +

= g

yr*g

yr;

g

yr =

0;

}

re

ad

_co

nfig

(in

pu

t, r

, v,

nu

mb

er,

tim

est

ep

, e

ne

rgy)

;

if

(co

mp

on

en

t =

= 0

) e

nd

+=

ab

s_sq

r(r[

0] −

r[n

um

be

r−1

]);

else

e

nd

+=

(r[

0] −

r[n

um

be

r−1

]).c

om

po

ne

nt(

com

po

ne

nt)

* (

r[0

] −

r[n

um

be

r−1

]).c

om

po

ne

nt(

com

po

ne

nt)

; g

yr +

= r

ad

ius(

r, n

um

be

r, c

om

po

ne

nt)

; }

// −

−−

−−

−−

−−

−−

−−

−−

−−

−−

−−

−−

−−

−−

−−

−−

−−

−−

−−

−−

−−

−−

−−

−−

−−

−−

−−

−−

−−

−−

−−

−−

−−

−−

−−

−−

−−

−−

−−

−−

−−

−−

−−

−−

−−

−−

F

ILE

* st

do

ut_

file

= fo

pe

n(a

rgv[

3], "

a");

fp

rin

tf(s

tdo

ut_

file

, "

delta

_t =

%lf,

zet

a_F

R*d

elta

_t =

%lf,

ste

ps/s

ampl

e =

%i,

epsi

lon_

SP

= %

lf\n

", d

elta

_t, z

eta

_F

R*d

elta

_t, s

tep

s/n

um

be

r_o

f_b

lock

s, e

psi

lon

_S

P);

fp

rin

tf(s

tdo

ut_

file

, "

step

s =

%i,

N =

%i,

T =

%lf,

zet

a_F

R =

%lf,

t =

%lf:

\n", st

ep

s, n

um

be

r, T

, ze

ta_

FR

, d

elta

_t*

ste

ps/

nu

mb

er_

of_

blo

cks)

;

do

ub

le e

xp_

en

d_

sqr

= (

sum

_e

nd

/nu

mb

er_

of_

blo

cks)

;

do

ub

le e

xp_

en

d_

sqr_

err

= s

qrt

((su

m_

en

d_

sqr−

sum

_e

nd

*su

m_

en

d/n

um

be

r_o

f_b

lock

s)/n

um

be

r_o

f_b

lock

s);

d

ou

ble

exp

_e

nd

= s

qrt

(exp

_e

nd

_sq

r);

d

ou

ble

exp

_e

nd

_e

rr =

exp

_e

nd

_sq

r_e

rr / 2

/ e

xp_

en

d;

d

ou

ble

exp

_g

yr_

sqr

= (

sum

_g

yr/n

um

be

r_o

f_b

lock

s);

d

ou

ble

exp

_g

yr_

sqr_

err

= s

qrt

((su

m_

gyr

_sq

r−su

m_

gyr

*su

m_

gyr

/nu

mb

er_

of_

blo

cks)

/nu

mb

er_

of_

blo

cks)

;

do

ub

le e

xp_

gyr

= s

qrt

(exp

_g

yr_

sqr)

;

do

ub

le e

xp_

gyr

_e

rr =

exp

_g

yr_

sqr_

err

/ 2

/ e

xp_

gyr

;

if

(co

mp

on

en

t =

= 0

) {

fp

rin

tf(s

tdo

ut_

file

, "

%i b

onds

: <R

_0 −

R_N

>

=

",

nu

mb

er−

1);

fp

rin

tf(s

tdo

ut_

file

, "

%lf

+−

%lf

\n",

exp

_e

nd

, e

xp_

en

d_

err

); fp

rin

tf(s

tdo

ut_

file

, "

%i b

onds

: <R

_g>

=

", n

um

be

r−1

); fp

rin

tf(s

tdo

ut_

file

, "

%lf

+−

%lf

\n",

exp

_g

yr, e

xp_

gyr

_e

rr);

fp

rin

tf(s

tdo

ut_

file

, "

%i b

onds

: <(R

_0−

R_N

)^2>

/ <

R_g

^2>

=

", n

um

be

r−1

); fp

rin

tf(s

tdo

ut_

file

, "

%lf

+−

%lf

\n",

exp

_e

nd

_sq

r / e

xp_

gyr

_sq

r,

e

xp_

en

d_

sqr

/ e

xp_

gyr

_sq

r *

sq

rt(p

ow

(exp

_e

nd

_sq

r_e

rr/e

xp_

en

d_

sqr,

2)

+ p

ow

(exp

_g

yr_

sqr_

err

/exp

_g

yr_

sqr,

2))

); }

else

{ fp

rin

tf(s

tdo

ut_

file

, "

%i b

onds

: <(R

_0 −

R_N

).%

i>

=

",

nu

mb

er−

1, co

mp

on

en

t);

fp

rin

tf(s

tdo

ut_

file

, "

%lf

+−

%lf

\n",

exp

_e

nd

, e

xp_

en

d_

err

); fp

rin

tf(s

tdo

ut_

file

, "

%i b

onds

: <R

_g.%

i>

=

", n

um

be

r−1

, co

mp

on

en

t);

fp

rin

tf(s

tdo

ut_

file

, "

%lf

+−

%lf

\n",

exp

_g

yr, e

xp_

gyr

_e

rr);

fp

rin

tf(s

tdo

ut_

file

, "

%i b

onds

: <(R

_0−

R_N

).%

i^2>

/ <

R_g

.%i^

2> =

",

nu

mb

er−

1, co

mp

on

en

t, c

om

po

ne

nt)

; fp

rin

tf(s

tdo

ut_

file

, "

%lf

+−

%lf

\n",

exp

_e

nd

_sq

r / e

xp_

gyr

_sq

r,

e

xp_

en

d_

sqr

/ e

xp_

gyr

_sq

r *

sq

rt(p

ow

(exp

_e

nd

_sq

r_e

rr/e

xp_

en

d_

sqr,

2)

+ p

ow

(exp

_g

yr_

sqr_

err

/exp

_g

yr_

sqr,

2))

); }

}

mai

n_a

nal

yse_

len

gth

.cc

8/17

mai

n_an

alys

e_le

ngth

.cc

#include

<io

stre

am

>#include

<io

ma

nip

>#include

<cm

ath

>#include

<fs

tre

am

>#include

<cs

tdlib

>#include

<st

dio

.h>

#include

"po

int.h

h"#include

"co

nsta

nts.

hh"

const

in

t n

um

be

r_o

f_b

lock

s =

10

;

po

int ro

use

_m

od

e(p

oin

t* r

, in

t N

, in

t m

od

e)

{// c

om

pu

tes

the

ro

use

−m

od

e o

f a

giv

en

co

nfig

ura

tion

p

oin

t h

elp

= p

oin

t(0

,0,0

);

for

(in

t i

= 0

; i <

N; i+

+)

h

elp

+=

co

s(m

od

e *

M_

PI *

(i+

0.5

) / N

) *

r[i];

h

elp

/=

(N

);

return

he

lp;

} int

de

lta(

int

i,

int

j) {

return

(i=

=j) ?

1 : 0

;

// s

imp

ly th

e K

ron

eck

er−

de

lta} vo

id r

ea

d_

con

fig(

FIL

E*&

inp

ut, p

oin

t* r

, p

oin

t* v

, in

t n

um

be

r,

int

tim

est

ep

s,

do

ub

le&

en

erg

y)// r

ea

ds

a c

om

ple

te c

on

figu

ratio

n fro

m a

file

{ in

t n

um

be

r2;

fs

can

f(in

pu

t, "

%*s

%i %

*s %

i", &

time

ste

ps,

&n

um

be

r2);

for

(in

t i

= 0

; i <

nu

mb

er;

i++

) {

do

ub

le x

, y,

z, vx

, vy

, vz

; fs

can

f(in

pu

t,"

%lf

%lf

%lf

%lf

%lf

%lf

", &

x, &

y, &

z, &

vx, &

vy, &

vz);

r[

i] =

po

int(

x, y

, z

);

v[

i] =

po

int(

vx, vy

, vz

); }

fs

can

f(in

pu

t, "

%*s

%lf"

,&

en

erg

y);

} int

ma

in(

int

arg

c,

cha

r*

arg

v[])

{

if

(a

rgc

!= 9

) {

ce

rr <

< "

Inva

lid n

umbe

r of

arg

umen

ts −

\n"; ce

rr <

< "

Cal

ling

sequ

ence

: <pr

ogra

m n

ame>

con

stan

t−fil

e ou

tput

−fil

e st

dout

−fil

e m

ode1

alp

ha

" <

< "

mod

e2 b

eta

times

teps

\n"; e

xit(

1);

}

in

t m

od

e1

, m

od

e2

, a

lph

a, b

eta

, tim

est

ep

s; m

od

e1

= a

toi(a

rgv[

4])

; a

lph

a =

ato

i(a

rgv[

5])

; m

od

e2

= a

toi(a

rgv[

6])

; b

eta

=

ato

i(a

rgv[

7])

; tim

est

ep

s =

ato

i(a

rgv[

8])

; co

ut <

< m

od

e1

<<

" "

<<

alp

ha

<<

" "

<<

mo

de

2 <

< "

" <

< b

eta

<<

" "

<<

tim

est

ep

s <

< e

nd

l; in

it_co

nst

an

ts(a

rgv[

1])

;

// r

ea

d c

on

sta

nts

F

ILE

* in

pu

t =

fo

pe

n(a

rgv[

2], "

r");

// o

pe

n file

co

nta

inin

g th

e c

on

figu

ratio

ns

if

(!in

pu

t) {

ce

rr <

< "

Inva

lid o

utpu

t−fil

e−na

me!

\n"; e

xit(

1);

}

in

t n

um

be

r;

cha

r c

;

do

ub

le e

ne

rgy;

in

t tim

est

ep

; fs

can

f(in

pu

t, "

# %

i # %

i", &

time

ste

p, &

nu

mb

er)

; p

oin

t *r

=

new

po

int[n

um

be

r];

p

oin

t *v

=

new

po

int[n

um

be

r];

fc

lose

(in

pu

t);

// r

ese

t re

ad

−b

uffe

r (n

ee

de

d fo

r su

cce

sfu

l re

ad

ing

of th

e

in

pu

t =

fo

pe

n(a

rgv[

2], "

r");

// first

co

nfig

ura

tion

) re

ad

_co

nfig

(in

pu

t, r

, v,

nu

mb

er,

tim

est

ep

, e

ne

rgy)

;

if

((m

od

e1

>=

nu

mb

er)

|| (m

od

e2

>=

nu

mb

er)

) {

// d

isca

rd a

ll in

valid

se

ttin

gs

of m

od

e ce

rr <

< "

Wro

ng m

odes

− n

ot e

noug

h m

onom

ers!

\n";

e

xit(

1);

}

if

((a

lph

a >

3)

|| (a

lph

a <

1)

|| (b

eta

> 3

) ||

(be

ta <

1))

{

// d

isca

rd a

ll in

valid

se

ttin

gs

of co

mp

on

en

t ce

rr <

< "

Wro

ng c

ompo

nent

s −

val

id a

re 1

to 3

!\n"; e

xit(

1);

}

d

ou

ble

ro

use

0 =

ro

use

_m

od

e(r

, n

um

be

r, m

od

e1

).co

mp

on

en

t(a

lph

a);

d

ou

ble

co

rre

l = 0

;

do

ub

le c

orr

el_

sqr

= 0

;

do

ub

le s

um

_ko

rre

l = 0

;

if

(tim

est

ep

s !=

0)

{

// d

isca

rd a

ll in

valid

se

ttin

gs

of tim

est

ep

if

(((

ste

ps/

nu

mb

er_

of_

blo

cks)

%tim

est

ep

s !=

0)

||

// (

blo

ck a

na

lysi

s d

oe

s n

ot g

et e

qu

al n

um

be

r o

f(t

ime

ste

ps%

(ste

ps/

nu

mb

er_

of_

ou

tpu

ts)

!= 0

)) {

// c

on

figu

ratio

ns

for

ea

ch b

lock

) ce

rr <

< "

Tim

este

ps d

o no

t fit

data

!\n"; e

xit(

1);

}

mai

n_a

nal

yse_

rou

se.c

c

for

(in

t i

= 1

; i <

= s

tep

s; i+

+)

{

if

(i%

(ste

ps/

nu

mb

er_

of_

blo

cks)

==

0)

{su

m_

korr

el /

=

do

ub

le(s

tep

s)/n

um

be

r_o

f_b

lock

s/tim

est

ep

s;co

rre

l +=

su

m_

korr

el;

corr

el_

sqr

+=

su

m_

korr

el*

sum

_ko

rre

l;su

m_

korr

el =

0;

}

if

(i%

(ste

ps/

nu

mb

er_

of_

ou

tpu

ts)

==

0)

rea

d_

con

fig(in

pu

t, r

, v,

nu

mb

er,

tim

est

ep

, e

ne

rgy)

;

if

(i %

tim

est

ep

s =

= 0

) {

sum

_ko

rre

l +=

ro

use

_m

od

e(r

, n

um

be

r, m

od

e2

).co

mp

on

en

t(b

eta

) *

rou

se0

;ro

use

0 =

ro

use

_m

od

e(r

, n

um

be

r, m

od

e1

).co

mp

on

en

t(a

lph

a);

}

}

}

else

{

for

(in

t i

= 1

; i <

= s

tep

s; i+

+)

{

if

(i%

(ste

ps/

nu

mb

er_

of_

blo

cks)

==

0)

{su

m_

korr

el /

=

do

ub

le(s

tep

s)/n

um

be

r_o

f_b

lock

s;co

rre

l +=

su

m_

korr

el;

corr

el_

sqr

+=

su

m_

korr

el*

sum

_ko

rre

l;su

m_

korr

el =

0;

}

if

(i%

(ste

ps/

nu

mb

er_

of_

ou

tpu

ts)

==

0)

{re

ad

_co

nfig

(in

pu

t, r

, v,

nu

mb

er,

tim

est

ep

, e

ne

rgy)

;su

m_

korr

el +

= r

ou

se_

mo

de

(r, n

um

be

r, m

od

e1

).co

mp

on

en

t(a

lph

a)

* ro

use

_m

od

e(r

, n

um

be

r, m

od

e2

).co

mp

on

en

t(b

eta

) *

do

ub

le(s

tep

s)/n

um

be

r_o

f_o

utp

uts

; }

}

}

// −

−−

−−

−−

−−

−−

−−

−−

−−

−−

−−

−−

−−

−−

−−

−−

−−

−−

−−

−−

−−

−−

−−

−−

−−

−−

−−

−−

−−

−−

−−

−−

−−

−−

−−

−−

−−

−−

−−

−−

−−

−−

−−

−−

−−

−−

−−

−−

−−

−−

−−

−−

−−

F

ILE

* st

do

ut_

file

= fo

pe

n(a

rgv[

3], "

a");

fp

rin

tf(s

tdo

ut_

file

, "

delta

_t =

%lf,

zet

a_F

R*d

elta

_t =

%lf,

ste

ps/s

ampl

e =

%i,

epsi

lon_

SP

= %

lf\n

", d

elta

_t, z

eta

_F

R*d

elta

_t, s

tep

s/n

um

be

r_o

f_b

lock

s, e

psi

lon

_S

P);

fp

rin

tf(s

tdo

ut_

file

, "

step

s =

%i,

N =

%i,

T =

%lf,

zet

a_F

R =

%lf,

t =

%lf:

\n", st

ep

s, n

um

be

r, T

, ze

ta_

FR

, d

elta

_t*

ste

ps/

nu

mb

er_

of_

blo

cks)

; fp

rin

tf(s

tdo

ut_

file

, "

<X

_%i,%

i(%lf)

,X_%

i,%i(0

)> =

", m

od

e1

, a

lph

a, d

elta

_t*

time

ste

ps,

mo

de

2, b

eta

);

do

ub

le k

orr

el =

co

rre

l/nu

mb

er_

of_

blo

cks;

d

ou

ble

err

_ko

rre

l = s

qrt

((co

rre

l_sq

r−co

rre

l*co

rre

l/(n

um

be

r_o

f_b

lock

s))/

(nu

mb

er_

of_

blo

cks−

1))

; fp

rin

tf(s

tdo

ut_

file

, "

%lf

+−

%lf"

, ko

rre

l, e

rr_

korr

el);

if

(m

od

e1

!=

0)

{

do

ub

le d

iscr

imin

an

t =

4*e

psi

lon

_S

P*s

in(m

od

e1

*M_

PI/2

/nu

mb

er)

*sin

(mo

de

1*M

_P

I/2

/nu

mb

er)

− z

eta

_F

R*z

eta

_F

R/4

;

do

ub

le X

_p

0 =

de

lta(m

od

e1

,mo

de

2)*

de

lta(a

lph

a, b

eta

)*e

psi

lon

_B

M*T

/ (8

*ep

silo

n_

SP

*nu

mb

er*

sin

(mo

de

1*M

_P

I/2

/nu

mb

er)

*sin

(mo

de

1*M

_P

I/2

/nu

mb

er)

);

do

ub

le t =

de

lta_

t*tim

est

ep

s;

if

(d

iscr

imin

an

t >

0)

{

do

ub

le o

me

ga

= s

qrt

(dis

crim

ina

nt)

;

do

ub

le e

xpe

cte

d =

X_

p0

* e

xp(−

zeta

_F

R/2

*t)

* (c

os(

om

eg

a*t

) +

ze

ta_

FR

/2/o

me

ga

*sin

(om

eg

a*t

));

fp

rin

tf(s

tdo

ut_

file

, "

, exp

ecte

d =

%lf,

osz

illan

ting", e

xpe

cte

d);

if

((k

orr

el−

err

_ko

rre

l < e

xpe

cte

d)

&&

(ko

rre

l+e

rr_

korr

el >

exp

ect

ed

))fp

rin

tf(s

tdo

ut_

file

, "

, with

in e

rror

bar

s.\n")

;

else

fp

rin

tf(s

tdo

ut_

file

, "

, out

of e

rror

bar

s.\n")

; fp

rin

tf(s

tdo

ut_

file

, "

X_p

0 =

%lf,

om

ega

= %

lf\n", X

_p

0, o

me

ga

); }

if

(d

iscr

imin

an

t =

= 0

) {

do

ub

le e

xpe

cte

d =

X_

p0

* e

xp(−

zeta

_F

R/2

*t)

* (1

+ z

eta

_F

R/2

*t);

fp

rin

tf(s

tdo

ut_

file

, "

, exp

ecte

d =

%lf,

ape

riodi

c", e

xpe

cte

d);

if

((k

orr

el−

err

_ko

rre

l < e

xpe

cte

d)

&&

(ko

rre

l+e

rr_

korr

el >

exp

ect

ed

))fp

rin

tf(s

tdo

ut_

file

, "

, with

in e

rror

bar

s.\n")

;

else

fp

rin

tf(s

tdo

ut_

file

, "

, out

of e

rror

bar

s.\n")

; fp

rin

tf(s

tdo

ut_

file

, "

X_p

0 =

%lf\

n", X

_p

0);

}

if

(d

iscr

imin

an

t <

0)

{

do

ub

le g

am

ma

= s

qrt

(−d

iscr

imin

an

t);

do

ub

le a

1 =

0.5

* (

1 +

ze

ta_

FR

/2/g

am

ma

);

do

ub

le a

2 =

0.5

* (

1 −

ze

ta_

FR

/2/g

am

ma

);

do

ub

le e

xpe

cte

d =

X_

p0

* e

xp(−

zeta

_F

R/2

*t)

* (a

1*e

xp(g

am

ma

*t)

+ a

2*e

xp(−

ga

mm

a*t

));

fp

rin

tf(s

tdo

ut_

file

, "

, exp

ecte

d =

%lf,

dam

ped

", e

xpe

cte

d);

if

((k

orr

el−

err

_ko

rre

l < e

xpe

cte

d)

&&

(ko

rre

l+e

rr_

korr

el >

exp

ect

ed

))fp

rin

tf(s

tdo

ut_

file

, "

, with

in e

rror

bar

s.\n")

;

else

fp

rin

tf(s

tdo

ut_

file

, "

, out

of e

rror

bar

s.\n")

; fp

rin

tf(s

tdo

ut_

file

, "

X_p

0 =

%lf,

gam

ma

= %

lf\n", X

_p

0, g

am

ma

); fp

rin

tf(s

tdo

ut_

file

, "

a1 =

%lf,

a2

= %

lf\n",

a1

, a

2);

}

}

fc

lose

(std

ou

t_fil

e);

}

mai

n_a

nal

yse_

rou

se.c

c

9/17

mai

n_an

alys

e_ro

use.

cc

#include

<io

stre

am

>#include

<io

ma

nip

>#include

<cm

ath

>#include

<fs

tre

am

>#include

<cs

tdlib

>#include

<st

dio

.h>

#include

<st

rin

g.h

>#include

"po

int.h

h"#include

"co

nsta

nts.

hh"

void

re

ad

_d

ata

(F

ILE

*& in

pu

t,

int

& tim

est

ep

s,

do

ub

le*

valu

e)

// r

ea

ds

the

da

ta fro

m th

e s

tdo

ut−

file

pro

du

ced

by

the

ma

in p

rog

ram

{ fs

can

f(in

pu

t, "

%i"

, &

time

ste

ps)

;

for

(in

t i

= 1

; i <

10

; i+

+)

fs

can

f(in

pu

t, "

%lf"

, &

valu

e[i]

);} vo

id c

rea

te_

tem

p_

file

s(ch

ar

* fil

en

am

e)

{

F

ILE

* tm

p =

fo

pe

n("

tmp.

out",

"r"

);

// if

ne

cce

ssa

ry, re

mo

ve o

ld te

mp

−fil

es

if

(tm

p)

{ fc

lose

(tm

p);

sy

ste

m("

rm −

f tm

p.ou

t");

}

tm

p =

fo

pe

n("

tmp.

num"

, "

r");

if

(tm

p)

{ fc

lose

(tm

p);

sy

ste

m("

rm −

f tm

p.nu

m");

}

tm

p =

fo

pe

n("

tmp.

step

s", "

r");

if

(tm

p)

{ fc

lose

(tm

p);

sy

ste

m("

rm −

f tm

p.st

eps")

; }

ch

ar

* co

mm

an

d =

new

ch

ar

[str

len

(file

na

me

) +

33

];

// c

rea

te n

ew

te

mp

−fil

es

sp

rin

tf(c

om

ma

nd

, "

cat %

s | g

rep

−v

’#’ >

tmp.

out

", file

na

me

); sy

ste

m(c

om

ma

nd

); sp

rin

tf(c

om

ma

nd

, "

cat %

s | g

rep

N >

tmp.

num", file

na

me

); sy

ste

m(c

om

ma

nd

); sp

rin

tf(c

om

ma

nd

, "

tail

−n

1 %

s >

tmp.

step

s", file

na

me

); sy

ste

m(c

om

ma

nd

);

delete

[] co

mm

an

d;

} int

ma

in(

int

arg

c,

cha

r*

arg

v[])

{

if

(a

rgc

!= 6

) {

ce

rr <

< "

Inva

lid n

umbe

r of

arg

umen

ts −

\n"; ce

rr <

< "

Cal

ling

sequ

ence

: <pr

ogra

m n

ame>

con

stan

t−fil

e ou

tput

−fil

e st

dout

_file

num

ber−

of−

bloc

ks"

<<

" c

ompo

nent

\n";

e

xit(

1);

}

in

it_co

nst

an

ts(a

rgv[

1])

;

// R

ea

d c

on

sta

nts

in

t s

ave

_st

ep

s =

ste

ps

/ n

um

be

r_o

f_o

utp

uts

;

int

co

mp

on

en

t =

ato

i(a

rgv[

5])

;

if

((c

om

po

ne

nt <

1)

|| (c

om

po

ne

nt >

9))

{ ce

rr <

< "

Inva

lid c

ompo

nent

− v

alid

are

1 to

9!\n";

e

xit(

1);

}

F

ILE

* in

pu

t =

fo

pe

n(a

rgv[

2], "

r");

// c

he

ck if

ou

tpu

t−fil

e e

xist

s

if

(!in

pu

t) {

ce

rr <

< "

Inva

lid o

utpu

t−fil

e−na

me!

\n"; e

xit(

1);

}

fc

lose

(in

pu

t);

cr

ea

te_

tem

p_

file

s(a

rgv[

2])

; in

pu

t =

fo

pe

n("

tmp.

num"

, "

r");

// R

ea

d n

um

be

r o

f m

on

om

ers

in

t n

um

be

r_o

f_m

on

om

ers

; fs

can

f(in

pu

t, "

# N

= %

i", &

nu

mb

er_

of_

mo

no

me

rs);

fc

lose

(in

pu

t);

in

pu

t =

fo

pe

n("

tmp.

step

s", "

r");

fs

can

f(in

pu

t, "

%i"

, &

ste

ps)

; fc

lose

(in

pu

t);

in

t n

um

be

r_o

f_b

lock

s =

ato

i(a

rgv[

4])

; n

um

be

r_o

f_o

utp

uts

= s

tep

s / sa

ve_

ste

ps;

n

um

be

r_o

f_o

utp

uts

−=

nu

mb

er_

of_

ou

tpu

ts %

nu

mb

er_

of_

blo

cks;

st

ep

s =

nu

mb

er_

of_

ou

tpu

ts *

sa

ve_

ste

ps;

in

t tim

est

ep

s =

0;

in

pu

t =

fo

pe

n("

tmp.

out",

"r"

);

// O

pe

n c

lea

ne

d o

utp

ut−

file

d

ou

ble

* d

ata

_se

t =

new

d

ou

ble

[10

];

do

ub

le s

um

_x

= 0

;

do

ub

le s

um

_x_

sqr

= 0

;

do

ub

le s

um

_y

= 0

;

do

ub

le s

um

_y_

sqr

= 0

;

do

ub

le s

um

_xy

= 0

;

do

ub

le s

um

_sl

op

e =

0;

mai

n_a

nal

yse_

spee

d.c

c

do

ub

le s

um

_sl

op

e_

sqr

= 0

;

for

(in

t i

= s

ave

_st

ep

s; i

<=

ste

ps;

i +

= s

ave

_st

ep

s) {

int

re

ad

_tim

est

ep

s =

−1

;

while

(re

ad

_tim

est

ep

s <

i) re

ad

_d

ata

(in

pu

t, r

ea

d_

time

ste

ps,

da

ta_

set)

; su

m_

x +

= i;

su

m_

x_sq

r +

=

do

ub

le(i)*

do

ub

le(i);

su

m_

y +

= d

ata

_se

t[co

mp

on

en

t];

su

m_

y_sq

r +

= d

ata

_se

t[co

mp

on

en

t] *

da

ta_

set[co

mp

on

en

t];

su

m_

xy +

= i

* d

ata

_se

t[co

mp

on

en

t];

if

(i %

(st

ep

s/n

um

be

r_o

f_b

lock

s) =

= 0

) {

do

ub

le s

lop

e =

((1

.0*n

um

be

r_o

f_o

utp

uts

/nu

mb

er_

of_

blo

cks)

*su

m_

xy −

su

m_

x*su

m_

y)/

((1

.0*n

um

be

r_o

f_o

utp

uts

/nu

mb

er_

of_

blo

cks)

*su

m_

x_sq

r−su

m_

x*su

m_

x);

su

m_

slo

pe

+=

slo

pe

; su

m_

slo

pe

_sq

r +

= s

lop

e*s

lop

e;

su

m_

x =

0;

su

m_

y =

0;

su

m_

x_sq

r =

0;

su

m_

y_sq

r =

0;

su

m_

xy =

0;

}

}

fc

lose

(in

pu

t);

d

ou

ble

slo

pe

= s

um

_sl

op

e / n

um

be

r_o

f_b

lock

s;

do

ub

le s

lop

e_

err

= s

qrt

((su

m_

slo

pe

_sq

r−su

m_

slo

pe

*su

m_

slo

pe

/nu

mb

er_

of_

blo

cks)

/(n

um

be

r_o

f_b

lock

s−1

));

F

ILE

* o

utp

ut =

fo

pe

n(a

rgv[

3], "

a");

if

((c

om

po

ne

nt =

= 1

) &

& e

psi

lon

_e

l_x

!= 0

) {

// w

he

n a

na

lysi

ng

th

e first

co

mp

on

en

t, c

om

pu

te th

e m

ob

ility

fp

rin

tf(o

utp

ut, "

# N

m

obili

ty

err

ste

ps

h

_top

h

_bot

l_

bot

l_to

p ti

lt f_

el\n

");

fp

rin

tf(o

utp

ut, "

%5i

%lf

%lf

%i

%lf

%lf

%lf

%lg

%lg

%lg

\n", n

um

be

r_o

f_m

on

om

ers

, −

slo

pe

/de

lta_

t/e

psi

lon

_e

l_x,

−sl

op

e_

err

/de

lta_

t/e

psi

lon

_e

l_x,

ste

ps,

h_

top

_tr

ap

, h

_b

otto

m_

tra

p, l_

bo

tto

m_

tra

p, l_

top

_tr

ap

, w

all_

tilt_

tra

p, −

ep

silo

n_

el_

x);

}

else

{

// o

the

rwis

e c

om

pu

te th

e s

lop

e fp

rin

tf(o

utp

ut, "

# N

sl

ope

err

step

s

h

_top

h

_bot

l_

bot

l_to

p ti

lt\n

");

fp

rin

tf(o

utp

ut, "

%5i

%lf

%lf

%i

%lf

%lf

%lf

%lg

%lg

\n", n

um

be

r_o

f_m

on

om

ers

, sl

op

e, sl

op

e_

err

, st

ep

s, h

_to

p_

tra

p,

h

_b

otto

m_

tra

p, l_

bo

tto

m_

tra

p, l_

top

_tr

ap

, w

all_

tilt_

tra

p);

}

fc

lose

(ou

tpu

t);

sy

ste

m("

rm −

f tm

p.ou

t");

// r

em

ove

te

mp

−fil

es

sy

ste

m("

rm −

f tm

p.nu

m");

sy

ste

m("

rm −

f tm

p.st

eps")

;}

mai

n_a

nal

yse_

spee

d.c

c

10/1

7m

ain_

anal

yse_

spee

d.cc

#include<iostream>

#include<iomanip>

#include<cmath>

#include<fstream>

#include<cstdlib>

#include<gsl/gsl_rng.h>

#include<gsl/gsl_randist.h>

#include<stdio.h>

#include "

poin

t.hh"

#include "

ipoi

nt.h

h"#include "

ipoi

nt_h

ash.

hh"#include "

cons

tant

s.hh"

#include "

pote

ntia

ls.h

h"#include "

verle

t.hh"

void

output(point* r, point* v, node**& box_ref, hash_ipoint idx_hash,

int number,

int total_timesteps,

cha

r* filename,

cha

r* open_mode) {

// T

his

fu

nct

ion

sa

ves

the

cu

rre

nt co

nfig

ura

tion

into

a file

. It is

po

ssib

le to

ap

pe

nd

or

ove

rwrite

// th

e c

he

ckp

oin

t.

FIL

E* out;

out = fopen(filename, open_mode);

fprintf(out, "

# %

i\n", total_timesteps);

fprintf(out, "

# %

i\n", number);

for (

int k = 0; k<number; k++) {

fprintf(out, "

%lf

%lf

%lf

", r[k].x_(), r[k].y_(), r[k].z_());

fprintf(out, "

%lf

%lf

%lf

\n", v[k].x_(), v[k].y_(), v[k].z_());

}

fprintf(out, "

# %

lf\n", energy_all(r, v, box_ref, number));

fprintf(out, "

\n");

fclose(out);

} point rouse_mode(point* r,

int N,

int mode) {

// C

om

pu

tes

the

Ro

use

−M

od

e o

f a

giv

en

co

nfig

ura

tion

. point help = point(0,0,0);

for (

int i = 0; i < N; i++)

help += cos(mode * M_PI * (i+0.5) / N) * r[i];

help /= (N);

return help;

} int main(

int argc,

cha

r* argv[]) {

if ((argc != 6) && (argc != 7)) {

cerr << "

Inva

lid n

umbe

r of

arg

umen

ts −

\n";

cerr << "

Cal

ling

sequ

ence

: <pr

ogra

m n

ame>

con

stan

t−fil

e st

art−

conf

ig−

file

num

ber−

of−

mon

omer

s ";

cerr << "

chec

kpoi

nt−

file

stdo

ut−

file

[app

end−

conf

ig−

file]

\n";

exit(1);

}

init_constants(argv[1]);

if (argc == 7)

append_configs = 1;

else

append_configs = 0;

if (display_constants)

displayconstants(argv[5]);

point *r;

point *v;

node **box_ref;

int number;

int old_timesteps;

hash_ipoint idx_hash;

init_polymer (r, v, box_ref, idx_hash, number, argv[2], argv[3], old_timesteps);

int enable_el_tmp = enable_el;

enable_el = 0;

// d

isa

ble

ele

ctric

field

du

rin

g in

itia

lisa

tion

for (

int i = 0; i < throw_away; i++) {

timestep(r, v, box_ref, idx_hash, number, i);

}

enable_el = enable_el_tmp;

FIL

E* stdout_file = fopen(argv[5], "

a");

fprintf(stdout_file, "

# N

= %

i\n", number);

fclose(stdout_file);

for (

int i = 0; i <= steps; i++) {

if (i % (steps/number_of_outputs) == 0) {

// s

ave

wh

en

th

e tim

e h

as

com

e point center = rouse_mode(r, number, 0);

do

ub

le min_x = r[0].x_();

do

ub

le max_x = r[0].x_();

do

ub

le min_y = r[0].y_();

do

ub

le max_y = r[0].y_();

do

ub

le min_z = r[0].z_();

do

ub

le max_z = r[0].z_();

for (

int j = 0; j < number; j++) {

if (min_x > r[j].x_()) min_x = r[j].x_();

if (max_x < r[j].x_()) max_x = r[j].x_();

if (min_y > r[j].y_()) min_y = r[j].y_();

if (max_y < r[j].y_()) max_y = r[j].y_();

if (min_z > r[j].z_()) min_z = r[j].z_();

mai

n_p

oly

mer

.cc

if (max_z < r[j].z_()) max_z = r[j].z_();

}

FIL

E* stdout_file = fopen(argv[5], "

a");

fprintf(stdout_file, "

%i

%lf

%lf

%lf

", i+old_timesteps, center.x_(), center.y_(), center.z_());

fprintf(stdout_file, "

%lf

%lf

%lf

", min_x, min_y, min_z);

fprintf(stdout_file, "

%lf

%lf

%lf

\n", max_x, max_y, max_z);

fclose(stdout_file);

output(r, v, box_ref, idx_hash, number, i + old_timesteps, argv[4], "

w");

if ((append_configs) && ((i != 0) || (throw_away != 0)))

output(r, v, box_ref, idx_hash, number, i + old_timesteps, argv[6], "

a");

}

if (i != steps) timestep(r, v, box_ref, idx_hash, number, i);

}

output(r, v, box_ref, idx_hash, number, steps + old_timesteps, argv[4], "

w");

}

mai

n_p

oly

mer

.cc

11/1

7m

ain_

poly

mer

.cc

po

lym

er : m

ain

_p

oly

me

r.o

co

nst

an

ts.o

ve

rle

t.o

ipo

int_

ha

sh.o

po

ten

tials

.og

++

−I/u

sr/lo

cal/i

ncl

ud

e −

L/u

sr/lo

cal/l

ib −

lgsl

−lg

slb

las

−lm

−O

3 −

mfp

−ro

un

din

g−

mo

de

=d

−o

po

lym

er.

ou

t m

ain

_p

oly

me

r.o

co

nst

an

ts.o

ve

rle

t.o

ipo

int_

ha

sh.o

po

ten

tials

.oan

alys

e_ro

use

: m

ain

_a

na

lyse

_ro

use

.o c

on

sta

nts

.o

g+

+ −

I/u

sr/lo

cal/i

ncl

ud

e −

L/u

sr/lo

cal/l

ib −

lgsl

−lg

slb

las

−lm

−O

3 −

mfp

−ro

un

din

g−

mo

de

=d

−o

an

aly

se_

rou

se

.ou

t m

ain

_a

na

lyse

_ro

use

.o c

on

sta

nts

.oan

alys

e_le

ng

th : m

ain

_a

na

lyse

_le

ng

th.o

co

nst

an

ts.o

g

++

−I/u

sr/lo

cal/i

ncl

ud

e −

L/u

sr/lo

cal/l

ib −

lgsl

−lg

slb

las

−lm

−O

3 −

mfp

−ro

un

din

g−

mo

de

=d

−o

an

aly

se_

len

gth

.ou

t m

ain

_a

na

lyse

_le

ng

th.o

co

nst

an

ts.o

anal

yse_

spee

d : m

ain

_a

na

lyse

_sp

ee

d.o

co

nst

an

ts.o

g

++

−I/u

sr/lo

cal/i

ncl

ud

e −

L/u

sr/lo

cal/l

ib −

lgsl

−lg

slb

las

−lm

−O

3 −

mfp

−ro

un

din

g−

mo

de

=d

−o

an

aly

se_

spe

ed

.ou

t m

ain

_a

na

lyse

_sp

ee

d.o

co

nst

an

ts.o

cl

ear

: rm

*.o

all :

ma

ke p

oly

me

r a

na

lyse

anal

yse

:m

ake

an

aly

se_

len

gth

an

aly

se_

rou

se a

na

lyse

_sp

ee

d

%.o

: %

.cc

g+

+ −

c −

O3

−m

fp−

rou

nd

ing

−m

od

e=

d −

I/u

sr/lo

cal/i

ncl

ud

e %

.cc

$<

mak

efile

#include

<io

stre

am

>#include

<io

ma

nip

>#include

<cm

ath

>

#ifndef

PO

INT

_IN

CL

UD

ED

#define

PO

INT

_IN

CL

UD

ED

// T

his

is a

cla

ss o

f 3

−d

ime

nsi

on

al p

oin

ts w

ith th

e c

om

mo

n a

lge

bra

ic r

ule

s.// T

he

fu

nct

ion

ab

s_sq

r ca

lcu

late

s th

e s

qa

re o

f th

e le

ng

th o

f a

po

int,

// it

avo

ide

s ta

kin

g th

e s

qu

are

−ro

ot.

class

po

int{

d

ou

ble

x;

do

ub

le y

; d

ou

ble

z;

public

: p

oin

t()

{}

p

oin

t(d

ou

ble

xx,

d

ou

ble

yy,

d

ou

ble

zz)

{

x

= x

x; y

= y

y; z

= z

z;

}

friend

d

ou

ble

ab

s(p

oin

t p

){

return

sq

rt(p

.x*p

.x +

p.y

*p.y

+ p

.z*p

.z);

}

friend

d

ou

ble

ab

s_sq

r(p

oin

t p

){

return

(p

.x*p

.x +

p.y

*p.y

+ p

.z*p

.z);

}

friend

ost

rea

m&

operator

<<

(ost

rea

m&

os,

po

int p

){

return

os

<<

"("

<<

se

tw(9

) <

< p

.x <

< "

," <

< s

etw

(9)

<<

p.y

<<

","

<

< s

etw

(9)

<<

p.z

<<

")"

; }

d

ou

ble

co

mp

on

en

t(in

t c

ho

ice

){

// a

llow

s d

ire

ct a

cce

ss to

a s

ing

le

return

(ch

oic

e=

=1

? x

: (

cho

ice

==

2 ?

y : z

));

// c

om

po

ne

nt, 1

.. 3

}

d

ou

ble

x_

(vo

id){

return

x;

}

d

ou

ble

y_

(vo

id){

return

y;

}

d

ou

ble

z_

(vo

id){

return

z;

}

friend

po

int

operator

+(p

oin

t p

1, p

oin

t p

2){

return

po

int(

p1

.x +

p2

.x, p

1.y

+ p

2.y

, p

1.z

+ p

2.z

); }

friend

po

int

operator

*(p

oin

t p

, d

ou

ble

d){

return

po

int(

p.x

* d

, p

.y *

d, p

.z *

d);

}

friend

po

int

operator

*(d

ou

ble

d, p

oin

t p

){

return

p *

d;

}

friend

d

ou

ble

operator

*(p

oin

t p

1, p

oin

t p

2){

return

p1

.x*p

2.x

+ p

1.y

*p2

.y +

p1

.z*p

2.z

; }

friend

po

int

operator

−(p

oin

t p

){

return

(−

1)*

p;

}

friend

po

int

operator

−(p

oin

t p

1, p

oin

t p

2){

return

p1

+(−

p2

); }

friend

po

int

operator

/(p

oin

t p

, d

ou

ble

d){

return

p*(

1/d

); }

p

oin

t

operator

*=(

do

ub

le d

){

x

*= d

; y*

=d

; z*

=d

;

return

po

int(

x,y,

z);

}

p

oin

t

operator

+=

(po

int p

){ x

+=

p.x

; y

+=

p.y

; z

+=

p.z

;

return

po

int(

x,y,

z);

}

p

oin

t

operator

−=

(po

int p

){ x

−=

p.x

; y

−=

p.y

; z

−=

p.z

;

return

po

int(

x,y,

z);

}

p

oin

t

operator

/=(

do

ub

le d

){

x

/= d

; y

/= d

; z

/= d

;

return

po

int(

x,y,

z);

}

friend

b

oo

l

operator

==

(po

int p

1, p

oin

t p

2){

return

(p

1.x

==

p2

.x)

&&

(p

1.y

==

p2

.y)

&&

(p

1.z

==

p2

.z);

}

friend

b

oo

l

operator

!=(p

oin

t p

1, p

oin

t p

2){

return

!(p

1 =

= p

2);

}

};

po

int.

hh

12/1

7m

akef

ile, p

oint

.hh

#endif

po

int.

hh

#include

<cm

ath

>#include

<st

dio

.h>

#include

<g

sl/g

sl_

rng

.h>

#include

<g

sl/g

sl_

ran

dis

t.h

>#include

"po

int.h

h"#include

"co

nsta

nts.

hh"

#ifndef

PO

TE

NT

IAL

S_

INC

LU

DE

D#define

PO

TE

NT

IAL

S_

INC

LU

DE

D

// T

his

file

co

nta

ins

the

inte

ract

ion

po

ten

tials

fo

r th

e m

on

om

ers

.// // −

Th

e L

en

na

rd−

Jon

es

po

ten

tial b

etw

ee

n tw

o m

on

om

ers

, w

hic

h is

cu

t o

ff a

t a

// le

ng

th o

f r

= 2

*sig

ma

// −

Th

e b

on

d−

len

gth

po

ten

tial b

etw

ee

n tw

o b

on

ds

on

th

e c

ha

in, lim

itin

g// th

e m

axi

mu

m b

on

d−

len

gth

// −

Th

e b

on

d−

an

gle

−p

ote

ntia

l, w

hic

h p

refe

rs b

on

d−

an

gle

= 0

// −

Th

e e

ntr

op

ic−

tra

p p

ote

ntia

l// −

Th

e s

prin

g−

po

ten

tial f

or

the

Ro

use

−m

od

el

// −

Th

e b

row

nia

n m

otio

n−

forc

e, w

itho

ut a

po

ten

tial

// −

Th

e e

xte

rna

l ele

ctric

field

, w

itho

ut a

po

ten

tial

do

ub

le V

_L

J(p

oin

t r)

{

do

ub

le r

_2

= a

bs_

sqr(

r);

if

(r_

2 >

l_cu

toff_

LJ_

sqr)

return

0;

d

ou

ble

r_

4 =

r_

2 *

r_

2;

return

ep

silo

n_

LJ*

(sig

ma

_6

*sig

ma

_6

/(r_

4*r

_4

*r_

4)

− s

igm

a_

6/(

r_4

*r_

2))

+

v_

cuto

ff_

LJ;

} po

int d

_V

_L

J_d

r_tim

es_

de

lta_

t_2

(po

int r)

{

do

ub

le r

_2

= a

bs_

sqr(

r);

if

(r_

2 >

l_cu

toff_

LJ_

sqr)

return

po

int(

0,0

,0);

d

ou

ble

r_

4 =

r_

2*r

_2

;

do

ub

le r

_8

= r

_4

*r_

4;

return

r *

(p

refa

cto

r_L

J *

(1/r

_8

− 2

*sig

ma

_6

/(r_

8*r

_4

*r_

2))

);} d

ou

ble

V_

BL

(d

ou

ble

d)

{

return

−e

psi

lon

_B

L/2

* d

_B

L_

sqr

* s

td::lo

g(1

− (

d−

d_

0_

BL

)*(d

−d

_0

_B

L)

/ d

_B

L_

sqr)

;} p

oin

t d

_V

_B

L_

dr_

time

s_d

elta

_t_

2(p

oin

t r,

d

ou

ble

d)

{

do

ub

le h

elp

= p

refa

cto

r_B

L / (

1 −

(d

−d

_0

_B

L)*

(d−

d_

0_

BL

)/d

_B

L_

sqr)

;

return

r *

((d

−d

_0

_B

L)/

d *

he

lp);

} do

ub

le V

_B

A(p

oin

t a

, p

oin

t b

, d

ou

ble

r_

a,

do

ub

le r

_b

) {

return

ep

silo

n_

BA

* (

1 −

a*b

/ (

r_a

*r_

b))

;} p

oin

t d

_V

_B

A_

da

_tim

es_

de

lta_

t_2

(po

int a

, p

oin

t b

, d

ou

ble

r_

a,

do

ub

le r

_b

) {

if

(a

==

b)

return

po

int(

0,0

,0);

p

oin

t h

elp

= (

a*b

) *

r_b

/r_

a *

a −

b *

r_

a *

r_

b;

return

he

lp *

(p

refa

cto

r_B

A / (

r_a

*r_

a *

r_

b*r

_b

));

} po

int d

_V

_B

A_

db

_tim

es_

de

lta_

t_2

(po

int a

, p

oin

t b

, d

ou

ble

r_

a,

do

ub

le r

_b

) {

if

(a

==

b)

return

po

int(

0,0

,0);

p

oin

t h

elp

= (

a*b

) *

r_a

/r_

b *

b −

a *

r_

a *

r_

b;

return

he

lp *

(p

refa

cto

r_B

A / (

r_a

*r_

a *

r_

b*r

_b

));

} do

ub

le V

_S

P(p

oin

t r)

{

return

ep

silo

n_

SP

/2 *

ab

s_sq

r(r)

;} p

oin

t d

_V

_S

P_

dr_

time

s_d

elta

_t_

2(p

oin

t r)

{

return

pre

fact

or_

SP

* r

;} p

oin

t d

_V

_B

M_

dr_

time

s_d

elta

_t(

void

) {

static

in

t in

it =

0;

static

gsl

_rn

g*

ran

do

m =

gsl

_rn

g_

allo

c(g

sl_

rng

_ta

us)

;

if

(!in

it) {

in

it =

1;

g

sl_

rng

_se

t(ra

nd

om

, ra

nd

_se

ed

); ra

nd

_se

ed

= g

sl_

rng

_g

et(

ran

do

m);

if

(sa

ve_

ne

w_

see

d)

{

FIL

E*

rse

ed

= fo

pe

n("

/hom

e/sc

hmid

/str

eek/

.ran

dom

_see

d", "

w")

; fp

rin

tf(r

see

d, "

%lu

\n",

ra

nd

_se

ed

); fc

lose

(rse

ed

); }

}

d

ou

ble

x =

gsl

_rn

g_

un

iform

(ra

nd

om

)−0

.5;

po

ten

tial

s.cc

13/1

7po

int.h

h, p

oten

tials

.cc

double

y =

gsl

_rn

g_

un

iform

(ra

nd

om

)−0

.5;

double

z =

gsl

_rn

g_

un

iform

(ra

nd

om

)−0

.5;

while

(w

hic

h_

BM

&&

(x*

x +

y*y

+ z

*z >

0.2

5))

{ x

= g

sl_

rng

_u

nifo

rm(r

an

do

m)−

0.5

; y

= g

sl_

rng

_u

nifo

rm(r

an

do

m)−

0.5

; z

= g

sl_

rng

_u

nifo

rm(r

an

do

m)−

0.5

; }

return

pre

fact

or_

BM

* p

oin

t(x,

y, z)

;} double

V_

LJ_

tra

p(

double

r_

2)

{

if

(r_

2 >

l_cu

toff_

LJ_

tra

p_

sqr)

return

0;

double

r_

4 =

r_

2 *

r_

2;

return

ep

silo

n_

LJ_

tra

p*(

sig

ma

_6

_tr

ap

*sig

ma

_6

_tr

ap

/(r_

4*r

_4

*r_

4)

− s

igm

a_

6_

tra

p/(

r_4

*r_

2))

+

v_

cuto

ff_

LJ_

tra

p;

} po

int d

_V

_L

J_tr

ap

_d

r_tim

es_

de

lta_

t_2

(po

int r)

{

double

r_

2 =

ab

s_sq

r(r)

;

if

(r_

2 >

l_cu

toff_

LJ_

tra

p_

sqr)

return

po

int(

0,0

,0);

double

r_

4 =

r_

2*r

_2

;

double

r_

8 =

r_

4*r

_4

;

return

r *

(p

refa

cto

r_L

J_tr

ap

* (

1/r

_8

− 2

*sig

ma

_6

_tr

ap

/(r_

8*r

_4

*r_

2))

);} double

V_

tra

p(p

oin

t r)

{

double

en

erg

y =

0;

double

x =

r.x

_()

;

double

z =

r.z

_()

; x

−=

std

::flo

or(

x / l_

surf

ace

_tr

ap

) *

l_su

rfa

ce_

tra

p;

if

(z

> (

h_

top

_tr

ap

−l_

cuto

ff_

LJ_

tra

p))

e

ne

rgy

+=

V_

LJ_

tra

p((

z−h

_to

p_

tra

p)*

(z−

h_

top

_tr

ap

));

if

(z

< l_

cuto

ff_

LJ_

tra

p)

{

if

(x

> l_

bo

tto

m_

tra

p)

e

ne

rgy

+=

V_

LJ_

tra

p(z

*z);

else

{

if

(z

< l_

cuto

ff_

LJ_

tra

p −

h_

bo

tto

m_

tra

p)

en

erg

y +

= V

_L

J_tr

ap

((z

+ h

_b

otto

m_

tra

p)*

(z +

h_

bo

tto

m_

tra

p))

;

if

(x

< l_

cuto

ff_

LJ_

tra

p +

h_

bo

tto

m_

tra

p *

wa

ll_til

t_tr

ap

) {

double

z_

= 1

/ (

1 +

wa

ll_til

t_tr

ap

*wa

ll_til

t_tr

ap

) *

(z −

wa

ll_til

t_tr

ap

* x

);if

(z_

> 0

) z_

= 0

;double

x_

= −

wa

ll_til

t_tr

ap

* z

_;

if

((x

− x

_)

< l_

cuto

ff_

LJ_

tra

p)

e

ne

rgy

+=

V_

LJ_

tra

p((

x−x_

)*(x

−x_

) +

(z−

z_)*

(z−

z_)

); }

if

(x

> l_

bo

tto

m_

tra

p −

l_cu

toff_

LJ_

tra

p −

h_

bo

tto

m_

tra

p *

wa

ll_til

t_tr

ap

) {

double

z_

= 1

/ (

1 +

wa

ll_til

t_tr

ap

*wa

ll_til

t_tr

ap

) *

(z −

wa

ll_til

t_tr

ap

* (

l_b

otto

m_

tra

p −

x))

;if

(z_

> 0

) z_

= 0

;double

x_

= l_

bo

tto

m_

tra

p +

wa

ll_til

t_tr

ap

* z

_;

if

((x

_ −

x)

< l_

cuto

ff_

LJ_

tra

p)

e

ne

rgy

+=

V_

LJ_

tra

p((

x−x_

)*(x

−x_

) +

(z−

z_)*

(z−

z_))

; }

}

}

return

en

erg

y;} p

oin

t d

_V

_tr

ap

_d

r_tim

es_

de

lta_

t_2

(po

int r)

{ p

oin

t fo

rce

= p

oin

t(0

,0,0

);

double

x =

r.x

_()

;

double

z =

r.z

_()

; x

−=

std

::flo

or(

x / l_

surf

ace

_tr

ap

) *

l_su

rfa

ce_

tra

p;

if

(z

> (

h_

top

_tr

ap

−l_

cuto

ff_

LJ_

tra

p))

fo

rce

+=

d_

V_

LJ_

tra

p_

dr_

time

s_d

elta

_t_

2(p

oin

t(0

, 0

, z−

h_

top

_tr

ap

));

if

(z

< l_

cuto

ff_

LJ_

tra

p)

{

if

(x

> l_

bo

tto

m_

tra

p)

fo

rce

+=

d_

V_

LJ_

tra

p_

dr_

time

s_d

elta

_t_

2(p

oin

t(0

, 0

, z)

);

else

{

if

(z

< l_

cuto

ff_

LJ_

tra

p −

h_

bo

tto

m_

tra

p)

forc

e +

= d

_V

_L

J_tr

ap

_d

r_tim

es_

de

lta_

t_2

(po

int(

0, 0

, z

+ h

_b

otto

m_

tra

p))

;

if

(x

< l_

cuto

ff_

LJ_

tra

p +

h_

bo

tto

m_

tra

p *

wa

ll_til

t_tr

ap

) {

double

z_

= 1

/ (

1 +

wa

ll_til

t_tr

ap

*wa

ll_til

t_tr

ap

) *

(z −

wa

ll_til

t_tr

ap

* x

);if

(z_

> 0

) z_

= 0

;double

x_

= −

wa

ll_til

t_tr

ap

* z

_;

if

((x

− x

_)

< l_

cuto

ff_

LJ_

tra

p)

fo

rce

+=

d_

V_

LJ_

tra

p_

dr_

time

s_d

elta

_t_

2(p

oin

t(x−

x_, 0

, z−

z_))

; }

if

(x

> l_

bo

tto

m_

tra

p −

l_cu

toff_

LJ_

tra

p −

h_

bo

tto

m_

tra

p *

wa

ll_til

t_tr

ap

) {

double

z_

= 1

/ (

1 +

wa

ll_til

t_tr

ap

*wa

ll_til

t_tr

ap

) *

(z −

wa

ll_til

t_tr

ap

* (

l_b

otto

m_

tra

p −

x))

;if

(z_

> 0

) z_

= 0

;double

x_

= l_

bo

tto

m_

tra

p +

wa

ll_til

t_tr

ap

* z

_;

po

ten

tial

s.cc

if

((x

_ −

x)

< l_

cuto

ff_

LJ_

tra

p)

fo

rce

+=

d_

V_

LJ_

tra

p_

dr_

time

s_d

elta

_t_

2(p

oin

t(x−

x_, 0

, z−

z_))

; }

}

}

return

fo

rce

;} p

oin

t d

_V

_e

l_d

r_tim

es_

de

lta_

t_2

(void

) {

return

po

int(

pre

fact

or_

el_

x, p

refa

cto

r_e

l_y,

pre

fact

or_

el_

z);

} #endif

po

ten

tial

s.cc

14/1

7po

tent

ials

.cc

#include

<cm

ath

>#include

<st

dio

.h>

#include

<g

sl/g

sl_

rng

.h>

#include

<g

sl/g

sl_

ran

dis

t.h

>#include

"po

int.h

h"#include

"co

nsta

nts.

hh"

#ifndef

PO

TE

NT

IAL

S_

INC

LU

DE

D#define

PO

TE

NT

IAL

S_

INC

LU

DE

D

// T

his

file

co

nta

ins

the

inte

ract

ion

po

ten

tials

fo

r th

e m

on

om

ers

.// // −

Th

e L

en

na

rd−

Jon

es

po

ten

tial b

etw

ee

n tw

o m

on

om

ers

, w

hic

h is

cu

t o

ff a

t a

// le

ng

th o

f r

= 2

*sig

ma

// −

Th

e b

on

d−

len

gth

po

ten

tial b

etw

ee

n tw

o b

on

ds

on

th

e c

ha

in, lim

itin

g// th

e m

axi

mu

m b

on

d−

len

gth

// −

Th

e b

on

d−

an

gle

−p

ote

ntia

l, w

hic

h p

refe

rs b

on

d−

an

gle

= 0

// −

Th

e e

ntr

op

ic−

tra

p p

ote

ntia

l// −

Th

e s

prin

g−

po

ten

tial f

or

the

Ro

use

−m

od

el

// −

Th

e b

row

nia

n m

otio

n−

forc

e, w

itho

ut a

po

ten

tial

// −

Th

e e

xte

rna

l ele

ctric

field

, w

itho

ut a

po

ten

tial

do

ub

le V

_L

J(p

oin

t r)

;

po

int d

_V

_L

J_d

r_tim

es_

de

lta_

t_2

(po

int r)

;d

ou

ble

V_

BL

(d

ou

ble

d);

po

int d

_V

_B

L_

dr_

time

s_d

elta

_t_

2(p

oin

t r,

d

ou

ble

d);

do

ub

le V

_B

A(p

oin

t a

, p

oin

t b

, d

ou

ble

r_

a,

do

ub

le r

_b

);

po

int d

_V

_B

A_

da

_tim

es_

de

lta_

t_2

(po

int a

, p

oin

t b

, d

ou

ble

r_

a,

do

ub

le r

_b

);p

oin

t d

_V

_B

A_

db

_tim

es_

de

lta_

t_2

(po

int a

, p

oin

t b

, d

ou

ble

r_

a,

do

ub

le r

_b

);

do

ub

le V

_S

P(p

oin

t r)

;

po

int d

_V

_S

P_

dr_

time

s_d

elta

_t_

2(p

oin

t r)

;p

oin

t d

_V

_B

M_

dr_

time

s_d

elta

_t(

void

);

do

ub

le V

_L

J_tr

ap

(d

ou

ble

r_

2);

do

ub

le V

_tr

ap

(po

int r)

;

po

int d

_V

_tr

ap

_d

r_tim

es_

de

lta_

t_2

(po

int r)

;p

oin

t d

_V

_e

l_d

r_tim

es_

de

lta_

t_2

(vo

id);

#endif

po

ten

tial

s.h

h#

0#

10

40

0 −

40

0

0 0

40

1 −

40

0

0 0

40

2 −

40

0

0 0

40

3 −

40

0

0 0

40

4 −

40

0

0 0

40

5 −

40

0

0 0

40

6 −

40

0

0 0

40

7 −

40

0

0 0

40

8 −

40

0

0 0

40

9 −

40

0

0 0

star

tfile

15/1

7po

tent

ials

.hh,

sta

rtfil

e

#include

<fs

tre

am

>#include

<io

stre

am

>#include

<cm

ath

>#include

<st

dio

.h>

#include

"po

int.h

h"#include

"ip

oint

.hh"

#include

"ip

oint

_has

h.hh"

#include

"co

nsta

nts.

hh"#include

"po

tent

ials

.hh"

#ifndef

VE

RL

ET

_IN

CL

UD

ED

#define

VE

RL

ET

_IN

CL

UD

ED

// T

his

file

incl

ud

es

the

Ve

rle

t−A

lgo

rith

m, a

fu

nct

ion

fo

r in

itia

lisa

tion

of th

e p

oly

me

r// a

nd

a fu

nct

ion

fo

r ca

lcu

latin

g th

e e

ne

rgy

of th

e p

oly

me

r.// // F

un

ctio

ns

use

d:

// −

init_

po

lym

er

: re

ad

s th

e p

oly

me

r−co

nfig

ura

tion

fro

m a

n in

pu

t−fil

e// −

LJ−

up

da

te : co

mp

ute

s th

e fo

rce

on

all

mo

no

me

rs c

rea

ted

by

the

Le

nn

ard

−Jo

ne

s−p

ote

ntia

l// −

v−

up

da

te : c

om

pu

tes

all

forc

es

act

ing

on

all

mo

no

me

rs c

rea

ted

by

the

po

ten

tials

,// d

oe

s n

ot co

mp

ute

th

e b

row

nia

n fo

rce

or

the

frict

ion

fo

rce

// −

tim

est

ep

: th

is is

th

e v

erle

t−a

lgo

rith

m u

sin

g a

ll e

na

ble

d fo

rce

s (in

clu

din

g

// frict

ion

an

d b

row

nia

n fo

rce

)// −

en

erg

y_a

ll : co

mp

ute

s th

e e

ne

rgy

of th

e p

oly

me

r, in

clu

din

g a

ll e

na

ble

d p

ote

ntia

lsvo

id in

it_p

oly

me

r(p

oin

t*&

r, p

oin

t*&

v, n

od

e**

& b

ox_

ref, h

ash

_ip

oin

t& id

x_h

ash

, in

t&

nu

mb

er,

ch

ar

* fil

en

am

e,

ch

ar

* n

um

be

r_o

f_m

on

om

ers

, in

t&

old

_tim

est

ep

s){

FIL

E*

inp

ut =

fo

pe

n(f

ilen

am

e, "

r");

if

(!in

pu

t) {

ce

rr <

< "

Inva

lid c

onfig

file

−na

me!

\n";

e

xit(

1);

}

n

um

be

r =

ato

i(n

um

be

r_o

f_m

on

om

ers

);

do

ub

le x

, y,

z;

in

t n

um

be

r_fil

e;

fs

can

f(in

pu

t, "

# %

i # %

i", &

old

_tim

est

ep

s, &

nu

mb

er_

file

);

if

(n

um

be

r >

nu

mb

er_

file

) {

ce

rr <

< "

Too

few

mon

omer

s in

sta

rt−

file!

\n"; e

xit(

1);

}

r

=

new

po

int[n

um

be

r];

v

=

new

po

int[n

um

be

r];

b

ox_

ref =

new

no

de

*[n

um

be

r];

for

(in

t i=

0; i <

nu

mb

er;

i++

) {

fs

can

f(in

pu

t,"

%lf"

, &

x);

fs

can

f(in

pu

t,"

%lf"

, &

y);

fs

can

f(in

pu

t,"

%lf"

, &

z);

r[

i] =

po

int(

x,y,

z);

fs

can

f(in

pu

t,"

%lf"

, &

x);

fs

can

f(in

pu

t,"

%lf"

, &

y);

fs

can

f(in

pu

t,"

%lf"

, &

z);

v[

i] =

po

int(

x,y,

z);

}

for

(in

t i

= 0

; i <

nu

mb

er;

i++

) h

ash

_in

sert

_in

de

x(id

x_h

ash

, r[

i]/l_

cuto

ff_

LJ,

bo

x_re

f[i],

i);

fc

lose

(in

pu

t);

} void

LJ_

up

da

te(p

oin

t* r

, p

oin

t* d

v, n

od

e**

& b

ox_

ref,

int

N)

{

for

(in

t i

= 0

; i <

N; i+

+)

{

int

nu

mb

er;

int

* n

eig

hb

ou

rs =

ne

igh

bo

urs

_o

f_cu

be

(bo

x_re

f[i],

nu

mb

er)

;

for

(in

t j

= 0

; j <

nu

mb

er;

j++

)

if

(n

eig

hb

ou

rs[j]

< i)

{p

oin

t d

_V

_L

J =

d_

V_

LJ_

dr_

time

s_d

elta

_t_

2(r

[ne

igh

bo

urs

[j]] −

r[i]

);d

v[i]

+

= d

_V

_L

J;d

v[n

eig

hb

ou

rs[j]

] +

= −

d_

V_

LJ;

}

delete

[]n

eig

hb

ou

rs;

}

} void

v_

up

da

te(p

oin

t* r

, p

oin

t* v

, p

oin

t* d

v, n

od

e**

& b

ox_

ref,

int

N,

int

tim

est

ep

s){

do

ub

le *

d =

new

d

ou

ble

[N];

if

(e

na

ble

_B

L || e

na

ble

_B

A)

for

(in

t i

= 0

; i <

= N

−2

; i+

+)

{ d

[i] =

ab

s(r[

i+1

] −

r[i]

); }

for

(in

t i

= 0

; i <

N; i+

+)

d

v[i]

= p

oin

t(0

,0,0

);

if

(e

na

ble

_L

J) fa

st_

LJ_

up

da

te(r

, d

v, b

ox_

ref, N

);

if

(e

na

ble

_B

L)

for

(in

t i

= 0

; i <

N−

1; i+

+)

{

verl

et.c

c

if

(d

[i] >

= d

_0

_B

L +

d_

BL

) {

cerr

<<

"R

uptu

re a

t t =

" <

< tim

est

ep

s*d

elta

_t <

< "

(" <

< tim

est

ep

s <

< "

tim

este

ps),

mon

omer

"

<

< i

<<

"!\n

";e

xit(

1);

}

p

oin

t V

_B

L_

dr

= d

_V

_B

L_

dr_

time

s_d

elta

_t_

2(r

[i+1

]−r[

i], d

[i]);

d

v[i]

+

= V

_B

L_

dr;

d

v[i+

1] +

= −

V_

BL

_d

r; }

if

(e

na

ble

_B

A)

for

(in

t i

= 1

; i <

N−

1; i+

+)

{ p

oin

t V

_B

A_

da

= d

_V

_B

A_

da

_tim

es_

de

lta_

t_2

(r[i−

1] −

r[i]

, r[

i] −

r[i+

1], d

[i−1

], d

[i]);

p

oin

t V

_B

A_

db

= d

_V

_B

A_

db

_tim

es_

de

lta_

t_2

(r[i−

1] −

r[i]

, r[

i] −

r[i+

1], d

[i−1

], d

[i]);

d

v[i−

1] +

= −

V_

BA

_d

a;

d

v[i]

+

= (

V_

BA

_d

a −

V_

BA

_d

b);

d

v[i+

1] +

= V

_B

A_

db

; }

if

(e

na

ble

_S

P)

for

(in

t i

= 0

; i <

N−

1; i+

+)

{ p

oin

t d

_V

_S

P =

d_

V_

SP

_d

r_tim

es_

de

lta_

t_2

(r[i+

1]−

r[i])

; d

v[i]

+

= d

_V

_S

P;

d

v[i+

1] +

= −

d_

V_

SP

; }

if

(e

na

ble

_tr

ap

)

for

(in

t i

= 0

; i <

N; i+

+)

d

v[i]

+=

−d

_V

_tr

ap

_d

r_tim

es_

de

lta_

t_2

(r[i]

);

if

(e

na

ble

_e

l) {

for

(in

t i

= 0

; i <

N; i+

+)

d

v[i]

+=

−d

_V

_e

l_d

r_tim

es_

de

lta_

t_2

();

}

delete

[]d

;} vo

id tim

est

ep

(po

int*

& r

, p

oin

t*&

v, n

od

e**

& b

ox_

ref, h

ash

_ip

oin

t& id

x_h

ash

, in

t N

, in

t tim

est

ep

s) {

static

ipo

int*

old

_cu

be

;

static

po

int *d

v;

static

in

t in

it =

0;

in

t c

ha

ng

ed

= 0

;

if

(!in

it) {

in

it =

1;

d

v =

new

po

int[N

]; v_

up

da

te(r

, v,

dv,

bo

x_re

f, N

, tim

est

ep

s);

o

ld_

cub

e =

new

ipo

int[N

];

for

(in

t i

= 0

; i <

N; i+

+)

o

ld_

cub

e[i]

= ip

oin

t(r[

i]/l_

cuto

ff_

LJ)

; }

for

(in

t i

= 0

; i <

N; i+

+)

{

if

(e

na

ble

_F

R)

v[

i] *=

de

cay_

v;

if

(e

na

ble

_B

M)

v[

i] +

= d

_V

_B

M_

dr_

time

s_d

elta

_t(

); v[

i] +

= d

v[i];

}

for

(in

t i

= 0

; i <

N; i+

+)

{ r[

i] +

= v

[i] *

de

lta_

t; ip

oin

t n

ew

_cu

be

= ip

oin

t(r[

i]/l_

cuto

ff_

LJ)

;

if

(o

ld_

cub

e[i]

!=

ne

w_

cub

e)

{ h

ash

_re

mo

ve_

ind

ex(

idx_

ha

sh, o

ld_

cub

e[i]

, b

ox_

ref[i],

i);

h

ash

_in

sert

_in

de

x(id

x_h

ash

, n

ew

_cu

be

, b

ox_

ref[i],

i);

o

ld_

cub

e[i]

= n

ew

_cu

be

; }

}

v_

up

da

te(r

, v,

dv,

bo

x_re

f, N

, tim

est

ep

s);

for

(in

t i

= 0

; i <

N; i+

+)

{ v[

i] +

= d

v[i];

}

} do

ub

le e

ne

rgy_

all(

po

int*

r, p

oin

t* v

, n

od

e**

& b

ox_

ref,

int

nu

mb

er)

{

do

ub

le e

ne

rgy

= 0

;

for

(in

t k

= 0

; k<

nu

mb

er;

k+

+)

e

ne

rgy

+=

v[k

]*v[

k]/2

;

if

(e

na

ble

_L

J)

for

(in

t k

= 0

; k<

nu

mb

er;

k+

+)

{

int

nu

mb

er;

int

* n

eig

hb

ou

rs =

ne

igh

bo

urs

_o

f_cu

be

(bo

x_re

f[k]

, n

um

be

r);

for

(in

t j

= 0

; j <

nu

mb

er;

j++

) if

(n

eig

hb

ou

rs[j]

< k

) e

ne

rgy

+=

V_

LJ(

r[k]

−r[

ne

igh

bo

urs

[j]])

;

delete

[]n

eig

hb

ou

rs;

}

if

(e

na

ble

_B

L)

for

(in

t k

= 0

; k<

nu

mb

er−

1; k+

+)

e

ne

rgy

+=

V_

BL

(ab

s(r[

k+1

]−r[

k]))

;

if

(e

na

ble

_B

A)

verl

et.c

c

16/1

7ve

rlet.c

c

for

(in

t k

= 1

; k<

nu

mb

er−

1; k+

+)

e

ne

rgy+

=V

_B

A(r

[k−

1]−

r[k]

,r[k

]−r[

k+1

],a

bs(

r[k−

1]−

r[k]

),a

bs(

r[k+

1]−

r[k]

));

if

(e

na

ble

_S

P)

for

(in

t k

= 0

; k<

nu

mb

er−

1; k+

+)

e

ne

rgy

+=

V_

SP

(r[k

+1

]−r[

k]);

if

(e

na

ble

_tr

ap

)

for

(in

t i

= 0

; i <

nu

mb

er;

i++

) e

ne

rgy

+=

V_

tra

p(r

[i]);

return

en

erg

y;} #endif

verl

et.c

c#include

<fs

tre

am

>#include

<io

stre

am

>#include

<cm

ath

>#include

<st

dio

.h>

#include

"po

int.h

h"#include

"ip

oint

.hh"

#include

"ip

oint

_has

h.hh

"#include

"co

nsta

nts.

hh"

#include

"po

tent

ials

.hh"

#ifndef

VE

RL

ET

_IN

CL

UD

ED

#define

VE

RL

ET

_IN

CL

UD

ED

// T

his

file

incl

ud

es

the

Ve

rle

t−A

lgo

rith

m, a

fu

nct

ion

fo

r in

itia

lisa

tion

of th

e p

oly

me

r// a

nd

a fu

nct

ion

fo

r ca

lcu

latin

g th

e e

ne

rgy

of th

e p

oly

me

r.// // F

un

ctio

ns

use

d:

// −

init_

po

lym

er

: re

ad

s th

e p

oly

me

r−co

nfig

ura

tion

fro

m a

n in

pu

t−fil

e// −

tim

est

ep

: th

is is

th

e v

erle

t−a

lgo

rith

m u

sin

g a

ll e

na

ble

d fo

rce

s (in

clu

din

g

// frict

ion

an

d b

row

nia

n fo

rce

)// −

en

erg

y_a

ll : co

mp

ute

s th

e e

ne

rgy

of th

e p

oly

me

r, in

clu

din

g a

ll e

na

ble

d p

ote

ntia

lsvo

id in

it_p

oly

me

r(p

oin

t*&

r, p

oin

t*&

v, n

od

e**

& b

ox_

ref, h

ash

_ip

oin

t& id

x_tr

ee

, in

t&

nu

mb

er,

ch

ar

* fil

e_

na

me

,

cha

r*

nu

mb

er_

of_

mo

no

me

rs,

int

& o

ld_

time

ste

ps)

;

void

tim

est

ep

(po

int*

& r

, p

oin

t*&

v, n

od

e**

& b

ox_

ref, h

ash

_ip

oin

t& id

x_tr

ee

, in

t N

, in

t tim

est

ep

s);

do

ub

le e

ne

rgy_

all(

po

int*

r, p

oin

t* v

, n

od

e**

& b

ox_

ref,

int

nu

mb

er)

;

#endif

verl

et.h

h

17/1

7ve

rlet.c

c, v

erle

t.hh

List of Figures

1.1 The recombination of double stranded DNA [19] . . . . . . . . . 31.2 Electrophoresis in a gel [19] . . . . . . . . . . . . . . . . . . . . . 51.3 The entropic trap . . . . . . . . . . . . . . . . . . . . . . . . . . . 61.4 Schematic drawing of a ratchet . . . . . . . . . . . . . . . . . . . 71.5 A ratchet device [13] . . . . . . . . . . . . . . . . . . . . . . . . . 7

3.1 The bead-spring model . . . . . . . . . . . . . . . . . . . . . . . . 143.2 The excluded volume . . . . . . . . . . . . . . . . . . . . . . . . . 173.3 The entanglement interaction . . . . . . . . . . . . . . . . . . . . 183.4 The persistence length . . . . . . . . . . . . . . . . . . . . . . . . 19

4.1 The binary search tree . . . . . . . . . . . . . . . . . . . . . . . . 224.2 End-to-end vector plotted against chain length . . . . . . . . . . 264.3 Correlation of mode 1 plotted against time . . . . . . . . . . . . 284.4 Correlation of mode 4 plotted against time . . . . . . . . . . . . 284.5 Correlation of mode 9 plotted against time . . . . . . . . . . . . 29

5.1 Snapshot of the polymer in free solution . . . . . . . . . . . . . . 315.2 Radius of gyration for different chain lengths . . . . . . . . . . . 325.3 Ratio of radius of gyration vs. end-to-end vector . . . . . . . . . 335.4 Movement in free solution . . . . . . . . . . . . . . . . . . . . . . 345.5 Mobility µ in free solution . . . . . . . . . . . . . . . . . . . . . . 355.6 Bond length potential with springs . . . . . . . . . . . . . . . . . 365.7 Two bonds crossing each other . . . . . . . . . . . . . . . . . . . 365.8 Distribution of the bond lengths . . . . . . . . . . . . . . . . . . 375.9 Folding of two neighboring bonds . . . . . . . . . . . . . . . . . . 375.10 Projection of the end-to-end vector on the first bond . . . . . . . 38

6.1 The corners of the walls . . . . . . . . . . . . . . . . . . . . . . . 406.2 The entropic trap . . . . . . . . . . . . . . . . . . . . . . . . . . . 406.3 1000 monomers in the initial trap - side view . . . . . . . . . . . 416.4 1000 monomers in the initial trap - view from above . . . . . . . 416.5 1000 monomers in the revised trap . . . . . . . . . . . . . . . . . 426.6 Movement of the center of mass in the trap . . . . . . . . . . . . 43

87

6.7 Movement of the center of mass in the trap with tilted force . . . 446.8 Trap with tilted walls . . . . . . . . . . . . . . . . . . . . . . . . 456.9 Movement of the center of mass in the trap with tilted walls . . . 46

7.1 Mobility plotted against the tilt of the wall for Ex = 0.01 . . . . 507.2 Mobility plotted against the tilt of the wall for Ex = 0.04 . . . . 517.3 Mobility plotted against the force for twall = 0 . . . . . . . . . . 517.4 Mobility plotted against the force for twall = 0.1 . . . . . . . . . 527.5 Mobility plotted against the force for twall = 0.2 . . . . . . . . . 527.6 Mobility plotted against the chain length for twall = 0. . . . . . . 537.7 Mobility plotted against the chain length for twall = 0.02 . . . . . 537.8 Mobility plotted against the chain length for twall = 0.1 . . . . . 547.9 Mobility plotted against the chain length for twall = 0.2 . . . . . 54

88

Bibliography

[1] R. Dean Astumian ”Thermodynamics and Kinetics of a Brownian Mo-tor”, Science Vol. 276, May 9, 1997, p. 917-922

[2] M. Doi and S. F. Edwards ”The Theory of Polymer Dynamics”, OxfordScience Publications 1986

[3] M. Doi ”Introduction to Polymer Physics”, Oxford University Press,reprint 1997

[4] B. Duenweg, W. Paul: ”Brownian Dynamics Simulation without Gaus-sian Random Numbers”, Int. J. Modern Physics C (Physics and Com-puters) 2, 817 (1991)

[5] D. Frenkel, B. Smit ”Understanding molecular Simulation: From Algo-rithms to Applications”, Academic Press 1996

[6] Homepage of the ”GSL - GNU Scientific Library - GNU Project - FreeSoftware Foundation (FSF)”,

http://www.gnu.org/software/gsl/gsl.html

[7] Le Guillou, J. C., and Zinn-Justin, J., ”Phys. Rev. Lett.” 39, 95 (1977)

[8] J. Han and H. G. Craighead ”Separation of Long DNA Molecules in aMicrofabricated Entropic Trap Array”, Science V. 288, 12 May 2000, p.1026 - 1029

[9] S. Hoffmann ”Conformations of Polymer Chains”, ”Soft Matter - Com-plex Materials on mesoscopic Scales”, Vorlesungsmanuskripte des 33. Fe-rienkurs des Instituts fur Festkorperforschung, Forschungszentrum JulichGmbH, 52425 Julich, Germany, Kapitel B.2

[10] McKenzie, D. S., ”Phys. Rep.” 27C, 35 (1976)

[11] A. Kolb, ”Molekular-Dynamik Untersuchungen steifer und flexibler Ket-ten” , Doktorarbeit am Fachbereich Physik der Johannes Gutenberg-Universitat in Mainz, Mai 1999

89

[12] A. Kopf, B. Dunweg, W. Paul ”Dynamics of polymer ’isotope’ mixtures:Molecular dynamics simulation and Rouse model analysis”, J. Chem.Physics 107(17), November 1997

[13] H. Linke ”Von Damonen und Elektronen”, Physikalische Blatter 56(2000) Nr. 5, p. 45-47

[14] W. Nolting ”Grundkurs: Theoretische Physik”, Band 1: ”KlassischeMechanik”, 3. Auflage 1993, Verlag Zimmermann-Neufang

[15] Science 2000, 288, 2294-2295 and 2304-2307

[16] Gary W. Slater et al. ”Theory of DNA electrophoresis: A look at somecurrent challenges”, Electrophoresis 2000, 21, 3873-3887

[17] M. N. Spiteri, F. Boue, A. Lapp, J. P. Cotton ”Polyelectrolyte Persistencelength in semidilute solution as a function of the ionic stregth”, PhysicaB 234-236 (1997) 303-305

[18] B. Tinland, A. Pluen, J. Sturm and G. Weill ”Persistence Length ofSingle-Stranded DNA”, Macromolecules 1997, 30, 5763-5765

[19] J.-L. Viovy ”Electrophoresis of DNA and other polyelectrolytes: Physicalmechanisms”, Reviews of Modern Physics, Vol. 72, No. 3, July 2000, 813-872

[20] R. Winter and F. Noll, ”Methoden der biophysikalischen Chemie”, Teub-ner 1998

[21] D. Wood ”Data Structures, Algorithms, and Performance”, Addison-Wesley 1993

90

Acknowledgements

At this point, I would like to express my gratitude to some people and organi-zations who supported me during my diploma thesis:

• Prof. Dr. Friederike Schmid for the helpful discussions I had with her andthe hints she gave me during my thesis.

• Dr. Alexandra Ros for having the idea of the topic for this thesis.

• Than-Tu Duong who examined entropic traps in experiments and suppliedme with hints on useful setups.

• Olaf Lenz for introducing me to the technique of the hash map as well asthe Standard Template Library.

• The Condensed Matter Group as well as the Mathematical Physics Groupof the University of Bielefeld for the good atmosphere I experienced whenwriting my thesis.

• My parents who provided me with the opportunity to study physics andsupported me wherever they could.

• Parts of this thesis were computed on a Linux-cluster running the jobqueuing system CONDOR. The Condor Software Program (Condor) wasdeveloped by the Condor Team at the Computer Sciences Departmentof the University of Wisconsin-Madison. All rights, title, and interest inCondor are owned by the Condor Team.

Hiermit versichere ich, daß ich außer den angegebenen Quellen keine weitereLiteratur verwendet und diese Arbeit selbst erstellt habe:

91