ICODE-ADC/lammps
4
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Merge pull request #1772 from ohenrich/user-cgdna

Browse Source 
USER-CGDNA oxRNA2 Implementation
This commit is contained in:
Richard Berger
2019-12-06 12:23:47 -07:00
committed by GitHub
parent b1668f28dd e64b395f5d
commit e68824a392
40 changed files with 8367 additions and 214 deletions
Download Patch File Download Diff File Expand all files Collapse all files

73
examples/USER/cgdna/README
View File

@ -1,8 +1,12 @@
This directory contains example data and input files
and utility scripts for the oxDNA coarse-grained model
for DNA.
as well as utility scripts for the oxDNA/oxDNA2/oxRNA2
coarse-grained model of DNA and RNA.
/******************************************************************************/
/examples/oxDNA/duplex1:
/examples/oxDNA2/duplex1:
/examples/duplex1:
Input, data and log files for a DNA duplex (double-stranded DNA)
consisiting of 5 base pairs. The duplex contains two strands with
complementary base pairs. The topology is
@ -11,7 +15,14 @@ A - C - G - T - A
| | | | |
T - G - C - A - T
/examples/duplex2:
Note that in this example the stacking and hydrogen-bonding interactions
are sequence-averaged (cf. keyword 'seqav' in according pair styles).
/******************************************************************************/
/examples/oxDNA/duplex2:
/examples/oxDNA2/duplex2:
Input, data and log files for a nicked DNA duplex (double-stranded DNA)
consisiting of 8 base pairs. The duplex contains strands with
complementary base pairs, but the backbone on one side is not continuous:
@ -22,16 +33,52 @@ A - C - G - T - A - C - G - T
| | | | | | | |
T - G - C - A T - G - C - A
/examples/duplex3:
This is basically the duplex1 run with sequence-dependent stacking
and hydrogen-bonding strengths enabled and both nucleotide mass and
moment of inertia set to the value of the standalone implementation
of oxDNA (M = I = 1). To achieve this, the masses can be set directly
in the input and data file, whereas the moment of inertia is set via
the diameter of the ellipsoid in the data file and has a value of 3.16227766.
Note that in this example the stacking and hydrogen-bonding interactions
are sequence-averaged (cf. keyword 'seqav' in according pair styles).
/******************************************************************************/
/examples/oxDNA2/duplex3:
This example uses the duplex1 with sequence-dependent stacking and
hydrogen-bonding interactions and both nucleotide mass and
moment of inertia set to the value used in the standalone implementation
of oxDNA (M = I = 1). The masses can be set directly in the input and
data file, whereas the moment of inertia is set via the diameter of the
ellipsoid in the data file and has a value of 3.16227766.
The change of mass and moment of inertia allows direct comparision of
e.g. trajectory data, energies or time-dependent observables on a per-timestep
basis until numerical noise causes deviations at later simulation times.
trajectory data or time-dependent observables on a per-timestep basis.
As mentioned above, the stacking and hydrogen-bonding interactions
are sequence-dependent (cf. keyword 'seqdep' in according pair styles).
/******************************************************************************/
/examples/oxDNA2/unique_bp:
This example uses atom types 1-8 to model a 13 base pair duplex.
The nucleotide types are assigned as follows:
A = 1,5; C = 2,6; G = 3,7; T = 4,8
The topology is
A C G T A C G T A C G T A
1 - 2 - 3 - 4 - 5 - 6 - 7 - 8 - 1 - 2 - 7 - 8 - 1
| | | | | | | | | | | | |
4 - 3 - 2 - 1 - 8 - 7 - 6 - 5 - 4 - 3 - 6 - 5 - 4
T G C A T G C A T G C A T
With a large (32 or 64) number of atom types quasi-unique base pairing
between two individual nucleotides can be established.
/******************************************************************************/
/examples/oxRNA2/duplex4
This example uses the duplex2 with the oxRNA2 force field instead of oxDNA or
oxDNA2 force field. Sequence-dependent stacking and hydrogen-bonding
strengths enabled (cf. keyword 'seqdep' in according pair styles).
/******************************************************************************/
/util:
This directory contains a simple python setup tool which creates

2
examples/USER/cgdna/examples/oxDNA2/duplex3/in.duplex3
View File

@ -1,4 +1,4 @@
variable number equal 1
variable number equal 3
variable ofreq equal 1000
variable efreq equal 1000
variable T equal 0.1

2
examples/USER/cgdna/examples/oxDNA2/duplex3/log.07Aug19.duplex3.g++.1
View File

@ -1,5 +1,5 @@
LAMMPS (7 Aug 2019)
variable number equal 1
variable number equal 3
variable ofreq equal 1000
variable efreq equal 1000
variable T equal 0.1

2
examples/USER/cgdna/examples/oxDNA2/duplex3/log.07Aug19.duplex3.g++.4
View File

@ -1,5 +1,5 @@
LAMMPS (7 Aug 2019)
variable number equal 1
variable number equal 3
variable ofreq equal 1000
variable efreq equal 1000
variable T equal 0.1

130
examples/USER/cgdna/examples/oxDNA2/unique_bp/data.duplex4.4type Normal file
View File

@ -0,0 +1,130 @@
# LAMMPS data file
26 atoms
26 ellipsoids
24 bonds
4 atom types
1 bond types
# System size
-20.000000 20.000000 xlo xhi
-20.000000 20.000000 ylo yhi
-20.000000 20.000000 zlo zhi
# Atom-ID, type, position, molecule-ID, ellipsoid flag, density
Atoms
1 1 -6.000000000000001e-01 0.000000000000000e+00 0.000000000000000e+00 1 1 1
2 2 -4.957432645895970e-01 -3.379920348381733e-01 3.897628551303122e-01 1 1 1
3 3 -2.192046146198370e-01 -5.585242491865227e-01 7.795257102606244e-01 1 1 1
4 4 1.335125603737887e-01 -5.849567473090943e-01 1.169288565390937e+00 1 1 1
5 1 4.398311230978960e-01 -4.081036426625517e-01 1.559051420521249e+00 1 1 1
6 2 5.932984957350773e-01 -8.942535970570469e-02 1.948814275651561e+00 1 1 1
7 3 5.405813207414517e-01 2.603302434705350e-01 2.338577130781873e+00 1 1 1
8 4 3.000000000000002e-01 5.196152422706634e-01 2.728339985912185e+00 1 1 1
9 1 -4.483805615185452e-02 5.983222783087083e-01 3.118102841042497e+00 1 1 1
10 2 -3.740938811152403e-01 4.690988894808181e-01 3.507865696172809e+00 1 1 1
11 3 -5.733436834716847e-01 1.768531046465427e-01 3.897628551303121e+00 1 1 1
12 4 -5.733436834716849e-01 -1.768531046465427e-01 4.287391406433434e+00 1 1 1
13 1 -3.740938811152403e-01 -4.690988894808182e-01 4.677154261563746e+00 1 1 1
14 4 3.740938811152403e-01 4.690988894808182e-01 4.677154261563746e+00 2 1 1
15 1 5.733436834716849e-01 1.768531046465427e-01 4.287391406433434e+00 2 1 1
16 2 5.733436834716849e-01 -1.768531046465426e-01 3.897628551303122e+00 2 1 1
17 3 3.740938811152403e-01 -4.690988894808181e-01 3.507865696172810e+00 2 1 1
18 4 4.483805615185462e-02 -5.983222783087085e-01 3.118102841042498e+00 2 1 1
19 1 -3.000000000000003e-01 -5.196152422706636e-01 2.728339985912186e+00 2 1 1
20 2 -5.405813207414519e-01 -2.603302434705351e-01 2.338577130781874e+00 2 1 1
21 3 -5.932984957350776e-01 8.942535970570474e-02 1.948814275651561e+00 2 1 1
22 4 -4.398311230978962e-01 4.081036426625520e-01 1.559051420521249e+00 2 1 1
23 1 -1.335125603737888e-01 5.849567473090947e-01 1.169288565390937e+00 2 1 1
24 2 2.192046146198373e-01 5.585242491865231e-01 7.795257102606246e-01 2 1 1
25 3 4.957432645895974e-01 3.379920348381736e-01 3.897628551303123e-01 2 1 1
26 4 6.000000000000006e-01 0.000000000000000e+00 1.110223024625157e-16 2 1 1
# Atom-ID, translational, rotational velocity
Velocities
1 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00
2 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00
3 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00
4 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00
5 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00
6 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00
7 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00
8 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00
9 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00
10 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00
11 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00
12 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00
13 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00
14 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00
15 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00
16 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00
17 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00
18 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00
19 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00
20 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00
21 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00
22 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00
23 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00
24 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00
25 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00
26 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00
# Atom-ID, shape, quaternion
Ellipsoids
1 1.173984503142341e+00 1.173984503142341e+00 1.173984503142341e+00 1.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00
2 1.173984503142341e+00 1.173984503142341e+00 1.173984503142341e+00 9.555728057861408e-01 0.000000000000000e+00 0.000000000000000e+00 2.947551744109042e-01
3 1.173984503142341e+00 1.173984503142341e+00 1.173984503142341e+00 8.262387743159949e-01 0.000000000000000e+00 0.000000000000000e+00 5.633200580636221e-01
4 1.173984503142341e+00 1.173984503142341e+00 1.173984503142341e+00 6.234898018587335e-01 0.000000000000000e+00 0.000000000000000e+00 7.818314824680299e-01
5 1.173984503142341e+00 1.173984503142341e+00 1.173984503142341e+00 3.653410243663950e-01 0.000000000000000e+00 0.000000000000000e+00 9.308737486442042e-01
6 1.173984503142341e+00 1.173984503142341e+00 1.173984503142341e+00 7.473009358642424e-02 0.000000000000000e+00 0.000000000000000e+00 9.972037971811802e-01
7 1.173984503142341e+00 1.173984503142341e+00 1.173984503142341e+00 -2.225209339563144e-01 0.000000000000000e+00 0.000000000000000e+00 9.749279121818237e-01
8 1.173984503142341e+00 1.173984503142341e+00 1.173984503142341e+00 -5.000000000000001e-01 0.000000000000000e+00 0.000000000000000e+00 8.660254037844387e-01
9 1.173984503142341e+00 1.173984503142341e+00 1.173984503142341e+00 7.330518718298263e-01 -0.000000000000000e+00 0.000000000000000e+00 -6.801727377709196e-01
10 1.173984503142341e+00 1.173984503142341e+00 1.173984503142341e+00 9.009688679024190e-01 -0.000000000000000e+00 0.000000000000000e+00 -4.338837391175581e-01
11 1.173984503142341e+00 1.173984503142341e+00 1.173984503142341e+00 9.888308262251286e-01 -0.000000000000000e+00 0.000000000000000e+00 -1.490422661761745e-01
12 1.173984503142341e+00 1.173984503142341e+00 1.173984503142341e+00 9.888308262251286e-01 0.000000000000000e+00 0.000000000000000e+00 1.490422661761745e-01
13 1.173984503142341e+00 1.173984503142341e+00 1.173984503142341e+00 9.009688679024190e-01 0.000000000000000e+00 0.000000000000000e+00 4.338837391175582e-01
14 1.173984503142341e+00 1.173984503142341e+00 1.173984503142341e+00 0.000000000000000e+00 -4.338837391175582e-01 9.009688679024190e-01 0.000000000000000e+00
15 1.173984503142341e+00 1.173984503142341e+00 1.173984503142341e+00 -0.000000000000000e+00 -1.490422661761746e-01 9.888308262251286e-01 0.000000000000000e+00
16 1.173984503142341e+00 1.173984503142341e+00 1.173984503142341e+00 -0.000000000000000e+00 1.490422661761745e-01 9.888308262251286e-01 -0.000000000000000e+00
17 1.173984503142341e+00 1.173984503142341e+00 1.173984503142341e+00 -0.000000000000000e+00 4.338837391175581e-01 9.009688679024190e-01 -0.000000000000000e+00
18 1.173984503142341e+00 1.173984503142341e+00 1.173984503142341e+00 0.000000000000000e+00 6.801727377709192e-01 7.330518718298267e-01 0.000000000000000e+00
19 1.173984503142341e+00 1.173984503142341e+00 1.173984503142341e+00 0.000000000000000e+00 8.660254037844386e-01 5.000000000000001e-01 0.000000000000000e+00
20 1.173984503142341e+00 1.173984503142341e+00 1.173984503142341e+00 0.000000000000000e+00 9.749279121818235e-01 2.225209339563145e-01 0.000000000000000e+00
21 1.173984503142341e+00 1.173984503142341e+00 1.173984503142341e+00 0.000000000000000e+00 9.972037971811801e-01 -7.473009358642428e-02 0.000000000000000e+00
22 1.173984503142341e+00 1.173984503142341e+00 1.173984503142341e+00 0.000000000000000e+00 9.308737486442041e-01 -3.653410243663952e-01 0.000000000000000e+00
23 1.173984503142341e+00 1.173984503142341e+00 1.173984503142341e+00 0.000000000000000e+00 7.818314824680296e-01 -6.234898018587339e-01 0.000000000000000e+00
24 1.173984503142341e+00 1.173984503142341e+00 1.173984503142341e+00 0.000000000000000e+00 5.633200580636215e-01 -8.262387743159952e-01 0.000000000000000e+00
25 1.173984503142341e+00 1.173984503142341e+00 1.173984503142341e+00 -0.000000000000000e+00 -2.947551744109044e-01 9.555728057861407e-01 0.000000000000000e+00
26 1.173984503142341e+00 1.173984503142341e+00 1.173984503142341e+00 -0.000000000000000e+00 0.000000000000000e+00 1.000000000000000e+00 -0.000000000000000e+00
# Bond topology
Bonds
1 1 1 2
2 1 2 3
3 1 3 4
4 1 4 5
5 1 5 6
6 1 6 7
7 1 7 8
8 1 8 9
9 1 9 10
10 1 10 11
11 1 11 12
12 1 12 13
13 1 14 15
14 1 15 16
15 1 16 17
16 1 17 18
17 1 18 19
18 1 19 20
19 1 20 21
20 1 21 22
21 1 22 23
22 1 23 24
23 1 24 25
24 1 25 26

130
examples/USER/cgdna/examples/oxDNA2/unique_bp/data.duplex4.8type Normal file
View File

