jianmu
/
lammps
forked from lijiext/lammps
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Merge pull request #1823 from dsbolin/rigid-gravity 

Adds an example using new options in fix rigid and fix gravity
This commit is contained in:
Axel Kohlmeyer 2020-01-08 09:15:55 -05:00 committed by GitHub
parent a88a00dbbb 7bc8c8e9d8
commit 1cd6d302bb
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
43 changed files with 8483 additions and 214 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

79
doc/src/bond_oxdna.rst
View File

@ -6,6 +6,9 @@ bond_style oxdna/fene command
bond_style oxdna2/fene command
==============================
bond_style oxrna2/fene command
==============================
Syntax
""""""
@ -16,6 +19,8 @@ Syntax
bond_style oxdna2/fene
bond_style oxrna2/fene
Examples
""""""""
@ -28,10 +33,13 @@ Examples
bond_style oxdna2/fene
bond_coeff * 2.0 0.25 0.7564
bond_style oxrna2/fene
bond_coeff \* 2.0 0.25 0.76107
Description
"""""""""""
The *oxdna/fene* and *oxdna2/fene* bond styles use the potential
The *oxdna/fene* , *oxdna2/fene* and *oxrna2/fene* bond styles use the potential
.. math::
@ -39,9 +47,9 @@ The *oxdna/fene* and *oxdna2/fene* bond styles use the potential
to define a modified finite extensible nonlinear elastic (FENE)
potential :ref:`(Ouldridge) <oxdna_fene>` to model the connectivity of the
phosphate backbone in the oxDNA force field for coarse-grained
modelling of DNA.
potential :ref:`(Ouldridge) <Ouldridge0>` to model the connectivity of the
phosphate backbone in the oxDNA/oxRNA force field for coarse-grained
modelling of DNA/RNA.
The following coefficients must be defined for the bond type via the
:doc:`bond\_coeff <bond_coeff>` command as given in the above example, or
@ -57,27 +65,36 @@ commands:
The oxDNA bond style has to be used together with the
corresponding oxDNA pair styles for excluded volume interaction
*oxdna/excv*\ , stacking *oxdna/stk*\ , cross-stacking *oxdna/xstk* and
*oxdna/excv* , stacking *oxdna/stk* , cross-stacking *oxdna/xstk* and
coaxial stacking interaction *oxdna/coaxstk* as well as
hydrogen-bonding interaction *oxdna/hbond* (see also documentation of
:doc:`pair\_style oxdna/excv <pair_oxdna>`). For the oxDNA2
:ref:`(Snodin) <oxdna2>` bond style the analogous pair styles and an
additional Debye-Hueckel pair style *oxdna2/dh* have to be defined.
:ref:`(Snodin) <Snodin0>` bond style the analogous pair styles
*oxdna2/excv* , *oxdna2/stk* , *oxdna2/xstk* , *oxdna2/coaxstk* ,
*oxdna2/hbond* and an additional Debye-Hueckel pair style
*oxdna2/dh* have to be defined. The same applies to the oxRNA2
:ref:`(Sulc1) <Sulc01>` styles.
The coefficients in the above example have to be kept fixed and cannot
be changed without reparameterizing the entire model.
Example input and data files for DNA duplexes can be found in
examples/USER/cgdna/examples/oxDNA/ and /oxDNA2/. A simple python
setup tool which creates single straight or helical DNA strands, DNA
duplexes or arrays of DNA duplexes can be found in
Example input and data files for DNA and RNA duplexes can be found in
examples/USER/cgdna/examples/oxDNA/ , /oxDNA2/ and /oxRNA2/. A simple python
setup tool which creates single straight or helical DNA strands, DNA/RNA
duplexes or arrays of DNA/RNA duplexes can be found in
examples/USER/cgdna/util/.
Please cite :ref:`(Henrich) <Henrich2>` and the relevant oxDNA articles in
any publication that uses this implementation. The article contains
more information on the model, the structure of the input file, the
setup tool and the performance of the LAMMPS-implementation of oxDNA.
The preprint version of the article can be found
Please cite :ref:`(Henrich) <Henrich0>` in any publication that uses
this implementation. The article contains general information
on the model, its implementation and performance as well as the structure of
the data and input file. The preprint version of the article can be found
`here <PDF/USER-CGDNA.pdf>`_.
Please cite also the relevant oxDNA/oxRNA publications. These are
:ref:`(Ouldridge) <Ouldridge0>` and
:ref:`(Ouldridge-DPhil) <Ouldridge-DPhil0>` for oxDNA,
:ref:`(Snodin) <Snodin0>` for oxDNA2,
:ref:`(Sulc1) <Sulc01>` for oxRNA2
and for sequence-specific hydrogen-bonding and stacking interactions
:ref:`(Sulc2) <Sulc02>`.
----------
@ -94,32 +111,36 @@ USER-CGDNA package and the MOLECULE and ASPHERE package. See the
Related commands
""""""""""""""""
:doc:`pair\_style oxdna/excv <pair_oxdna>`, :doc:`pair\_style oxdna2/excv <pair_oxdna2>`, :doc:`fix nve/dotc/langevin <fix_nve_dotc_langevin>`,
:doc:`bond\_coeff <bond_coeff>`
:doc:`pair\_style oxdna/excv <pair_oxdna>`, :doc:`pair\_style oxdna2/excv <pair_oxdna2>`, :doc:`pair\_style oxrna2/excv <pair_oxrna2>`,
:doc:`bond\_coeff <bond_coeff>`, :doc:`fix nve/dotc/langevin <fix_nve_dotc_langevin>`
**Default:** none
**Default:**
none
----------
.. _Henrich0:
.. _Henrich2:
**(Henrich)** O. Henrich, Y. A. Gutierrez-Fosado, T. Curk, T. E. Ouldridge, Eur. Phys. J. E 41, 57 (2018).
.. _Ouldridge-DPhil0:
**(Ouldridge-DPhil)** T.E. Ouldridge, Coarse-grained modelling of DNA and DNA self-assembly, DPhil. University of Oxford (2011).
**(Henrich)** O. Henrich, Y. A. Gutierrez-Fosado, T. Curk,
T. E. Ouldridge, Eur. Phys. J. E 41, 57 (2018).
.. _Ouldridge0:
.. _oxdna\_fene:
**(Ouldridge)** T.E. Ouldridge, A.A. Louis, J.P.K. Doye, J. Chem. Phys. 134, 085101 (2011).
.. _Snodin0:
**(Snodin)** B.E. Snodin, F. Randisi, M. Mosayebi, et al., J. Chem. Phys. 142, 234901 (2015).
**(Ouldridge)** T.E. Ouldridge, A.A. Louis, J.P.K. Doye,
J. Chem. Phys. 134, 085101 (2011).
.. _Sulc01:
.. _oxdna2:
**(Sulc1)** P. Sulc, F. Romano, T. E. Ouldridge, et al., J. Chem. Phys. 140, 235102 (2014).
.. _Sulc02:
**(Snodin)** B.E. Snodin, F. Randisi, M. Mosayebi, et al.,
J. Chem. Phys. 142, 234901 (2015).
**(Sulc2)** P. Sulc, F. Romano, T.E. Ouldridge, L. Rovigatti, J.P.K. Doye, A.A. Louis, J. Chem. Phys. 137, 135101 (2012).

16
doc/src/fix_nve_dot.rst
View File

