Added unique base pairing
This commit is contained in:
Oliver Henrich
@ -55,6 +55,23 @@ 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
|
||||
When a large number of atom types is used, this feature allows
|
||||
quasi-unique base pairing between two individual nucleotides.
|
||||
|
||||
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
|
||||
|
||||
/******************************************************************************/
|
||||
|
||||
/examples/oxRNA2/duplex4
|
||||
|
||||
This is the duplex2 run with the oxRNA2 force field instead of the oxDNA or
|
||||
|
||||
130
examples/USER/cgdna/examples/oxDNA2/unique_bp/data.duplex4.4type
Normal file
130
examples/USER/cgdna/examples/oxDNA2/unique_bp/data.duplex4.4type
Normal 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
130
examples/USER/cgdna/examples/oxDNA2/unique_bp/data.duplex4.8type
Normal 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
828
examples/USER/cgdna/examples/oxDNA2/unique_bp/generate_unique.py
Normal 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)
|
||||
@ -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}.*
|
||||
|
||||
|
||||
@ -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}.*
|
||||
|
||||
|
||||
@ -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
|
||||
@ -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
|
||||
@ -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
|
||||
@ -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
|
||||
Reference in New Issue
Block a user