@ -0,0 +1,130 @@
# LAMMPS data file
26 atoms
26 ellipsoids
24 bonds
8 atom types
1 bond types
# System size
-20.000000 20.000000 xlo xhi
-20.000000 20.000000 ylo yhi
-20.000000 20.000000 zlo zhi
# Atom-ID, type, position, molecule-ID, ellipsoid flag, density
Atoms
1 1 -6.000000000000001e-01 0.000000000000000e+00 0.000000000000000e+00 1 1 1
2 2 -4.957432645895970e-01 -3.379920348381733e-01 3.897628551303122e-01 1 1 1
3 3 -2.192046146198370e-01 -5.585242491865227e-01 7.795257102606244e-01 1 1 1
4 4 1.335125603737887e-01 -5.849567473090943e-01 1.169288565390937e+00 1 1 1
5 5 4.398311230978960e-01 -4.081036426625517e-01 1.559051420521249e+00 1 1 1
6 6 5.932984957350773e-01 -8.942535970570469e-02 1.948814275651561e+00 1 1 1
7 7 5.405813207414517e-01 2.603302434705350e-01 2.338577130781873e+00 1 1 1
8 8 3.000000000000002e-01 5.196152422706634e-01 2.728339985912185e+00 1 1 1
9 1 -4.483805615185452e-02 5.983222783087083e-01 3.118102841042497e+00 1 1 1
10 2 -3.740938811152403e-01 4.690988894808181e-01 3.507865696172809e+00 1 1 1
11 7 -5.733436834716847e-01 1.768531046465427e-01 3.897628551303121e+00 1 1 1
12 8 -5.733436834716849e-01 -1.768531046465427e-01 4.287391406433434e+00 1 1 1
13 1 -3.740938811152403e-01 -4.690988894808182e-01 4.677154261563746e+00 1 1 1
14 4 3.740938811152403e-01 4.690988894808182e-01 4.677154261563746e+00 2 1 1
15 5 5.733436834716849e-01 1.768531046465427e-01 4.287391406433434e+00 2 1 1
16 6 5.733436834716849e-01 -1.768531046465426e-01 3.897628551303122e+00 2 1 1
17 3 3.740938811152403e-01 -4.690988894808181e-01 3.507865696172810e+00 2 1 1
18 4 4.483805615185462e-02 -5.983222783087085e-01 3.118102841042498e+00 2 1 1
19 5 -3.000000000000003e-01 -5.196152422706636e-01 2.728339985912186e+00 2 1 1
20 6 -5.405813207414519e-01 -2.603302434705351e-01 2.338577130781874e+00 2 1 1
21 7 -5.932984957350776e-01 8.942535970570474e-02 1.948814275651561e+00 2 1 1
22 8 -4.398311230978962e-01 4.081036426625520e-01 1.559051420521249e+00 2 1 1
23 1 -1.335125603737888e-01 5.849567473090947e-01 1.169288565390937e+00 2 1 1
24 2 2.192046146198373e-01 5.585242491865231e-01 7.795257102606246e-01 2 1 1
25 3 4.957432645895974e-01 3.379920348381736e-01 3.897628551303123e-01 2 1 1
26 4 6.000000000000006e-01 0.000000000000000e+00 1.110223024625157e-16 2 1 1
# Atom-ID, translational, rotational velocity
Velocities
1 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00
2 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00
3 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00
4 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00
5 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00
6 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00
7 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00
8 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00
9 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00
10 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00
11 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00
12 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00
13 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00
14 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00
15 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00
16 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00
17 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00
18 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00
19 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00
20 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00
21 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00
22 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00
23 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00
24 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00
25 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00
26 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00
# Atom-ID, shape, quaternion
Ellipsoids
1 1.173984503142341e+00 1.173984503142341e+00 1.173984503142341e+00 1.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00
2 1.173984503142341e+00 1.173984503142341e+00 1.173984503142341e+00 9.555728057861408e-01 0.000000000000000e+00 0.000000000000000e+00 2.947551744109042e-01
3 1.173984503142341e+00 1.173984503142341e+00 1.173984503142341e+00 8.262387743159949e-01 0.000000000000000e+00 0.000000000000000e+00 5.633200580636221e-01
4 1.173984503142341e+00 1.173984503142341e+00 1.173984503142341e+00 6.234898018587335e-01 0.000000000000000e+00 0.000000000000000e+00 7.818314824680299e-01
5 1.173984503142341e+00 1.173984503142341e+00 1.173984503142341e+00 3.653410243663950e-01 0.000000000000000e+00 0.000000000000000e+00 9.308737486442042e-01
6 1.173984503142341e+00 1.173984503142341e+00 1.173984503142341e+00 7.473009358642424e-02 0.000000000000000e+00 0.000000000000000e+00 9.972037971811802e-01
7 1.173984503142341e+00 1.173984503142341e+00 1.173984503142341e+00 -2.225209339563144e-01 0.000000000000000e+00 0.000000000000000e+00 9.749279121818237e-01
8 1.173984503142341e+00 1.173984503142341e+00 1.173984503142341e+00 -5.000000000000001e-01 0.000000000000000e+00 0.000000000000000e+00 8.660254037844387e-01
9 1.173984503142341e+00 1.173984503142341e+00 1.173984503142341e+00 7.330518718298263e-01 -0.000000000000000e+00 0.000000000000000e+00 -6.801727377709196e-01
10 1.173984503142341e+00 1.173984503142341e+00 1.173984503142341e+00 9.009688679024190e-01 -0.000000000000000e+00 0.000000000000000e+00 -4.338837391175581e-01
11 1.173984503142341e+00 1.173984503142341e+00 1.173984503142341e+00 9.888308262251286e-01 -0.000000000000000e+00 0.000000000000000e+00 -1.490422661761745e-01
12 1.173984503142341e+00 1.173984503142341e+00 1.173984503142341e+00 9.888308262251286e-01 0.000000000000000e+00 0.000000000000000e+00 1.490422661761745e-01
13 1.173984503142341e+00 1.173984503142341e+00 1.173984503142341e+00 9.009688679024190e-01 0.000000000000000e+00 0.000000000000000e+00 4.338837391175582e-01
14 1.173984503142341e+00 1.173984503142341e+00 1.173984503142341e+00 0.000000000000000e+00 -4.338837391175582e-01 9.009688679024190e-01 0.000000000000000e+00
15 1.173984503142341e+00 1.173984503142341e+00 1.173984503142341e+00 -0.000000000000000e+00 -1.490422661761746e-01 9.888308262251286e-01 0.000000000000000e+00
16 1.173984503142341e+00 1.173984503142341e+00 1.173984503142341e+00 -0.000000000000000e+00 1.490422661761745e-01 9.888308262251286e-01 -0.000000000000000e+00
17 1.173984503142341e+00 1.173984503142341e+00 1.173984503142341e+00 -0.000000000000000e+00 4.338837391175581e-01 9.009688679024190e-01 -0.000000000000000e+00
18 1.173984503142341e+00 1.173984503142341e+00 1.173984503142341e+00 0.000000000000000e+00 6.801727377709192e-01 7.330518718298267e-01 0.000000000000000e+00
19 1.173984503142341e+00 1.173984503142341e+00 1.173984503142341e+00 0.000000000000000e+00 8.660254037844386e-01 5.000000000000001e-01 0.000000000000000e+00
20 1.173984503142341e+00 1.173984503142341e+00 1.173984503142341e+00 0.000000000000000e+00 9.749279121818235e-01 2.225209339563145e-01 0.000000000000000e+00
21 1.173984503142341e+00 1.173984503142341e+00 1.173984503142341e+00 0.000000000000000e+00 9.972037971811801e-01 -7.473009358642428e-02 0.000000000000000e+00
22 1.173984503142341e+00 1.173984503142341e+00 1.173984503142341e+00 0.000000000000000e+00 9.308737486442041e-01 -3.653410243663952e-01 0.000000000000000e+00
23 1.173984503142341e+00 1.173984503142341e+00 1.173984503142341e+00 0.000000000000000e+00 7.818314824680296e-01 -6.234898018587339e-01 0.000000000000000e+00
24 1.173984503142341e+00 1.173984503142341e+00 1.173984503142341e+00 0.000000000000000e+00 5.633200580636215e-01 -8.262387743159952e-01 0.000000000000000e+00
25 1.173984503142341e+00 1.173984503142341e+00 1.173984503142341e+00 -0.000000000000000e+00 -2.947551744109044e-01 9.555728057861407e-01 0.000000000000000e+00
26 1.173984503142341e+00 1.173984503142341e+00 1.173984503142341e+00 -0.000000000000000e+00 0.000000000000000e+00 1.000000000000000e+00 -0.000000000000000e+00
# Bond topology
Bonds
1 1 1 2
2 1 2 3
3 1 3 4
4 1 4 5
5 1 5 6
6 1 6 7
7 1 7 8
8 1 8 9
9 1 9 10
10 1 10 11
11 1 11 12
12 1 12 13
13 1 14 15
14 1 15 16
15 1 16 17
16 1 17 18
17 1 18 19
18 1 19 20
19 1 20 21
20 1 21 22
21 1 22 23
22 1 23 24
23 1 24 25
24 1 25 26

828
examples/USER/cgdna/examples/oxDNA2/unique_bp/generate_unique.py Normal file
View File