@ -26,11 +26,11 @@ Examples
Description
"""""""""""
Apply a rigid-body integrator as described in :ref:`(Davidchack) <Davidchack1>`
Apply a rigid-body integrator as described in :ref:`(Davidchack) <Davidchack4>`
to a group of atoms, but without Langevin dynamics.
This command performs Molecular dynamics (MD)
via a velocity-Verlet algorithm and an evolution operator that rotates
the quaternion degrees of freedom, similar to the scheme outlined in :ref:`(Miller) <Miller1>`.
the quaternion degrees of freedom, similar to the scheme outlined in :ref:`(Miller) <Miller4>`.
This command is the equivalent of the :doc:`fix nve/dotc/langevin <fix_nve_dotc_langevin>`
without damping and noise and can be used to determine the stability range
@ -40,7 +40,7 @@ The command is equivalent to the :doc:`fix nve <fix_nve>`.
The particles are always considered to have a finite size.
An example input file can be found in /examples/USER/cgdna/examples/duplex1/.
Further details of the implementation and stability of the integrator are contained in :ref:`(Henrich) <Henrich3>`.
Further details of the implementation and stability of the integrator are contained in :ref:`(Henrich) <Henrich4>`.
The preprint version of the article can be found `here <PDF/USER-CGDNA.pdf>`_.
@ -66,19 +66,15 @@ Related commands
----------
.. _Davidchack1:
.. _Miller1:
.. _Davidchack4:
**(Davidchack)** R.L Davidchack, T.E. Ouldridge, and M.V. Tretyakov. J. Chem. Phys. 142, 144114 (2015).
.. _Henrich3:
.. _Miller4:
**(Miller)** T. F. Miller III, M. Eleftheriou, P. Pattnaik, A. Ndirango, G. J. Martyna, J. Chem. Phys., 116, 8649-8659 (2002).
.. _Henrich4:
**(Henrich)** O. Henrich, Y. A. Gutierrez-Fosado, T. Curk, T. E. Ouldridge, Eur. Phys. J. E 41, 57 (2018).

21
doc/src/fix_nve_dotc_langevin.rst
View File

@ -38,14 +38,14 @@ Description
"""""""""""
Apply a rigid-body Langevin-type integrator of the kind "Langevin C"
as described in :ref:`(Davidchack) <Davidchack2>`
as described in :ref:`(Davidchack) <Davidchack5>`
to a group of atoms, which models an interaction with an implicit background
solvent. This command performs Brownian dynamics (BD)
via a technique that splits the integration into a deterministic Hamiltonian
part and the Ornstein-Uhlenbeck process for noise and damping.
The quaternion degrees of freedom are updated though an evolution
operator which performs a rotation in quaternion space, preserves
the quaternion norm and is akin to :ref:`(Miller) <Miller2>`.
the quaternion norm and is akin to :ref:`(Miller) <Miller5>`.
In terms of syntax this command has been closely modelled on the
:doc:`fix langevin <fix_langevin>` and its *angmom* option. But it combines
@ -86,7 +86,7 @@ dt damp), where Kb is the Boltzmann constant, T is the desired
temperature, m is the mass of the particle, dt is the timestep size,
and damp is the damping factor. Random numbers are used to randomize
the direction and magnitude of this force as described in
:ref:`(Dunweg) <Dunweg3>`, where a uniform random number is used (instead of
:ref:`(Dunweg) <Dunweg5>`, where a uniform random number is used (instead of
a Gaussian random number) for speed.
@ -128,7 +128,7 @@ The scale factor after the *angmom* keyword gives the ratio of the rotational to
the translational friction coefficient.
An example input file can be found in /examples/USER/cgdna/examples/duplex2/.
Further details of the implementation and stability of the integrators are contained in :ref:`(Henrich) <Henrich4>`.
Further details of the implementation and stability of the integrators are contained in :ref:`(Henrich) <Henrich5>`.
The preprint version of the article can be found `here <PDF/USER-CGDNA.pdf>`_.
@ -154,24 +154,19 @@ Related commands
----------
.. _Davidchack2:
.. _Miller2:
.. _Davidchack5:
**(Davidchack)** R.L Davidchack, T.E. Ouldridge, M.V. Tretyakov. J. Chem. Phys. 142, 144114 (2015).
.. _Dunweg3:
.. _Miller5:
**(Miller)** T. F. Miller III, M. Eleftheriou, P. Pattnaik, A. Ndirango, G. J. Martyna, J. Chem. Phys., 116, 8649-8659 (2002).
.. _Henrich4:
.. _Dunweg5:
**(Dunweg)** B. Dunweg, W. Paul, Int. J. Mod. Phys. C, 2, 817-27 (1991).
.. _Henrich5:
**(Henrich)** O. Henrich, Y. A. Gutierrez-Fosado, T. Curk, T. E. Ouldridge, Eur. Phys. J. E 41, 57 (2018).

45
doc/src/pair_oxdna.rst
View File

@ -36,8 +36,8 @@ Syntax
*oxdna/stk* args = seq T xi kappa 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
seq = seqav (for average sequence stacking strength) or seqdep (for sequence-dependent stacking strength)
T = temperature (oxDNA units, 0.1 = 300 K)
xi = temperature-independent coefficient in stacking strength
kappa = coefficient of linear temperature dependence in stacking strength
xi = 1.3448 (temperature-independent coefficient in stacking strength)
kappa = 2.6568 (coefficient of linear temperature dependence in stacking strength)
*oxdna/hbond* args = seq eps 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
seq = seqav (for average sequence base-pairing strength) or seqdep (for sequence-dependent base-pairing strength)
eps = 1.077 (between base pairs A-T and C-G) or 0 (all other pairs)
@ -94,11 +94,15 @@ Example input and data files for DNA duplexes can be found in examples/USER/cgdn
A simple python setup tool which creates single straight or helical DNA strands,
DNA duplexes or arrays of DNA duplexes can be found in examples/USER/cgdna/util/.
Please cite :ref:`(Henrich) <Henrich1>` and the relevant oxDNA articles in any publication that uses this implementation.
The article contains more information on the model, the structure of the input file, the setup tool
and the performance of the LAMMPS-implementation of oxDNA.
The preprint version of the article can be found `here <PDF/USER-CGDNA.pdf>`_.
Please cite :ref:`(Henrich) <Henrich1>` in any publication that uses
this implementation. The article contains general information
on the model, its implementation and performance as well as the structure of
the data and input file. The preprint version of the article can be found
`here <PDF/USER-CGDNA.pdf>`_.
Please cite also the relevant oxDNA publications
:ref:`(Ouldridge) <Ouldridge1>`,
:ref:`(Ouldridge-DPhil) <Ouldridge-DPhil1>`
and :ref:`(Sulc) <Sulc1>`.
----------
@ -114,39 +118,32 @@ USER-CGDNA package and the MOLECULE and ASPHERE package. See the
Related commands
""""""""""""""""
:doc:`bond\_style oxdna/fene <bond_oxdna>`, :doc:`fix nve/dotc/langevin <fix_nve_dotc_langevin>`, :doc:`pair\_coeff <pair_coeff>`,
:doc:`bond\_style oxdna2/fene <bond_oxdna>`, :doc:`pair\_style oxdna2/excv <pair_oxdna2>`
:doc:`bond\_style oxdna/fene <bond_oxdna>`, :doc:`pair\_coeff <pair_coeff>`,
:doc:`bond\_style oxdna2/fene <bond_oxdna>`, :doc:`pair\_style oxdna2/excv <pair_oxdna2>`,
:doc:`bond\_style oxrna2/fene <bond_oxdna>`, :doc:`pair\_style oxrna2/excv <pair_oxrna2>`,
:doc:`fix nve/dotc/langevin <fix_nve_dotc_langevin>`
**Default:** none
----------
.. _Henrich1:
**(Henrich)** O. Henrich, Y. A. Gutierrez-Fosado, T. Curk, T. E. Ouldridge, Eur. Phys. J. E 41, 57 (2018).
.. _Sulc1:
**(Sulc)** P. Sulc, F. Romano, T.E. Ouldridge, L. Rovigatti, J.P.K. Doye, A.A. Louis, J. Chem. Phys. 137, 135101 (2012).
.. _Ouldridge-DPhil1:
**(Ouldrigde-DPhil)** T.E. Ouldridge, Coarse-grained modelling of DNA and DNA self-assembly, DPhil. University of Oxford (2011).
**(Ouldridge-DPhil)** T.E. Ouldridge, Coarse-grained modelling of DNA and DNA self-assembly, DPhil. University of Oxford (2011).
.. _Ouldridge1:
**(Ouldridge)** T.E. Ouldridge, A.A. Louis, J.P.K. Doye, J. Chem. Phys. 134, 085101 (2011).
.. _Sulc1:
**(Sulc)** P. Sulc, F. Romano, T.E. Ouldridge, L. Rovigatti, J.P.K. Doye, A.A. Louis, J. Chem. Phys. 137, 135101 (2012).
.. _lws: http://lammps.sandia.gov
.. _ld: Manual.html

55
doc/src/pair_oxdna2.rst
View File

@ -39,15 +39,15 @@ Syntax
*oxdna2/stk* args = seq T xi kappa 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
seq = seqav (for average sequence stacking strength) or seqdep (for sequence-dependent stacking strength)
T = temperature (oxDNA units, 0.1 = 300 K)
xi = temperature-independent coefficient in stacking strength
kappa = coefficient of linear temperature dependence in stacking strength
xi = 1.3523 (temperature-independent coefficient in stacking strength)
kappa = 2.6717 (coefficient of linear temperature dependence in stacking strength)
*oxdna2/hbond* args = seq eps 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
seq = seqav (for average sequence base-pairing strength) or seqdep (for sequence-dependent base-pairing strength)
eps = 1.0678 (between base pairs A-T and C-G) or 0 (all other pairs)
*oxdna2/dh* args = T rhos qeff
T = temperature (oxDNA units, 0.1 = 300 K)
rhos = salt concentration (mole per litre)
qeff = effective charge (elementary charges)
qeff = 0.815 (effective charge in elementary charges)
Examples
""""""""
@ -63,7 +63,7 @@ Examples
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
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 0.1 1.0 0.815
pair_coeff \* \* oxdna2/dh 0.1 0.5 0.815
Description
"""""""""""
@ -83,7 +83,7 @@ The exact functional form of the pair styles is rather complex.
The individual potentials consist of products of modulation factors,
which themselves are constructed from a number of more basic potentials
(Morse, Lennard-Jones, harmonic angle and distance) as well as quadratic smoothing and modulation terms.
We refer to :ref:`(Snodin) <Snodin>` and the original oxDNA publications :ref:`(Ouldridge-DPhil) <Ouldridge-DPhil2>`
We refer to :ref:`(Snodin) <Snodin2>` and the original oxDNA publications :ref:`(Ouldridge-DPhil) <Ouldridge-DPhil2>`
and :ref:`(Ouldridge) <Ouldridge2>` for a detailed description of the oxDNA2 force field.
.. note::
@ -94,7 +94,7 @@ and :ref:`(Ouldridge) <Ouldridge2>` for a detailed description of the oxDNA2 fo
in the above example have to be kept fixed and cannot be changed without reparameterizing the entire model.
Exceptions are the first four coefficients after *oxdna2/stk* (seq=seqdep, T=0.1, xi=1.3523 and kappa=2.6717 in the above example),
the first coefficient after *oxdna2/hbond* (seq=seqdep in the above example) and the three coefficients
after *oxdna2/dh* (T=0.1, rhos=1.0, qeff=0.815 in the above example). When using a Langevin thermostat
after *oxdna2/dh* (T=0.1, rhos=0.5, qeff=0.815 in the above example). When using a Langevin thermostat
e.g. through :doc:`fix langevin <fix_langevin>` or :doc:`fix nve/dotc/langevin <fix_nve_dotc_langevin>`
the temperature coefficients have to be matched to the one used in the fix.
@ -102,11 +102,13 @@ Example input and data files for DNA duplexes can be found in examples/USER/cgdn
A simple python setup tool which creates single straight or helical DNA strands,
DNA duplexes or arrays of DNA duplexes can be found in examples/USER/cgdna/util/.
Please cite :ref:`(Henrich) <Henrich>` and the relevant oxDNA articles in any publication that uses this implementation.
The article contains more information on the model, the structure of the input file, the setup tool
and the performance of the LAMMPS-implementation of oxDNA.
The preprint version of the article can be found `here <PDF/USER-CGDNA.pdf>`_.
Please cite :ref:`(Henrich) <Henrich2>` in any publication that uses
this implementation. The article contains general information
on the model, its implementation and performance as well as the structure of
the data and input file. The preprint version of the article can be found
`here <PDF/USER-CGDNA.pdf>`_.
Please cite also the relevant oxDNA2 publications
:ref:`(Snodin) <Snodin2>` and :ref:`(Sulc) <Sulc2>`.
----------
@ -122,43 +124,34 @@ USER-CGDNA package and the MOLECULE and ASPHERE package. See the
Related commands
""""""""""""""""
:doc:`bond\_style oxdna2/fene <bond_oxdna>`, :doc:`fix nve/dotc/langevin <fix_nve_dotc_langevin>`, :doc:`pair\_coeff <pair_coeff>`,
:doc:`bond\_style oxdna/fene <bond_oxdna>`, :doc:`pair\_style oxdna/excv <pair_oxdna>`
:doc:`bond\_style oxdna2/fene <bond_oxdna>`, :doc:`pair\_coeff <pair_coeff>`,
:doc:`bond\_style oxdna/fene <bond_oxdna>`, :doc:`pair\_style oxdna/excv <pair_oxdna>`,
:doc:`bond\_style oxrna2/fene <bond_oxdna>`, :doc:`pair\_style oxrna2/excv <pair_oxrna2>`,
:doc:`fix nve/dotc/langevin <fix_nve_dotc_langevin>`
**Default:** none
----------
.. _Henrich:
.. _Henrich2:
**(Henrich)** O. Henrich, Y. A. Gutierrez-Fosado, T. Curk, T. E. Ouldridge, Eur. Phys. J. E 41, 57 (2018).
.. _Sulc2:
**(Sulc)** P. Sulc, F. Romano, T.E. Ouldridge, L. Rovigatti, J.P.K. Doye, A.A. Louis, J. Chem. Phys. 137, 135101 (2012).
.. _Snodin:
.. _Snodin2:
**(Snodin)** B.E. Snodin, F. Randisi, M. Mosayebi, et al., J. Chem. Phys. 142, 234901 (2015).
.. _Sulc2:
**(Sulc)** P. Sulc, F. Romano, T.E. Ouldridge, L. Rovigatti, J.P.K. Doye, A.A. Louis, J. Chem. Phys. 137, 135101 (2012).
.. _Ouldridge-DPhil2:
**(Ouldrigde-DPhil)** T.E. Ouldridge, Coarse-grained modelling of DNA and DNA self-assembly, DPhil. University of Oxford (2011).
**(Ouldridge-DPhil)** T.E. Ouldridge, Coarse-grained modelling of DNA and DNA self-assembly, DPhil. University of Oxford (2011).
.. _Ouldridge2:
**(Ouldridge)** T.E. Ouldridge, A.A. Louis, J.P.K. Doye, J. Chem. Phys. 134, 085101 (2011).

158
doc/src/pair_oxrna2.rst Normal file
View File

@ -0,0 +1,158 @@
.. index:: pair_style oxrna2/excv
pair_style oxrna2/excv command
==============================
pair_style oxrna2/stk command
=============================
pair_style oxrna2/hbond command
===============================
pair_style oxrna2/xstk command
==============================
pair_style oxrna2/coaxstk command
=================================
pair_style oxrna2/dh command
============================
Syntax
""""""
.. code-block:: LAMMPS
pair_style style1
pair_coeff * * style2 args
* style1 = *hybrid/overlay oxrna2/excv oxrna2/stk oxrna2/hbond oxrna2/xstk oxrna2/coaxstk oxrna2/dh*
* style2 = *oxrna2/excv* or *oxrna2/stk* or *oxrna2/hbond* or *oxrna2/xstk* or *oxrna2/coaxstk* or *oxrna2/dh*
* args = list of arguments for these particular styles
.. parsed-literal::
*oxrna2/stk* args = seq T xi kappa 6.0 0.43 0.93 0.35 0.78 0.9 0 0.95 0.9 0 0.95 1.3 0 0.8 1.3 0 0.8 2.0 0.65 2.0 0.65
seq = seqav (for average sequence stacking strength) or seqdep (for sequence-dependent stacking strength)
T = temperature (oxDNA units, 0.1 = 300 K)
xi = 1.40206 (temperature-independent coefficient in stacking strength)
kappa = 2.77 (coefficient of linear temperature dependence in stacking strength)
*oxrna2/hbond* args = seq eps 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
seq = seqav (for average sequence base-pairing strength) or seqdep (for sequence-dependent base-pairing strength)
eps = 0.870439 (between base pairs A-T, C-G and G-T) or 0 (all other pairs)
*oxrna2/dh* args = T rhos qeff
T = temperature (oxDNA units, 0.1 = 300 K)
rhos = salt concentration (mole per litre)
qeff = 1.02455 (effective charge in elementary charges)
Examples
""""""""
.. code-block:: LAMMPS
pair_style hybrid/overlay oxrna2/excv oxrna2/stk oxrna2/hbond oxrna2/xstk oxrna2/coaxstk oxrna2/dh
pair_coeff * * oxrna2/excv 2.0 0.7 0.675 2.0 0.515 0.5 2.0 0.33 0.32
pair_coeff * * oxrna2/stk seqdep 0.1 1.40206 2.77 6.0 0.43 0.93 0.35 0.78 0.9 0 0.95 0.9 0 0.95 1.3 0 0.8 1.3 0 0.8 2.0 0.65 2.0 0.65
pair_coeff * * oxrna2/hbond seqdep 0.0 8.0 0.4 0.75 0.34 0.7 1.5 0 0.7 1.5 0 0.7 1.5 0 0.7 0.46 3.141592653589793 0.7 4.0 1.5707963267948966 0.45 4.0 1.5707963267948966 0.45
pair_coeff 1 4 oxrna2/hbond seqdep 0.870439 8.0 0.4 0.75 0.34 0.7 1.5 0 0.7 1.5 0 0.7 1.5 0 0.7 0.46 3.141592653589793 0.7 4.0 1.5707963267948966 0.45 4.0 1.5707963267948966 0.45
pair_coeff 2 3 oxrna2/hbond seqdep 0.870439 8.0 0.4 0.75 0.34 0.7 1.5 0 0.7 1.5 0 0.7 1.5 0 0.7 0.46 3.141592653589793 0.7 4.0 1.5707963267948966 0.45 4.0 1.5707963267948966 0.45
pair_coeff 3 4 oxrna2/hbond seqdep 0.870439 8.0 0.4 0.75 0.34 0.7 1.5 0 0.7 1.5 0 0.7 1.5 0 0.7 0.46 3.141592653589793 0.7 4.0 1.5707963267948966 0.45 4.0 1.5707963267948966 0.45
pair_coeff * * oxrna2/xstk 59.9626 0.5 0.6 0.42 0.58 2.25 0.505 0.58 1.7 1.266 0.68 1.7 1.266 0.68 1.7 0.309 0.68 1.7 0.309 0.68
pair_coeff * * oxrna2/coaxstk 80 0.5 0.6 0.42 0.58 2.0 2.592 0.65 1.3 0.151 0.8 0.9 0.685 0.95 0.9 0.685 0.95 2.0 -0.65 2.0 -0.65
pair_coeff * * oxrna2/dh 0.1 0.5 1.02455
Description
"""""""""""
The *oxrna2* pair styles compute the pairwise-additive parts of the oxDNA force field
for coarse-grained modelling of DNA. The effective interaction between the nucleotides consists of potentials for the
excluded volume interaction *oxrna2/excv*\ , the stacking *oxrna2/stk*\ , cross-stacking *oxrna2/xstk*
and coaxial stacking interaction *oxrna2/coaxstk*\ , electrostatic Debye-Hueckel interaction *oxrna2/dh*
as well as the hydrogen-bonding interaction *oxrna2/hbond* between complementary pairs of nucleotides on
opposite strands. Average sequence or sequence-dependent stacking and base-pairing strengths
are supported :ref:`(Sulc2) <Sulc32>`. Quasi-unique base-pairing between nucleotides can be achieved by using
more complementary pairs of atom types like 5-8 and 6-7, 9-12 and 10-11, 13-16 and 14-15, etc.
This prevents the hybridization of in principle complementary bases within Ntypes/4 bases
up and down along the backbone.
The exact functional form of the pair styles is rather complex.
The individual potentials consist of products of modulation factors,
which themselves are constructed from a number of more basic potentials
(Morse, Lennard-Jones, harmonic angle and distance) as well as quadratic smoothing and modulation terms.
We refer to :ref:`(Sulc1) <Sulc31>` and the original oxDNA publications :ref:`(Ouldridge-DPhil) <Ouldridge-DPhil3>`
and :ref:`(Ouldridge) <Ouldridge3>` for a detailed description of the oxRNA2 force field.
.. note::
These pair styles have to be used together with the related oxDNA2 bond style
*oxrna2/fene* for the connectivity of the phosphate backbone (see also documentation of
:doc:`bond\_style oxrna2/fene <bond_oxdna>`). Most of the coefficients
in the above example have to be kept fixed and cannot be changed without reparameterizing the entire model.
Exceptions are the first four coefficients after *oxrna2/stk* (seq=seqdep, T=0.1, xi=1.40206 and kappa=2.77 in the above example),
the first coefficient after *oxrna2/hbond* (seq=seqdep in the above example) and the three coefficients
after *oxrna2/dh* (T=0.1, rhos=0.5, qeff=1.02455 in the above example). When using a Langevin thermostat
e.g. through :doc:`fix langevin <fix_langevin>` or :doc:`fix nve/dotc/langevin <fix_nve_dotc_langevin>`
the temperature coefficients have to be matched to the one used in the fix.
Example input and data files for DNA duplexes can be found in examples/USER/cgdna/examples/oxDNA/ and /oxDNA2/.
A simple python setup tool which creates single straight or helical DNA strands,
DNA duplexes or arrays of DNA duplexes can be found in examples/USER/cgdna/util/.
Please cite :ref:`(Henrich) <Henrich3>` in any publication that uses
this implementation. The article contains general information
on the model, its implementation and performance as well as the structure of
the data and input file. The preprint version of the article can be found
`here <PDF/USER-CGDNA.pdf>`_.
Please cite also the relevant oxRNA2 publications
:ref:`(Sulc1) <Sulc31>` and :ref:`(Sulc2) <Sulc32>`.
----------
Restrictions
""""""""""""
These pair styles can only be used if LAMMPS was built with the
USER-CGDNA package and the MOLECULE and ASPHERE package. See the
:doc:`Build package <Build_package>` doc page for more info.
Related commands
""""""""""""""""
:doc:`bond\_style oxrna2/fene <bond_oxdna>`, :doc:`pair\_coeff <pair_coeff>`,
:doc:`bond\_style oxdna/fene <bond_oxdna>`, :doc:`pair\_style oxdna/excv <pair_oxdna>`,
:doc:`bond\_style oxdna2/fene <bond_oxdna>`, :doc:`pair\_style oxdna2/excv <pair_oxdna2>`,
:doc:`fix nve/dotc/langevin <fix_nve_dotc_langevin>`
**Default:**
none
----------
.. _Henrich3:
**(Henrich)** O. Henrich, Y. A. Gutierrez-Fosado, T. Curk, T. E. Ouldridge, Eur. Phys. J. E 41, 57 (2018).
.. _Sulc31:
**(Sulc1)** P. Sulc, F. Romano, T. E. Ouldridge, et al., J. Chem. Phys. 140, 235102 (2014).
.. _Sulc32:
**(Sulc2)** P. Sulc, F. Romano, T.E. Ouldridge, L. Rovigatti, J.P.K. Doye, A.A. Louis, J. Chem. Phys. 137, 135101 (2012).
.. _Ouldridge-DPhil3:
**(Ouldridge-DPhil)** T.E. Ouldridge, Coarse-grained modelling of DNA and DNA self-assembly, DPhil. University of Oxford (2011).
.. _Ouldridge3:
**(Ouldridge)** T.E. Ouldridge, A.A. Louis, J.P.K. Doye, J. Chem. Phys. 134, 085101 (2011).

73
examples/USER/cgdna/README
View File

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

File diff suppressed because it is too large Load Diff

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

File diff suppressed because it is too large Load Diff

63
examples/rigid/in.rigid.gravity Normal file
View File

@ -0,0 +1,63 @@
#Pour composite granular particles on flat wall
newton on
atom_style sphere
atom_modify map array sort 0 0
thermo_modify flush yes
units si
variable minrad equal 0.5
variable maxrad equal 1.4
variable skin equal 0.3*${maxrad}
boundary p p f
region reg block 0 20 0 20 0 200 units box
create_box 1 reg
fix prop all property/atom mol ghost yes
variable dumpfreq equal 1000
variable logfreq equal 1000
pair_style gran/hooke/history 4e5 NULL 1e2 NULL 0.5 0
pair_coeff * *
timestep 0.0001
group particles type 1
atom_modify first particles
neighbor ${skin} bin
group rigid type 1
neigh_modify every 1 delay 0 check yes exclude molecule/intra all
thermo ${logfreq}
thermo_style custom step cpu atoms ke
thermo_modify flush yes lost warn
comm_modify vel yes cutoff 3
molecule mymol molecule.data
region pourreg block 5 15 5 15 80 100 side in units box
#Note: in versions prior to 1/2020, the 'disable' keyword to fix/gravity
# and the 'gravity' keyword to fix rigid/small were not available.
# These settings produce undesirable behavior, where gravity can induce
# torque on rigid bodies.
#fix gravfix all gravity 9.8 vector 0 0 -1 #disable
#fix rigidfix all rigid/small molecule mol mymol #gravity gravfix
#The correct behavior is recovered with the following settings:
fix gravfix all gravity 9.8 vector 0 0 -1 disable
fix rigidfix all rigid/small molecule mol mymol gravity gravfix
fix pourfix all pour 5 0 1234 region pourreg mol mymol rigid rigidfix
fix zwall all wall/gran hooke/history 4000.0 NULL 100.0 NULL 0.5 0 zplane 0.1 NULL
#dump 1 all custom 1000 molecule_pour.dump id type mass radius x y z fx fy fz
run 100000

41
examples/rigid/molecule.data Normal file
View File

@ -0,0 +1,41 @@
LAMMPS data file created for rigid body molecule template
5 atoms
2.3388800000000005 mass
6.002239704473936 4.99 4.989999999999999 com
116.79265620480001 144.26721336320003 144.26721336320006 -70.05220681600004 -70.05220681600002 -58.238345888000005 inertia
Coords
1 5 5 5
2 5.1 5.0 5.0
3 5.2 5.0 5.0
4 6.2 5.0 5.0
5 7.2 5.0 5.0
Types
1 1
2 1
3 1
4 1
5 1
Diameters
1 1.0
2 0.9
3 1.2
4 1.2
5 1.0
Masses
1 0.5235987755982988
2 0.3817035074111599
3 0.9047786842338602
4 0.9047786842338602
5 0.5235987755982988

12
src/USER-CGDNA/Install.sh
View File

@ -30,8 +30,10 @@ action () {
action bond_oxdna_fene.cpp bond_fene.h
action bond_oxdna2_fene.cpp bond_fene.h
action bond_oxrna2_fene.cpp bond_fene.h
action bond_oxdna_fene.h bond_fene.h
action bond_oxdna2_fene.h bond_fene.h
action bond_oxrna2_fene.h bond_fene.h
action fix_nve_dotc_langevin.cpp atom_vec_ellipsoid.h
action fix_nve_dotc_langevin.h atom_vec_ellipsoid.h
action fix_nve_dot.cpp atom_vec_ellipsoid.h
@ -43,13 +45,23 @@ action pair_oxdna_coaxstk.h atom_vec_ellipsoid.h
action pair_oxdna2_coaxstk.h atom_vec_ellipsoid.h
action pair_oxdna_excv.cpp atom_vec_ellipsoid.h
action pair_oxdna2_excv.cpp atom_vec_ellipsoid.h
action pair_oxrna2_excv.cpp atom_vec_ellipsoid.h
action pair_oxdna_excv.h atom_vec_ellipsoid.h
action pair_oxdna2_excv.h atom_vec_ellipsoid.h
action pair_oxrna2_excv.h atom_vec_ellipsoid.h
action pair_oxdna_hbond.cpp atom_vec_ellipsoid.h
action pair_oxdna_hbond.h atom_vec_ellipsoid.h
action pair_oxrna2_hbond.cpp atom_vec_ellipsoid.h
action pair_oxrna2_hbond.h atom_vec_ellipsoid.h
action pair_oxdna_stk.cpp atom_vec_ellipsoid.h
action pair_oxdna_stk.h atom_vec_ellipsoid.h
action pair_oxrna2_stk.cpp atom_vec_ellipsoid.h
action pair_oxrna2_stk.h atom_vec_ellipsoid.h
action pair_oxdna_xstk.cpp atom_vec_ellipsoid.h
action pair_oxdna_xstk.h atom_vec_ellipsoid.h
action pair_oxrna2_xstk.cpp atom_vec_ellipsoid.h
action pair_oxrna2_xstk.h atom_vec_ellipsoid.h
action pair_oxdna2_dh.cpp atom_vec_ellipsoid.h
action pair_oxdna2_dh.h atom_vec_ellipsoid.h
action pair_oxrna2_dh.cpp atom_vec_ellipsoid.h
action pair_oxrna2_dh.h atom_vec_ellipsoid.h

87
src/USER-CGDNA/README
View File

@ -1,87 +0,0 @@
This package contains a LAMMPS implementation of coarse-grained
models of DNA, which can be used to model sequence-specific
DNA strands.
Please cite [1] and the relevant oxDNA articles in any publication
that uses this package.
See the doc pages and [2,3,4] for the individual bond and pair styles.
The packages contains also a new Langevin-type rigid-body integrator,
which has also its own doc page and is explained in [5].
[1] O. Henrich, Y. A. Gutierrez-Fosado, T. Curk, T. E. Ouldridge,
"Coarse-grained simulation of DNA using LAMMPS",
Eur. Phys. J. E 41, 57 (2018).
[2] T. Ouldridge, A. Louis, J. Doye, "Structural, mechanical,
and thermodynamic properties of a coarse-grained DNA model",
J. Chem. Phys. 134, 085101 (2011).
[3] T.E. Ouldridge, Coarse-grained modelling of DNA and DNA
self-assembly, DPhil. University of Oxford (2011).
[4] B.E. Snodin, F. Randisi, M. Mosayebi, et al., Introducing
Improved Structural Properties and Salt Dependence into a Coarse-Grained
Model of DNA, J. Chem. Phys. 142, 234901 (2015).
[5] R. Davidchack, T. Ouldridge, M. Tretyakov, "New Langevin and
gradient thermostats for rigid body dynamics", J. Chem. Phys. 142,
144114 (2015).
Example input and data files can be found in
/examples/USER/cgdna/examples/oxDNA/ and /oxDNA2/. Python setup
tools which create single straight or helical DNA strands as
well as DNA duplexes or arrays of duplexes can be found in
/examples/USER/cgdna/util/. A technical report with more information
on the models, the structure of the input and data file, the setup tool
and the performance of the LAMMPS-implementation of oxDNA can be found
in /doc/src/PDF/USER-CGDNA.pdf.
IMPORTANT NOTE: This package can only be used if LAMMPS is compiled
with the MOLECULE and ASPHERE packages. These should be included in
the LAMMPS build by typing "make yes-asphere yes-molecule" prior to
the usual compilation (see the "Including/excluding packages" section
of the LAMMPS manual).
The creator of this package is:
Dr Oliver Henrich
University of Strathclyde, Glasgow, UK
oliver d o t henrich a t strath d o t ac d o t uk
--------------------------------------------------------------------------
** Bond styles provided by this package:
bond_oxdna_fene.cpp: backbone connectivity, a modified FENE potential
bond_oxdna2_fene.cpp: corresponding bond style in oxDNA2 (see [3])
** Pair styles provided by this package:
pair_oxdna_excv.cpp: excluded volume interaction between the nucleotides
pair_oxdna_stk.cpp: stacking interaction between consecutive nucleotides
on the same strand
pair_oxdna_hbond.cpp: hydrogen-bonding interaction between complementary
nucleotides on different strands, e.g. A-T and C-G
pair_oxdna_xstk.cpp: cross-stacking interaction between nucleotides
pair_oxdna_coaxstk.cpp: coaxial stacking interaction between nucleotides
pair_oxdna2_excv.cpp, pair_oxdna2_coaxstk.cpp:
corresponding pair styles in oxDNA2 (see [3])
pair_oxdna2_dh.cpp: Debye-Hueckel electrostatic interaction between backbone
sites
** Fixes provided by this package:
fix_nve_dotc_langevin.cpp: fix for Langevin-type rigid body integrator "C"
in above Ref. [3]
fix_nve_dot.cpp: NVE-type rigid body integrator without noise

105
src/USER-CGDNA/README.md Normal file
View File

@ -0,0 +1,105 @@
This package contains a LAMMPS implementation of coarse-grained
models of DNA, which can be used to model sequence-specific
DNA strands.
Please cite
[![DOI](https://zenodo.org/badge/132764768.svg)](https://zenodo.org/badge/latestdoi/132764768)
as well as [1] and the relevant oxDNA, oxDNA2 and oxRNA2 articles
in any publication that uses this package.
See the doc pages and [2,3,4,5,6] for the individual bond and pair styles.
The packages contains also a new Langevin-type rigid-body integrator,
which has also its own doc page and is explained in [7].
[1] O. Henrich, Y. A. Gutierrez-Fosado, T. Curk, T. E. Ouldridge,
"Coarse-grained simulation of DNA using LAMMPS",
Eur. Phys. J. E 41, 57 (2018).
[2] T. Ouldridge, A. Louis, J. Doye, "Structural, mechanical,
and thermodynamic properties of a coarse-grained DNA model",
J. Chem. Phys. 134, 085101 (2011).
[3] T.E. Ouldridge, "Coarse-grained modelling of DNA and DNA
self-assembly", DPhil. University of Oxford (2011).
[4] B.E. Snodin, F. Randisi, M. Mosayebi, et al., "Introducing
Improved structural properties and salt dependence into a coarse-grained
model of DNA", J. Chem. Phys. 142, 234901 (2015).
[5] P. Sulc, F. Romano, T.E. Ouldridge, et al., "A nucleotide-level
coarse-grained model of RNA", J. Chem. Phys. 140, 235102 (2014).
[6] P. Sulc, F. Romano, T.E. Ouldridge, et al., "Sequence-dependent
thermodynamics of a coarse-grained DNA model",
J. Chem. Phys. 137, 135101 (2012).
[7] R. Davidchack, T. Ouldridge, M. Tretyakov, "New Langevin and
gradient thermostats for rigid body dynamics", J. Chem. Phys. 142,
144114 (2015).
Example input and data files can be found in
/examples/USER/cgdna/examples/oxDNA/, /oxDNA2/ and /oxRNA2/.
Python setup tools which create single straight or helical DNA or RNA
strands as well as DNA or RNA duplexes or arrays of duplexes can be
found in /examples/USER/cgdna/util/. A technical report with more
general information on the model, its implementation and performance
as well as the structure of the data and input file can be found
in /doc/src/PDF/USER-CGDNA.pdf.
IMPORTANT NOTE: This package can only be used if LAMMPS is compiled
with the MOLECULE and ASPHERE packages. These should be included in
the LAMMPS build by typing "make yes-asphere yes-molecule" prior to
the usual compilation (see the "Including/excluding packages" section
of the LAMMPS manual).
The creator of this package is:
Dr Oliver Henrich
University of Strathclyde, Glasgow, UK
oliver d o t henrich a t strath d o t ac d o t uk
--------------------------------------------------------------------------
** Bond styles provided by this package:
bond_oxdna_fene.cpp: backbone connectivity,
a modified FENE potential (see [2,3])
bond_oxdna2_fene.cpp: corresponding bond style in oxDNA2 (see [4])
bond_oxrna2_fene.cpp: corresponding bond style in oxRNA2 (see [5])
** Pair styles provided by this package:
pair_oxdna_excv.cpp: excluded volume interaction between the nucleotides
pair_oxdna_stk.cpp: stacking interaction between consecutive nucleotides
on the same strand
pair_oxdna_hbond.cpp: hydrogen-bonding interaction between complementary
nucleotides on different strands, e.g. A-T and C-G
pair_oxdna_xstk.cpp: cross-stacking interaction between nucleotides
pair_oxdna_coaxstk.cpp: coaxial stacking interaction between nucleotides
pair_oxdna2_excv.cpp, pair_oxdna2_coaxstk.cpp:
corresponding pair styles in oxDNA2 (see [4])
pair_oxrna2_excv.cpp, pair_oxrna2_stk.cpp, pair_oxrna2_hbond.cpp,
pair_oxrna2_xstk.cpp:
corresponding pair styles in oxDNA2 (see [5])
pair_oxdna2_dh.cpp, pair_oxrna2_dh.cpp:
Debye-Hueckel electrostatic interaction between backbone sites
** Fixes provided by this package:
fix_nve_dotc_langevin.cpp: fix for Langevin-type rigid body integrator "C"
in above Ref. [7]
fix_nve_dot.cpp: NVE-type rigid body integrator without noise

49
src/USER-CGDNA/bond_oxrna2_fene.cpp Normal file
View File

@ -0,0 +1,49 @@
/* ----------------------------------------------------------------------
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)
------------------------------------------------------------------------- */
#include <cmath>
#include <cstdlib>
#include "bond_oxrna2_fene.h"
using namespace LAMMPS_NS;
/* ---------------------------------------------------------------------- */
BondOxrna2Fene::BondOxrna2Fene(LAMMPS *lmp) : BondOxdnaFene(lmp)
{
}
/* ---------------------------------------------------------------------- */
BondOxrna2Fene::~BondOxrna2Fene()
{
}
/* ----------------------------------------------------------------------
compute vector COM-sugar-phosphate backbone interaction site in oxRNA2
------------------------------------------------------------------------- */
void BondOxrna2Fene::compute_interaction_sites(double e1[3], double /*e2*/[3],
double e3[3], double r[3])
{
double d_cs_x=-0.4, d_cs_z=+0.2;
r[0] = d_cs_x*e1[0] + d_cs_z*e3[0];
r[1] = d_cs_x*e1[1] + d_cs_z*e3[1];
r[2] = d_cs_x*e1[2] + d_cs_z*e3[2];
}

65
src/USER-CGDNA/bond_oxrna2_fene.h Normal file
View File

@ -0,0 +1,65 @@
/* -*- c++ -*- ----------------------------------------------------------
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.
------------------------------------------------------------------------- */
#ifdef BOND_CLASS
BondStyle(oxrna2/fene,BondOxrna2Fene)
#else
#ifndef LMP_BOND_OXRNA2_FENE_H
#define LMP_BOND_OXRNA2_FENE_H
#include "bond_oxdna_fene.h"
namespace LAMMPS_NS {
class BondOxrna2Fene : public BondOxdnaFene {
public:
BondOxrna2Fene(class LAMMPS *);
virtual ~BondOxrna2Fene();
virtual void compute_interaction_sites(double *, double *, double *,
double *);
};
}
#endif
#endif
/* ERROR/WARNING messages:
W: FENE bond too long: %ld %d %d %g
A FENE bond has stretched dangerously far. It's interaction strength
will be truncated to attempt to prevent the bond from blowing up.
E: Bad FENE bond
Two atoms in a FENE bond have become so far apart that the bond cannot
be computed.
E: Incorrect args for bond coefficients
Self-explanatory. Check the input script or data file.
W: Use special bonds = 0,1,1 with bond style oxrna
Most FENE models need this setting for the special_bonds command.
W: FENE bond too long: %ld %g
A FENE bond has stretched dangerously far. It's interaction strength
will be truncated to attempt to prevent the bond from blowing up.
*/

6
src/USER-CGDNA/pair_oxdna2_dh.cpp
View File

@ -64,7 +64,7 @@ PairOxdna2Dh::~PairOxdna2Dh()
compute vector COM-sugar-phosphate backbone interaction site in oxDNA2
------------------------------------------------------------------------- */
void PairOxdna2Dh::compute_interaction_sites(double e1[3],
double e2[3], double r[3])
double e2[3], double /*e3*/[3], double r[3])
{
double d_cs_x=-0.34, d_cs_y=+0.3408;
@ -125,7 +125,7 @@ void PairOxdna2Dh::compute(int eflag, int vflag)
MathExtra::q_to_exyz(qa,ax,ay,az);
// vector COM-backbone site a
compute_interaction_sites(ax,ay,ra_cs);
compute_interaction_sites(ax,ay,az,ra_cs);
rtmp_s[0] = x[a][0] + ra_cs[0];
rtmp_s[1] = x[a][1] + ra_cs[1];
@ -145,7 +145,7 @@ void PairOxdna2Dh::compute(int eflag, int vflag)
MathExtra::q_to_exyz(qb,bx,by,bz);
// vector COM-backbone site b
compute_interaction_sites(bx,by,rb_cs);
compute_interaction_sites(bx,by,bz,rb_cs);
// vector backbone site b to a
delr[0] = rtmp_s[0] - x[b][0] - rb_cs[0];

3
src/USER-CGDNA/pair_oxdna2_dh.h
View File

@ -28,7 +28,8 @@ class PairOxdna2Dh : public Pair {
public:
PairOxdna2Dh(class LAMMPS *);
virtual ~PairOxdna2Dh();
virtual void compute_interaction_sites(double *, double *, double *);
virtual void compute_interaction_sites(double *, double *, double *,
double *);
virtual void compute(int, int);
void settings(int, char **);
void coeff(int, char **);

1
src/USER-CGDNA/pair_oxdna_coaxstk.h
View File

@ -14,6 +14,7 @@
#ifdef PAIR_CLASS
PairStyle(oxdna/coaxstk,PairOxdnaCoaxstk)
PairStyle(oxrna2/coaxstk,PairOxdnaCoaxstk)
#else

47
src/USER-CGDNA/pair_oxrna2_dh.cpp Normal file
View File

@ -0,0 +1,47 @@
/* ----------------------------------------------------------------------
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)
------------------------------------------------------------------------- */
#include "pair_oxrna2_dh.h"
using namespace LAMMPS_NS;
/* ---------------------------------------------------------------------- */
PairOxrna2Dh::PairOxrna2Dh(LAMMPS *lmp) : PairOxdna2Dh(lmp)
{
}
/* ---------------------------------------------------------------------- */
PairOxrna2Dh::~PairOxrna2Dh()
{
}
/* ----------------------------------------------------------------------
compute vector COM-sugar-phosphate backbone interaction site in oxRNA2
------------------------------------------------------------------------- */
void PairOxrna2Dh::compute_interaction_sites(double e1[3], double /*e2*/[3],
double e3[3], double r[3])
{
double d_cs_x=-0.4, d_cs_z=+0.2;
r[0] = d_cs_x*e1[0] + d_cs_z*e3[0];
r[1] = d_cs_x*e1[1] + d_cs_z*e3[1];
r[2] = d_cs_x*e1[2] + d_cs_z*e3[2];
}

53
src/USER-CGDNA/pair_oxrna2_dh.h Normal file
View File

@ -0,0 +1,53 @@
/* -*- c++ -*- ----------------------------------------------------------
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.
------------------------------------------------------------------------- */
#ifdef PAIR_CLASS
PairStyle(oxrna2/dh,PairOxrna2Dh)
#else
#ifndef LMP_PAIR_OXRNA2_DH_H
#define LMP_PAIR_OXRNA2_DH_H
#include "pair_oxdna2_dh.h"
namespace LAMMPS_NS {
class PairOxrna2Dh : public PairOxdna2Dh {
public:
PairOxrna2Dh(class LAMMPS *);
virtual ~PairOxrna2Dh();
virtual void compute_interaction_sites(double *, double *, double *,
double *);
};
}
#endif
#endif
/* ERROR/WARNING messages:
E: Illegal ... command
Self-explanatory. Check the input script syntax and compare to the
documentation for the command. You can use -echo screen as a
command-line option when running LAMMPS to see the offending line.
E: Incorrect args for pair coefficients
Self-explanatory. Check the input script or data file.
*/

55
src/USER-CGDNA/pair_oxrna2_excv.cpp Normal file
View File

@ -0,0 +1,55 @@
/* ----------------------------------------------------------------------
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)
------------------------------------------------------------------------- */
#include <cmath>
#include <cstdio>
#include <cstdlib>
#include <cstring>
#include "pair_oxrna2_excv.h"
using namespace LAMMPS_NS;
/* ---------------------------------------------------------------------- */
PairOxrna2Excv::PairOxrna2Excv(LAMMPS *lmp) : PairOxdnaExcv(lmp)
{
}
/* ---------------------------------------------------------------------- */
PairOxrna2Excv::~PairOxrna2Excv()
{
}
/* ----------------------------------------------------------------------
compute vector COM-excluded volume interaction sites in oxRNA2
------------------------------------------------------------------------- */
void PairOxrna2Excv::compute_interaction_sites(double e1[3], double /*e2*/[3],
double e3[3], double rs[3], double rb[3])
{
double d_cs_x=-0.4, d_cs_z=+0.2, d_cb=+0.4;
rs[0] = d_cs_x*e1[0] + d_cs_z*e3[0];
rs[1] = d_cs_x*e1[1] + d_cs_z*e3[1];
rs[2] = d_cs_x*e1[2] + d_cs_z*e3[2];
rb[0] = d_cb*e1[0];
rb[1] = d_cb*e1[1];
rb[2] = d_cb*e1[2];
}

52
src/USER-CGDNA/pair_oxrna2_excv.h Normal file
View File

@ -0,0 +1,52 @@
/* -*- c++ -*- ----------------------------------------------------------
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.
------------------------------------------------------------------------- */
#ifdef PAIR_CLASS
PairStyle(oxrna2/excv,PairOxrna2Excv)
#else
#ifndef LMP_PAIR_OXRNA2_EXCV_H
#define LMP_PAIR_OXRNA2_EXCV_H
#include "pair_oxdna_excv.h"
namespace LAMMPS_NS {
class PairOxrna2Excv : public PairOxdnaExcv {
public:
PairOxrna2Excv(class LAMMPS *);
virtual ~PairOxrna2Excv();
virtual void compute_interaction_sites(double *, double *,
double *, double *, double *);
};
}
#endif
#endif
/* ERROR/WARNING messages:
E: Illegal ... command
Self-explanatory. Check the input script syntax and compare to the
documentation for the command. You can use -echo screen as a
command-line option when running LAMMPS to see the offending line.
E: Incorrect args for pair coefficients
Self-explanatory. Check the input script or data file.
*/

58
src/USER-CGDNA/pair_oxrna2_hbond.cpp Normal file
View File

@ -0,0 +1,58 @@
/* ----------------------------------------------------------------------
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)
------------------------------------------------------------------------- */
#include "pair_oxrna2_hbond.h"
using namespace LAMMPS_NS;
/* ---------------------------------------------------------------------- */
PairOxrna2Hbond::PairOxrna2Hbond(LAMMPS *lmp) : PairOxdnaHbond(lmp)
{
single_enable = 0;
writedata = 1;
// sequence-specific base-pairing strength
// A:0 C:1 G:2 U:3, 5'- [i][j] -3'
alpha_hb[0][0] = 1.00000;
alpha_hb[0][1] = 1.00000;
alpha_hb[0][2] = 1.00000;
alpha_hb[0][3] = 0.94253;
alpha_hb[1][0] = 1.00000;
alpha_hb[1][1] = 1.00000;
alpha_hb[1][2] = 1.22288;
alpha_hb[1][3] = 1.00000;
alpha_hb[2][0] = 1.00000;
alpha_hb[2][1] = 1.22288;
alpha_hb[2][2] = 1.00000;
alpha_hb[2][3] = 0.58655;
alpha_hb[3][0] = 0.94253;
alpha_hb[3][1] = 1.00000;
alpha_hb[3][2] = 0.58655;
alpha_hb[3][3] = 1.00000;
}
/* ---------------------------------------------------------------------- */
PairOxrna2Hbond::~PairOxrna2Hbond()
{
}

50
src/USER-CGDNA/pair_oxrna2_hbond.h Normal file
View File

@ -0,0 +1,50 @@
/* -*- c++ -*- ----------------------------------------------------------
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.
------------------------------------------------------------------------- */
#ifdef PAIR_CLASS
PairStyle(oxrna2/hbond,PairOxrna2Hbond)
#else
#ifndef LMP_PAIR_OXRNA2_HBOND_H
#define LMP_PAIR_OXRNA2_HBOND_H
#include "pair_oxdna_hbond.h"
namespace LAMMPS_NS {
class PairOxrna2Hbond : public PairOxdnaHbond {
public:
PairOxrna2Hbond(class LAMMPS *);
virtual ~PairOxrna2Hbond();
};
}
#endif
#endif
/* ERROR/WARNING messages:
E: Illegal ... command
Self-explanatory. Check the input script syntax and compare to the
documentation for the command. You can use -echo screen as a
command-line option when running LAMMPS to see the offending line.
E: Incorrect args for pair coefficients
Self-explanatory. Check the input script or data file.
*/

1435
src/USER-CGDNA/pair_oxrna2_stk.cpp Normal file
View File

File diff suppressed because it is too large Load Diff

86
src/USER-CGDNA/pair_oxrna2_stk.h Normal file
View File

@ -0,0 +1,86 @@
/* -*- c++ -*- ----------------------------------------------------------
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.
------------------------------------------------------------------------- */
#ifdef PAIR_CLASS
PairStyle(oxrna2/stk,PairOxrna2Stk)
#else
#ifndef LMP_PAIR_OXRNA2_STK_H
#define LMP_PAIR_OXRNA2_STK_H
#include "pair.h"
namespace LAMMPS_NS {
class PairOxrna2Stk : public Pair {
public:
PairOxrna2Stk(class LAMMPS *);
virtual ~PairOxrna2Stk();
virtual void compute(int, int);
void settings(int, char **);
void coeff(int, char **);
void init_list(int, class NeighList *);
double init_one(int, int);
void write_restart(FILE *);
void read_restart(FILE *);
void write_restart_settings(FILE *);
void read_restart_settings(FILE *);
void write_data(FILE *);
void write_data_all(FILE *);
void *extract(const char *, int &);
protected:
// stacking interaction
double eta_st[4][4];
double stacking_strength(double, double, double);
double **epsilon_st, **a_st, **cut_st_0, **cut_st_c;
double **cut_st_lo, **cut_st_hi;
double **cut_st_lc, **cut_st_hc, **b_st_lo, **b_st_hi, **shift_st;
double **cutsq_st_hc;
double **a_st5, **theta_st5_0, **dtheta_st5_ast;
double **b_st5, **dtheta_st5_c;
double **a_st6, **theta_st6_0, **dtheta_st6_ast;
double **b_st6, **dtheta_st6_c;
double **a_st9, **theta_st9_0, **dtheta_st9_ast;
double **b_st9, **dtheta_st9_c;
double **a_st10, **theta_st10_0, **dtheta_st10_ast;
double **b_st10, **dtheta_st10_c;
double **a_st1, **cosphi_st1_ast, **b_st1, **cosphi_st1_c;
double **a_st2, **cosphi_st2_ast, **b_st2, **cosphi_st2_c;
int seqdepflag;
virtual void allocate();
void ev_tally_xyz(int, int, int, int, double, double, double, double, double, double, double);
};
}
#endif
#endif
/* ERROR/WARNING messages:
E: Illegal ... command
Self-explanatory. Check the input script syntax and compare to the
documentation for the command. You can use -echo screen as a
command-line option when running LAMMPS to see the offending line.
E: Incorrect args for pair coefficients
Self-explanatory. Check the input script or data file.
*/

1079
src/USER-CGDNA/pair_oxrna2_xstk.cpp Normal file
View File

File diff suppressed because it is too large Load Diff

85
src/USER-CGDNA/pair_oxrna2_xstk.h Normal file
View File

@ -0,0 +1,85 @@
/* -*- c++ -*- ----------------------------------------------------------
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.
------------------------------------------------------------------------- */
#ifdef PAIR_CLASS
PairStyle(oxrna2/xstk,PairOxrna2Xstk)
#else
#ifndef LMP_PAIR_OXRNA2_XSTK_H
#define LMP_PAIR_OXRNA2_XSTK_H
#include "pair.h"
namespace LAMMPS_NS {
class PairOxrna2Xstk : public Pair {
public:
PairOxrna2Xstk(class LAMMPS *);
virtual ~PairOxrna2Xstk();
virtual void compute(int, int);
void settings(int, char **);
void coeff(int, char **);
void init_list(int, class NeighList *);
double init_one(int, int);
void write_restart(FILE *);
void read_restart(FILE *);
void write_restart_settings(FILE *);
void read_restart_settings(FILE *);
void write_data(FILE *);
void write_data_all(FILE *);
void *extract(const char *, int &);
protected:
// cross-stacking interaction
double **k_xst, **cut_xst_0, **cut_xst_c, **cut_xst_lo, **cut_xst_hi;
double **cut_xst_lc, **cut_xst_hc, **b_xst_lo, **b_xst_hi;
double **cutsq_xst_hc;
double **a_xst1, **theta_xst1_0, **dtheta_xst1_ast;
double **b_xst1, **dtheta_xst1_c;
double **a_xst2, **theta_xst2_0, **dtheta_xst2_ast;
double **b_xst2, **dtheta_xst2_c;
double **a_xst3, **theta_xst3_0, **dtheta_xst3_ast;
double **b_xst3, **dtheta_xst3_c;
double **a_xst7, **theta_xst7_0, **dtheta_xst7_ast;
double **b_xst7, **dtheta_xst7_c;
double **a_xst8, **theta_xst8_0, **dtheta_xst8_ast;
double **b_xst8, **dtheta_xst8_c;
virtual void allocate();
};
}
#endif
#endif
/* ERROR/WARNING messages:
E: Illegal ... command
Self-explanatory. Check the input script syntax and compare to the
documentation for the command. You can use -echo screen as a
command-line option when running LAMMPS to see the offending line.
E: Incorrect args for pair coefficients
Self-explanatory. Check the input script or data file.
*/