@ -0,0 +1,828 @@
#!/usr/bin/env python
"""
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
Contributing author: Oliver Henrich (University of Strathclyde, Glasgow)
------------------------------------------------------------------------- */
"""
"""
Creates unique base-pairings to avoid asymmetrical H-bonds.
Modified to create the bead wall setup.
N_BEADS is the number of beads along one direction, the final system will have N_BEADS^2 beads in the wall. N_BEADS should be set to be odd number.
"""
#Define number of base-pairs per turn for B-form DNA
N = 10.5
#Define distance between the big bead and the centre of mass of the last base-pair
BEAD_OFFSET = 2.0
WALL_PARTICLE_SIZE = 2.0
N_BEADS = 11
#Number of unique base type groups (1-4) ACGT counts as one group
N_BASE_TYPES = 20
"""
Import basic modules
"""
import sys, os, timeit
from timeit import default_timer as timer
start_time = timer()
"""
Try to import numpy; if failed, import a local version mynumpy
which needs to be provided
"""
try:
import numpy as np
except:
print("numpy not found. Exiting.", file=sys.stderr)
sys.exit(1)
"""
Check that the required arguments (box offset and size in simulation units
and the sequence file were provided
"""
try:
box_offset = float(sys.argv[1])
box_length = float(sys.argv[2])
infile = sys.argv[3]
if len(sys.argv) == 4:
topo = 'strand'
lk = 0
elif len(sys.argv) == 5:
topo = 'strand'
lk = int(sys.argv[4])
except:
print("Usage: %s <%s> <%s> <%s> <%s> " % (sys.argv[0], \
"box offset", "box length", "file with sequences", "[Lk]"), file=sys.stderr)
sys.exit(1)
box = np.array ([box_length, box_length, box_length])
"""
Try to open the file and fail gracefully if file cannot be opened
"""
try:
inp = open (infile, 'r')
inp.close()
except:
print("Could not open file '%s' for reading. \
Aborting." % infile, file=sys.stderr)
sys.exit(2)
# return parts of a string
def partition(s, d):
if d in s:
sp = s.split(d, 1)
return sp[0], d, sp[1]
else:
return s, "", ""
"""
Define the model constants
"""
# set model constants
PI = np.pi
POS_BASE = 0.4
POS_BACK = -0.4
EXCL_RC1 = 0.711879214356
EXCL_RC2 = 0.335388426126
EXCL_RC3 = 0.52329943261
"""
Define auxillary variables for the construction of a helix
"""
# center of the double strand
COM_CENTRE_DS = POS_BASE + 0.2
# ideal rise between two consecutive nucleotides on the
# same strand which are to be base paired in a duplex
BASE_BASE = 0.3897628551303122
# cutoff distance for overlap check
RC2 = 16
# squares of the excluded volume distances for overlap check
RC2_BACK = EXCL_RC1**2
RC2_BASE = EXCL_RC2**2
RC2_BACK_BASE = EXCL_RC3**2
# enumeration to translate from letters to numbers and vice versa
number_to_base = {1 : 'A', 2 : 'C', 3 : 'G', 4 : 'T'}
base_to_number = {'A' : 1, 'a' : 1, 'C' : 2, 'c' : 2,
'G' : 3, 'g' : 3, 'T' : 4, 't' : 4}
# auxillary arrays
positions = []
a1s = []
a3s = []
quaternions = []
newpositions = []
newa1s = []
newa3s = []
basetype = []
strandnum = []
bonds = []
"""
Convert local body frame to quaternion DOF
"""
def exyz_to_quat (mya1, mya3):
mya2 = np.cross(mya3, mya1)
myquat = [1,0,0,0]
q0sq = 0.25 * (mya1[0] + mya2[1] + mya3[2] + 1.0)
q1sq = q0sq - 0.5 * (mya2[1] + mya3[2])
q2sq = q0sq - 0.5 * (mya1[0] + mya3[2])
q3sq = q0sq - 0.5 * (mya1[0] + mya2[1])
# some component must be greater than 1/4 since they sum to 1
# compute other components from it
if q0sq >= 0.25:
myquat[0] = np.sqrt(q0sq)
myquat[1] = (mya2[2] - mya3[1]) / (4.0*myquat[0])
myquat[2] = (mya3[0] - mya1[2]) / (4.0*myquat[0])
myquat[3] = (mya1[1] - mya2[0]) / (4.0*myquat[0])
elif q1sq >= 0.25:
myquat[1] = np.sqrt(q1sq)
myquat[0] = (mya2[2] - mya3[1]) / (4.0*myquat[1])
myquat[2] = (mya2[0] + mya1[1]) / (4.0*myquat[1])
myquat[3] = (mya1[2] + mya3[0]) / (4.0*myquat[1])
elif q2sq >= 0.25:
myquat[2] = np.sqrt(q2sq)
myquat[0] = (mya3[0] - mya1[2]) / (4.0*myquat[2])
myquat[1] = (mya2[0] + mya1[1]) / (4.0*myquat[2])
myquat[3] = (mya3[1] + mya2[2]) / (4.0*myquat[2])
elif q3sq >= 0.25:
myquat[3] = np.sqrt(q3sq)
myquat[0] = (mya1[1] - mya2[0]) / (4.0*myquat[3])
myquat[1] = (mya3[0] + mya1[2]) / (4.0*myquat[3])
myquat[2] = (mya3[1] + mya2[2]) / (4.0*myquat[3])
norm = 1.0/np.sqrt(myquat[0]*myquat[0] + myquat[1]*myquat[1] + \
myquat[2]*myquat[2] + myquat[3]*myquat[3])
myquat[0] *= norm
myquat[1] *= norm
myquat[2] *= norm
myquat[3] *= norm
return np.array([myquat[0],myquat[1],myquat[2],myquat[3]])
"""
Adds a strand to the system by appending it to the array of previous strands
"""
def add_strands (mynewpositions, mynewa1s, mynewa3s):
overlap = False
# This is a simple check for each of the particles where for previously
# placed particles i we check whether it overlaps with any of the
# newly created particles j
print("## Checking for overlaps", file=sys.stdout)
for i in range(len(positions)):
p = positions[i]
pa1 = a1s[i]
for j in range (len(mynewpositions)):
q = mynewpositions[j]
qa1 = mynewa1s[j]
# skip particles that are anyway too far away
dr = p - q
dr -= box * np.rint (dr / box)
if np.dot(dr, dr) > RC2:
continue
# base site and backbone site of the two particles
p_pos_back = p + pa1 * POS_BACK
p_pos_base = p + pa1 * POS_BASE
q_pos_back = q + qa1 * POS_BACK
q_pos_base = q + qa1 * POS_BASE
# check for no overlap between the two backbone sites
dr = p_pos_back - q_pos_back
dr -= box * np.rint (dr / box)
if np.dot(dr, dr) < RC2_BACK:
overlap = True
# check for no overlap between the two base sites
dr = p_pos_base - q_pos_base
dr -= box * np.rint (dr / box)
if np.dot(dr, dr) < RC2_BASE:
overlap = True
# check for no overlap between backbone site of particle p
# with base site of particle q
dr = p_pos_back - q_pos_base
dr -= box * np.rint (dr / box)
if np.dot(dr, dr) < RC2_BACK_BASE:
overlap = True
# check for no overlap between base site of particle p and
# backbone site of particle q
dr = p_pos_base - q_pos_back
dr -= box * np.rint (dr / box)
if np.dot(dr, dr) < RC2_BACK_BASE:
overlap = True
# exit if there is an overlap
if overlap:
return False
# append to the existing list if no overlap is found
if not overlap:
for p in mynewpositions:
positions.append(p)
for p in mynewa1s:
a1s.append (p)
for p in mynewa3s:
a3s.append (p)
# calculate quaternion from local body frame and append
for ia in range(len(mynewpositions)):
mynewquaternions = exyz_to_quat(mynewa1s[ia],mynewa3s[ia])
quaternions.append(mynewquaternions)
return True
"""
Calculate angle of rotation site to site
"""
def get_angle(bp):
#n, minimal number of bases per turn
n = 10.5
found = False
while found == False:
turns = bp/n
diff = abs( turns - round(turns))
if diff < 0.03:
found = True
turns = round(turns)+lk
angle = (360*turns)/bp
angle = round (angle,2)
#angle =round( 360/n,2)
elif n > 11.5:
angle = 35.9
found = True
else:
n += 0.02
return angle
def get_angle2(bp):
turns = bp/N + lk
angle = (360*turns)/bp
return angle
"""
Returns the rotation matrix defined by an axis and angle
"""
def get_rotation_matrix(axis, anglest, nbp=0):
# The argument anglest can be either an angle in radiants
# (accepted types are float, int or np.float64 or np.float64)
# or a tuple [angle, units] where angle is a number and
# units is a string. It tells the routine whether to use degrees,
# radiants (the default) or base pairs turns.
if not isinstance (anglest, (np.float64, np.float32, float, int)):
if len(anglest) > 1:
if anglest[1] in ["degrees", "deg", "o"]:
angle = (np.pi / 180.) * (anglest[0])
elif anglest[1] in ["bp"]:
if nbp == 0:
angle = int(anglest[0]) * (np.pi / 180.) * (35.9)
else:
ang = get_angle2(nbp)
angle = int(anglest[0]) * (np.pi / 180.) * (ang)
else:
angle = float(anglest[0])
else:
angle = float(anglest[0])
else:
angle = float(anglest) # in degrees (?)
axis = np.array(axis)
axis /= np.sqrt(np.dot(axis, axis))
ct = np.cos(angle)
st = np.sin(angle)
olc = 1. - ct
x, y, z = axis
return np.array([[olc*x*x+ct, olc*x*y-st*z, olc*x*z+st*y],
[olc*x*y+st*z, olc*y*y+ct, olc*y*z-st*x],
[olc*x*z-st*y, olc*y*z+st*x, olc*z*z+ct]])
"""
Generates the position and orientation vectors of a
(single or double) strand from a sequence string
"""
def generate_strand(bp, sequence=None, start_pos=np.array([0, 0, 0]), \
dir=np.array([0, 0, 1]), perp=False, double=True, rot=0.):
# generate empty arrays
mynewpositions, mynewa1s, mynewa3s = [], [], []
# cast the provided start_pos array into a numpy array
start_pos = np.array(start_pos, dtype=float)
# overall direction of the helix
dir = np.array(dir, dtype=float)
#if sequence == None:
# sequence = np.random.randint(1, 5, bp)
# the elseif here is most likely redundant
#elif len(sequence) != bp:
# n = bp - len(sequence)
# sequence += np.random.randint(1, 5, n)
# print("sequence is too short, adding %d random bases" % n, file=sys.stderr)
# normalize direction
dir_norm = np.sqrt(np.dot(dir,dir))
if dir_norm < 1e-10:
print("direction must be a valid vector,\
defaulting to (0, 0, 1)", file=sys.stderr)
dir = np.array([0, 0, 1])
else: dir /= dir_norm
# find a vector orthogonal to dir to act as helix direction,
# if not provided switch off random orientation
if perp is None or perp is False:
v1 = np.random.random_sample(3)
# comment in to suppress randomised base vector
v1 = [1,0,0]
v1 -= dir * (np.dot(dir, v1))
v1 /= np.sqrt(sum(v1*v1))
else:
v1 = perp;
# generate rotational matrix representing the overall rotation of the helix
R0 = get_rotation_matrix(dir, rot)
# rotation matrix corresponding to one step along the helix
R = get_rotation_matrix(dir, [1, "bp"],bp)
# set the vector a1 (backbone to base) to v1
a1 = v1
# apply the global rotation to a1
a1 = np.dot(R0, a1)
# set the position of the fist backbone site to start_pos
rb = np.array(start_pos)
# set a3 to the direction of the helix
a3 = dir
for i in range(bp):
# work out the position of the centre of mass of the nucleotide
rcom = rb - COM_CENTRE_DS * a1
# append to newpositions
mynewpositions.append(rcom)
mynewa1s.append(a1)
mynewa3s.append(a3)
# if we are not at the end of the helix, we work out a1 and rb for the
# next nucleotide along the helix
if i != bp - 1:
a1 = np.dot(R, a1)
rb += a3 * BASE_BASE
# if we are working on a double strand, we do a cycle similar
# to the previous one but backwards
if double == True:
a1 = -a1
a3 = -dir
R = R.transpose()
for i in range(bp):
rcom = rb - COM_CENTRE_DS * a1
mynewpositions.append (rcom)
mynewa1s.append (a1)
mynewa3s.append (a3)
a1 = np.dot(R, a1)
rb += a3 * BASE_BASE
#Calculate the positions of the bead wall
last_base1 = mynewpositions[int( len(mynewpositions)/2 - 1) ]
last_base2 = mynewpositions[int( len(mynewpositions)/2) ]
mid_point = (last_base1 + last_base2) / 2
NN = N_BEADS**2
p1 = [mid_point[0] - (N_BEADS-1)*WALL_PARTICLE_SIZE, mid_point[1] - (N_BEADS-1)*WALL_PARTICLE_SIZE, mid_point[2] + BEAD_OFFSET ]
for i in range(N_BEADS):
for j in range(N_BEADS):
position = [ p1[0] + 2*i*WALL_PARTICLE_SIZE, p1[1] + 2*j*WALL_PARTICLE_SIZE, p1[2]]
mynewa1s.append([1,0,0])
mynewa3s.append([1,0,0])
mynewpositions.append(position)
assert (len (mynewpositions) > 0)
return [mynewpositions, mynewa1s, mynewa3s]
"""
Main function for this script.
Reads a text file with the following format:
- Each line contains the sequence for a single strand (A,C,G,T)
- Lines beginning with the keyword 'DOUBLE' produce double-stranded DNA
Ex: Two ssDNA (single stranded DNA)
ATATATA
GCGCGCG
Ex: Two strands, one double stranded, the other single stranded.
DOUBLE AGGGCT
CCTGTA
"""
def read_strands(filename):
try:
infile = open (filename)
except:
print("Could not open file '%s'. Aborting." % filename, file=sys.stderr)
sys.exit(2)
# This block works out the number of nucleotides and strands by reading
# the number of non-empty lines in the input file and the number of letters,
# taking the possible DOUBLE keyword into account.
nstrands, nnucl, nbonds = 0, 0, 0
lines = infile.readlines()
for line in lines:
line = line.upper().strip()
if len(line) == 0:
continue
if line[:6] == 'DOUBLE':
line = line.split()[1]
length = len(line)
print("## Found duplex of %i base pairs" % length, file=sys.stdout)
nnucl += 2*length
nstrands += 2
nbonds+= 2*length
else:
line = line.split()[0]
length = len(line)
print("## Found single strand of %i bases" % length, file=sys.stdout)
nnucl += length
nstrands += 1
if topo == 'ring':
nbonds =+ length
else:
nbonds += length+1
# rewind the sequence input file
infile.seek(0)
print("## nstrands, nnucl = ", nstrands, nnucl, file=sys.stdout)
# generate the data file in LAMMPS format
try:
out = open ("data.oxdna", "w")
except:
print("Could not open data file for writing. Aborting.", file=sys.stderr)
sys.exit(2)
lines = infile.readlines()
nlines = len(lines)
i = 1
myns = 0
noffset = 1
for line in lines:
line = line.upper().strip()
# skip empty lines
if len(line) == 0:
i += 1
continue
# block for duplexes: last argument of the generate function
# is set to 'True'
if line[:6] == 'DOUBLE':
line = line.split()[1]
length = len(line)
seq = [(base_to_number[x]) for x in line]
seq = np.array(seq,dtype=int)
n_a, n_c, n_g, n_t = 0, 0, 0, 0
for s in range(seq.size):
if seq[s] == 1:
n_a += 1
elif seq[s] == 2:
n_c += 1
elif seq[s] ==3:
n_g += 1
elif seq[s] == 4:
n_t += 1
smallest_n_bases = n_c
if n_a < n_c:
smallest_n_bases = n_a
if smallest_n_bases > n_t:
smallest_n_bases = n_t
if smallest_n_bases > n_g:
smallest_n_bases = n_g
if smallest_n_bases < N_BASE_TYPES:
print('## Not enough occurances of base types in the sequence for ' + str(N_BASE_TYPES))
print('## unique base types, switching to ' + str(smallest_n_bases) + ' unique types')
else:
smallest_n_bases = N_BASE_TYPES
a, c, g, t = -3, -2, -1, 0
for s in range(seq.size):
if seq[s] == 1:
if a < (smallest_n_bases*4-3):
a += 4
else:
a = 1
seq[s] = a
elif seq[s] == 2:
if c < (smallest_n_bases*4-2):
c += 4
else:
c = 2
seq[s] = c
elif seq[s] == 3:
if g < (smallest_n_bases*4-1):
g += 4
else:
g = 3
seq[s] = g
elif seq[s] == 4:
if t < (smallest_n_bases*4):
t += 4
else:
t = 4
seq[s] = t
myns += 1
for b in range(length):
basetype.append(seq[b])
strandnum.append(myns)
for b in range(length-1):
bondpair = [noffset + b, noffset + b + 1]
bonds.append(bondpair)
noffset += length
# create the sequence of the second strand as made of
# complementary bases
#seq2 = [5-s for s in seq]
seq2 = seq
for s in range(seq2.size):
if seq2[s]%4 == 1:
seq2[s] += 3
elif seq2[s]%4 == 2:
seq2[s] += 1
elif seq2[s]%4 == 3:
seq2[s] -= 1
elif seq2[s]%4 == 0:
seq2[s] -= 3
#seq2.reverse()
myns += 1
for b in range(length):
basetype.append(seq2[b])
strandnum.append(myns)
for b in range(length-1):
bondpair = [noffset + b, noffset + b + 1]
bonds.append(bondpair)
#create wall bead types
bead_type = 4*smallest_n_bases + 1
for i in range(N_BEADS**2):
basetype.append(bead_type)
basetype.append(bead_type)
strandnum.append(bead_type)
strandnum.append(bead_type)
#bonds.append([length, noffset + length])
#bonds.append([length+1, noffset + length])
noffset += length
print("## Created duplex of %i bases" % (2*length), file=sys.stdout)
# generate random position of the first nucleotide
com = box_offset + np.random.random_sample(3) * box
# comment out to randomise
com = [0,0,0]
# generate the random direction of the helix
axis = np.random.random_sample(3)
# comment out to randomise
axis = [0,0,1]
axis /= np.sqrt(np.dot(axis, axis))
# use the generate function defined above to create
# the position and orientation vector of the strand
if topo == 'ring':
newpositions, newa1s, newa3s = generate_ring(len(line), \
sequence=seq, dir=axis, start_pos=com, double=True)
else:
newpositions, newa1s, newa3s = generate_strand(len(line), \
sequence=seq, dir=axis, start_pos=com, double=True)
# generate a new position for the strand until it does not overlap
# with anything already present
start = timer()
while not add_strands(newpositions, newa1s, newa3s):
com = box_offset + np.random.random_sample(3) * box
axis = np.random.random_sample(3)
axis /= np.sqrt(np.dot(axis, axis))
if topo == 'ring':
newpositions, newa1s, newa3s = generate_ring(len(line), \
sequence=seq, dir=axis, start_pos=com, double=True)
else:
newpositions, newa1s, newa3s = generate_strand(len(line), \
sequence=seq, dir=axis, start_pos=com, double=True)
print("## Trying %i" % i, file=sys.stdout)
end = timer()
print("## Added duplex of %i bases (line %i/%i) in %.2fs, now at %i/%i" % \
(2*length, i, nlines, end-start, len(positions), nnucl), file=sys.stdout)
# block for single strands: last argument of the generate function
# is set to 'False'
else:
length = len(line)
seq = [(base_to_number[x]) for x in line]
myns += 1
for b in range(length):
basetype.append(seq[b])
strandnum.append(myns)
for b in range(length-1):
bondpair = [noffset + b, noffset + b + 1]
bonds.append(bondpair)
if topo == 'ring':
bondpair = [noffset, noffset + length-1]
bonds.append(bondpair)
noffset += length
# generate random position of the first nucleotide
com = box_offset + np.random.random_sample(3) * box
# comment out to randomise
com = [-30,0,0]
# generate the random direction of the helix
axis = np.random.random_sample(3)
# comment out to randomise
axis = [0,0,1]
axis /= np.sqrt(np.dot(axis, axis))
print("## Created single strand of %i bases" % length, file=sys.stdout)
if topo == 'ring':
newpositions, newa1s, newa3s = generate_ring(length, \
sequence=seq, dir=axis, start_pos=com, double=False)
else:
newpositions, newa1s, newa3s = generate_strand(length, \
sequence=seq, dir=axis, start_pos=com, double=False)
start = timer()
while not add_strands(newpositions, newa1s, newa3s):
com = box_offset + np.random.random_sample(3) * box
axis = np.random.random_sample(3)
axis /= np.sqrt(np.dot(axis, axis))
if topo == 'ring':
newpositions, newa1s, newa3s = generate_ring(length, \
sequence=seq, dir=axis, start_pos=com, double=False)
else:
newpositions, newa1s, newa3s = generate_strand(length, \
sequence=seq, dir=axis, start_pos=com, double=False)
print("## Trying %i" % (i), file=sys.stdout)
end = timer()
print("## Added single strand of %i bases (line %i/%i) in %.2fs, now at %i/%i" % \
(length, i, nlines, end-start,len(positions), nnucl), file=sys.stdout)
i += 1
# sanity check
#if not len(positions) == nnucl:
# print(len(positions), nnucl)
# raise AssertionError
nnucl = nnucl + (N_BEADS**2)
nbonds -= 4
out.write('# LAMMPS data file\n')
out.write('%d atoms\n' % nnucl)
out.write('%d ellipsoids\n' % nnucl)
out.write('%d bonds\n' % nbonds)
out.write('\n')
out.write('%d atom types\n' %bead_type )
out.write('1 bond types\n')
out.write('\n')
out.write('# System size\n')
out.write('%f %f xlo xhi\n' % (box_offset,box_offset+box_length))
out.write('%f %f ylo yhi\n' % (box_offset,box_offset+box_length))
out.write('%f %f zlo zhi\n' % (0,box_length))
#out.write('\n')
#out.write('Masses\n')
#out.write('\n')
#out.write('1 3.1575\n')
#out.write('2 3.1575\n')
#out.write('3 3.1575\n')
#out.write('4 3.1575\n')
#out.write('5 3.1575\n')
# for each nucleotide print a line under the headers
# Atoms, Velocities, Ellipsoids and Bonds
out.write('\n')
out.write(\
'# Atom-ID, type, position, molecule-ID, ellipsoid flag, density\n')
out.write('Atoms\n')
out.write('\n')
for i in range(nnucl):
out.write('%d %d %22.15le %22.15le %22.15le %d 1 1\n' \
% (i+1, basetype[i], \
positions[i][0], positions[i][1], positions[i][2], \
strandnum[i]))
out.write('\n')
out.write('# Atom-ID, translational, rotational velocity\n')
out.write('Velocities\n')
out.write('\n')
for i in range(nnucl):
out.write("%d %22.15le %22.15le %22.15le %22.15le %22.15le %22.15le\n" \
% (i+1,0.0,0.0,0.0,0.0,0.0,0.0))
out.write('\n')
out.write('# Atom-ID, shape, quaternion\n')
out.write('Ellipsoids\n')
out.write('\n')
for i in range(nnucl):
out.write(\
"%d %22.15le %22.15le %22.15le %22.15le %22.15le %22.15le %22.15le\n" \
% (i+1,1.1739845031423408,1.1739845031423408,1.1739845031423408, \
quaternions[i][0],quaternions[i][1], quaternions[i][2],quaternions[i][3]))
out.write('\n')
out.write('# Bond topology\n')
out.write('Bonds\n')
out.write('\n')
for i in range(nbonds):
if i < nbonds-2:
out.write("%d %d %d %d\n" % (i+1,1,bonds[i][0],bonds[i][1]))
#else:
#out.write("%d %d %d %d\n" % (i+1,2,bonds[i][0],bonds[i][1]))
out.close()
print("## Wrote data to 'data.oxdna'", file=sys.stdout)
print("## DONE", file=sys.stdout)
# call the above main() function, which executes the program
read_strands (infile)
end_time=timer()
runtime = end_time-start_time
hours = runtime/3600
minutes = (runtime-np.rint(hours)*3600)/60
seconds = (runtime-np.rint(hours)*3600-np.rint(minutes)*60)%60
print("## Total runtime %ih:%im:%.2fs" % (hours,minutes,seconds), file=sys.stdout)

94
examples/USER/cgdna/examples/oxDNA2/unique_bp/in.duplex4.4type Normal file
View File

@ -0,0 +1,94 @@
variable number equal 1
variable ofreq equal 10000
variable efreq equal 10000
variable ntype equal 4
variable T equal 0.1
units lj
dimension 3
newton off
boundary p p p
atom_style hybrid bond ellipsoid
atom_modify sort 0 1.0
# Pair interactions require lists of neighbours to be calculated
neighbor 2.0 bin
neigh_modify every 10 delay 0 check yes
read_data data.duplex4.4type
mass * 3.1575
group all type 1 4
# oxDNA bond interactions - FENE backbone
bond_style oxdna2/fene
bond_coeff * 2.0 0.25 0.7564
# oxDNA pair interactions
pair_style hybrid/overlay oxdna2/excv oxdna2/stk oxdna2/hbond oxdna2/xstk oxdna2/coaxstk oxdna2/dh
pair_coeff * * oxdna2/excv 2.0 0.7 0.675 2.0 0.515 0.5 2.0 0.33 0.32
pair_coeff * * oxdna2/stk seqdep ${T} 1.3523 2.6717 6.0 0.4 0.9 0.32 0.75 1.3 0 0.8 0.9 0 0.95 0.9 0 0.95 2.0 0.65 2.0 0.65
pair_coeff * * oxdna2/hbond seqdep 0.0 8.0 0.4 0.75 0.34 0.7 1.5 0 0.7 1.5 0 0.7 1.5 0 0.7 0.46 3.141592653589793 0.7 4.0 1.5707963267948966 0.45 4.0 1.5707963267948966 0.45
label loop
variable base loop ${ntype}
variable basemod equal ${base}%4
if "${basemod} == 1" then &
"variable comp equal ${base}+3" &
"pair_coeff ${base} ${comp} oxdna2/hbond seqdep 1.0678 8.0 0.4 0.75 0.34 0.7 1.5 0 0.7 1.5 0 0.7 1.5 0 0.7 0.46 3.141592653589793 0.7 4.0 1.5707963267948966 0.45 4.0 1.5707963267948966 0.45"
if "${basemod} == 2" then &
"variable comp equal ${base}+1" &
"pair_coeff ${base} ${comp} oxdna2/hbond seqdep 1.0678 8.0 0.4 0.75 0.34 0.7 1.5 0 0.7 1.5 0 0.7 1.5 0 0.7 0.46 3.141592653589793 0.7 4.0 1.5707963267948966 0.45 4.0 1.5707963267948966 0.45"
next base
jump in.duplex4.4type loop
pair_coeff * * oxdna2/xstk 47.5 0.575 0.675 0.495 0.655 2.25 0.791592653589793 0.58 1.7 1.0 0.68 1.7 1.0 0.68 1.5 0 0.65 1.7 0.875 0.68 1.7 0.875 0.68
pair_coeff * * oxdna2/coaxstk 58.5 0.4 0.6 0.22 0.58 2.0 2.891592653589793 0.65 1.3 0 0.8 0.9 0 0.95 0.9 0 0.95 40.0 3.116592653589793
pair_coeff * * oxdna2/dh ${T} 0.2 0.815
# Langevin dynamics
fix 1 all nve/asphere
fix 2 all langevin ${T} ${T} 25.0 457145 angmom 10
timestep 1e-4
#comm_style tiled
#fix 3 all balance 10000 1.1 rcb
#compute mol all chunk/atom molecule
#compute mychunk all vcm/chunk mol
#fix 4 all ave/time 10000 1 10000 c_mychunk[1] c_mychunk[2] c_mychunk[3] file vcm.txt mode vector
#dump pos all xyz ${ofreq} traj.${number}.xyz
compute quat all property/atom quatw quati quatj quatk
#dump quat all custom ${ofreq} quat.${number}.txt id c_quat[1] c_quat[2] c_quat[3] c_quat[4]
#dump_modify quat sort id
#dump_modify quat format line "%d %13.6le %13.6le %13.6le %13.6le"
compute erot all erotate/asphere
compute ekin all ke
compute epot all pe
variable erot equal c_erot
variable ekin equal c_ekin
variable epot equal c_epot
variable etot equal c_erot+c_ekin+c_epot
fix 5 all print ${efreq} "$(step) ekin = ${ekin} | erot = ${erot} | epot = ${epot} | etot = ${etot}" screen yes
dump out all custom ${ofreq} out.${number}.txt id mol type x y z ix iy iz c_quat[1] c_quat[2] c_quat[3] c_quat[4] vx vy vz
dump_modify out sort id
dump_modify out format line "%d %d %d %13.6le %13.6le %13.6le %d %d %d %13.6le %13.6le %13.6le %13.6le %13.6le %13.6le %13.6le "
#restart 10000 config0_restart config1_restart
run 100000
#write_restart config.${number}.*

94
examples/USER/cgdna/examples/oxDNA2/unique_bp/in.duplex4.8type Normal file
View File

@ -0,0 +1,94 @@
variable number equal 1
variable ofreq equal 10000
variable efreq equal 10000
variable ntype equal 8
variable T equal 0.1
units lj
dimension 3
newton on
boundary p p p
atom_style hybrid bond ellipsoid
atom_modify sort 0 1.0
# Pair interactions require lists of neighbours to be calculated
neighbor 1.0 bin
neigh_modify every 10 delay 0 check yes
read_data data.duplex4.8type
mass * 3.1575
group all type 1 8
# oxDNA bond interactions - FENE backbone
bond_style oxdna2/fene
bond_coeff * 2.0 0.25 0.7564
# oxDNA pair interactions
pair_style hybrid/overlay oxdna2/excv oxdna2/stk oxdna2/hbond oxdna2/xstk oxdna2/coaxstk oxdna2/dh
pair_coeff * * oxdna2/excv 2.0 0.7 0.675 2.0 0.515 0.5 2.0 0.33 0.32
pair_coeff * * oxdna2/stk seqdep ${T} 1.3523 2.6717 6.0 0.4 0.9 0.32 0.75 1.3 0 0.8 0.9 0 0.95 0.9 0 0.95 2.0 0.65 2.0 0.65
pair_coeff * * oxdna2/hbond seqdep 0.0 8.0 0.4 0.75 0.34 0.7 1.5 0 0.7 1.5 0 0.7 1.5 0 0.7 0.46 3.141592653589793 0.7 4.0 1.5707963267948966 0.45 4.0 1.5707963267948966 0.45
label loop
variable base loop ${ntype}
variable basemod equal ${base}%4
if "${basemod} == 1" then &
"variable comp equal ${base}+3" &
"pair_coeff ${base} ${comp} oxdna2/hbond seqdep 1.0678 8.0 0.4 0.75 0.34 0.7 1.5 0 0.7 1.5 0 0.7 1.5 0 0.7 0.46 3.141592653589793 0.7 4.0 1.5707963267948966 0.45 4.0 1.5707963267948966 0.45"
if "${basemod} == 2" then &
"variable comp equal ${base}+1" &
"pair_coeff ${base} ${comp} oxdna2/hbond seqdep 1.0678 8.0 0.4 0.75 0.34 0.7 1.5 0 0.7 1.5 0 0.7 1.5 0 0.7 0.46 3.141592653589793 0.7 4.0 1.5707963267948966 0.45 4.0 1.5707963267948966 0.45"
next base
jump in.duplex4.8type loop
pair_coeff * * oxdna2/xstk 47.5 0.575 0.675 0.495 0.655 2.25 0.791592653589793 0.58 1.7 1.0 0.68 1.7 1.0 0.68 1.5 0 0.65 1.7 0.875 0.68 1.7 0.875 0.68
pair_coeff * * oxdna2/coaxstk 58.5 0.4 0.6 0.22 0.58 2.0 2.891592653589793 0.65 1.3 0 0.8 0.9 0 0.95 0.9 0 0.95 40.0 3.116592653589793
pair_coeff * * oxdna2/dh ${T} 0.2 0.815
# Langevin dynamics
fix 1 all nve/asphere
fix 2 all langevin ${T} ${T} 25.0 457145 angmom 10
timestep 1e-4
#comm_style tiled
#fix 3 all balance 10000 1.1 rcb
#compute mol all chunk/atom molecule
#compute mychunk all vcm/chunk mol
#fix 4 all ave/time 10000 1 10000 c_mychunk[1] c_mychunk[2] c_mychunk[3] file vcm.txt mode vector
#dump pos all xyz ${ofreq} traj.${number}.xyz
compute quat all property/atom quatw quati quatj quatk
#dump quat all custom ${ofreq} quat.${number}.txt id c_quat[1] c_quat[2] c_quat[3] c_quat[4]
#dump_modify quat sort id
#dump_modify quat format line "%d %13.6le %13.6le %13.6le %13.6le"
compute erot all erotate/asphere
compute ekin all ke
compute epot all pe
variable erot equal c_erot
variable ekin equal c_ekin
variable epot equal c_epot
variable etot equal c_erot+c_ekin+c_epot
fix 5 all print ${efreq} "$(step) ekin = ${ekin} | erot = ${erot} | epot = ${epot} | etot = ${etot}" screen yes
dump out all custom ${ofreq} out.${number}.txt id mol type x y z ix iy iz c_quat[1] c_quat[2] c_quat[3] c_quat[4] vx vy vz
dump_modify out sort id
dump_modify out format line "%d %d %d %13.6le %13.6le %13.6le %d %d %d %13.6le %13.6le %13.6le %13.6le %13.6le %13.6le %13.6le "
#restart 10000 config0_restart config1_restart
run 100000
#write_restart config.${number}.*

232
examples/USER/cgdna/examples/oxDNA2/unique_bp/log.07Aug19.duplex4.4type.g++1 Normal file
View File

@ -0,0 +1,232 @@
LAMMPS (7 Aug 2019)
variable number equal 1
variable ofreq equal 10000
variable efreq equal 10000
variable ntype equal 4
variable T equal 0.1
units lj
dimension 3
newton off
boundary p p p
atom_style hybrid bond ellipsoid
atom_modify sort 0 1.0
# Pair interactions require lists of neighbours to be calculated
neighbor 2.0 bin
neigh_modify every 10 delay 0 check yes
read_data data.duplex4.4type
orthogonal box = (-20 -20 -20) to (20 20 20)
1 by 1 by 1 MPI processor grid
reading atoms ...
26 atoms
reading velocities ...
26 velocities
26 ellipsoids
scanning bonds ...
2 = max bonds/atom
reading bonds ...
24 bonds
2 = max # of 1-2 neighbors
2 = max # of 1-3 neighbors
4 = max # of 1-4 neighbors
6 = max # of special neighbors
special bonds CPU = 6.2e-05 secs
read_data CPU = 0.001124 secs
mass * 3.1575
group all type 1 4
26 atoms in group all
# oxDNA bond interactions - FENE backbone
bond_style oxdna2/fene
bond_coeff * 2.0 0.25 0.7564
# oxDNA pair interactions
pair_style hybrid/overlay oxdna2/excv oxdna2/stk oxdna2/hbond oxdna2/xstk oxdna2/coaxstk oxdna2/dh
pair_coeff * * oxdna2/excv 2.0 0.7 0.675 2.0 0.515 0.5 2.0 0.33 0.32
pair_coeff * * oxdna2/stk seqdep ${T} 1.3523 2.6717 6.0 0.4 0.9 0.32 0.75 1.3 0 0.8 0.9 0 0.95 0.9 0 0.95 2.0 0.65 2.0 0.65
pair_coeff * * oxdna2/stk seqdep 0.1 1.3523 2.6717 6.0 0.4 0.9 0.32 0.75 1.3 0 0.8 0.9 0 0.95 0.9 0 0.95 2.0 0.65 2.0 0.65
pair_coeff * * oxdna2/hbond seqdep 0.0 8.0 0.4 0.75 0.34 0.7 1.5 0 0.7 1.5 0 0.7 1.5 0 0.7 0.46 3.141592653589793 0.7 4.0 1.5707963267948966 0.45 4.0 1.5707963267948966 0.45
label loop
variable base loop ${ntype}
variable base loop 4
variable basemod equal ${base}%4
variable basemod equal 1%4
if "${basemod} == 1" then "variable comp equal ${base}+3" "pair_coeff ${base} ${comp} oxdna2/hbond seqdep 1.0678 8.0 0.4 0.75 0.34 0.7 1.5 0 0.7 1.5 0 0.7 1.5 0 0.7 0.46 3.141592653589793 0.7 4.0 1.5707963267948966 0.45 4.0 1.5707963267948966 0.45"
variable comp equal ${base}+3
variable comp equal 1+3
pair_coeff ${base} ${comp} oxdna2/hbond seqdep 1.0678 8.0 0.4 0.75 0.34 0.7 1.5 0 0.7 1.5 0 0.7 1.5 0 0.7 0.46 3.141592653589793 0.7 4.0 1.5707963267948966 0.45 4.0 1.5707963267948966 0.45
pair_coeff 1 ${comp} oxdna2/hbond seqdep 1.0678 8.0 0.4 0.75 0.34 0.7 1.5 0 0.7 1.5 0 0.7 1.5 0 0.7 0.46 3.141592653589793 0.7 4.0 1.5707963267948966 0.45 4.0 1.5707963267948966 0.45
pair_coeff 1 4 oxdna2/hbond seqdep 1.0678 8.0 0.4 0.75 0.34 0.7 1.5 0 0.7 1.5 0 0.7 1.5 0 0.7 0.46 3.141592653589793 0.7 4.0 1.5707963267948966 0.45 4.0 1.5707963267948966 0.45
if "${basemod} == 2" then "variable comp equal ${base}+1" "pair_coeff ${base} ${comp} oxdna2/hbond seqdep 1.0678 8.0 0.4 0.75 0.34 0.7 1.5 0 0.7 1.5 0 0.7 1.5 0 0.7 0.46 3.141592653589793 0.7 4.0 1.5707963267948966 0.45 4.0 1.5707963267948966 0.45"
next base
jump in.duplex4.4type loop
variable base loop ${ntype}
variable base loop 4
variable basemod equal ${base}%4
variable basemod equal 2%4
if "${basemod} == 1" then "variable comp equal ${base}+3" "pair_coeff ${base} ${comp} oxdna2/hbond seqdep 1.0678 8.0 0.4 0.75 0.34 0.7 1.5 0 0.7 1.5 0 0.7 1.5 0 0.7 0.46 3.141592653589793 0.7 4.0 1.5707963267948966 0.45 4.0 1.5707963267948966 0.45"
if "${basemod} == 2" then "variable comp equal ${base}+1" "pair_coeff ${base} ${comp} oxdna2/hbond seqdep 1.0678 8.0 0.4 0.75 0.34 0.7 1.5 0 0.7 1.5 0 0.7 1.5 0 0.7 0.46 3.141592653589793 0.7 4.0 1.5707963267948966 0.45 4.0 1.5707963267948966 0.45"
variable comp equal ${base}+1
variable comp equal 2+1
pair_coeff ${base} ${comp} oxdna2/hbond seqdep 1.0678 8.0 0.4 0.75 0.34 0.7 1.5 0 0.7 1.5 0 0.7 1.5 0 0.7 0.46 3.141592653589793 0.7 4.0 1.5707963267948966 0.45 4.0 1.5707963267948966 0.45
pair_coeff 2 ${comp} oxdna2/hbond seqdep 1.0678 8.0 0.4 0.75 0.34 0.7 1.5 0 0.7 1.5 0 0.7 1.5 0 0.7 0.46 3.141592653589793 0.7 4.0 1.5707963267948966 0.45 4.0 1.5707963267948966 0.45
pair_coeff 2 3 oxdna2/hbond seqdep 1.0678 8.0 0.4 0.75 0.34 0.7 1.5 0 0.7 1.5 0 0.7 1.5 0 0.7 0.46 3.141592653589793 0.7 4.0 1.5707963267948966 0.45 4.0 1.5707963267948966 0.45
next base
jump in.duplex4.4type loop
variable base loop ${ntype}
variable base loop 4
variable basemod equal ${base}%4
variable basemod equal 3%4
if "${basemod} == 1" then "variable comp equal ${base}+3" "pair_coeff ${base} ${comp} oxdna2/hbond seqdep 1.0678 8.0 0.4 0.75 0.34 0.7 1.5 0 0.7 1.5 0 0.7 1.5 0 0.7 0.46 3.141592653589793 0.7 4.0 1.5707963267948966 0.45 4.0 1.5707963267948966 0.45"
if "${basemod} == 2" then "variable comp equal ${base}+1" "pair_coeff ${base} ${comp} oxdna2/hbond seqdep 1.0678 8.0 0.4 0.75 0.34 0.7 1.5 0 0.7 1.5 0 0.7 1.5 0 0.7 0.46 3.141592653589793 0.7 4.0 1.5707963267948966 0.45 4.0 1.5707963267948966 0.45"
next base
jump in.duplex4.4type loop
variable base loop ${ntype}
variable base loop 4
variable basemod equal ${base}%4
variable basemod equal 4%4
if "${basemod} == 1" then "variable comp equal ${base}+3" "pair_coeff ${base} ${comp} oxdna2/hbond seqdep 1.0678 8.0 0.4 0.75 0.34 0.7 1.5 0 0.7 1.5 0 0.7 1.5 0 0.7 0.46 3.141592653589793 0.7 4.0 1.5707963267948966 0.45 4.0 1.5707963267948966 0.45"
if "${basemod} == 2" then "variable comp equal ${base}+1" "pair_coeff ${base} ${comp} oxdna2/hbond seqdep 1.0678 8.0 0.4 0.75 0.34 0.7 1.5 0 0.7 1.5 0 0.7 1.5 0 0.7 0.46 3.141592653589793 0.7 4.0 1.5707963267948966 0.45 4.0 1.5707963267948966 0.45"
next base
jump in.duplex4.4type loop
pair_coeff * * oxdna2/xstk 47.5 0.575 0.675 0.495 0.655 2.25 0.791592653589793 0.58 1.7 1.0 0.68 1.7 1.0 0.68 1.5 0 0.65 1.7 0.875 0.68 1.7 0.875 0.68
pair_coeff * * oxdna2/coaxstk 58.5 0.4 0.6 0.22 0.58 2.0 2.891592653589793 0.65 1.3 0 0.8 0.9 0 0.95 0.9 0 0.95 40.0 3.116592653589793
pair_coeff * * oxdna2/dh ${T} 0.2 0.815
pair_coeff * * oxdna2/dh 0.1 0.2 0.815
# Langevin dynamics
fix 1 all nve/asphere
fix 2 all langevin ${T} ${T} 25.0 457145 angmom 10
fix 2 all langevin 0.1 ${T} 25.0 457145 angmom 10
fix 2 all langevin 0.1 0.1 25.0 457145 angmom 10
timestep 1e-4
#comm_style tiled
#fix 3 all balance 10000 1.1 rcb
#compute mol all chunk/atom molecule
#compute mychunk all vcm/chunk mol
#fix 4 all ave/time 10000 1 10000 c_mychunk[1] c_mychunk[2] c_mychunk[3] file vcm.txt mode vector
#dump pos all xyz ${ofreq} traj.${number}.xyz
compute quat all property/atom quatw quati quatj quatk
#dump quat all custom ${ofreq} quat.${number}.txt id c_quat[1] c_quat[2] c_quat[3] c_quat[4]
#dump_modify quat sort id
#dump_modify quat format line "%d %13.6le %13.6le %13.6le %13.6le"
compute erot all erotate/asphere
compute ekin all ke
compute epot all pe
variable erot equal c_erot
variable ekin equal c_ekin
variable epot equal c_epot
variable etot equal c_erot+c_ekin+c_epot
fix 5 all print ${efreq} "$(step) ekin = ${ekin} | erot = ${erot} | epot = ${epot} | etot = ${etot}" screen yes
fix 5 all print 10000 "$(step) ekin = ${ekin} | erot = ${erot} | epot = ${epot} | etot = ${etot}" screen yes
dump out all custom ${ofreq} out.${number}.txt id mol type x y z ix iy iz c_quat[1] c_quat[2] c_quat[3] c_quat[4] vx vy vz
dump out all custom 10000 out.${number}.txt id mol type x y z ix iy iz c_quat[1] c_quat[2] c_quat[3] c_quat[4] vx vy vz
dump out all custom 10000 out.1.txt id mol type x y z ix iy iz c_quat[1] c_quat[2] c_quat[3] c_quat[4] vx vy vz
dump_modify out sort id
dump_modify out format line "%d %d %d %13.6le %13.6le %13.6le %d %d %d %13.6le %13.6le %13.6le %13.6le %13.6le %13.6le %13.6le "
#restart 10000 config0_restart config1_restart
run 100000
Neighbor list info ...
update every 10 steps, delay 0 steps, check yes
max neighbors/atom: 2000, page size: 100000
master list distance cutoff = 5.63899
ghost atom cutoff = 5.63899
binsize = 2.81949, bins = 15 15 15
6 neighbor lists, perpetual/occasional/extra = 6 0 0
(1) pair oxdna2/excv, perpetual
attributes: half, newton off
pair build: half/bin/newtoff
stencil: half/bin/3d/newtoff
bin: standard
(2) pair oxdna2/stk, perpetual, copy from (1)
attributes: half, newton off
pair build: copy
stencil: none
bin: none
(3) pair oxdna2/hbond, perpetual, copy from (1)
attributes: half, newton off
pair build: copy
stencil: none
bin: none
(4) pair oxdna2/xstk, perpetual, copy from (1)
attributes: half, newton off
pair build: copy
stencil: none
bin: none
(5) pair oxdna2/coaxstk, perpetual, copy from (1)
attributes: half, newton off
pair build: copy
stencil: none
bin: none
(6) pair oxdna2/dh, perpetual, copy from (1)
attributes: half, newton off
pair build: copy
stencil: none
bin: none
Per MPI rank memory allocation (min/avg/max) = 3.89 | 3.89 | 3.89 Mbytes
Step Temp E_pair E_mol TotEng Press
0 0 -1.6384018 0.037304147 -1.6010976 6.7769766e-05
10000 ekin = 1.01296379553452 | erot = 1.38476360245491 | epot = -41.9785360036035 | etot = -39.5808086056141
20000 ekin = 1.75607695499755 | erot = 1.86620102096623 | epot = -41.9558496477913 | etot = -38.3335716718275
30000 ekin = 2.25878673688255 | erot = 2.60035907744906 | epot = -41.2869827431736 | etot = -36.427836928842
40000 ekin = 2.11105717434601 | erot = 3.16611217122795 | epot = -40.2762731426449 | etot = -34.999103797071
50000 ekin = 3.05379083545185 | erot = 2.45050275456902 | epot = -40.3032957287999 | etot = -34.799002138779
60000 ekin = 3.05154113920288 | erot = 2.52539576708931 | epot = -39.6800099891302 | etot = -34.103073082838
70000 ekin = 3.23212209366503 | erot = 3.28914763888578 | epot = -40.1240667487288 | etot = -33.602797016178
80000 ekin = 2.82623224207568 | erot = 2.91292948677805 | epot = -39.0983962904507 | etot = -33.3592345615969
90000 ekin = 3.05995314905566 | erot = 3.5429982411034 | epot = -39.1646070343971 | etot = -32.561655644238
100000 ekin = 3.46139546969836 | erot = 4.11704803295073 | epot = -38.5124679837256 | etot = -30.9340244810765
100000 0.092303879 -1.5243833 0.043134544 -1.3481182 6.862374e-06
Loop time of 7.99505 on 1 procs for 100000 steps with 26 atoms
Performance: 108066.893 tau/day, 12507.742 timesteps/s
99.8% CPU use with 1 MPI tasks x no OpenMP threads
MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total
---------------------------------------------------------------
Pair | 7.2072 | 7.2072 | 7.2072 | 0.0 | 90.15
Bond | 0.14292 | 0.14292 | 0.14292 | 0.0 | 1.79
Neigh | 1.9e-05 | 1.9e-05 | 1.9e-05 | 0.0 | 0.00
Comm | 0.015676 | 0.015676 | 0.015676 | 0.0 | 0.20
Output | 0.001372 | 0.001372 | 0.001372 | 0.0 | 0.02
Modify | 0.61071 | 0.61071 | 0.61071 | 0.0 | 7.64
Other | | 0.01714 | | | 0.21
Nlocal: 26 ave 26 max 26 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Nghost: 0 ave 0 max 0 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Neighs: 325 ave 325 max 325 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Total # of neighbors = 325
Ave neighs/atom = 12.5
Ave special neighs/atom = 5.07692
Neighbor list builds = 1
Dangerous builds = 0
#write_restart config.${number}.*
Total wall time: 0:00:07

232
examples/USER/cgdna/examples/oxDNA2/unique_bp/log.07Aug19.duplex4.4type.g++4 Normal file
View File

@ -0,0 +1,232 @@
LAMMPS (7 Aug 2019)
variable number equal 1
variable ofreq equal 10000
variable efreq equal 10000
variable ntype equal 4
variable T equal 0.1
units lj
dimension 3
newton off
boundary p p p
atom_style hybrid bond ellipsoid
atom_modify sort 0 1.0
# Pair interactions require lists of neighbours to be calculated
neighbor 2.0 bin
neigh_modify every 10 delay 0 check yes
read_data data.duplex4.4type
orthogonal box = (-20 -20 -20) to (20 20 20)
1 by 2 by 2 MPI processor grid
reading atoms ...
26 atoms
reading velocities ...
26 velocities
26 ellipsoids
scanning bonds ...
2 = max bonds/atom
reading bonds ...
24 bonds
2 = max # of 1-2 neighbors
2 = max # of 1-3 neighbors
4 = max # of 1-4 neighbors
6 = max # of special neighbors
special bonds CPU = 0.000173 secs
read_data CPU = 0.00242 secs
mass * 3.1575
group all type 1 4
26 atoms in group all
# oxDNA bond interactions - FENE backbone
bond_style oxdna2/fene
bond_coeff * 2.0 0.25 0.7564
# oxDNA pair interactions
pair_style hybrid/overlay oxdna2/excv oxdna2/stk oxdna2/hbond oxdna2/xstk oxdna2/coaxstk oxdna2/dh
pair_coeff * * oxdna2/excv 2.0 0.7 0.675 2.0 0.515 0.5 2.0 0.33 0.32
pair_coeff * * oxdna2/stk seqdep ${T} 1.3523 2.6717 6.0 0.4 0.9 0.32 0.75 1.3 0 0.8 0.9 0 0.95 0.9 0 0.95 2.0 0.65 2.0 0.65
pair_coeff * * oxdna2/stk seqdep 0.1 1.3523 2.6717 6.0 0.4 0.9 0.32 0.75 1.3 0 0.8 0.9 0 0.95 0.9 0 0.95 2.0 0.65 2.0 0.65
pair_coeff * * oxdna2/hbond seqdep 0.0 8.0 0.4 0.75 0.34 0.7 1.5 0 0.7 1.5 0 0.7 1.5 0 0.7 0.46 3.141592653589793 0.7 4.0 1.5707963267948966 0.45 4.0 1.5707963267948966 0.45
label loop
variable base loop ${ntype}
variable base loop 4
variable basemod equal ${base}%4
variable basemod equal 1%4
if "${basemod} == 1" then "variable comp equal ${base}+3" "pair_coeff ${base} ${comp} oxdna2/hbond seqdep 1.0678 8.0 0.4 0.75 0.34 0.7 1.5 0 0.7 1.5 0 0.7 1.5 0 0.7 0.46 3.141592653589793 0.7 4.0 1.5707963267948966 0.45 4.0 1.5707963267948966 0.45"
variable comp equal ${base}+3
variable comp equal 1+3
pair_coeff ${base} ${comp} oxdna2/hbond seqdep 1.0678 8.0 0.4 0.75 0.34 0.7 1.5 0 0.7 1.5 0 0.7 1.5 0 0.7 0.46 3.141592653589793 0.7 4.0 1.5707963267948966 0.45 4.0 1.5707963267948966 0.45
pair_coeff 1 ${comp} oxdna2/hbond seqdep 1.0678 8.0 0.4 0.75 0.34 0.7 1.5 0 0.7 1.5 0 0.7 1.5 0 0.7 0.46 3.141592653589793 0.7 4.0 1.5707963267948966 0.45 4.0 1.5707963267948966 0.45
pair_coeff 1 4 oxdna2/hbond seqdep 1.0678 8.0 0.4 0.75 0.34 0.7 1.5 0 0.7 1.5 0 0.7 1.5 0 0.7 0.46 3.141592653589793 0.7 4.0 1.5707963267948966 0.45 4.0 1.5707963267948966 0.45
if "${basemod} == 2" then "variable comp equal ${base}+1" "pair_coeff ${base} ${comp} oxdna2/hbond seqdep 1.0678 8.0 0.4 0.75 0.34 0.7 1.5 0 0.7 1.5 0 0.7 1.5 0 0.7 0.46 3.141592653589793 0.7 4.0 1.5707963267948966 0.45 4.0 1.5707963267948966 0.45"
next base
jump in.duplex4.4type loop
variable base loop ${ntype}
variable base loop 4
variable basemod equal ${base}%4
variable basemod equal 2%4
if "${basemod} == 1" then "variable comp equal ${base}+3" "pair_coeff ${base} ${comp} oxdna2/hbond seqdep 1.0678 8.0 0.4 0.75 0.34 0.7 1.5 0 0.7 1.5 0 0.7 1.5 0 0.7 0.46 3.141592653589793 0.7 4.0 1.5707963267948966 0.45 4.0 1.5707963267948966 0.45"
if "${basemod} == 2" then "variable comp equal ${base}+1" "pair_coeff ${base} ${comp} oxdna2/hbond seqdep 1.0678 8.0 0.4 0.75 0.34 0.7 1.5 0 0.7 1.5 0 0.7 1.5 0 0.7 0.46 3.141592653589793 0.7 4.0 1.5707963267948966 0.45 4.0 1.5707963267948966 0.45"
variable comp equal ${base}+1
variable comp equal 2+1
pair_coeff ${base} ${comp} oxdna2/hbond seqdep 1.0678 8.0 0.4 0.75 0.34 0.7 1.5 0 0.7 1.5 0 0.7 1.5 0 0.7 0.46 3.141592653589793 0.7 4.0 1.5707963267948966 0.45 4.0 1.5707963267948966 0.45
pair_coeff 2 ${comp} oxdna2/hbond seqdep 1.0678 8.0 0.4 0.75 0.34 0.7 1.5 0 0.7 1.5 0 0.7 1.5 0 0.7 0.46 3.141592653589793 0.7 4.0 1.5707963267948966 0.45 4.0 1.5707963267948966 0.45
pair_coeff 2 3 oxdna2/hbond seqdep 1.0678 8.0 0.4 0.75 0.34 0.7 1.5 0 0.7 1.5 0 0.7 1.5 0 0.7 0.46 3.141592653589793 0.7 4.0 1.5707963267948966 0.45 4.0 1.5707963267948966 0.45
next base
jump in.duplex4.4type loop
variable base loop ${ntype}
variable base loop 4
variable basemod equal ${base}%4
variable basemod equal 3%4
if "${basemod} == 1" then "variable comp equal ${base}+3" "pair_coeff ${base} ${comp} oxdna2/hbond seqdep 1.0678 8.0 0.4 0.75 0.34 0.7 1.5 0 0.7 1.5 0 0.7 1.5 0 0.7 0.46 3.141592653589793 0.7 4.0 1.5707963267948966 0.45 4.0 1.5707963267948966 0.45"
if "${basemod} == 2" then "variable comp equal ${base}+1" "pair_coeff ${base} ${comp} oxdna2/hbond seqdep 1.0678 8.0 0.4 0.75 0.34 0.7 1.5 0 0.7 1.5 0 0.7 1.5 0 0.7 0.46 3.141592653589793 0.7 4.0 1.5707963267948966 0.45 4.0 1.5707963267948966 0.45"
next base
jump in.duplex4.4type loop
variable base loop ${ntype}
variable base loop 4
variable basemod equal ${base}%4
variable basemod equal 4%4
if "${basemod} == 1" then "variable comp equal ${base}+3" "pair_coeff ${base} ${comp} oxdna2/hbond seqdep 1.0678 8.0 0.4 0.75 0.34 0.7 1.5 0 0.7 1.5 0 0.7 1.5 0 0.7 0.46 3.141592653589793 0.7 4.0 1.5707963267948966 0.45 4.0 1.5707963267948966 0.45"
if "${basemod} == 2" then "variable comp equal ${base}+1" "pair_coeff ${base} ${comp} oxdna2/hbond seqdep 1.0678 8.0 0.4 0.75 0.34 0.7 1.5 0 0.7 1.5 0 0.7 1.5 0 0.7 0.46 3.141592653589793 0.7 4.0 1.5707963267948966 0.45 4.0 1.5707963267948966 0.45"
next base
jump in.duplex4.4type loop
pair_coeff * * oxdna2/xstk 47.5 0.575 0.675 0.495 0.655 2.25 0.791592653589793 0.58 1.7 1.0 0.68 1.7 1.0 0.68 1.5 0 0.65 1.7 0.875 0.68 1.7 0.875 0.68
pair_coeff * * oxdna2/coaxstk 58.5 0.4 0.6 0.22 0.58 2.0 2.891592653589793 0.65 1.3 0 0.8 0.9 0 0.95 0.9 0 0.95 40.0 3.116592653589793
pair_coeff * * oxdna2/dh ${T} 0.2 0.815
pair_coeff * * oxdna2/dh 0.1 0.2 0.815
# Langevin dynamics
fix 1 all nve/asphere
fix 2 all langevin ${T} ${T} 25.0 457145 angmom 10
fix 2 all langevin 0.1 ${T} 25.0 457145 angmom 10
fix 2 all langevin 0.1 0.1 25.0 457145 angmom 10
timestep 1e-4
#comm_style tiled
#fix 3 all balance 10000 1.1 rcb
#compute mol all chunk/atom molecule
#compute mychunk all vcm/chunk mol
#fix 4 all ave/time 10000 1 10000 c_mychunk[1] c_mychunk[2] c_mychunk[3] file vcm.txt mode vector
#dump pos all xyz ${ofreq} traj.${number}.xyz
compute quat all property/atom quatw quati quatj quatk
#dump quat all custom ${ofreq} quat.${number}.txt id c_quat[1] c_quat[2] c_quat[3] c_quat[4]
#dump_modify quat sort id
#dump_modify quat format line "%d %13.6le %13.6le %13.6le %13.6le"
compute erot all erotate/asphere
compute ekin all ke
compute epot all pe
variable erot equal c_erot
variable ekin equal c_ekin
variable epot equal c_epot
variable etot equal c_erot+c_ekin+c_epot
fix 5 all print ${efreq} "$(step) ekin = ${ekin} | erot = ${erot} | epot = ${epot} | etot = ${etot}" screen yes
fix 5 all print 10000 "$(step) ekin = ${ekin} | erot = ${erot} | epot = ${epot} | etot = ${etot}" screen yes
dump out all custom ${ofreq} out.${number}.txt id mol type x y z ix iy iz c_quat[1] c_quat[2] c_quat[3] c_quat[4] vx vy vz
dump out all custom 10000 out.${number}.txt id mol type x y z ix iy iz c_quat[1] c_quat[2] c_quat[3] c_quat[4] vx vy vz
dump out all custom 10000 out.1.txt id mol type x y z ix iy iz c_quat[1] c_quat[2] c_quat[3] c_quat[4] vx vy vz
dump_modify out sort id
dump_modify out format line "%d %d %d %13.6le %13.6le %13.6le %d %d %d %13.6le %13.6le %13.6le %13.6le %13.6le %13.6le %13.6le "
#restart 10000 config0_restart config1_restart
run 100000
Neighbor list info ...
update every 10 steps, delay 0 steps, check yes
max neighbors/atom: 2000, page size: 100000
master list distance cutoff = 5.63899
ghost atom cutoff = 5.63899
binsize = 2.81949, bins = 15 15 15
6 neighbor lists, perpetual/occasional/extra = 6 0 0
(1) pair oxdna2/excv, perpetual
attributes: half, newton off
pair build: half/bin/newtoff
stencil: half/bin/3d/newtoff
bin: standard
(2) pair oxdna2/stk, perpetual, copy from (1)
attributes: half, newton off
pair build: copy
stencil: none
bin: none
(3) pair oxdna2/hbond, perpetual, copy from (1)
attributes: half, newton off
pair build: copy
stencil: none
bin: none
(4) pair oxdna2/xstk, perpetual, copy from (1)
attributes: half, newton off
pair build: copy
stencil: none
bin: none
(5) pair oxdna2/coaxstk, perpetual, copy from (1)
attributes: half, newton off
pair build: copy
stencil: none
bin: none
(6) pair oxdna2/dh, perpetual, copy from (1)
attributes: half, newton off
pair build: copy
stencil: none
bin: none
Per MPI rank memory allocation (min/avg/max) = 9.499 | 9.682 | 9.864 Mbytes
Step Temp E_pair E_mol TotEng Press
0 0 -1.6384018 0.037304147 -1.6010976 6.7769766e-05
10000 ekin = 1.20524984175817 | erot = 1.2679122962101 | epot = -41.9183873815797 | etot = -39.4452252436114
20000 ekin = 1.60856540355485 | erot = 2.64140880508964 | epot = -41.0106168994545 | etot = -36.76064269081
30000 ekin = 2.821101638199 | erot = 3.06373823252434 | epot = -40.9348041675027 | etot = -35.0499642967793
40000 ekin = 2.36296917392486 | erot = 3.53353555114673 | epot = -39.5591676362755 | etot = -33.6626629112039
50000 ekin = 2.93383938780818 | erot = 3.39920211088179 | epot = -40.3602175538497 | etot = -34.0271760551597
60000 ekin = 2.42811021532829 | erot = 2.75130045540831 | epot = -39.3178531324554 | etot = -34.1384424617188
70000 ekin = 2.7282200818863 | erot = 3.76502037022117 | epot = -40.5673105068195 | etot = -34.074070054712
80000 ekin = 2.90877851656705 | erot = 2.43617899356904 | epot = -38.4099618426175 | etot = -33.0650043324814
90000 ekin = 3.89203816080263 | erot = 2.67194811393274 | epot = -39.1880466354802 | etot = -32.6240603607448
100000 ekin = 3.98445325397832 | erot = 2.77079124049636 | epot = -39.2805505658488 | etot = -32.5253060713741
100000 0.10625209 -1.5603519 0.049561514 -1.3575422 -1.7755557e-05
Loop time of 7.14827 on 4 procs for 100000 steps with 26 atoms
Performance: 120868.376 tau/day, 13989.395 timesteps/s
99.7% CPU use with 4 MPI tasks x no OpenMP threads
MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total
---------------------------------------------------------------
Pair | 0.047776 | 3.0726 | 6.3143 | 172.8 | 42.98
Bond | 0.00736 | 0.056846 | 0.10864 | 20.7 | 0.80
Neigh | 1.9e-05 | 1.9e-05 | 1.9e-05 | 0.0 | 0.00
Comm | 0.30925 | 3.451 | 6.3739 | 157.5 | 48.28
Output | 0.000896 | 0.0010197 | 0.001301 | 0.5 | 0.01
Modify | 0.015512 | 0.18417 | 0.37845 | 39.4 | 2.58
Other | | 0.3826 | | | 5.35
Nlocal: 6.5 ave 13 max 0 min
Histogram: 2 0 0 0 0 0 0 0 0 2
Nghost: 19.5 ave 26 max 13 min
Histogram: 2 0 0 0 0 0 0 0 0 2
Neighs: 123.5 ave 247 max 0 min
Histogram: 2 0 0 0 0 0 0 0 0 2
Total # of neighbors = 494
Ave neighs/atom = 19
Ave special neighs/atom = 5.07692
Neighbor list builds = 1
Dangerous builds = 0
#write_restart config.${number}.*
Total wall time: 0:00:07

274
examples/USER/cgdna/examples/oxDNA2/unique_bp/log.07Aug19.duplex4.8type.g++1 Normal file
View File

@ -0,0 +1,274 @@
LAMMPS (7 Aug 2019)
variable number equal 1
variable ofreq equal 10000
variable efreq equal 10000
variable ntype equal 8
variable T equal 0.1
units lj
dimension 3
newton on
boundary p p p
atom_style hybrid bond ellipsoid
atom_modify sort 0 1.0
# Pair interactions require lists of neighbours to be calculated
neighbor 1.0 bin
neigh_modify every 10 delay 0 check yes
read_data data.duplex4.8type
orthogonal box = (-20 -20 -20) to (20 20 20)
1 by 1 by 1 MPI processor grid
reading atoms ...
26 atoms
reading velocities ...
26 velocities
26 ellipsoids
scanning bonds ...
1 = max bonds/atom
reading bonds ...
24 bonds
2 = max # of 1-2 neighbors
2 = max # of 1-3 neighbors
4 = max # of 1-4 neighbors
6 = max # of special neighbors
special bonds CPU = 6.5e-05 secs
read_data CPU = 0.001139 secs
mass * 3.1575
group all type 1 8
26 atoms in group all
# oxDNA bond interactions - FENE backbone
bond_style oxdna2/fene
bond_coeff * 2.0 0.25 0.7564
# oxDNA pair interactions
pair_style hybrid/overlay oxdna2/excv oxdna2/stk oxdna2/hbond oxdna2/xstk oxdna2/coaxstk oxdna2/dh
pair_coeff * * oxdna2/excv 2.0 0.7 0.675 2.0 0.515 0.5 2.0 0.33 0.32
pair_coeff * * oxdna2/stk seqdep ${T} 1.3523 2.6717 6.0 0.4 0.9 0.32 0.75 1.3 0 0.8 0.9 0 0.95 0.9 0 0.95 2.0 0.65 2.0 0.65
pair_coeff * * oxdna2/stk seqdep 0.1 1.3523 2.6717 6.0 0.4 0.9 0.32 0.75 1.3 0 0.8 0.9 0 0.95 0.9 0 0.95 2.0 0.65 2.0 0.65
pair_coeff * * oxdna2/hbond seqdep 0.0 8.0 0.4 0.75 0.34 0.7 1.5 0 0.7 1.5 0 0.7 1.5 0 0.7 0.46 3.141592653589793 0.7 4.0 1.5707963267948966 0.45 4.0 1.5707963267948966 0.45
label loop
variable base loop ${ntype}
variable base loop 8
variable basemod equal ${base}%4
variable basemod equal 1%4
if "${basemod} == 1" then "variable comp equal ${base}+3" "pair_coeff ${base} ${comp} oxdna2/hbond seqdep 1.0678 8.0 0.4 0.75 0.34 0.7 1.5 0 0.7 1.5 0 0.7 1.5 0 0.7 0.46 3.141592653589793 0.7 4.0 1.5707963267948966 0.45 4.0 1.5707963267948966 0.45"
variable comp equal ${base}+3
variable comp equal 1+3
pair_coeff ${base} ${comp} oxdna2/hbond seqdep 1.0678 8.0 0.4 0.75 0.34 0.7 1.5 0 0.7 1.5 0 0.7 1.5 0 0.7 0.46 3.141592653589793 0.7 4.0 1.5707963267948966 0.45 4.0 1.5707963267948966 0.45
pair_coeff 1 ${comp} oxdna2/hbond seqdep 1.0678 8.0 0.4 0.75 0.34 0.7 1.5 0 0.7 1.5 0 0.7 1.5 0 0.7 0.46 3.141592653589793 0.7 4.0 1.5707963267948966 0.45 4.0 1.5707963267948966 0.45
pair_coeff 1 4 oxdna2/hbond seqdep 1.0678 8.0 0.4 0.75 0.34 0.7 1.5 0 0.7 1.5 0 0.7 1.5 0 0.7 0.46 3.141592653589793 0.7 4.0 1.5707963267948966 0.45 4.0 1.5707963267948966 0.45
if "${basemod} == 2" then "variable comp equal ${base}+1" "pair_coeff ${base} ${comp} oxdna2/hbond seqdep 1.0678 8.0 0.4 0.75 0.34 0.7 1.5 0 0.7 1.5 0 0.7 1.5 0 0.7 0.46 3.141592653589793 0.7 4.0 1.5707963267948966 0.45 4.0 1.5707963267948966 0.45"
next base
jump in.duplex4.8type loop
variable base loop ${ntype}
variable base loop 8
variable basemod equal ${base}%4
variable basemod equal 2%4
if "${basemod} == 1" then "variable comp equal ${base}+3" "pair_coeff ${base} ${comp} oxdna2/hbond seqdep 1.0678 8.0 0.4 0.75 0.34 0.7 1.5 0 0.7 1.5 0 0.7 1.5 0 0.7 0.46 3.141592653589793 0.7 4.0 1.5707963267948966 0.45 4.0 1.5707963267948966 0.45"
if "${basemod} == 2" then "variable comp equal ${base}+1" "pair_coeff ${base} ${comp} oxdna2/hbond seqdep 1.0678 8.0 0.4 0.75 0.34 0.7 1.5 0 0.7 1.5 0 0.7 1.5 0 0.7 0.46 3.141592653589793 0.7 4.0 1.5707963267948966 0.45 4.0 1.5707963267948966 0.45"
variable comp equal ${base}+1
variable comp equal 2+1
pair_coeff ${base} ${comp} oxdna2/hbond seqdep 1.0678 8.0 0.4 0.75 0.34 0.7 1.5 0 0.7 1.5 0 0.7 1.5 0 0.7 0.46 3.141592653589793 0.7 4.0 1.5707963267948966 0.45 4.0 1.5707963267948966 0.45
pair_coeff 2 ${comp} oxdna2/hbond seqdep 1.0678 8.0 0.4 0.75 0.34 0.7 1.5 0 0.7 1.5 0 0.7 1.5 0 0.7 0.46 3.141592653589793 0.7 4.0 1.5707963267948966 0.45 4.0 1.5707963267948966 0.45
pair_coeff 2 3 oxdna2/hbond seqdep 1.0678 8.0 0.4 0.75 0.34 0.7 1.5 0 0.7 1.5 0 0.7 1.5 0 0.7 0.46 3.141592653589793 0.7 4.0 1.5707963267948966 0.45 4.0 1.5707963267948966 0.45
next base
jump in.duplex4.8type loop
variable base loop ${ntype}
variable base loop 8
variable basemod equal ${base}%4
variable basemod equal 3%4
if "${basemod} == 1" then "variable comp equal ${base}+3" "pair_coeff ${base} ${comp} oxdna2/hbond seqdep 1.0678 8.0 0.4 0.75 0.34 0.7 1.5 0 0.7 1.5 0 0.7 1.5 0 0.7 0.46 3.141592653589793 0.7 4.0 1.5707963267948966 0.45 4.0 1.5707963267948966 0.45"
if "${basemod} == 2" then "variable comp equal ${base}+1" "pair_coeff ${base} ${comp} oxdna2/hbond seqdep 1.0678 8.0 0.4 0.75 0.34 0.7 1.5 0 0.7 1.5 0 0.7 1.5 0 0.7 0.46 3.141592653589793 0.7 4.0 1.5707963267948966 0.45 4.0 1.5707963267948966 0.45"
next base
jump in.duplex4.8type loop
variable base loop ${ntype}
variable base loop 8
variable basemod equal ${base}%4
variable basemod equal 4%4
if "${basemod} == 1" then "variable comp equal ${base}+3" "pair_coeff ${base} ${comp} oxdna2/hbond seqdep 1.0678 8.0 0.4 0.75 0.34 0.7 1.5 0 0.7 1.5 0 0.7 1.5 0 0.7 0.46 3.141592653589793 0.7 4.0 1.5707963267948966 0.45 4.0 1.5707963267948966 0.45"
if "${basemod} == 2" then "variable comp equal ${base}+1" "pair_coeff ${base} ${comp} oxdna2/hbond seqdep 1.0678 8.0 0.4 0.75 0.34 0.7 1.5 0 0.7 1.5 0 0.7 1.5 0 0.7 0.46 3.141592653589793 0.7 4.0 1.5707963267948966 0.45 4.0 1.5707963267948966 0.45"
next base
jump in.duplex4.8type loop
variable base loop ${ntype}
variable base loop 8
variable basemod equal ${base}%4
variable basemod equal 5%4
if "${basemod} == 1" then "variable comp equal ${base}+3" "pair_coeff ${base} ${comp} oxdna2/hbond seqdep 1.0678 8.0 0.4 0.75 0.34 0.7 1.5 0 0.7 1.5 0 0.7 1.5 0 0.7 0.46 3.141592653589793 0.7 4.0 1.5707963267948966 0.45 4.0 1.5707963267948966 0.45"
variable comp equal ${base}+3
variable comp equal 5+3
pair_coeff ${base} ${comp} oxdna2/hbond seqdep 1.0678 8.0 0.4 0.75 0.34 0.7 1.5 0 0.7 1.5 0 0.7 1.5 0 0.7 0.46 3.141592653589793 0.7 4.0 1.5707963267948966 0.45 4.0 1.5707963267948966 0.45
pair_coeff 5 ${comp} oxdna2/hbond seqdep 1.0678 8.0 0.4 0.75 0.34 0.7 1.5 0 0.7 1.5 0 0.7 1.5 0 0.7 0.46 3.141592653589793 0.7 4.0 1.5707963267948966 0.45 4.0 1.5707963267948966 0.45
pair_coeff 5 8 oxdna2/hbond seqdep 1.0678 8.0 0.4 0.75 0.34 0.7 1.5 0 0.7 1.5 0 0.7 1.5 0 0.7 0.46 3.141592653589793 0.7 4.0 1.5707963267948966 0.45 4.0 1.5707963267948966 0.45
if "${basemod} == 2" then "variable comp equal ${base}+1" "pair_coeff ${base} ${comp} oxdna2/hbond seqdep 1.0678 8.0 0.4 0.75 0.34 0.7 1.5 0 0.7 1.5 0 0.7 1.5 0 0.7 0.46 3.141592653589793 0.7 4.0 1.5707963267948966 0.45 4.0 1.5707963267948966 0.45"
next base
jump in.duplex4.8type loop
variable base loop ${ntype}
variable base loop 8
variable basemod equal ${base}%4
variable basemod equal 6%4
if "${basemod} == 1" then "variable comp equal ${base}+3" "pair_coeff ${base} ${comp} oxdna2/hbond seqdep 1.0678 8.0 0.4 0.75 0.34 0.7 1.5 0 0.7 1.5 0 0.7 1.5 0 0.7 0.46 3.141592653589793 0.7 4.0 1.5707963267948966 0.45 4.0 1.5707963267948966 0.45"
if "${basemod} == 2" then "variable comp equal ${base}+1" "pair_coeff ${base} ${comp} oxdna2/hbond seqdep 1.0678 8.0 0.4 0.75 0.34 0.7 1.5 0 0.7 1.5 0 0.7 1.5 0 0.7 0.46 3.141592653589793 0.7 4.0 1.5707963267948966 0.45 4.0 1.5707963267948966 0.45"
variable comp equal ${base}+1
variable comp equal 6+1
pair_coeff ${base} ${comp} oxdna2/hbond seqdep 1.0678 8.0 0.4 0.75 0.34 0.7 1.5 0 0.7 1.5 0 0.7 1.5 0 0.7 0.46 3.141592653589793 0.7 4.0 1.5707963267948966 0.45 4.0 1.5707963267948966 0.45
pair_coeff 6 ${comp} oxdna2/hbond seqdep 1.0678 8.0 0.4 0.75 0.34 0.7 1.5 0 0.7 1.5 0 0.7 1.5 0 0.7 0.46 3.141592653589793 0.7 4.0 1.5707963267948966 0.45 4.0 1.5707963267948966 0.45
pair_coeff 6 7 oxdna2/hbond seqdep 1.0678 8.0 0.4 0.75 0.34 0.7 1.5 0 0.7 1.5 0 0.7 1.5 0 0.7 0.46 3.141592653589793 0.7 4.0 1.5707963267948966 0.45 4.0 1.5707963267948966 0.45
next base
jump in.duplex4.8type loop
variable base loop ${ntype}
variable base loop 8
variable basemod equal ${base}%4
variable basemod equal 7%4
if "${basemod} == 1" then "variable comp equal ${base}+3" "pair_coeff ${base} ${comp} oxdna2/hbond seqdep 1.0678 8.0 0.4 0.75 0.34 0.7 1.5 0 0.7 1.5 0 0.7 1.5 0 0.7 0.46 3.141592653589793 0.7 4.0 1.5707963267948966 0.45 4.0 1.5707963267948966 0.45"
if "${basemod} == 2" then "variable comp equal ${base}+1" "pair_coeff ${base} ${comp} oxdna2/hbond seqdep 1.0678 8.0 0.4 0.75 0.34 0.7 1.5 0 0.7 1.5 0 0.7 1.5 0 0.7 0.46 3.141592653589793 0.7 4.0 1.5707963267948966 0.45 4.0 1.5707963267948966 0.45"
next base
jump in.duplex4.8type loop
variable base loop ${ntype}
variable base loop 8
variable basemod equal ${base}%4
variable basemod equal 8%4
if "${basemod} == 1" then "variable comp equal ${base}+3" "pair_coeff ${base} ${comp} oxdna2/hbond seqdep 1.0678 8.0 0.4 0.75 0.34 0.7 1.5 0 0.7 1.5 0 0.7 1.5 0 0.7 0.46 3.141592653589793 0.7 4.0 1.5707963267948966 0.45 4.0 1.5707963267948966 0.45"
if "${basemod} == 2" then "variable comp equal ${base}+1" "pair_coeff ${base} ${comp} oxdna2/hbond seqdep 1.0678 8.0 0.4 0.75 0.34 0.7 1.5 0 0.7 1.5 0 0.7 1.5 0 0.7 0.46 3.141592653589793 0.7 4.0 1.5707963267948966 0.45 4.0 1.5707963267948966 0.45"
next base
jump in.duplex4.8type loop
pair_coeff * * oxdna2/xstk 47.5 0.575 0.675 0.495 0.655 2.25 0.791592653589793 0.58 1.7 1.0 0.68 1.7 1.0 0.68 1.5 0 0.65 1.7 0.875 0.68 1.7 0.875 0.68
pair_coeff * * oxdna2/coaxstk 58.5 0.4 0.6 0.22 0.58 2.0 2.891592653589793 0.65 1.3 0 0.8 0.9 0 0.95 0.9 0 0.95 40.0 3.116592653589793
pair_coeff * * oxdna2/dh ${T} 0.2 0.815
pair_coeff * * oxdna2/dh 0.1 0.2 0.815
# Langevin dynamics
fix 1 all nve/asphere
fix 2 all langevin ${T} ${T} 25.0 457145 angmom 10
fix 2 all langevin 0.1 ${T} 25.0 457145 angmom 10
fix 2 all langevin 0.1 0.1 25.0 457145 angmom 10
timestep 1e-4
#comm_style tiled
#fix 3 all balance 10000 1.1 rcb
#compute mol all chunk/atom molecule
#compute mychunk all vcm/chunk mol
#fix 4 all ave/time 10000 1 10000 c_mychunk[1] c_mychunk[2] c_mychunk[3] file vcm.txt mode vector
#dump pos all xyz ${ofreq} traj.${number}.xyz
compute quat all property/atom quatw quati quatj quatk
#dump quat all custom ${ofreq} quat.${number}.txt id c_quat[1] c_quat[2] c_quat[3] c_quat[4]
#dump_modify quat sort id
#dump_modify quat format line "%d %13.6le %13.6le %13.6le %13.6le"
compute erot all erotate/asphere
compute ekin all ke
compute epot all pe
variable erot equal c_erot
variable ekin equal c_ekin
variable epot equal c_epot
variable etot equal c_erot+c_ekin+c_epot
fix 5 all print ${efreq} "$(step) ekin = ${ekin} | erot = ${erot} | epot = ${epot} | etot = ${etot}" screen yes
fix 5 all print 10000 "$(step) ekin = ${ekin} | erot = ${erot} | epot = ${epot} | etot = ${etot}" screen yes
dump out all custom ${ofreq} out.${number}.txt id mol type x y z ix iy iz c_quat[1] c_quat[2] c_quat[3] c_quat[4] vx vy vz
dump out all custom 10000 out.${number}.txt id mol type x y z ix iy iz c_quat[1] c_quat[2] c_quat[3] c_quat[4] vx vy vz
dump out all custom 10000 out.1.txt id mol type x y z ix iy iz c_quat[1] c_quat[2] c_quat[3] c_quat[4] vx vy vz
dump_modify out sort id
dump_modify out format line "%d %d %d %13.6le %13.6le %13.6le %d %d %d %13.6le %13.6le %13.6le %13.6le %13.6le %13.6le %13.6le "
#restart 10000 config0_restart config1_restart
run 100000
Neighbor list info ...
update every 10 steps, delay 0 steps, check yes
max neighbors/atom: 2000, page size: 100000
master list distance cutoff = 4.63899
ghost atom cutoff = 4.63899
binsize = 2.31949, bins = 18 18 18
6 neighbor lists, perpetual/occasional/extra = 6 0 0
(1) pair oxdna2/excv, perpetual
attributes: half, newton on
pair build: half/bin/newton
stencil: half/bin/3d/newton
bin: standard
(2) pair oxdna2/stk, perpetual, copy from (1)
attributes: half, newton on
pair build: copy
stencil: none
bin: none
(3) pair oxdna2/hbond, perpetual, copy from (1)
attributes: half, newton on
pair build: copy
stencil: none
bin: none
(4) pair oxdna2/xstk, perpetual, copy from (1)
attributes: half, newton on
pair build: copy
stencil: none
bin: none
(5) pair oxdna2/coaxstk, perpetual, copy from (1)
attributes: half, newton on
pair build: copy
stencil: none
bin: none
(6) pair oxdna2/dh, perpetual, copy from (1)
attributes: half, newton on
pair build: copy
stencil: none
bin: none
Per MPI rank memory allocation (min/avg/max) = 3.905 | 3.905 | 3.905 Mbytes
Step Temp E_pair E_mol TotEng Press
0 0 -1.6384018 0.037304147 -1.6010976 6.7769766e-05
10000 ekin = 1.0129637955345 | erot = 1.38476360245492 | epot = -41.9785360036034 | etot = -39.580808605614
20000 ekin = 1.75607695499755 | erot = 1.86620102096618 | epot = -41.9558496477911 | etot = -38.3335716718274
30000 ekin = 2.25878673688259 | erot = 2.60035907744898 | epot = -41.2869827431736 | etot = -36.427836928842
40000 ekin = 2.11105717434584 | erot = 3.16611217122799 | epot = -40.2762731426449 | etot = -34.9991037970711
50000 ekin = 3.05379083545187 | erot = 2.45050275456888 | epot = -40.3032957287998 | etot = -34.7990021387791
60000 ekin = 3.05154113920291 | erot = 2.5253957670889 | epot = -39.6800099891299 | etot = -34.1030730828381
70000 ekin = 3.23212209366464 | erot = 3.28914763888543 | epot = -40.1240667487278 | etot = -33.6027970161777
80000 ekin = 2.82623224207591 | erot = 2.91292948677732 | epot = -39.0983962904496 | etot = -33.3592345615963
90000 ekin = 3.05995314905694 | erot = 3.54299824110274 | epot = -39.164607034397 | etot = -32.5616556442373
100000 ekin = 3.46139546969892 | erot = 4.11704803295143 | epot = -38.512467983727 | etot = -30.9340244810767
100000 0.092303879 -1.5243833 0.043134544 -1.3481182 6.862374e-06
Loop time of 7.94086 on 1 procs for 100000 steps with 26 atoms
Performance: 108804.336 tau/day, 12593.094 timesteps/s
99.8% CPU use with 1 MPI tasks x no OpenMP threads
MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total
---------------------------------------------------------------
Pair | 7.1507 | 7.1507 | 7.1507 | 0.0 | 90.05
Bond | 0.14487 | 0.14487 | 0.14487 | 0.0 | 1.82
Neigh | 6.4e-05 | 6.4e-05 | 6.4e-05 | 0.0 | 0.00
Comm | 0.032259 | 0.032259 | 0.032259 | 0.0 | 0.41
Output | 0.003484 | 0.003484 | 0.003484 | 0.0 | 0.04
Modify | 0.59013 | 0.59013 | 0.59013 | 0.0 | 7.43
Other | | 0.01934 | | | 0.24
Nlocal: 26 ave 26 max 26 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Nghost: 0 ave 0 max 0 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Neighs: 323 ave 323 max 323 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Total # of neighbors = 323
Ave neighs/atom = 12.4231
Ave special neighs/atom = 5.07692
Neighbor list builds = 4
Dangerous builds = 0
#write_restart config.${number}.*
Total wall time: 0:00:07

274
examples/USER/cgdna/examples/oxDNA2/unique_bp/log.07Aug19.duplex4.8type.g++4 Normal file
View File

@ -0,0 +1,274 @@
LAMMPS (7 Aug 2019)
variable number equal 1
variable ofreq equal 10000
variable efreq equal 10000
variable ntype equal 8
variable T equal 0.1
units lj
dimension 3
newton on
boundary p p p
atom_style hybrid bond ellipsoid
atom_modify sort 0 1.0
# Pair interactions require lists of neighbours to be calculated
neighbor 1.0 bin
neigh_modify every 10 delay 0 check yes
read_data data.duplex4.8type
orthogonal box = (-20 -20 -20) to (20 20 20)
1 by 2 by 2 MPI processor grid
reading atoms ...
26 atoms
reading velocities ...
26 velocities
26 ellipsoids
scanning bonds ...
1 = max bonds/atom
reading bonds ...
24 bonds
2 = max # of 1-2 neighbors
2 = max # of 1-3 neighbors
4 = max # of 1-4 neighbors
6 = max # of special neighbors
special bonds CPU = 0.000232 secs
read_data CPU = 0.002447 secs
mass * 3.1575
group all type 1 8
26 atoms in group all
# oxDNA bond interactions - FENE backbone
bond_style oxdna2/fene
bond_coeff * 2.0 0.25 0.7564
# oxDNA pair interactions
pair_style hybrid/overlay oxdna2/excv oxdna2/stk oxdna2/hbond oxdna2/xstk oxdna2/coaxstk oxdna2/dh
pair_coeff * * oxdna2/excv 2.0 0.7 0.675 2.0 0.515 0.5 2.0 0.33 0.32
pair_coeff * * oxdna2/stk seqdep ${T} 1.3523 2.6717 6.0 0.4 0.9 0.32 0.75 1.3 0 0.8 0.9 0 0.95 0.9 0 0.95 2.0 0.65 2.0 0.65
pair_coeff * * oxdna2/stk seqdep 0.1 1.3523 2.6717 6.0 0.4 0.9 0.32 0.75 1.3 0 0.8 0.9 0 0.95 0.9 0 0.95 2.0 0.65 2.0 0.65
pair_coeff * * oxdna2/hbond seqdep 0.0 8.0 0.4 0.75 0.34 0.7 1.5 0 0.7 1.5 0 0.7 1.5 0 0.7 0.46 3.141592653589793 0.7 4.0 1.5707963267948966 0.45 4.0 1.5707963267948966 0.45
label loop
variable base loop ${ntype}
variable base loop 8
variable basemod equal ${base}%4
variable basemod equal 1%4
if "${basemod} == 1" then "variable comp equal ${base}+3" "pair_coeff ${base} ${comp} oxdna2/hbond seqdep 1.0678 8.0 0.4 0.75 0.34 0.7 1.5 0 0.7 1.5 0 0.7 1.5 0 0.7 0.46 3.141592653589793 0.7 4.0 1.5707963267948966 0.45 4.0 1.5707963267948966 0.45"
variable comp equal ${base}+3
variable comp equal 1+3
pair_coeff ${base} ${comp} oxdna2/hbond seqdep 1.0678 8.0 0.4 0.75 0.34 0.7 1.5 0 0.7 1.5 0 0.7 1.5 0 0.7 0.46 3.141592653589793 0.7 4.0 1.5707963267948966 0.45 4.0 1.5707963267948966 0.45
pair_coeff 1 ${comp} oxdna2/hbond seqdep 1.0678 8.0 0.4 0.75 0.34 0.7 1.5 0 0.7 1.5 0 0.7 1.5 0 0.7 0.46 3.141592653589793 0.7 4.0 1.5707963267948966 0.45 4.0 1.5707963267948966 0.45
pair_coeff 1 4 oxdna2/hbond seqdep 1.0678 8.0 0.4 0.75 0.34 0.7 1.5 0 0.7 1.5 0 0.7 1.5 0 0.7 0.46 3.141592653589793 0.7 4.0 1.5707963267948966 0.45 4.0 1.5707963267948966 0.45
if "${basemod} == 2" then "variable comp equal ${base}+1" "pair_coeff ${base} ${comp} oxdna2/hbond seqdep 1.0678 8.0 0.4 0.75 0.34 0.7 1.5 0 0.7 1.5 0 0.7 1.5 0 0.7 0.46 3.141592653589793 0.7 4.0 1.5707963267948966 0.45 4.0 1.5707963267948966 0.45"
next base
jump in.duplex4.8type loop
variable base loop ${ntype}
variable base loop 8
variable basemod equal ${base}%4
variable basemod equal 2%4
if "${basemod} == 1" then "variable comp equal ${base}+3" "pair_coeff ${base} ${comp} oxdna2/hbond seqdep 1.0678 8.0 0.4 0.75 0.34 0.7 1.5 0 0.7 1.5 0 0.7 1.5 0 0.7 0.46 3.141592653589793 0.7 4.0 1.5707963267948966 0.45 4.0 1.5707963267948966 0.45"
if "${basemod} == 2" then "variable comp equal ${base}+1" "pair_coeff ${base} ${comp} oxdna2/hbond seqdep 1.0678 8.0 0.4 0.75 0.34 0.7 1.5 0 0.7 1.5 0 0.7 1.5 0 0.7 0.46 3.141592653589793 0.7 4.0 1.5707963267948966 0.45 4.0 1.5707963267948966 0.45"
variable comp equal ${base}+1
variable comp equal 2+1
pair_coeff ${base} ${comp} oxdna2/hbond seqdep 1.0678 8.0 0.4 0.75 0.34 0.7 1.5 0 0.7 1.5 0 0.7 1.5 0 0.7 0.46 3.141592653589793 0.7 4.0 1.5707963267948966 0.45 4.0 1.5707963267948966 0.45
pair_coeff 2 ${comp} oxdna2/hbond seqdep 1.0678 8.0 0.4 0.75 0.34 0.7 1.5 0 0.7 1.5 0 0.7 1.5 0 0.7 0.46 3.141592653589793 0.7 4.0 1.5707963267948966 0.45 4.0 1.5707963267948966 0.45
pair_coeff 2 3 oxdna2/hbond seqdep 1.0678 8.0 0.4 0.75 0.34 0.7 1.5 0 0.7 1.5 0 0.7 1.5 0 0.7 0.46 3.141592653589793 0.7 4.0 1.5707963267948966 0.45 4.0 1.5707963267948966 0.45
next base
jump in.duplex4.8type loop
variable base loop ${ntype}
variable base loop 8
variable basemod equal ${base}%4
variable basemod equal 3%4
if "${basemod} == 1" then "variable comp equal ${base}+3" "pair_coeff ${base} ${comp} oxdna2/hbond seqdep 1.0678 8.0 0.4 0.75 0.34 0.7 1.5 0 0.7 1.5 0 0.7 1.5 0 0.7 0.46 3.141592653589793 0.7 4.0 1.5707963267948966 0.45 4.0 1.5707963267948966 0.45"
if "${basemod} == 2" then "variable comp equal ${base}+1" "pair_coeff ${base} ${comp} oxdna2/hbond seqdep 1.0678 8.0 0.4 0.75 0.34 0.7 1.5 0 0.7 1.5 0 0.7 1.5 0 0.7 0.46 3.141592653589793 0.7 4.0 1.5707963267948966 0.45 4.0 1.5707963267948966 0.45"
next base
jump in.duplex4.8type loop
variable base loop ${ntype}
variable base loop 8
variable basemod equal ${base}%4
variable basemod equal 4%4
if "${basemod} == 1" then "variable comp equal ${base}+3" "pair_coeff ${base} ${comp} oxdna2/hbond seqdep 1.0678 8.0 0.4 0.75 0.34 0.7 1.5 0 0.7 1.5 0 0.7 1.5 0 0.7 0.46 3.141592653589793 0.7 4.0 1.5707963267948966 0.45 4.0 1.5707963267948966 0.45"
if "${basemod} == 2" then "variable comp equal ${base}+1" "pair_coeff ${base} ${comp} oxdna2/hbond seqdep 1.0678 8.0 0.4 0.75 0.34 0.7 1.5 0 0.7 1.5 0 0.7 1.5 0 0.7 0.46 3.141592653589793 0.7 4.0 1.5707963267948966 0.45 4.0 1.5707963267948966 0.45"
next base
jump in.duplex4.8type loop
variable base loop ${ntype}
variable base loop 8
variable basemod equal ${base}%4
variable basemod equal 5%4
if "${basemod} == 1" then "variable comp equal ${base}+3" "pair_coeff ${base} ${comp} oxdna2/hbond seqdep 1.0678 8.0 0.4 0.75 0.34 0.7 1.5 0 0.7 1.5 0 0.7 1.5 0 0.7 0.46 3.141592653589793 0.7 4.0 1.5707963267948966 0.45 4.0 1.5707963267948966 0.45"
variable comp equal ${base}+3
variable comp equal 5+3
pair_coeff ${base} ${comp} oxdna2/hbond seqdep 1.0678 8.0 0.4 0.75 0.34 0.7 1.5 0 0.7 1.5 0 0.7 1.5 0 0.7 0.46 3.141592653589793 0.7 4.0 1.5707963267948966 0.45 4.0 1.5707963267948966 0.45
pair_coeff 5 ${comp} oxdna2/hbond seqdep 1.0678 8.0 0.4 0.75 0.34 0.7 1.5 0 0.7 1.5 0 0.7 1.5 0 0.7 0.46 3.141592653589793 0.7 4.0 1.5707963267948966 0.45 4.0 1.5707963267948966 0.45
pair_coeff 5 8 oxdna2/hbond seqdep 1.0678 8.0 0.4 0.75 0.34 0.7 1.5 0 0.7 1.5 0 0.7 1.5 0 0.7 0.46 3.141592653589793 0.7 4.0 1.5707963267948966 0.45 4.0 1.5707963267948966 0.45
if "${basemod} == 2" then "variable comp equal ${base}+1" "pair_coeff ${base} ${comp} oxdna2/hbond seqdep 1.0678 8.0 0.4 0.75 0.34 0.7 1.5 0 0.7 1.5 0 0.7 1.5 0 0.7 0.46 3.141592653589793 0.7 4.0 1.5707963267948966 0.45 4.0 1.5707963267948966 0.45"
next base
jump in.duplex4.8type loop
variable base loop ${ntype}
variable base loop 8
variable basemod equal ${base}%4
variable basemod equal 6%4
if "${basemod} == 1" then "variable comp equal ${base}+3" "pair_coeff ${base} ${comp} oxdna2/hbond seqdep 1.0678 8.0 0.4 0.75 0.34 0.7 1.5 0 0.7 1.5 0 0.7 1.5 0 0.7 0.46 3.141592653589793 0.7 4.0 1.5707963267948966 0.45 4.0 1.5707963267948966 0.45"
if "${basemod} == 2" then "variable comp equal ${base}+1" "pair_coeff ${base} ${comp} oxdna2/hbond seqdep 1.0678 8.0 0.4 0.75 0.34 0.7 1.5 0 0.7 1.5 0 0.7 1.5 0 0.7 0.46 3.141592653589793 0.7 4.0 1.5707963267948966 0.45 4.0 1.5707963267948966 0.45"
variable comp equal ${base}+1
variable comp equal 6+1
pair_coeff ${base} ${comp} oxdna2/hbond seqdep 1.0678 8.0 0.4 0.75 0.34 0.7 1.5 0 0.7 1.5 0 0.7 1.5 0 0.7 0.46 3.141592653589793 0.7 4.0 1.5707963267948966 0.45 4.0 1.5707963267948966 0.45
pair_coeff 6 ${comp} oxdna2/hbond seqdep 1.0678 8.0 0.4 0.75 0.34 0.7 1.5 0 0.7 1.5 0 0.7 1.5 0 0.7 0.46 3.141592653589793 0.7 4.0 1.5707963267948966 0.45 4.0 1.5707963267948966 0.45
pair_coeff 6 7 oxdna2/hbond seqdep 1.0678 8.0 0.4 0.75 0.34 0.7 1.5 0 0.7 1.5 0 0.7 1.5 0 0.7 0.46 3.141592653589793 0.7 4.0 1.5707963267948966 0.45 4.0 1.5707963267948966 0.45
next base
jump in.duplex4.8type loop
variable base loop ${ntype}
variable base loop 8
variable basemod equal ${base}%4
variable basemod equal 7%4
if "${basemod} == 1" then "variable comp equal ${base}+3" "pair_coeff ${base} ${comp} oxdna2/hbond seqdep 1.0678 8.0 0.4 0.75 0.34 0.7 1.5 0 0.7 1.5 0 0.7 1.5 0 0.7 0.46 3.141592653589793 0.7 4.0 1.5707963267948966 0.45 4.0 1.5707963267948966 0.45"
if "${basemod} == 2" then "variable comp equal ${base}+1" "pair_coeff ${base} ${comp} oxdna2/hbond seqdep 1.0678 8.0 0.4 0.75 0.34 0.7 1.5 0 0.7 1.5 0 0.7 1.5 0 0.7 0.46 3.141592653589793 0.7 4.0 1.5707963267948966 0.45 4.0 1.5707963267948966 0.45"
next base
jump in.duplex4.8type loop
variable base loop ${ntype}
variable base loop 8
variable basemod equal ${base}%4
variable basemod equal 8%4
if "${basemod} == 1" then "variable comp equal ${base}+3" "pair_coeff ${base} ${comp} oxdna2/hbond seqdep 1.0678 8.0 0.4 0.75 0.34 0.7 1.5 0 0.7 1.5 0 0.7 1.5 0 0.7 0.46 3.141592653589793 0.7 4.0 1.5707963267948966 0.45 4.0 1.5707963267948966 0.45"
if "${basemod} == 2" then "variable comp equal ${base}+1" "pair_coeff ${base} ${comp} oxdna2/hbond seqdep 1.0678 8.0 0.4 0.75 0.34 0.7 1.5 0 0.7 1.5 0 0.7 1.5 0 0.7 0.46 3.141592653589793 0.7 4.0 1.5707963267948966 0.45 4.0 1.5707963267948966 0.45"
next base
jump in.duplex4.8type loop
pair_coeff * * oxdna2/xstk 47.5 0.575 0.675 0.495 0.655 2.25 0.791592653589793 0.58 1.7 1.0 0.68 1.7 1.0 0.68 1.5 0 0.65 1.7 0.875 0.68 1.7 0.875 0.68
pair_coeff * * oxdna2/coaxstk 58.5 0.4 0.6 0.22 0.58 2.0 2.891592653589793 0.65 1.3 0 0.8 0.9 0 0.95 0.9 0 0.95 40.0 3.116592653589793
pair_coeff * * oxdna2/dh ${T} 0.2 0.815
pair_coeff * * oxdna2/dh 0.1 0.2 0.815
# Langevin dynamics
fix 1 all nve/asphere
fix 2 all langevin ${T} ${T} 25.0 457145 angmom 10
fix 2 all langevin 0.1 ${T} 25.0 457145 angmom 10
fix 2 all langevin 0.1 0.1 25.0 457145 angmom 10
timestep 1e-4
#comm_style tiled
#fix 3 all balance 10000 1.1 rcb
#compute mol all chunk/atom molecule
#compute mychunk all vcm/chunk mol
#fix 4 all ave/time 10000 1 10000 c_mychunk[1] c_mychunk[2] c_mychunk[3] file vcm.txt mode vector
#dump pos all xyz ${ofreq} traj.${number}.xyz
compute quat all property/atom quatw quati quatj quatk
#dump quat all custom ${ofreq} quat.${number}.txt id c_quat[1] c_quat[2] c_quat[3] c_quat[4]
#dump_modify quat sort id
#dump_modify quat format line "%d %13.6le %13.6le %13.6le %13.6le"
compute erot all erotate/asphere
compute ekin all ke
compute epot all pe
variable erot equal c_erot
variable ekin equal c_ekin
variable epot equal c_epot
variable etot equal c_erot+c_ekin+c_epot
fix 5 all print ${efreq} "$(step) ekin = ${ekin} | erot = ${erot} | epot = ${epot} | etot = ${etot}" screen yes
fix 5 all print 10000 "$(step) ekin = ${ekin} | erot = ${erot} | epot = ${epot} | etot = ${etot}" screen yes
dump out all custom ${ofreq} out.${number}.txt id mol type x y z ix iy iz c_quat[1] c_quat[2] c_quat[3] c_quat[4] vx vy vz
dump out all custom 10000 out.${number}.txt id mol type x y z ix iy iz c_quat[1] c_quat[2] c_quat[3] c_quat[4] vx vy vz
dump out all custom 10000 out.1.txt id mol type x y z ix iy iz c_quat[1] c_quat[2] c_quat[3] c_quat[4] vx vy vz
dump_modify out sort id
dump_modify out format line "%d %d %d %13.6le %13.6le %13.6le %d %d %d %13.6le %13.6le %13.6le %13.6le %13.6le %13.6le %13.6le "
#restart 10000 config0_restart config1_restart
run 100000
Neighbor list info ...
update every 10 steps, delay 0 steps, check yes
max neighbors/atom: 2000, page size: 100000
master list distance cutoff = 4.63899
ghost atom cutoff = 4.63899
binsize = 2.31949, bins = 18 18 18
6 neighbor lists, perpetual/occasional/extra = 6 0 0
(1) pair oxdna2/excv, perpetual
attributes: half, newton on
pair build: half/bin/newton
stencil: half/bin/3d/newton
bin: standard
(2) pair oxdna2/stk, perpetual, copy from (1)
attributes: half, newton on
pair build: copy
stencil: none
bin: none
(3) pair oxdna2/hbond, perpetual, copy from (1)
attributes: half, newton on
pair build: copy
stencil: none
bin: none
(4) pair oxdna2/xstk, perpetual, copy from (1)
attributes: half, newton on
pair build: copy
stencil: none
bin: none
(5) pair oxdna2/coaxstk, perpetual, copy from (1)
attributes: half, newton on
pair build: copy
stencil: none
bin: none
(6) pair oxdna2/dh, perpetual, copy from (1)
attributes: half, newton on
pair build: copy
stencil: none
bin: none
Per MPI rank memory allocation (min/avg/max) = 9.381 | 9.563 | 9.746 Mbytes
Step Temp E_pair E_mol TotEng Press
0 0 -1.6384018 0.037304147 -1.6010976 6.7769766e-05
10000 ekin = 1.20524984175816 | erot = 1.26791229621008 | epot = -41.9183873815797 | etot = -39.4452252436115
20000 ekin = 1.6085654035548 | erot = 2.6414088050897 | epot = -41.0106168994546 | etot = -36.7606426908101
30000 ekin = 2.82110163819891 | erot = 3.06373823252419 | epot = -40.9348041675026 | etot = -35.0499642967795
40000 ekin = 2.35553603700794 | erot = 3.60405945883761 | epot = -39.5418687902286 | etot = -33.5822732943831
50000 ekin = 3.06061403125833 | erot = 3.35548669087246 | epot = -40.3236585928169 | etot = -33.9075578706861
60000 ekin = 2.7347216991605 | erot = 4.23926242244084 | epot = -39.5379959923263 | etot = -32.564011870725
70000 ekin = 3.46562604991873 | erot = 2.8521307045909 | epot = -38.7237000550356 | etot = -32.4059433005259
80000 ekin = 3.17815543267806 | erot = 3.11078099027451 | epot = -39.7525036398366 | etot = -33.463567216884
90000 ekin = 3.51635117668857 | erot = 2.58384069017333 | epot = -39.4762533092413 | etot = -33.3760614423795
100000 ekin = 3.64218174280962 | erot = 2.88607181210736 | epot = -37.6002449519119 | etot = -31.0719913969949
100000 0.097124846 -1.5164679 0.070304678 -1.3060794 -2.9111933e-05
Loop time of 5.64462 on 4 procs for 100000 steps with 26 atoms
Performance: 153066.031 tau/day, 17715.976 timesteps/s
99.0% CPU use with 4 MPI tasks x no OpenMP threads
MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total
---------------------------------------------------------------
Pair | 0.069668 | 2.1601 | 4.2189 | 138.7 | 38.27
Bond | 0.008542 | 0.052544 | 0.095929 | 18.7 | 0.93
Neigh | 4.8e-05 | 5.6e-05 | 6.4e-05 | 0.0 | 0.00
Comm | 0.90402 | 3.1443 | 5.4155 | 124.9 | 55.70
Output | 0.001003 | 0.0011317 | 0.001468 | 0.6 | 0.02
Modify | 0.021098 | 0.20024 | 0.38154 | 39.2 | 3.55
Other | | 0.0863 | | | 1.53
Nlocal: 6.5 ave 13 max 0 min
Histogram: 2 0 0 0 0 0 0 0 0 2
Nghost: 17.5 ave 22 max 13 min
Histogram: 2 0 0 0 0 0 0 0 0 2
Neighs: 80.25 ave 167 max 0 min
Histogram: 2 0 0 0 0 0 0 0 0 2
Total # of neighbors = 321
Ave neighs/atom = 12.3462
Ave special neighs/atom = 5.07692
Neighbor list builds = 3
Dangerous builds = 0
#write_restart config.${number}.*
Total wall time: 0:00:05

96
examples/USER/cgdna/examples/oxRNA2/duplex4/data.duplex4 Normal file
View File

@ -0,0 +1,96 @@
# LAMMPS data file
16 atoms
16 ellipsoids
13 bonds
4 atom types
1 bond types
# System size
-20.000000 20.000000 xlo xhi
-20.000000 20.000000 ylo yhi
-20.000000 20.000000 zlo zhi
Masses
1 3.1575
2 3.1575
3 3.1575
4 3.1575
# Atom-ID, type, position, molecule-ID, ellipsoid flag, density
Atoms
1 1 -6.000000000000001e-01 0.000000000000000e+00 0.000000000000000e+00 1 1 1
2 2 -4.860249842674776e-01 -3.518234140414736e-01 3.897628551303122e-01 1 1 1
3 3 -1.874009511073395e-01 -5.699832309147915e-01 7.795257102606244e-01 1 1 1
4 4 1.824198365552941e-01 -5.715968887521518e-01 1.169288565390937e+00 1 1 1
5 1 4.829362784135484e-01 -3.560513319622209e-01 1.559051420521249e+00 1 1 1
6 2 5.999771538385027e-01 -5.235921299024461e-03 1.948814275651561e+00 1 1 1
7 3 4.890766774371325e-01 3.475687034056071e-01 2.338577130781873e+00 1 1 1
8 4 1.923677943514057e-01 5.683261666476170e-01 2.728339985912185e+00 1 1 1
9 1 -1.923677943514057e-01 -5.683261666476170e-01 2.728339985912185e+00 2 1 1
10 2 -4.890766774371324e-01 -3.475687034056071e-01 2.338577130781873e+00 2 1 1
11 3 -5.999771538385025e-01 5.235921299024461e-03 1.948814275651561e+00 2 1 1
12 4 -4.829362784135481e-01 3.560513319622207e-01 1.559051420521249e+00 2 1 1
13 1 -1.824198365552940e-01 5.715968887521514e-01 1.169288565390936e+00 2 1 1
14 2 1.874009511073395e-01 5.699832309147912e-01 7.795257102606241e-01 2 1 1
15 3 4.860249842674773e-01 3.518234140414733e-01 3.897628551303119e-01 2 1 1
16 4 5.999999999999995e-01 -3.330669073875470e-17 -3.330669073875470e-16 2 1 1
# Atom-ID, translational velocity, angular momentum
Velocities
1 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00
2 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00
3 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00
4 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00
5 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00
6 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00
7 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00
8 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00
9 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00
10 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00
11 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00
12 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00
13 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00
14 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00
15 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00
16 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00
# Atom-ID, shape, quaternion
Ellipsoids
1 1.173984503142341e+00 1.173984503142341e+00 1.173984503142341e+00 1.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00
2 1.173984503142341e+00 1.173984503142341e+00 1.173984503142341e+00 9.513258223252946e-01 0.000000000000000e+00 0.000000000000000e+00 3.081869234362515e-01
3 1.173984503142341e+00 1.173984503142341e+00 1.173984503142341e+00 8.100416404457962e-01 0.000000000000000e+00 0.000000000000000e+00 5.863723567357894e-01
4 1.173984503142341e+00 1.173984503142341e+00 1.173984503142341e+00 5.899012371043606e-01 0.000000000000000e+00 0.000000000000000e+00 8.074754054847398e-01
5 1.173984503142341e+00 1.173984503142341e+00 1.173984503142341e+00 3.123349185122326e-01 0.000000000000000e+00 0.000000000000000e+00 9.499720515246527e-01
6 1.173984503142341e+00 1.173984503142341e+00 1.173984503142341e+00 4.363309284746654e-03 0.000000000000000e+00 0.000000000000000e+00 9.999904807207346e-01
7 1.173984503142341e+00 1.173984503142341e+00 1.173984503142341e+00 -3.040330609254902e-01 0.000000000000000e+00 0.000000000000000e+00 9.526614812535865e-01
8 1.173984503142341e+00 1.173984503142341e+00 1.173984503142341e+00 5.828323126827837e-01 0.000000000000000e+00 0.000000000000000e+00 -8.125924533816677e-01
9 1.173984503142341e+00 1.173984503142341e+00 1.173984503142341e+00 0.000000000000000e+00 8.125924533816681e-01 5.828323126827832e-01 -0.000000000000000e+00
10 1.173984503142341e+00 1.173984503142341e+00 1.173984503142341e+00 0.000000000000000e+00 9.526614812535864e-01 3.040330609254902e-01 0.000000000000000e+00
11 1.173984503142341e+00 1.173984503142341e+00 1.173984503142341e+00 0.000000000000000e+00 9.999904807207346e-01 -4.363309284746654e-03 0.000000000000000e+00
12 1.173984503142341e+00 1.173984503142341e+00 1.173984503142341e+00 0.000000000000000e+00 9.499720515246526e-01 -3.123349185122325e-01 0.000000000000000e+00
13 1.173984503142341e+00 1.173984503142341e+00 1.173984503142341e+00 0.000000000000000e+00 8.074754054847402e-01 -5.899012371043603e-01 0.000000000000000e+00
14 1.173984503142341e+00 1.173984503142341e+00 1.173984503142341e+00 0.000000000000000e+00 5.863723567357898e-01 -8.100416404457959e-01 0.000000000000000e+00
15 1.173984503142341e+00 1.173984503142341e+00 1.173984503142341e+00 -0.000000000000000e+00 -3.081869234362514e-01 9.513258223252948e-01 0.000000000000000e+00
16 1.173984503142341e+00 1.173984503142341e+00 1.173984503142341e+00 -0.000000000000000e+00 2.775557561562893e-17 1.000000000000000e+00 -0.000000000000000e+00
# Bond topology
Bonds
1 1 1 2
2 1 2 3
3 1 3 4
4 1 4 5
5 1 5 6
6 1 6 7
7 1 7 8
8 1 9 10
9 1 10 11
10 1 11 12
11 1 13 14
12 1 14 15
13 1 15 16

80
examples/USER/cgdna/examples/oxRNA2/duplex4/in.duplex4 Normal file
View File

@ -0,0 +1,80 @@
variable number equal 4
variable ofreq equal 1000
variable efreq equal 1000
variable T equal 0.1
units lj
dimension 3
newton off
boundary p p p
atom_style hybrid bond ellipsoid
atom_modify sort 0 1.0
# Pair interactions require lists of neighbours to be calculated
neighbor 1.0 bin
neigh_modify every 1 delay 0 check yes
read_data data.duplex4
set atom * mass 3.1575
group all type 1 4
# oxRNA2 bond interactions - FENE backbone
bond_style oxrna2/fene
bond_coeff * 2.0 0.25 0.761070781051
# oxRNA2 pair interactions
pair_style hybrid/overlay oxrna2/excv oxrna2/stk oxrna2/hbond oxrna2/xstk oxrna2/coaxstk oxrna2/dh
pair_coeff * * oxrna2/excv 2.0 0.7 0.675 2.0 0.515 0.5 2.0 0.33 0.32
pair_coeff * * oxrna2/stk seqdep ${T} 1.40206 2.77 6.0 0.43 0.93 0.35 0.78 0.9 0 0.95 0.9 0 0.95 1.3 0 0.8 1.3 0 0.8 2.0 0.65 2.0 0.65
pair_coeff * * oxrna2/hbond seqdep 0.0 8.0 0.4 0.75 0.34 0.7 1.5 0 0.7 1.5 0 0.7 1.5 0 0.7 0.46 3.141592653589793 0.7 4.0 1.5707963267948966 0.45 4.0 1.5707963267948966 0.45
pair_coeff 1 4 oxrna2/hbond seqdep 0.870439 8.0 0.4 0.75 0.34 0.7 1.5 0 0.7 1.5 0 0.7 1.5 0 0.7 0.46 3.141592653589793 0.7 4.0 1.5707963267948966 0.45 4.0 1.5707963267948966 0.45
pair_coeff 2 3 oxrna2/hbond seqdep 0.870439 8.0 0.4 0.75 0.34 0.7 1.5 0 0.7 1.5 0 0.7 1.5 0 0.7 0.46 3.141592653589793 0.7 4.0 1.5707963267948966 0.45 4.0 1.5707963267948966 0.45
pair_coeff 3 4 oxrna2/hbond seqdep 0.870439 8.0 0.4 0.75 0.34 0.7 1.5 0 0.7 1.5 0 0.7 1.5 0 0.7 0.46 3.141592653589793 0.7 4.0 1.5707963267948966 0.45 4.0 1.5707963267948966 0.45
pair_coeff * * oxrna2/xstk 59.9626 0.5 0.6 0.42 0.58 2.25 0.505 0.58 1.7 1.266 0.68 1.7 1.266 0.68 1.7 0.309 0.68 1.7 0.309 0.68
pair_coeff * * oxrna2/coaxstk 80 0.5 0.6 0.42 0.58 2.0 2.592 0.65 1.3 0.151 0.8 0.9 0.685 0.95 0.9 0.685 0.95 2.0 -0.65 2.0 -0.65
pair_coeff * * oxrna2/dh ${T} 0.5 1.02455
# NVE ensemble
#fix 1 all nve/dotc/langevin 0.1 0.1 0.03 457145 angmom 10
#fix 1 all nve/dot
fix 1 all nve/asphere
timestep 1e-5
#comm_style tiled
#fix 3 all balance 10000 1.1 rcb
#compute mol all chunk/atom molecule
#compute mychunk all vcm/chunk mol
#fix 4 all ave/time 10000 1 10000 c_mychunk[1] c_mychunk[2] c_mychunk[3] file vcm.txt mode vector
#dump pos all xyz ${ofreq} traj.${number}.xyz
#compute quat all property/atom quatw quati quatj quatk
#dump quat all custom ${ofreq} quat.${number}.txt id c_quat[1] c_quat[2] c_quat[3] c_quat[4]
#dump_modify quat sort id
#dump_modify quat format line "%d %13.6le %13.6le %13.6le %13.6le"
compute erot all erotate/asphere
compute ekin all ke
compute epot all pe
variable erot equal c_erot
variable ekin equal c_ekin
variable epot equal c_epot
variable etot equal c_erot+c_ekin+c_epot
fix 5 all print ${efreq} "$(step) ekin = ${ekin} | erot = ${erot} | epot = ${epot} | etot = ${etot}" screen yes
#dump out all custom ${ofreq} out.${number}.txt id x y z vx vy vz fx fy fz tqx tqy tqz
#dump_modify out sort id
#dump_modify out format line "%d %13.9f %13.9f %13.9f %13.9f %13.9f %13.9f %13.9f %13.9f %13.9f %13.9f %13.9f %13.9f"
run 1000000
#write_restart config.${number}.*

1174
examples/USER/cgdna/examples/oxRNA2/duplex4/log.30Oct19.duplex4.g++.1 Normal file
View File

File diff suppressed because it is too large Load Diff

1174
examples/USER/cgdna/examples/oxRNA2/duplex4/log.30Oct19.duplex4.g++.4 Normal file
View File

File diff suppressed because it is too large Load Diff