diff --git a/doc/src/bond_oxdna.rst b/doc/src/bond_oxdna.rst index 1bce606da7..8c2ddfa5a0 100644 --- a/doc/src/bond_oxdna.rst +++ b/doc/src/bond_oxdna.rst @@ -27,6 +27,7 @@ Examples .. code-block:: LAMMPS + # LJ units bond_style oxdna/fene bond_coeff * 2.0 0.25 0.7525 @@ -36,6 +37,30 @@ Examples bond_style oxrna2/fene bond_coeff * 2.0 0.25 0.76107 + bond_style oxdna/fene + bond_coeff * oxdna.lj + + # Real units + bond_style oxdna/fene + bond_coeff * 11.92337812042065 2.1295 6.409795 + + bond_style oxdna2/fene + bond_coeff * 11.92337812042065 2.1295 6.4430152 + + bond_style oxrna2/fene + bond_coeff * 11.92337812042065 2.1295 6.482800913 + + bond_style oxrna2/fene + bond_coeff * oxrna2.real + +.. note:: + The coefficients in the above examples have to be kept fixed and cannot + be changed without reparameterizing the entire model. They are provided in forms + compatible with both *units lj* and *units real* (see documentation of :doc:`units `). + These can also be read from a potential file with correct unit style by specifying the name + of the file. Several potential files for each unit style are included in the + /potentials/ directory of the LAMMPS distribution. + Description """"""""""" @@ -73,8 +98,6 @@ commands: *oxdna2/hbond* and an additional Debye-Hueckel pair style *oxdna2/dh* have to be defined. The same applies to the oxRNA2 :ref:`(Sulc1) ` styles. - The coefficients in the above example have to be kept fixed and cannot - be changed without reparameterizing the entire model. .. note:: @@ -113,6 +136,29 @@ and for sequence-specific hydrogen-bonding and stacking interactions ---------- +Potential file reading +"""""""""""""""""""""" + +For each style oxdna, oxdna2 and oxrna2, the first parameter argument can be a filename, and if it is, no further arguments should be supplied. Therefore the following command: + +.. code-block:: LAMMPS + + bond_style oxdna/fene + bond_coeff * oxdna.lj + +will be interpreted as a request to read the (FENE) potential :ref:`(Ouldridge) ` parameters from the file with the given name. +The file can define multiple potential parameters for both bonded and pair interactions, but for the above bonded interactions there must exist in the file a line of the form: + +.. code-block:: LAMMPS + + * fene epsilon delta r0 + +There are sample potential files for each unit style in the /potentials/ directory of the LAMMPS distribution. The potential file unit system must align with +the units defined via the :doc:`units ` command. For conversion between different *LJ* and *real* unit systems for oxDNA, the python tool *lj2real.py* located in the examples/PACKAGES/cgdna/util/ +directory can be used. This tool assumes similar file structure to the examples found in examples/PACKAGES/cgdna/examples/. + +---------- + Restrictions """""""""""" diff --git a/doc/src/pair_oxdna.rst b/doc/src/pair_oxdna.rst index dab1c2a230..2781cffab8 100644 --- a/doc/src/pair_oxdna.rst +++ b/doc/src/pair_oxdna.rst @@ -49,6 +49,7 @@ Examples .. code-block:: LAMMPS + # LJ units pair_style hybrid/overlay oxdna/excv oxdna/stk oxdna/hbond oxdna/xstk oxdna/coaxstk pair_coeff * * oxdna/excv 2.0 0.7 0.675 2.0 0.515 0.5 2.0 0.33 0.32 pair_coeff * * oxdna/stk seqdep 0.1 1.3448 2.6568 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 @@ -58,6 +59,40 @@ Examples pair_coeff * * oxdna/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 * * oxdna/coaxstk 46.0 0.4 0.6 0.22 0.58 2.0 2.541592653589793 0.65 1.3 0 0.8 0.9 0 0.95 0.9 0 0.95 2.0 -0.65 2.0 -0.65 + pair_style hybrid/overlay oxdna/excv oxdna/stk oxdna/hbond oxdna/xstk oxdna/coaxstk + pair_coeff * * oxdna/excv oxdna.lj + pair_coeff * * oxdna/stk seqav 0.1 1.3448 2.6568 oxdna.lj + pair_coeff * * oxdna/hbond seqav oxdna.lj + pair_coeff 1 4 oxdna/hbond seqav oxdna.lj + pair_coeff 2 3 oxdna/hbond seqav oxdna.lj + pair_coeff * * oxdna/xstk oxdna.lj + pair_coeff * * oxdna/coaxstk oxdna.lj + + # Real units + pair_style hybrid/overlay oxdna/excv oxdna/stk oxdna/hbond oxdna/xstk oxdna/coaxstk + pair_coeff * * oxdna/excv 11.92337812042065 5.9626 5.74965 11.92337812042065 4.38677 4.259 11.92337812042065 2.81094 2.72576 + pair_coeff * * oxdna/stk seqdep 300.0 8.01727944817084 0.005279604 0.70439070204273 3.4072 7.6662 2.72576 6.3885 1.3 0.0 0.8 0.9 0.0 0.95 0.9 0.0 0.95 2.0 0.65 2.0 0.65 + pair_coeff * * oxdna/hbond seqdep 0.0 0.93918760272364 3.4072 6.3885 2.89612 5.9626 1.5 0.0 0.7 1.5 0.0 0.7 1.5 0.0 0.7 0.46 3.141592654 0.7 4.0 1.570796327 0.45 4.0 1.570796327 0.45 + pair_coeff 1 4 oxdna/hbond seqdep 6.42073911784652 0.93918760272364 3.4072 6.3885 2.89612 5.9626 1.5 0.0 0.7 1.5 0.0 0.7 1.5 0.0 0.7 0.46 3.141592654 0.7 4.0 1.570796327 0.45 4.0 1.570796327 0.45 + pair_coeff 2 3 oxdna/hbond seqdep 6.42073911784652 0.93918760272364 3.4072 6.3885 2.89612 5.9626 1.5 0.0 0.7 1.5 0.0 0.7 1.5 0.0 0.7 0.46 3.141592654 0.7 4.0 1.570796327 0.45 4.0 1.570796327 0.45 + pair_coeff * * oxdna/xstk 3.9029021145006 4.89785 5.74965 4.21641 5.57929 2.25 0.791592654 0.58 1.7 1.0 0.68 1.7 1.0 0.68 1.5 0.0 0.65 1.7 0.875 0.68 1.7 0.875 0.68 + pair_coeff * * oxdna/coaxstk 3.77965257404268 3.4072 5.1108 1.87396 4.94044 2.0 2.541592654 0.65 1.3 0.0 0.8 0.9 0.0 0.95 0.9 0.0 0.95 2.0 -0.65 2.0 -0.65 + + pair_style hybrid/overlay oxdna/excv oxdna/stk oxdna/hbond oxdna/xstk oxdna/coaxstk + pair_coeff * * oxdna/excv oxdna.real + pair_coeff * * oxdna/stk seqav 300.0 8.01727944817084 0.005279604 oxdna.real + pair_coeff * * oxdna/hbond seqav oxdna.real + pair_coeff 1 4 oxdna/hbond seqav oxdna.real + pair_coeff 2 3 oxdna/hbond seqav oxdna.real + pair_coeff * * oxdna/xstk oxdna.real + pair_coeff * * oxdna/coaxstk oxdna.real + +.. note:: + + The coefficients in the above examples are provided in forms compatible with both *units lj* and *units real* (see documentation of :doc:`units `). + These can also be read from a potential file with correct unit style by specifying the name of the file. Several potential files for each unit style are included in the + /potentials/ directory of the LAMMPS distribution. + Description """"""""""" @@ -85,7 +120,7 @@ for a detailed description of the oxDNA force field. *oxdna/fene* for the connectivity of the phosphate backbone (see also documentation of :doc:`bond_style oxdna/fene `). 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 *oxdna/stk* (seq=seqdep, T=0.1, xi=1.3448 and kappa=2.6568 in the above example) + Exceptions are the first four coefficients after *oxdna/stk* (seq=seqdep, T=0.1, xi=1.3448 and kappa=2.6568 and corresponding *real unit* equivalents in the above examples) and the first coefficient after *oxdna/hbond* (seq=seqdep in the above example). When using a Langevin thermostat, e.g. through :doc:`fix langevin ` or :doc:`fix nve/dotc/langevin ` @@ -115,6 +150,43 @@ and :ref:`(Sulc) `. ---------- +Potential file reading +"""""""""""""""""""""" + +For each pair style above the first non-modifiable argument can be a filename, and if it is, no further arguments should be supplied. Therefore the following command: + +.. code-block:: LAMMPS + + pair_coeff 1 4 oxdna/hbond seqav oxdna.lj + +will be interpreted as a request to read the corresponding hydrogen bonding potential parameters from the file with the given name. +The file can define multiple potential parameters for both bonded and pair interactions, but for the example pair interaction above there must exist in the file a line of the form: + +.. code-block:: LAMMPS + + 1 4 hbond + +If potential customization is required, the potential file reading can be mixed with the manual specification of the potential parameters. For example, the following command: + +.. code-block:: LAMMPS + + pair_style hybrid/overlay oxdna/excv oxdna/stk oxdna/hbond oxdna/xstk oxdna/coaxstk + pair_coeff * * oxdna/excv oxdna.lj + pair_coeff * * oxdna/stk seqav 0.1 1.3448 2.6568 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 * * oxdna/hbond seqav oxdna.lj + pair_coeff 1 4 oxdna/hbond seqav oxdna.lj + pair_coeff 2 3 oxdna/hbond seqav oxdna.lj + pair_coeff * * oxdna/xstk oxdna.lj + pair_coeff * * oxdna/coaxstk 46.0 0.4 0.6 0.22 0.58 2.0 2.541592653589793 0.65 1.3 0 0.8 0.9 0 0.95 0.9 0 0.95 2.0 -0.65 2.0 -0.65 + +will read the stacking and coaxial stacking potential parameters from the manual specification and all others from the potential file *oxdna.lj*. + +There are sample potential files for each unit style in the /potentials/ directory of the LAMMPS distribution. The potential file unit system must align with +the units defined via the :doc:`units ` command. For conversion between different *LJ* and *real* unit systems for oxDNA, the python tool *lj2real.py* located in the examples/PACKAGES/cgdna/util/ +directory can be used. This tool assumes similar file structure to the examples found in examples/PACKAGES/cgdna/examples/. + +---------- + Restrictions """""""""""" diff --git a/doc/src/pair_oxdna2.rst b/doc/src/pair_oxdna2.rst index 5cac7d8f77..d1d9d741e2 100644 --- a/doc/src/pair_oxdna2.rst +++ b/doc/src/pair_oxdna2.rst @@ -57,6 +57,7 @@ Examples .. code-block:: LAMMPS + # LJ units 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 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 @@ -67,6 +68,43 @@ Examples 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 0.5 0.815 + pair_style hybrid/overlay oxdna2/excv oxdna2/stk oxdna2/hbond oxdna2/xstk oxdna2/coaxstk oxdna2/dh + pair_coeff * * oxdna2/excv oxdna2.lj + pair_coeff * * oxdna2/stk seqdep 0.1 1.3523 2.6717 oxdna2.lj + pair_coeff * * oxdna2/hbond seqdep oxdna2.lj + pair_coeff 1 4 oxdna2/hbond seqdep oxdna2.lj + pair_coeff 2 3 oxdna2/hbond seqdep oxdna2.lj + pair_coeff * * oxdna2/xstk oxdna2.lj + pair_coeff * * oxdna2/coaxstk oxdna2.lj + pair_coeff * * oxdna2/dh 0.1 0.5 oxdna2.lj + + # Real units + pair_style hybrid/overlay oxdna2/excv oxdna2/stk oxdna2/hbond oxdna2/xstk oxdna2/coaxstk oxdna2/dh + pair_coeff * * oxdna2/excv 11.92337812042065 5.9626 5.74965 11.92337812042065 4.38677 4.259 11.92337812042065 2.81094 2.72576 + pair_coeff * * oxdna2/stk seqdep 300.0 8.06199211612242 0.005309213 0.70439070204273 3.4072 7.6662 2.72576 6.3885 1.3 0.0 0.8 0.9 0.0 0.95 0.9 0.0 0.95 2.0 0.65 2.0 0.65 + pair_coeff * * oxdna2/hbond seqdep 0.0 0.93918760272364 3.4072 6.3885 2.89612 5.9626 1.5 0.0 0.7 1.5 0.0 0.7 1.5 0.0 0.7 0.46 3.141592654 0.7 4.0 1.570796327 0.45 4.0 1.570796327 0.45 + pair_coeff 1 4 oxdna2/hbond seqdep 6.36589157849259 0.93918760272364 3.4072 6.3885 2.89612 5.9626 1.5 0.0 0.7 1.5 0.0 0.7 1.5 0.0 0.7 0.46 3.141592654 0.7 4.0 1.570796327 0.45 4.0 1.570796327 0.45 + pair_coeff 2 3 oxdna2/hbond seqdep 6.36589157849259 0.93918760272364 3.4072 6.3885 2.89612 5.9626 1.5 0.0 0.7 1.5 0.0 0.7 1.5 0.0 0.7 0.46 3.141592654 0.7 4.0 1.570796327 0.45 4.0 1.570796327 0.45 + pair_coeff * * oxdna2/xstk 3.9029021145006 4.89785 5.74965 4.21641 5.57929 2.25 0.791592654 0.58 1.7 1.0 0.68 1.7 1.0 0.68 1.5 0.0 0.65 1.7 0.875 0.68 1.7 0.875 0.68 + pair_coeff * * oxdna2/coaxstk 4.80673207785863 3.4072 5.1108 1.87396 4.94044 2.0 2.891592653589793 0.65 1.3 0.0 0.8 0.9 0.0 0.95 0.9 0.0 0.95 40.0 3.116592653589793 + pair_coeff * * oxdna2/dh 300.0 0.5 0.815 + + pair_style hybrid/overlay oxdna2/excv oxdna2/stk oxdna2/hbond oxdna2/xstk oxdna2/coaxstk oxdna2/dh + pair_coeff * * oxdna2/excv oxdna2.real + pair_coeff * * oxdna2/stk seqdep 300.0 8.06199211612242 0.005309213 oxdna2.real + pair_coeff * * oxdna2/hbond seqdep oxdna2.real + pair_coeff 1 4 oxdna2/hbond seqdep oxdna2.real + pair_coeff 2 3 oxdna2/hbond seqdep oxdna2.real + pair_coeff * * oxdna2/xstk oxdna2.real + pair_coeff * * oxdna2/coaxstk oxdna2.real + pair_coeff * * oxdna2/dh 300.0 0.5 oxdna2.real + +.. note:: + + The coefficients in the above examples are provided in forms compatible with both *units lj* and *units real* (see documentation of :doc:`units `). + These can also be read from a potential file with correct unit style by specifying the name of the file. Several potential files for each unit style are included in the + /potentials/ directory of the LAMMPS distribution. + Description """"""""""" @@ -94,7 +132,7 @@ and :ref:`(Ouldridge) ` for a detailed description of the oxDNA2 fo *oxdna2/fene* for the connectivity of the phosphate backbone (see also documentation of :doc:`bond_style oxdna2/fene `). 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 *oxdna2/stk* (seq=seqdep, T=0.1, xi=1.3523 and kappa=2.6717 in the above example), + Exceptions are the first four coefficients after *oxdna2/stk* (seq=seqdep, T=0.1, xi=1.3523 and kappa=2.6717 and corresponding *real unit* equivalents in the above examples). the first coefficient after *oxdna2/hbond* (seq=seqdep in the above example) and the three coefficients 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 ` or :doc:`fix nve/dotc/langevin ` @@ -122,6 +160,45 @@ Please cite also the relevant oxDNA2 publications ---------- +Potential file reading +"""""""""""""""""""""" + +For each pair style above the first non-modifiable argument can be a filename (with exception of Debye-Hueckel, for which the effective charge argument can be a filename), and if it is, no further arguments should be supplied. +Therefore the following command: + +.. code-block:: LAMMPS + + pair_coeff 1 4 oxdna2/hbond seqdep oxdna.real + +will be interpreted as a request to read the corresponding hydrogen bonding potential parameters from the file with the given name. +The file can define multiple potential parameters for both bonded and pair interactions, but for the example pair interaction above there must exist in the file a line of the form: + +.. code-block:: LAMMPS + + 1 4 hbond + +If potential customization is required, the potential file reading can be mixed with the manual specification of the potential parameters. For example, the following command: + +.. code-block:: LAMMPS + + 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 0.1 1.3523 2.6717 oxdna2.lj + pair_coeff * * oxdna2/hbond seqdep oxdna2.lj + pair_coeff 1 4 oxdna2/hbond seqdep oxdna2.lj + pair_coeff 2 3 oxdna2/hbond seqdep oxdna2.lj + pair_coeff * * oxdna2/xstk oxdna2.lj + pair_coeff * * oxdna2/coaxstk oxdna2.lj + pair_coeff * * oxdna2/dh 0.1 0.5 0.815 + +will read the excluded volume and Debye-Hueckel effective charge *qeff* parameters from the manual specification and all others from the potential file *oxdna2.lj*. + +There are sample potential files for each unit style in the /potentials/ directory of the LAMMPS distribution. The potential file unit system must align with +the units defined via the :doc:`units ` command. For conversion between different *LJ* and *real* unit systems for oxDNA, the python tool *lj2real.py* located in the examples/PACKAGES/cgdna/util/ +directory can be used. This tool assumes similar file structure to the examples found in examples/PACKAGES/cgdna/examples/. + +---------- + Restrictions """""""""""" diff --git a/doc/src/pair_oxrna2.rst b/doc/src/pair_oxrna2.rst index 7b8220740c..1726e3cbff 100644 --- a/doc/src/pair_oxrna2.rst +++ b/doc/src/pair_oxrna2.rst @@ -57,6 +57,7 @@ Examples .. code-block:: LAMMPS + # LJ units 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 @@ -68,11 +69,51 @@ Examples 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 + pair_style hybrid/overlay oxrna2/excv oxrna2/stk oxrna2/hbond oxrna2/xstk oxrna2/coaxstk oxrna2/dh + pair_coeff * * oxrna2/excv oxrna2.lj + pair_coeff * * oxrna2/stk seqdep 0.1 1.40206 2.77 oxrna2.lj + pair_coeff * * oxrna2/hbond seqdep oxrna2.lj + pair_coeff 1 4 oxrna2/hbond seqdep oxrna2.lj + pair_coeff 2 3 oxrna2/hbond seqdep oxrna2.lj + pair_coeff 3 4 oxrna2/hbond seqdep oxrna2.lj + pair_coeff * * oxrna2/xstk oxrna2.lj + pair_coeff * * oxrna2/coaxstk oxrna2.lj + pair_coeff * * oxrna2/dh 0.1 0.5 oxrna2.lj + + # Real units + pair_style hybrid/overlay oxrna2/excv oxrna2/stk oxrna2/hbond oxrna2/xstk oxrna2/coaxstk oxrna2/dh + pair_coeff * * oxrna2/excv 11.92337812042065 5.9626 5.74965 11.92337812042065 4.38677 4.259 11.92337812042065 2.81094 2.72576 + pair_coeff * * oxrna2/stk seqdep 300.0 8.35864576375849 0.005504556 0.70439070204273 3.66274 7.92174 2.9813 6.64404 0.9 0.0 0.95 0.9 0.0 0.95 1.3 0.0 0.8 1.3 0.0 0.8 2.0 0.65 2.0 0.65 + pair_coeff * * oxrna2/hbond seqdep 0.0 0.93918760272364 3.4072 6.3885 2.89612 5.9626 1.5 0.0 0.7 1.5 0.0 0.7 1.5 0.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 5.18928666388042 0.93918760272364 3.4072 6.3885 2.89612 5.9626 1.5 0.0 0.7 1.5 0.0 0.7 1.5 0.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 5.18928666388042 0.93918760272364 3.4072 6.3885 2.89612 5.9626 1.5 0.0 0.7 1.5 0.0 0.7 1.5 0.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 5.18928666388042 0.93918760272364 3.4072 6.3885 2.89612 5.9626 1.5 0.0 0.7 1.5 0.0 0.7 1.5 0.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 4.92690859644113 4.259 5.1108 3.57756 4.94044 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 6.57330882442206 4.259 5.1108 3.57756 4.94044 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 300.0 0.5 1.02455 + + pair_style hybrid/overlay oxrna2/excv oxrna2/stk oxrna2/hbond oxrna2/xstk oxrna2/coaxstk oxrna2/dh + pair_coeff * * oxrna2/excv oxrna2.real + pair_coeff * * oxrna2/stk seqdep 300.0 8.35864576375849 0.005504556 oxrna2.real + pair_coeff * * oxrna2/hbond seqdep oxrna2.real + pair_coeff 1 4 oxrna2/hbond seqdep oxrna2.real + pair_coeff 2 3 oxrna2/hbond seqdep oxrna2.real + pair_coeff 3 4 oxrna2/hbond seqdep oxrna2.real + pair_coeff * * oxrna2/xstk oxrna2.real + pair_coeff * * oxrna2/coaxstk oxrna2.real + pair_coeff * * oxrna2/dh 300.0 0.5 oxrna2.real + +.. note:: + + The coefficients in the above examples are provided in forms compatible with both *units lj* and *units real* (see documentation of :doc:`units `). + These can also be read from a potential file with correct unit style by specifying the name of the file. Several potential files for each unit style are included in the + /potentials/ directory of the LAMMPS distribution. + 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 +for coarse-grained modelling of RNA. 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 @@ -95,7 +136,7 @@ and :ref:`(Ouldridge) ` for a detailed description of the oxRNA2 fo *oxrna2/fene* for the connectivity of the phosphate backbone (see also documentation of :doc:`bond_style oxrna2/fene `). 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), + Exceptions are the first four coefficients after *oxrna2/stk* (seq=seqdep, T=0.1, xi=1.40206 and kappa=2.77 and corresponding *real unit* equivalents in the above examples), 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 ` or :doc:`fix nve/dotc/langevin ` @@ -123,6 +164,46 @@ Please cite also the relevant oxRNA2 publications ---------- +Potential file reading +"""""""""""""""""""""" + +For each pair style above the first non-modifiable argument can be a filename (with exception of Debye-Hueckel, for which the effective charge argument can be a filename), and if it is, no further arguments should be supplied. +Therefore the following command: + +.. code-block:: LAMMPS + + pair_coeff 3 4 oxrna2/hbond seqdep oxrna2.lj + +will be interpreted as a request to read the corresponding hydrogen bonding potential parameters from the file with the given name. +The file can define multiple potential parameters for both bonded and pair interactions, but for the example pair interaction above there must exist in the file a line of the form: + +.. code-block:: LAMMPS + + 3 4 hbond + +If potential customization is required, the potential file reading can be mixed with the manual specification of the potential parameters. For example, the following command: + +.. 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 oxrna2.lj + pair_coeff * * oxrna2/hbond seqdep oxrna2.lj + pair_coeff 1 4 oxrna2/hbond seqdep oxrna2.lj + pair_coeff 2 3 oxrna2/hbond seqdep oxrna2.lj + pair_coeff 3 4 oxrna2/hbond seqdep oxrna2.lj + pair_coeff * * oxrna2/xstk oxrna2.lj + pair_coeff * * oxrna2/coaxstk oxrna2.lj + pair_coeff * * oxrna2/dh 0.1 0.5 1.02455 + +will read the excluded volume and Debye-Hueckel effective charge *qeff* parameters from the manual specification and all others from the potential file *oxrna2.lj*. + +There are sample potential files for each unit style in the /potentials/ directory of the LAMMPS distribution. The potential file unit system must align with +the units defined via the :doc:`units ` command. For conversion between different *LJ* and *real* unit systems for oxDNA, the python tool *lj2real.py* located in the examples/PACKAGES/cgdna/util/ +directory can be used. This tool assumes similar file structure to the examples found in examples/PACKAGES/cgdna/examples/. + +---------- + Restrictions """""""""""" diff --git a/examples/PACKAGES/cgdna/util/lj2real.py b/examples/PACKAGES/cgdna/util/lj2real.py new file mode 100644 index 0000000000..7b4f3bac28 --- /dev/null +++ b/examples/PACKAGES/cgdna/util/lj2real.py @@ -0,0 +1,604 @@ +""" +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + https://www.lammps.org/ 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 authors: Kierran Falloon (University of Strathclyde, Glasgow) + Oliver Henrich (University of Strathclyde, Glasgow) +------------------------------------------------------------------------- */ + +Program: + +Usage: +$$ python lj2real.py [-i] +$$ [-i] flag is optional and is used to convert real -> LJ units. + +Requirements: +LAMMPS data file and input file. + +This script assumes a input and data file structure similar to those found in examples/PACKAGES/cgdna/examples/. +""" + +import datetime +import os +import sys + + +class Sections: + """Sections of the data file""" + + def __init__( + self, + bounds: bool, + masses: bool, + atoms: bool, + velocities: bool, + ellipsoids: bool, + ): + self.bounds = bounds # xlo, xhi, ylo, yhi, zlo, zhi + self.masses = masses # Masses + self.atoms = atoms # Atoms + self.velocities = velocities # Velocities + self.ellipsoids = ellipsoids # Ellipsoids + + +# Conversion factors +class ConversionFactors: + """Conversion factors for LJ to real units""" + + def __init__(self, invert: bool = False): + self.inverted = False + self.temp_conv_factor = 3000.038822 + self.energy_conv_factor = 5.961689060210325 + self.kT_conv_factor = 0.001987204155 + self.mass_conv_factor = 100.0277580236 + self.length_conv_factor = 8.518 + self.time_conv_factor = 1706.0 + self.vel_conv_factor = 0.004992966002344666 + self.angular_mom_conv_factor = 4.254188991883894 + self.density_conv_factor = 0.2687551067436886 + + self.oxdna_fene_string = "11.92337812042065 2.1295 6.409795" + self.oxdna_excv_string = "11.92337812042065 5.9626 5.74965 11.92337812042065 4.38677 4.259 11.92337812042065 2.81094 2.72576" + self.oxdna_stk_string = "0.70439070204273 3.4072 7.6662 2.72576 6.3885 1.3 0.0 0.8 0.9 0.0 0.95 0.9 0.0 0.95 2.0 0.65 2.0 0.65" + self.oxdna_hbond_string = "0.0 0.93918760272364 3.4072 6.3885 2.89612 5.9626 1.5 0.0 0.7 1.5 0 0.7 1.5 0 0.7 0.46 3.141592654 0.7 4.0 1.570796327 0.45 4.0 1.570796327 0.45" + self.oxdna_hbond_1_4_2_3_string = "6.42073911784652 0.93918760272364 3.4072 6.3885 2.89612 5.9626 1.5 0 0.7 1.5 0.0 0.7 1.5 0 0.7 0.46 3.141592654 0.7 4.0 1.570796327 0.45 4.0 1.570796327 0.45" + self.oxdna_xstk_string = "3.9029021145006 4.89785 5.74965 4.21641 5.57929 2.25 0.791592654 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" + self.oxdna_coaxstk_string = "3.77965257404268 3.4072 5.1108 1.87396 4.94044 2.0 2.541592654 0.65 1.3 0 0.8 0.9 0 0.95 0.9 0 0.95 2 -0.65 2 -0.65" + + self.oxdna2_fene_string = "11.92337812042065 2.1295 6.4430152" + self.oxdna2_excv_string = "11.92337812042065 5.9626 5.74965 11.92337812042065 4.38677 4.259 11.92337812042065 2.81094 2.72576" + self.oxdna2_stk_string = "0.70439070204273 3.4072 7.6662 2.72576 6.3885 1.3 0.0 0.8 0.9 0.0 0.95 0.9 0.0 0.95 2.0 0.65 2.0 0.65" + self.oxdna2_hbond_string = "0.0 0.93918760272364 3.4072 6.3885 2.89612 5.9626 1.5 0.0 0.7 1.5 0 0.7 1.5 0 0.7 0.46 3.141592654 0.7 4.0 1.570796327 0.45 4.0 1.570796327 0.45" + self.oxdna2_hbond_1_4_2_3_string = "6.36589157849259 0.93918760272364 3.4072 6.3885 2.89612 5.9626 1.5 0 0.7 1.5 0.0 0.7 1.5 0 0.7 0.46 3.141592654 0.7 4.0 1.570796327 0.45 4.0 1.570796327 0.45" + self.oxdna2_xstk_string = "3.9029021145006 4.89785 5.74965 4.21641 5.57929 2.25 0.791592654 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" + self.oxdna2_coaxstk_string = "4.80673207785863 3.4072 5.1108 1.87396 4.94044 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" + + self.oxrna2_fene_string = "11.92337812042065 2.1295 6.482800913" + self.oxrna2_excv_string = "11.92337812042065 5.9626 5.74965 11.92337812042065 4.38677 4.259 11.92337812042065 2.81094 2.72576" + self.oxrna2_stk_string = "0.70439070204273 3.66274 7.92174 2.9813 6.64404 0.9 0.0 0.95 0.9 0.0 0.95 1.3 0.0 0.8 1.3 0.0 0.8 2.0 0.65 2.0 0.65" + self.oxrna2_hbond_string = "0.0 0.93918760272364 3.4072 6.3885 2.89612 5.9626 1.5 0.0 0.7 1.5 0.0 0.7 1.5 0.0 0.7 0.46 3.141592653589793 0.7 4.0 1.5707963267948966 0.45 4.0 1.5707963267948966 0.45" + self.oxrna2_hbond_1_4_2_3_3_4_string = "5.18928666388042 0.93918760272364 3.4072 6.3885 2.89612 5.9626 1.5 0.0 0.7 1.5 0.0 0.7 1.5 0.0 0.7 0.46 3.141592653589793 0.7 4.0 1.5707963267948966 0.45 4.0 1.5707963267948966 0.45" + self.oxrna2_xstk_string = "4.92690859644113 4.259 5.1108 3.57756 4.94044 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" + self.oxrna2_coaxstk_string = "6.57330882442206 4.259 5.1108 3.57756 4.94044 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" + + if invert: + self.invert() + + def invert(self): + """Inverts the conversion factors for real -> LJ""" + self.inverted = True + self.temp_conv_factor = 1.0 / self.temp_conv_factor + self.energy_conv_factor = 1.0 / self.energy_conv_factor + self.kT_conv_factor = 1.0 / self.kT_conv_factor + self.mass_conv_factor = 1.0 / self.mass_conv_factor + self.length_conv_factor = 1.0 / self.length_conv_factor + self.time_conv_factor = 1.0 / self.time_conv_factor + self.vel_conv_factor = 1.0 / self.vel_conv_factor + self.angular_mom_conv_factor = 1.0 / self.angular_mom_conv_factor + self.density_conv_factor = 1.0 / self.density_conv_factor + + self.oxdna_fene_string = "2.0 0.25 0.7525" + self.oxdna_excv_string = "2.0 0.7 0.675 2.0 0.515 0.5 2.0 0.33 0.32" + self.oxdna_stk_string = ( + "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" + ) + self.oxdna_hbond_string = "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" + self.oxdna_hbond_1_4_2_3_string = "1.077 8.0 0.4 0.75 0.34 0.7 1.5 0 0.7 1.5 0 0.7 1.5 0 0.7 0.46 3.141592653589793 0.7 4.0 1.5707963267948966 0.45 4.0 1.5707963267948966 0.45" + self.oxdna_xstk_string = "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" + self.oxdna_coaxstk_string = "46.0 0.4 0.6 0.22 0.58 2.0 2.541592653589793 0.65 1.3 0 0.8 0.9 0 0.95 0.9 0 0.95 2.0 -0.65 2.0 -0.65" + + self.oxdna2_fene_string = "2.0 0.25 0.7564" + self.oxdna2_excv_string = "2.0 0.7 0.675 2.0 0.515 0.5 2.0 0.33 0.32" + self.oxdna2_stk_string = ( + "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" + ) + self.oxdna2_hbond_string = "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" + self.oxdna2_hbond_1_4_2_3_string = "1.0678 8.0 0.4 0.75 0.34 0.7 1.5 0 0.7 1.5 0 0.7 1.5 0 0.7 0.46 3.141592653589793 0.7 4.0 1.5707963267948966 0.45 4.0 1.5707963267948966 0.45" + self.oxdna2_xstk_string = "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" + self.oxdna2_coaxstk_string = "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" + + self.oxrna2_fene_string = "2.0 0.25 0.761070781051" + self.oxrna2_excv_string = "2.0 0.7 0.675 2.0 0.515 0.5 2.0 0.33 0.32" + self.oxrna2_stk_string = "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" + self.oxrna2_hbond_string = "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" + self.oxrna2_hbond_1_4_2_3_3_4_string = "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" + self.oxrna2_xstk_string = "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" + self.oxrna2_coaxstk_string = "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" + + +def check_datafile_header(line: str, sections: Sections): + """Checks for headers to modify corresponding data, since datafile is split into headers. + Modifies the Sections object to keep track of the current section. + + Args: + line (str): The line to check + masses_section (bool): If the current section is the masses section + atoms_section (bool): If the current section is the atoms section + velocities_section (bool): If the current section is the velocities section + ellipsoids_section (bool): If the current section is the ellipsoids section + """ + + if any(header in line for header in ["xlo", "xhi", "ylo", "yhi", "zlo", "zhi"]): + sections.bounds = True + sections.masses = False + sections.atoms = False + sections.velocities = False + sections.ellipsoids = False + elif "Masses" in line: + sections.bounds = False + sections.masses = True + sections.atoms = False + sections.velocities = False + sections.ellipsoids = False + elif "Atoms" in line: + sections.bounds = False + sections.masses = False + sections.atoms = True + sections.velocities = False + sections.ellipsoids = False + elif "Velocities" in line: + sections.bounds = False + sections.masses = False + sections.atoms = False + sections.velocities = True + sections.ellipsoids = False + elif "Ellipsoids" in line: + sections.bounds = False + sections.masses = False + sections.atoms = False + sections.velocities = False + sections.ellipsoids = True + elif "Bonds" in line: + sections.bounds = False + sections.masses = False + sections.atoms = False + sections.velocities = False + sections.ellipsoids = False + + +def modify_datafile(datafile_path: str, conversion_factors: ConversionFactors): + """Modifies the file by header to use real units. + + Args: + datafile_path (str): The path to the file to modify + """ + lines_changed = 0 + current_section = Sections(False, False, False, False, False) + + with open(datafile_path, "r", encoding="UTF-8") as file: + lines = file.readlines() + if conversion_factors.inverted: + lines[0] = ( + "LAMMPS data file in LJ units via oxdna lj2real.py, date " + + str(datetime.date.today()) + + "\n" + ) + else: + lines[0] = ( + "LAMMPS data file in real units via oxdna lj2real.py, date " + + str(datetime.date.today()) + + "\n" + ) + + for i, line in enumerate(lines): + check_datafile_header(line, current_section) # check for headers + + elements = line.split() + if ( + not elements + or elements[0] == "#" + or any( + header in line + for header in ["Masses", "Atoms", "Velocities", "Ellipsoids", "Bonds"] + ) + ): + continue + + # modify the line based on the current section it is in + if current_section.bounds: + elements[0:2] = [ + str(int(float(x) * conversion_factors.length_conv_factor)) + for x in elements[0:2] + ] + lines[i] = " ".join(elements) + "\n" + lines_changed += 1 + if current_section.masses: + elements[1] = str( + round(float(elements[1]) * conversion_factors.mass_conv_factor, 4) + ) + lines[i] = " ".join(elements) + "\n" + lines_changed += 1 + elif current_section.atoms: + elements[2:5] = [ + str(float(x) * conversion_factors.length_conv_factor) + for x in elements[2:5] + ] + elements[7] = str( + float(elements[7]) * conversion_factors.density_conv_factor + ) + lines[i] = " ".join(elements) + "\n" + lines_changed += 1 + elif current_section.velocities: + elements[1:4] = [ + str(float(x) * conversion_factors.vel_conv_factor) + for x in elements[1:4] + ] + elements[4:7] = [ + str(float(x) * conversion_factors.angular_mom_conv_factor) + for x in elements[4:7] + ] + lines[i] = " ".join(elements) + "\n" + lines_changed += 1 + elif current_section.ellipsoids: + elements[1:4] = [ + str(float(x) * conversion_factors.length_conv_factor) + for x in elements[1:4] + ] + lines[i] = " ".join(elements) + "\n" + lines_changed += 1 + + if conversion_factors.inverted: + new_datafile_path = datafile_path + "_lj" + else: + new_datafile_path = datafile_path + "_real" + + with open(new_datafile_path, "w", encoding="UTF-8") as file: + file.writelines(lines) + if lines_changed == 0: + print( + "Warning: No lines changed in data file. Ensure correct usage: python lj2real.py [-i]" + ) + else: + print(f"Data file lines changed: {lines_changed}") + + return new_datafile_path + + +def modify_inputfile(inputfile_path: str, conversion_factors: ConversionFactors): + """Modifies the input file line by line to use real units. + + Args: + inputfile_path (str): The path to the input file to modify + """ + + lines_changed = 0 + oxdna2_flag, oxrna2_flag = False, False + + with open(inputfile_path, "r", encoding="UTF-8") as file: + lines = file.readlines() + + for i, line in enumerate(lines): + if "oxdna2" in line and not oxdna2_flag: + oxdna2_flag = True + print("Note: oxdna2 found in input file. Using oxdna2 conversion factors.") + if "oxrna2" in line and not oxrna2_flag: + oxrna2_flag = True + print("Note: oxrna2 found in input file. Using oxrna2 conversion factors.") + if oxdna2_flag and oxrna2_flag: + print( + "Warning: Both oxdna2 and oxrna2 found in input file. Output will likely be incorrect." + ) + + if "variable T" in line: + old_value = line.split()[3] + + new_value = str( + round(float(old_value) * conversion_factors.temp_conv_factor, 1) + ) + lines[i] = line.replace(old_value, new_value) + lines_changed += 1 + + elif "units" in line: + if conversion_factors.inverted: + lines[i] = "units lj\n" + else: + lines[i] = "units real\n" + lines_changed += 1 + + elif "atom_modify" in line: + elements = line.split() + elements[3] = str( + round(float(elements[3]) * conversion_factors.length_conv_factor, 3) + ) + lines[i] = " ".join(elements) + "\n" + lines_changed += 1 + + elif "neighbor" in line: + elements = line.split() + elements[1] = str( + round(float(elements[1]) * conversion_factors.length_conv_factor, 3) + ) + lines[i] = " ".join(elements) + "\n" + lines_changed += 1 + + elif "read_data" in line: + elements = line.split() + if conversion_factors.inverted: + elements[1] = elements[1] + "_lj" + else: + elements[1] = ( + elements[1] + "_real" + ) # naming convention of datafile after conversion + lines[i] = " ".join(elements) + "\n" + lines_changed += 1 + + elif "mass" in line: + elements = line.split() + elements[4] = str( + round(float(elements[4]) * conversion_factors.mass_conv_factor, 4) + ) + lines[i] = " ".join(elements) + "\n" + lines_changed += 1 + + elif "bond_coeff" in line or "pair_coeff" in line: + if ".lj" in line or ".real" in line: + if conversion_factors.inverted: + line = line.replace(".real", ".lj") + else: + line = line.replace(".lj", ".real") + lines[i] = line + lines_changed += 1 + + if "stk" in line and "xstk" not in line and "coaxstk" not in line: + elements = line.split() + elements[6] = str( # convert xi + round( + float(elements[6]) * conversion_factors.energy_conv_factor, + 14, + ) + ) + elements[7] = str( # convert kappa + round(float(elements[7]) * conversion_factors.kT_conv_factor, 9) + ) + lines[i] = " ".join(elements) + "\n" + + else: + elements = line.split() + + if "bond_coeff" in line: + if oxdna2_flag: + elements[2:] = conversion_factors.oxdna2_fene_string.split() + elif oxrna2_flag: + elements[2:] = conversion_factors.oxrna2_fene_string.split() + else: + elements[2:] = conversion_factors.oxdna_fene_string.split() + + elif "excv" in line: + if oxdna2_flag: + elements[4:] = conversion_factors.oxdna2_excv_string.split() + elif oxrna2_flag: + elements[4:] = conversion_factors.oxrna2_excv_string.split() + else: + elements[4:] = conversion_factors.oxdna_excv_string.split() + + elif "stk" in line: + + if "coaxstk" in line: + if oxdna2_flag: + elements[4:] = ( + conversion_factors.oxdna2_coaxstk_string.split() + ) + elif oxrna2_flag: + elements[4:] = ( + conversion_factors.oxrna2_coaxstk_string.split() + ) + else: + elements[4:] = ( + conversion_factors.oxdna_coaxstk_string.split() + ) + + elif "xstk" in line: + if oxdna2_flag: + elements[4:] = conversion_factors.oxdna2_xstk_string.split() + elif oxrna2_flag: + elements[4:] = conversion_factors.oxrna2_xstk_string.split() + else: + elements[4:] = conversion_factors.oxdna_xstk_string.split() + + else: # stk + elements[6] = str( # convert xi + round( + float(elements[6]) + * conversion_factors.energy_conv_factor, + 14, + ) + ) + elements[7] = str( # convert kappa + round( + float(elements[7]) * conversion_factors.kT_conv_factor, + 9, + ) + ) + if oxdna2_flag: + elements[8:] = conversion_factors.oxdna2_stk_string.split() + elif oxrna2_flag: + elements[8:] = conversion_factors.oxrna2_stk_string.split() + else: + elements[8:] = conversion_factors.oxdna_stk_string.split() + + elif "hbond" in line: + if elements[1] == "*" and elements[2] == "*": + if oxdna2_flag: + elements[5:] = ( + conversion_factors.oxdna2_hbond_string.split() + ) + elif oxrna2_flag: + elements[5:] = ( + conversion_factors.oxrna2_hbond_string.split() + ) + else: + elements[5:] = conversion_factors.oxdna_hbond_string.split() + else: + if oxdna2_flag: + elements[5:] = ( + conversion_factors.oxdna2_hbond_1_4_2_3_string.split() + ) + elif oxrna2_flag: + elements[5:] = ( + conversion_factors.oxrna2_hbond_1_4_2_3_3_4_string.split() + ) + else: + elements[5:] = ( + conversion_factors.oxdna_hbond_1_4_2_3_string.split() + ) + lines[i] = " ".join(elements) + "\n" + lines_changed += 1 + + elif "langevin" in line: + elements = line.split() + elements[6] = str( + round(float(elements[6]) * conversion_factors.time_conv_factor, 2) + ) + lines[i] = " ".join(elements) + "\n" + lines_changed += 1 + + elif "timestep" in line: + elements = line.split() + elements[1] = str( + round(float(elements[1]) * conversion_factors.time_conv_factor, 5) + ) + lines[i] = " ".join(elements) + "\n" + lines_changed += 1 + + elif "comm_modify" in line: + elements = line.split() + elements[2] = str( + round(float(elements[2]) * conversion_factors.length_conv_factor, 1) + ) + lines[i] = " ".join(elements) + "\n" + lines_changed += 1 + + else: + continue + + if conversion_factors.inverted: + new_inputfile_path = inputfile_path + "_lj" + else: + new_inputfile_path = inputfile_path + "_real" + + with open(new_inputfile_path, "w", encoding="UTF-8") as file: + if conversion_factors.inverted: + file.write( + "# LAMMPS input file in LJ units via oxdna lj2real.py, date " + + str(datetime.date.today()) + + "\n" + ) + else: + file.write( + "# LAMMPS input file in real units via oxdna lj2real.py, date " + + str(datetime.date.today()) + + "\n" + ) + file.writelines(lines) + if lines_changed == 0: + print( + "Warning: No lines changed in input file. Ensure correct usage: python lj2real.py [-i]" + ) + else: + print(f"Input file lines changed: {lines_changed}") + + return new_inputfile_path + + +def main(): + """Main function""" + + print( + "\nLAMMPS data and input file conversion to real units via oxdna convert_data.py\n" + "Note: This script assumes a input and data file structure similar to those found in examples/PACKAGES/cgdna/examples/.\n" + "Ensure output is checked for correctness." + ) + + if len(sys.argv) > 2: + datafile_path = sys.argv[1] + inputfile_path = sys.argv[2] + invert = False + if len(sys.argv) > 3: + if sys.argv[3] == "-i": + invert = True + print("Performing real -> LJ conversion.") + else: + print( + "Invalid flag. Usage: python lj2real.py [-i]" + ) + sys.exit(1) + + if invert: + conversion_factors = ConversionFactors(invert=True) + print( + "\nUsing conversion factors (real T, m, l, t, v, ρ -> LJ T*, m*, l*, t*, v*, ρ*):" + ) + + else: + conversion_factors = ConversionFactors(invert=False) + print( + "\nUsing conversion factors (LJ T*, m*, l*, t*, v*, ρ* -> real T, m, l, t, v, ρ):" + ) + + else: + print("\nUsage: python lj2real.py [-i]") + print("\t[-i] flag is optional and is used to convert real -> LJ units.\n") + sys.exit(1) + + conversion_factors_string = ( + f"Temperature T\t≈ {round(conversion_factors.temp_conv_factor, 6)} T* (K)\n" + f"Energy ε\t≈ {round(conversion_factors.energy_conv_factor, 6)} ε* (kcal/mol)\n" + f"Mass m\t\t≈ {round(conversion_factors.mass_conv_factor, 6)} m* (g/mol)\n" + f"Length l\t≈ {round(conversion_factors.length_conv_factor, 6)} l* (Å)\n" + f"Time t\t\t≈ {round(conversion_factors.time_conv_factor, 6)} t* (fs)\n" + f"Velocity v\t≈ {round(conversion_factors.vel_conv_factor, 6)} v* (Å/fs)\n" + f"AngMom l\t≈ {round(conversion_factors.angular_mom_conv_factor, 6)} (g/mol Å^2/fs)\n" + f"Density ρ\t≈ {round(conversion_factors.density_conv_factor, 6)} ρ* (g/cm^3)\n" + f"Calculated using Sengar, Ouldridge, Henrich, Rovigatti, & Šulc. Front Mol Biosci 8 (2021). & https://docs.lammps.org/units.html.\n" + f"See examples/PACKAGES/cgdna/util/lj2real.py for exact conversion factors.\n" + ) + print(conversion_factors_string) + + print("Current directory: ", os.getcwd()) + + try: + new_datafile_path = modify_datafile(datafile_path, conversion_factors) + print(f"New data file: {new_datafile_path}") + except Exception as e: + print(f"Error modifying data file: {e}") + + try: + new_inputfile_path = modify_inputfile(inputfile_path, conversion_factors) + print(f"New input file: {new_inputfile_path}\n") + except Exception as e: + print(f"Error modifying input file: {e}") + + +if __name__ == "__main__": + main() diff --git a/potentials/oxdna.lj b/potentials/oxdna.lj new file mode 100644 index 0000000000..e737d68e15 --- /dev/null +++ b/potentials/oxdna.lj @@ -0,0 +1,10 @@ +# DATE: 2024-04-21 UNITS: lj CONTRIBUTOR: Oliver Henrich, oliver.henrich@strath.ac.uk CITATION: Ouldridge, Louis, and Doye, J Chem Phys, 134, 8 (2011) +# +* fene 2.0 0.25 0.7525 +* * excv 2.0 0.7 0.675 2.0 0.515 0.5 2.0 0.33 0.32 +* * stk 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 +* * hbond 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 +1 4 hbond 1.077 8.0 0.4 0.75 0.34 0.7 1.5 0 0.7 1.5 0 0.7 1.5 0 0.7 0.46 3.141592653589793 0.7 4.0 1.5707963267948966 0.45 4.0 1.5707963267948966 0.45 +2 3 hbond 1.077 8.0 0.4 0.75 0.34 0.7 1.5 0 0.7 1.5 0 0.7 1.5 0 0.7 0.46 3.141592653589793 0.7 4.0 1.5707963267948966 0.45 4.0 1.5707963267948966 0.45 +* * 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 +* * coaxstk 46.0 0.4 0.6 0.22 0.58 2.0 2.541592653589793 0.65 1.3 0 0.8 0.9 0 0.95 0.9 0 0.95 2.0 -0.65 2.0 -0.65 diff --git a/potentials/oxdna.real b/potentials/oxdna.real new file mode 100644 index 0000000000..d687cbe1f4 --- /dev/null +++ b/potentials/oxdna.real @@ -0,0 +1,10 @@ +# DATE: 2024-04-26 UNITS: real CONTRIBUTOR: Oliver Henrich, oliver.henrich@strath.ac.uk +# +* fene 11.92337812042065 2.1295 6.409795 +* * excv 11.92337812042065 5.9626 5.74965 11.92337812042065 4.38677 4.259 11.92337812042065 2.81094 2.72576 +* * stk 0.70439070204273 3.4072 7.6662 2.72576 6.3885 1.3 0.0 0.8 0.9 0.0 0.95 0.9 0.0 0.95 2.0 0.65 2.0 0.65 +* * hbond 0.0 0.93918760272364 3.4072 6.3885 2.89612 5.9626 1.5 0.0 0.7 1.5 0.0 0.7 1.5 0.0 0.7 0.46 3.141592654 0.7 4.0 1.570796327 0.45 4.0 1.570796327 0.45 +1 4 hbond 6.42073911784652 0.93918760272364 3.4072 6.3885 2.89612 5.9626 1.5 0.0 0.7 1.5 0.0 0.7 1.5 0.0 0.7 0.46 3.141592654 0.7 4.0 1.570796327 0.45 4.0 1.570796327 0.45 +2 3 hbond 6.42073911784652 0.93918760272364 3.4072 6.3885 2.89612 5.9626 1.5 0.0 0.7 1.5 0.0 0.7 1.5 0.0 0.7 0.46 3.141592654 0.7 4.0 1.570796327 0.45 4.0 1.570796327 0.45 +* * xstk 3.902852174 4.89785 5.74965 4.21641 5.57929 2.25 0.791592654 0.58 1.7 1.0 0.68 1.7 1.0 0.68 1.5 0.0 0.65 1.7 0.875 0.68 1.7 0.875 0.68 +* * coaxstk 3.77965257404268 3.4072 5.1108 1.87396 4.94044 2.0 2.541592654 0.65 1.3 0.0 0.8 0.9 0.0 0.95 0.9 0.0 0.95 2.0 -0.65 2.0 -0.65 diff --git a/potentials/oxdna2.lj b/potentials/oxdna2.lj new file mode 100644 index 0000000000..15c82f52c6 --- /dev/null +++ b/potentials/oxdna2.lj @@ -0,0 +1,11 @@ +# DATE: 2024-04-21 UNITS: lj CONTRIBUTOR: Oliver Henrich, oliver.henrich@strath.ac.uk CITATION: Snodin, Randisi, Mosayebi, Šulc et. al., J Chem Phys, 142, 23 (2015) +# +* fene 2.0 0.25 0.7564 +* * excv 2.0 0.7 0.675 2.0 0.515 0.5 2.0 0.33 0.32 +* * stk 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 +* * hbond 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 +1 4 hbond 1.0678 8.0 0.4 0.75 0.34 0.7 1.5 0 0.7 1.5 0 0.7 1.5 0 0.7 0.46 3.141592653589793 0.7 4.0 1.5707963267948966 0.45 4.0 1.5707963267948966 0.45 +2 3 hbond 1.0678 8.0 0.4 0.75 0.34 0.7 1.5 0 0.7 1.5 0 0.7 1.5 0 0.7 0.46 3.141592653589793 0.7 4.0 1.5707963267948966 0.45 4.0 1.5707963267948966 0.45 +* * 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 +* * 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 +* * dh 0.815 diff --git a/potentials/oxdna2.real b/potentials/oxdna2.real new file mode 100644 index 0000000000..fa42412f2b --- /dev/null +++ b/potentials/oxdna2.real @@ -0,0 +1,11 @@ +# DATE: 2024-04-26 UNITS: real CONTRIBUTOR: Oliver Henrich, oliver.henrich@strath.ac.uk +# +* fene 11.92337812042065 2.1295 6.4430152 +* * excv 11.92337812042065 5.9626 5.74965 11.92337812042065 4.38677 4.259 11.92337812042065 2.81094 2.72576 +* * stk 0.70439070204273 3.4072 7.6662 2.72576 6.3885 1.3 0.0 0.8 0.9 0.0 0.95 0.9 0.0 0.95 2.0 0.65 2.0 0.65 +* * hbond 0.0 0.93918760272364 3.4072 6.3885 2.89612 5.9626 1.5 0.0 0.7 1.5 0.0 0.7 1.5 0.0 0.7 0.46 3.141592654 0.7 4.0 1.570796327 0.45 4.0 1.570796327 0.45 +1 4 hbond 6.36589157849259 0.93918760272364 3.4072 6.3885 2.89612 5.9626 1.5 0.0 0.7 1.5 0.0 0.7 1.5 0.0 0.7 0.46 3.141592654 0.7 4.0 1.570796327 0.45 4.0 1.570796327 0.45 +2 3 hbond 6.36589157849259 0.93918760272364 3.4072 6.3885 2.89612 5.9626 1.5 0.0 0.7 1.5 0.0 0.7 1.5 0.0 0.7 0.46 3.141592654 0.7 4.0 1.570796327 0.45 4.0 1.570796327 0.45 +* * xstk 3.9029021145006 4.89785 5.74965 4.21641 5.57929 2.25 0.791592654 0.58 1.7 1.0 0.68 1.7 1.0 0.68 1.5 0.0 0.65 1.7 0.875 0.68 1.7 0.875 0.68 +* * coaxstk 4.80673207785863 3.4072 5.1108 1.87396 4.94044 2.0 2.891592653589793 0.65 1.3 0.0 0.8 0.9 0.0 0.95 0.9 0.0 0.95 40.0 3.116592653589793 +* * dh 0.815 diff --git a/potentials/oxrna2.lj b/potentials/oxrna2.lj new file mode 100644 index 0000000000..f4e7c59dbb --- /dev/null +++ b/potentials/oxrna2.lj @@ -0,0 +1,12 @@ +# DATE: 2024-04-19 UNITS: lj CONTRIBUTOR: Oliver Henrich, oliver.henrich@strath.ac.uk CITATION: Šulc, Romano, Ouldridge, Rovigatti, Doye, Louis, J Chem Phys, 140, 23 (2014) +# +* fene 2.0 0.25 0.761070781051 +* * excv 2.0 0.7 0.675 2.0 0.515 0.5 2.0 0.33 0.32 +* * stk 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 +* * hbond 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 +1 4 hbond 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 +2 3 hbond 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 +3 4 hbond 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 +* * 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 +* * 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 +* * dh 1.02455 diff --git a/potentials/oxrna2.real b/potentials/oxrna2.real new file mode 100644 index 0000000000..10a8d8a8b5 --- /dev/null +++ b/potentials/oxrna2.real @@ -0,0 +1,12 @@ +# DATE: 2024-04-26 UNITS: real CONTRIBUTOR: Oliver Henrich, oliver.henrich@strath.ac.uk +# +* fene 11.92337812042065 2.1295 6.482800913 +* * excv 11.92337812042065 5.9626 5.74965 11.92337812042065 4.38677 4.259 11.92337812042065 2.81094 2.72576 +* * stk 0.70439070204273 3.66274 7.92174 2.9813 6.64404 0.9 0.0 0.95 0.9 0.0 0.95 1.3 0.0 0.8 1.3 0.0 0.8 2.0 0.65 2.0 0.65 +* * hbond 0.0 0.93918760272364 3.4072 6.3885 2.89612 5.9626 1.5 0.0 0.7 1.5 0.0 0.7 1.5 0.0 0.7 0.46 3.141592653589793 0.7 4.0 1.5707963267948966 0.45 4.0 1.5707963267948966 0.45 +1 4 hbond 5.18928666388042 0.93918760272364 3.4072 6.3885 2.89612 5.9626 1.5 0.0 0.7 1.5 0.0 0.7 1.5 0.0 0.7 0.46 3.141592653589793 0.7 4.0 1.5707963267948966 0.45 4.0 1.5707963267948966 0.45 +2 3 hbond 5.18928666388042 0.93918760272364 3.4072 6.3885 2.89612 5.9626 1.5 0.0 0.7 1.5 0.0 0.7 1.5 0.0 0.7 0.46 3.141592653589793 0.7 4.0 1.5707963267948966 0.45 4.0 1.5707963267948966 0.45 +3 4 hbond 5.18928666388042 0.93918760272364 3.4072 6.3885 2.89612 5.9626 1.5 0.0 0.7 1.5 0.0 0.7 1.5 0.0 0.7 0.46 3.141592653589793 0.7 4.0 1.5707963267948966 0.45 4.0 1.5707963267948966 0.45 +* * xstk 4.92690859644113 4.259 5.1108 3.57756 4.94044 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 +* * coaxstk 6.57330882442206 4.259 5.1108 3.57756 4.94044 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 +* * dh 1.02455 diff --git a/src/CG-DNA/atom_vec_oxdna.cpp b/src/CG-DNA/atom_vec_oxdna.cpp index d7aa7a3d01..38f78f94bf 100644 --- a/src/CG-DNA/atom_vec_oxdna.cpp +++ b/src/CG-DNA/atom_vec_oxdna.cpp @@ -12,6 +12,7 @@ ------------------------------------------------------------------------- */ #include "atom_vec_oxdna.h" +#include "constants_oxdna.h" #include "atom.h" #include "error.h" @@ -45,6 +46,9 @@ AtomVecOxdna::AtomVecOxdna(LAMMPS *lmp) : AtomVec(lmp) if (!force->newton_bond) error->warning(FLERR, "Write_data command requires newton on to preserve 3'->5' bond polarity"); + + // initialize oxDNA units + ConstantsOxdna constants(lmp); } /* ---------------------------------------------------------------------- diff --git a/src/CG-DNA/bond_oxdna2_fene.cpp b/src/CG-DNA/bond_oxdna2_fene.cpp index ace0742963..73ad827d72 100644 --- a/src/CG-DNA/bond_oxdna2_fene.cpp +++ b/src/CG-DNA/bond_oxdna2_fene.cpp @@ -15,6 +15,7 @@ ------------------------------------------------------------------------- */ #include "bond_oxdna2_fene.h" +#include "constants_oxdna.h" using namespace LAMMPS_NS; @@ -24,8 +25,8 @@ using namespace LAMMPS_NS; void BondOxdna2Fene::compute_interaction_sites(double e1[3], double e2[3], double /*e3*/[3], double r[3]) const { - constexpr double d_cs_x = -0.34; - constexpr double d_cs_y = +0.3408; + double d_cs_x = ConstantsOxdna::get_d_cs_x(); + double d_cs_y = ConstantsOxdna::get_d_cs_y(); r[0] = d_cs_x * e1[0] + d_cs_y * e2[0]; r[1] = d_cs_x * e1[1] + d_cs_y * e2[1]; diff --git a/src/CG-DNA/bond_oxdna_fene.cpp b/src/CG-DNA/bond_oxdna_fene.cpp index 780b71e44c..39444a738a 100644 --- a/src/CG-DNA/bond_oxdna_fene.cpp +++ b/src/CG-DNA/bond_oxdna_fene.cpp @@ -18,10 +18,12 @@ #include "atom.h" #include "comm.h" +#include "constants_oxdna.h" #include "error.h" #include "force.h" #include "memory.h" #include "neighbor.h" +#include "potential_file_reader.h" #include "update.h" #include "math_extra.h" @@ -49,7 +51,7 @@ BondOxdnaFene::~BondOxdnaFene() void BondOxdnaFene::compute_interaction_sites(double e1[3], double /*e2*/[3], double /*e3*/[3], double r[3]) const { - constexpr double d_cs = -0.4; + double d_cs = ConstantsOxdna::get_d_cs(); r[0] = d_cs * e1[0]; r[1] = d_cs * e1[1]; @@ -295,15 +297,49 @@ void BondOxdnaFene::allocate() void BondOxdnaFene::coeff(int narg, char **arg) { - if (narg != 4) error->all(FLERR, "Incorrect args for bond coefficients in oxdna/fene"); + if (narg != 2 && narg != 4) error->all(FLERR, "Incorrect args for bond coefficients in oxdna/fene"); if (!allocated) allocate(); int ilo, ihi; utils::bounds(FLERR, arg[0], 1, atom->nbondtypes, ilo, ihi, error); - double k_one = utils::numeric(FLERR, arg[1], false, lmp); - double Delta_one = utils::numeric(FLERR, arg[2], false, lmp); - double r0_one = utils::numeric(FLERR, arg[3], false, lmp); + double k_one; + double Delta_one; + double r0_one; + + if (narg == 4) { + k_one = utils::numeric(FLERR, arg[1], false, lmp); + Delta_one = utils::numeric(FLERR, arg[2], false, lmp); + r0_one = utils::numeric(FLERR, arg[3], false, lmp); + } else { + if (comm->me == 0) { // read values from potential file + PotentialFileReader reader(lmp, arg[1], "oxdna potential", " (fene)"); + char * line; + std::string iloc, potential_name; + + while(line = reader.next_line()) { + try { + ValueTokenizer values(line); + iloc = values.next_string(); + potential_name = values.next_string(); + if (iloc == arg[0] && potential_name == "fene") { + k_one = values.next_double(); + Delta_one = values.next_double(); + r0_one = values.next_double(); + + break; + } else continue; + } catch (std::exception &e) { + error->one(FLERR, "Problem parsing oxDNA potential file: {}", e.what()); + } + } + if (iloc != arg[0] || potential_name != "fene") error->one(FLERR, "No corresponding fene potential found in file {} for bond type {}", arg[1], arg[0]); + } + + MPI_Bcast(&k_one, 1, MPI_DOUBLE, 0, world); + MPI_Bcast(&Delta_one, 1, MPI_DOUBLE, 0, world); + MPI_Bcast(&r0_one, 1, MPI_DOUBLE, 0, world); + } int count = 0; diff --git a/src/CG-DNA/bond_oxrna2_fene.cpp b/src/CG-DNA/bond_oxrna2_fene.cpp index 5d28f744a8..4faa95ef6f 100644 --- a/src/CG-DNA/bond_oxrna2_fene.cpp +++ b/src/CG-DNA/bond_oxrna2_fene.cpp @@ -16,6 +16,8 @@ #include "bond_oxrna2_fene.h" +#include "constants_oxdna.h" + using namespace LAMMPS_NS; /* ---------------------------------------------------------------------- @@ -25,8 +27,8 @@ using namespace LAMMPS_NS; void BondOxrna2Fene::compute_interaction_sites(double e1[3], double /*e2*/[3], double e3[3], double r[3]) const { - constexpr double d_cs_x = -0.4; - constexpr double d_cs_z = +0.2; + double d_cs_x = ConstantsOxdna::get_d_cs(); + double d_cs_z = ConstantsOxdna::get_d_cs_z(); r[0] = d_cs_x * e1[0] + d_cs_z * e3[0]; r[1] = d_cs_x * e1[1] + d_cs_z * e3[1]; diff --git a/src/CG-DNA/constants_oxdna.cpp b/src/CG-DNA/constants_oxdna.cpp new file mode 100644 index 0000000000..f3623f4dab --- /dev/null +++ b/src/CG-DNA/constants_oxdna.cpp @@ -0,0 +1,69 @@ +/* -*- c++ -*- ---------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + https://www.lammps.org/, Sandia National Laboratories + LAMMPS development team: developers@lammps.org + + 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 authors: Oliver Henrich (University of Strathclyde, Glasgow) + Kierran Falloon (University of Strathclyde, Glasgow) +------------------------------------------------------------------------- */ + +#include "constants_oxdna.h" + +namespace LAMMPS_NS { + +ConstantsOxdna::ConstantsOxdna(class LAMMPS *lmp) : Pointers(lmp) +{ + // set oxDNA units + units = update->unit_style; + real_flag = utils::strmatch(units.c_str(), "^real"); + if (real_flag) set_real_units(); +} + +// default to lj units +// oxDNA 1 parameters +double ConstantsOxdna::d_cs = -0.4; +double ConstantsOxdna::d_cst = +0.34; +double ConstantsOxdna::d_chb = +0.4; +double ConstantsOxdna::d_cb = +0.4; +// oxDNA 2 parameters +double ConstantsOxdna::d_cs_x = -0.34; +double ConstantsOxdna::d_cs_y = +0.3408; +double ConstantsOxdna::lambda_dh_one_prefactor = +0.3616455075438555; // = C1 +double ConstantsOxdna::qeff_dh_pf_one_prefactor = +0.08173808693529228; // = C2 +// oxRNA 2 parameters +double ConstantsOxdna::d_cs_z = +0.2; +double ConstantsOxdna::d_cst_x_3p = +0.4; +double ConstantsOxdna::d_cst_y_3p = +0.1; +double ConstantsOxdna::d_cst_x_5p = +0.124906078525; +double ConstantsOxdna::d_cst_y_5p = -0.00866274917473; + +void ConstantsOxdna::set_real_units() +{ + // oxDNA 1 parameters in real units + d_cs = -3.4072; + d_cst = +2.89612; + d_chb = +3.4072; + d_cb = +3.4072; + // oxDNA 2 parameters in real units + d_cs_x = -2.89612; + d_cs_y = +2.9029344; + lambda_dh_one_prefactor = +0.05624154892; // = C1 * 8.518 * sqrt(k_B/4.142e-20) + qeff_dh_pf_one_prefactor = +4.15079634587587; // = C2 * 5.961689060210325 * 8.518 + // oxRNA 2 parameters in real units + // d_cs_x = -3.4072 = d_cs for RNA + d_cs_z = +1.7036; + d_cst_x_3p = +3.4072, + d_cst_y_3p = +0.8518; + d_cst_x_5p = +1.063949977, + d_cst_y_5p = -0.07378929747; +}; + +} // namespace LAMMPS_NS diff --git a/src/CG-DNA/constants_oxdna.h b/src/CG-DNA/constants_oxdna.h new file mode 100644 index 0000000000..69defbc0ea --- /dev/null +++ b/src/CG-DNA/constants_oxdna.h @@ -0,0 +1,65 @@ +/* -*- c++ -*- ---------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + https://www.lammps.org/, Sandia National Laboratories + LAMMPS development team: developers@lammps.org + + 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. +------------------------------------------------------------------------- */ + +#ifndef CONSTANTS_OXDNA_H +#define CONSTANTS_OXDNA_H + +#include "update.h" + +namespace LAMMPS_NS { + +class ConstantsOxdna : protected Pointers { + public: + ConstantsOxdna(class LAMMPS *lmp); + virtual ~ConstantsOxdna(){}; + + // oxDNA 1 getters + static double get_d_cs() { return d_cs; } + static double get_d_cst() { return d_cst; } + static double get_d_chb() { return d_chb; } + static double get_d_cb() { return d_cb; } + + // oxDNA 2 getters + static double get_d_cs_x() { return d_cs_x; } + static double get_d_cs_y() { return d_cs_y; } + static double get_lambda_dh_one_prefactor() { return lambda_dh_one_prefactor; } + static double get_qeff_dh_pf_one_prefactor() { return qeff_dh_pf_one_prefactor; } + + // oxRNA 2 getters + static double get_d_cs_z() { return d_cs_z; } + static double get_d_cst_x_3p() { return d_cst_x_3p; } + static double get_d_cst_y_3p() { return d_cst_y_3p; } + static double get_d_cst_x_5p() { return d_cst_x_5p; } + static double get_d_cst_y_5p() { return d_cst_y_5p; } + + private: + std::string units; + bool real_flag; + void set_real_units(); + + // oxDNA 1 parameters + static double d_cs, d_cst, d_chb, d_cb; + + // oxDNA 2 parameters + static double d_cs_x, d_cs_y; + static double lambda_dh_one_prefactor, qeff_dh_pf_one_prefactor; + + // oxRNA 2 parameters + static double d_cs_z; + static double d_cst_x_3p, d_cst_y_3p; + static double d_cst_x_5p, d_cst_y_5p; +}; + +} // namespace LAMMPS_NS + +#endif diff --git a/src/CG-DNA/pair_oxdna2_coaxstk.cpp b/src/CG-DNA/pair_oxdna2_coaxstk.cpp index b2666c2f0f..bf8c8e545a 100644 --- a/src/CG-DNA/pair_oxdna2_coaxstk.cpp +++ b/src/CG-DNA/pair_oxdna2_coaxstk.cpp @@ -19,6 +19,7 @@ #include "atom.h" #include "comm.h" +#include "constants_oxdna.h" #include "error.h" #include "force.h" #include "math_const.h" @@ -26,6 +27,7 @@ #include "memory.h" #include "mf_oxdna.h" #include "neigh_list.h" +#include "potential_file_reader.h" #include #include @@ -113,7 +115,8 @@ void PairOxdna2Coaxstk::compute(int eflag, int vflag) double cosphi3; // distances COM-backbone site, COM-stacking site - double d_cs=-0.4, d_cst=+0.34; + double d_cs = ConstantsOxdna::get_d_cs(); + double d_cst = ConstantsOxdna::get_d_cst(); // vectors COM-backbone site, COM-stacking site in lab frame double ra_cs[3],ra_cst[3]; double rb_cs[3],rb_cst[3]; @@ -557,7 +560,7 @@ void PairOxdna2Coaxstk::coeff(int narg, char **arg) { int count; - if (narg != 21) error->all(FLERR,"Incorrect args for pair coefficients in oxdna2/coaxstk"); + if (narg != 3 && narg != 21) error->all(FLERR,"Incorrect args for pair coefficients in oxdna2/coaxstk"); if (!allocated) allocate(); int ilo,ihi,jlo,jhi; @@ -584,30 +587,103 @@ void PairOxdna2Coaxstk::coeff(int narg, char **arg) double AA_cxst1_one, BB_cxst1_one; - k_cxst_one = utils::numeric(FLERR,arg[2],false,lmp); - cut_cxst_0_one = utils::numeric(FLERR,arg[3],false,lmp); - cut_cxst_c_one = utils::numeric(FLERR,arg[4],false,lmp); - cut_cxst_lo_one = utils::numeric(FLERR,arg[5],false,lmp); - cut_cxst_hi_one = utils::numeric(FLERR,arg[6],false,lmp); + if (narg == 21) { + k_cxst_one = utils::numeric(FLERR,arg[2],false,lmp); + cut_cxst_0_one = utils::numeric(FLERR,arg[3],false,lmp); + cut_cxst_c_one = utils::numeric(FLERR,arg[4],false,lmp); + cut_cxst_lo_one = utils::numeric(FLERR,arg[5],false,lmp); + cut_cxst_hi_one = utils::numeric(FLERR,arg[6],false,lmp); - a_cxst1_one = utils::numeric(FLERR,arg[7],false,lmp); - theta_cxst1_0_one = utils::numeric(FLERR,arg[8],false,lmp); - dtheta_cxst1_ast_one = utils::numeric(FLERR,arg[9],false,lmp); + a_cxst1_one = utils::numeric(FLERR,arg[7],false,lmp); + theta_cxst1_0_one = utils::numeric(FLERR,arg[8],false,lmp); + dtheta_cxst1_ast_one = utils::numeric(FLERR,arg[9],false,lmp); - a_cxst4_one = utils::numeric(FLERR,arg[10],false,lmp); - theta_cxst4_0_one = utils::numeric(FLERR,arg[11],false,lmp); - dtheta_cxst4_ast_one = utils::numeric(FLERR,arg[12],false,lmp); + a_cxst4_one = utils::numeric(FLERR,arg[10],false,lmp); + theta_cxst4_0_one = utils::numeric(FLERR,arg[11],false,lmp); + dtheta_cxst4_ast_one = utils::numeric(FLERR,arg[12],false,lmp); - a_cxst5_one = utils::numeric(FLERR,arg[13],false,lmp); - theta_cxst5_0_one = utils::numeric(FLERR,arg[14],false,lmp); - dtheta_cxst5_ast_one = utils::numeric(FLERR,arg[15],false,lmp); + a_cxst5_one = utils::numeric(FLERR,arg[13],false,lmp); + theta_cxst5_0_one = utils::numeric(FLERR,arg[14],false,lmp); + dtheta_cxst5_ast_one = utils::numeric(FLERR,arg[15],false,lmp); - a_cxst6_one = utils::numeric(FLERR,arg[16],false,lmp); - theta_cxst6_0_one = utils::numeric(FLERR,arg[17],false,lmp); - dtheta_cxst6_ast_one = utils::numeric(FLERR,arg[18],false,lmp); + a_cxst6_one = utils::numeric(FLERR,arg[16],false,lmp); + theta_cxst6_0_one = utils::numeric(FLERR,arg[17],false,lmp); + dtheta_cxst6_ast_one = utils::numeric(FLERR,arg[18],false,lmp); - AA_cxst1_one = utils::numeric(FLERR,arg[19],false,lmp); - BB_cxst1_one = utils::numeric(FLERR,arg[20],false,lmp); + AA_cxst1_one = utils::numeric(FLERR,arg[19],false,lmp); + BB_cxst1_one = utils::numeric(FLERR,arg[20],false,lmp); + } else { + if (comm->me == 0) { // read values from potential file + PotentialFileReader reader(lmp, arg[2], "oxdna potential", " (coaxstk)"); + char * line; + std::string iloc, jloc, potential_name; + + while(line = reader.next_line()) { + try { + ValueTokenizer values(line); + iloc = values.next_string(); + jloc = values.next_string(); + potential_name = values.next_string(); + if (iloc == arg[0] && jloc == arg[1] && potential_name == "coaxstk") { + k_cxst_one = values.next_double(); + cut_cxst_0_one = values.next_double(); + cut_cxst_c_one = values.next_double(); + cut_cxst_lo_one = values.next_double(); + cut_cxst_hi_one = values.next_double(); + + a_cxst1_one = values.next_double(); + theta_cxst1_0_one = values.next_double(); + dtheta_cxst1_ast_one = values.next_double(); + + a_cxst4_one = values.next_double(); + theta_cxst4_0_one = values.next_double(); + dtheta_cxst4_ast_one = values.next_double(); + + a_cxst5_one = values.next_double(); + theta_cxst5_0_one = values.next_double(); + dtheta_cxst5_ast_one = values.next_double(); + + a_cxst6_one = values.next_double(); + theta_cxst6_0_one = values.next_double(); + dtheta_cxst6_ast_one = values.next_double(); + + AA_cxst1_one = values.next_double(); + BB_cxst1_one = values.next_double(); + + break; + } else continue; + } catch (std::exception &e) { + error->one(FLERR, "Problem parsing oxDNA2 potential file: {}", e.what()); + } + } + if (iloc != arg[0] || jloc != arg[1] || potential_name != "coaxstk") error->one(FLERR, "No corresponding coaxstk potential found in file {} for pair type {} {}", arg[2], arg[0], arg[1]); + } + + MPI_Bcast(&k_cxst_one, 1, MPI_DOUBLE, 0, world); + MPI_Bcast(&cut_cxst_0_one, 1, MPI_DOUBLE, 0, world); + MPI_Bcast(&cut_cxst_c_one, 1, MPI_DOUBLE, 0, world); + MPI_Bcast(&cut_cxst_lo_one, 1, MPI_DOUBLE, 0, world); + MPI_Bcast(&cut_cxst_hi_one, 1, MPI_DOUBLE, 0, world); + + MPI_Bcast(&a_cxst1_one, 1, MPI_DOUBLE, 0, world); + MPI_Bcast(&theta_cxst1_0_one, 1, MPI_DOUBLE, 0, world); + MPI_Bcast(&dtheta_cxst1_ast_one, 1, MPI_DOUBLE, 0, world); + + MPI_Bcast(&a_cxst4_one, 1, MPI_DOUBLE, 0, world); + MPI_Bcast(&theta_cxst4_0_one, 1, MPI_DOUBLE, 0, world); + MPI_Bcast(&dtheta_cxst4_ast_one, 1, MPI_DOUBLE, 0, world); + + MPI_Bcast(&a_cxst5_one, 1, MPI_DOUBLE, 0, world); + MPI_Bcast(&theta_cxst5_0_one, 1, MPI_DOUBLE, 0, world); + MPI_Bcast(&dtheta_cxst5_ast_one, 1, MPI_DOUBLE, 0, world); + + MPI_Bcast(&a_cxst6_one, 1, MPI_DOUBLE, 0, world); + MPI_Bcast(&theta_cxst6_0_one, 1, MPI_DOUBLE, 0, world); + MPI_Bcast(&dtheta_cxst6_ast_one, 1, MPI_DOUBLE, 0, world); + + MPI_Bcast(&AA_cxst1_one, 1, MPI_DOUBLE, 0, world); + MPI_Bcast(&BB_cxst1_one, 1, MPI_DOUBLE, 0, world); + } b_cxst_lo_one = 0.25 * (cut_cxst_lo_one - cut_cxst_0_one) * (cut_cxst_lo_one - cut_cxst_0_one)/ (0.5 * (cut_cxst_lo_one - cut_cxst_0_one) * (cut_cxst_lo_one - cut_cxst_0_one) - diff --git a/src/CG-DNA/pair_oxdna2_dh.cpp b/src/CG-DNA/pair_oxdna2_dh.cpp index d60342e5e2..7f38c4b96f 100644 --- a/src/CG-DNA/pair_oxdna2_dh.cpp +++ b/src/CG-DNA/pair_oxdna2_dh.cpp @@ -19,11 +19,13 @@ #include "atom.h" #include "comm.h" +#include "constants_oxdna.h" #include "error.h" #include "force.h" #include "math_extra.h" #include "memory.h" #include "neigh_list.h" +#include "potential_file_reader.h" #include #include @@ -66,7 +68,8 @@ PairOxdna2Dh::~PairOxdna2Dh() void PairOxdna2Dh::compute_interaction_sites(double e1[3], double e2[3], double /*e3*/[3], double r[3]) { - double d_cs_x=-0.34, d_cs_y=+0.3408; + double d_cs_x = ConstantsOxdna::get_d_cs_x(); + double d_cs_y = ConstantsOxdna::get_d_cs_y(); r[0] = d_cs_x*e1[0] + d_cs_y*e2[0]; r[1] = d_cs_x*e1[1] + d_cs_y*e2[1]; @@ -303,7 +306,30 @@ void PairOxdna2Dh::coeff(int narg, char **arg) T = utils::numeric(FLERR,arg[2],false,lmp); rhos_dh_one = utils::numeric(FLERR,arg[3],false,lmp); - qeff_dh_one = utils::numeric(FLERR,arg[4],false,lmp); + + if (utils::strmatch(arg[4], "^[a-zA-Z0-9]*\\.[a-zA-Z]+$") == true) { // if last arg is a potential file + if (comm->me == 0) { // read value from potential file + PotentialFileReader reader(lmp, arg[4], "oxdna potential", " (dh)"); + char * line; + std::string iloc, jloc, potential_name; + while(line = reader.next_line()) { + try { + ValueTokenizer values(line); + iloc = values.next_string(); + jloc = values.next_string(); + potential_name = values.next_string(); + if (iloc == arg[0] && jloc == arg[1] && potential_name == "dh") { + qeff_dh_one = values.next_double(); + break; + } else continue; + } catch (std::exception &e) { + error->one(FLERR, "Problem parsing oxDNA2 potential file: {}", e.what()); + } + } + if (iloc != arg[0] || jloc != arg[1] || potential_name != "dh") error->one(FLERR, "No corresponding dh potential found in file {} for pair type {} {}", arg[4], arg[0], arg[1]); + } + MPI_Bcast(&qeff_dh_one, 1, MPI_DOUBLE, 0, world); + } else qeff_dh_one = utils::numeric(FLERR,arg[4],false,lmp); // else, it is effective charge double lambda_dh_one, kappa_dh_one, qeff_dh_pf_one; double b_dh_one, cut_dh_ast_one, cut_dh_c_one; @@ -315,7 +341,8 @@ void PairOxdna2Dh::coeff(int narg, char **arg) The numerical factor is the Debye length in s.u. lambda(T = 300 K = 0.1) = sqrt(eps_0 * eps_r * k_B * T/(2 * N_A * e^2 * 1000 mol/m^3)) - * 1/oxDNA_energy_unit + * 1/oxDNA_length_unit for LJ units, or; + * [(8.518 * sqrt(k_B / 4.142e-20))/oxDNA_length_unit] for real units (see B. Snodin et al., J. Chem. Phys. 142, 234901 (2015).) We use @@ -328,7 +355,7 @@ void PairOxdna2Dh::coeff(int narg, char **arg) oxDNA_length_unit = 8.518e-10 m */ - lambda_dh_one = 0.3616455075438555*sqrt(T/0.1/rhos_dh_one); + lambda_dh_one = ConstantsOxdna::get_lambda_dh_one_prefactor()*sqrt(T/0.1/rhos_dh_one); kappa_dh_one = 1.0/lambda_dh_one; // prefactor in DH interaction containing qeff^2 @@ -337,14 +364,15 @@ void PairOxdna2Dh::coeff(int narg, char **arg) NOTE: The numerical factor is qeff_dh_pf = e^2/(4 * pi * eps_0 * eps_r) - * 1/(oxDNA_energy_unit * oxDNA_length_unit) + * 1/(oxDNA_energy_unit * oxDNA_length_unit) for LJ units, or; + * [(~5.96169* 8.518)/(oxDNA_energy_unit * oxDNA_length_unit)] for real units (see B. Snodin et al., J. Chem. Phys. 142, 234901 (2015).) In addition to the above units we use oxDNA_energy_unit = 4.142e-20 J */ - qeff_dh_pf_one = 0.08173808693529228*qeff_dh_one*qeff_dh_one; + qeff_dh_pf_one = ConstantsOxdna::get_qeff_dh_pf_one_prefactor()*qeff_dh_one*qeff_dh_one; // smoothing parameters - determined through continuity and differentiability diff --git a/src/CG-DNA/pair_oxdna2_excv.cpp b/src/CG-DNA/pair_oxdna2_excv.cpp index f0f9ca904e..2b047dae8e 100644 --- a/src/CG-DNA/pair_oxdna2_excv.cpp +++ b/src/CG-DNA/pair_oxdna2_excv.cpp @@ -14,6 +14,7 @@ Contributing author: Oliver Henrich (University of Strathclyde, Glasgow) ------------------------------------------------------------------------- */ +#include "constants_oxdna.h" #include "pair_oxdna2_excv.h" using namespace LAMMPS_NS; @@ -24,7 +25,9 @@ using namespace LAMMPS_NS; void PairOxdna2Excv::compute_interaction_sites(double e1[3], double e2[3], double /*e3*/[3], double rs[3], double rb[3]) { - double d_cs_x = -0.34, d_cs_y = +0.3408, d_cb = +0.4; + double d_cs_x = ConstantsOxdna::get_d_cs_x(); + double d_cs_y = ConstantsOxdna::get_d_cs_y(); + double d_cb = ConstantsOxdna::get_d_cb(); rs[0] = d_cs_x * e1[0] + d_cs_y * e2[0]; rs[1] = d_cs_x * e1[1] + d_cs_y * e2[1]; diff --git a/src/CG-DNA/pair_oxdna_coaxstk.cpp b/src/CG-DNA/pair_oxdna_coaxstk.cpp index 1fb0ad9b00..3b955a7db4 100644 --- a/src/CG-DNA/pair_oxdna_coaxstk.cpp +++ b/src/CG-DNA/pair_oxdna_coaxstk.cpp @@ -19,6 +19,7 @@ #include "atom.h" #include "comm.h" +#include "constants_oxdna.h" #include "error.h" #include "force.h" #include "math_const.h" @@ -26,6 +27,7 @@ #include "memory.h" #include "mf_oxdna.h" #include "neigh_list.h" +#include "potential_file_reader.h" #include #include @@ -124,7 +126,8 @@ void PairOxdnaCoaxstk::compute(int eflag, int vflag) double dcdrax,dcdray,dcdraz; // distances COM-backbone site, COM-stacking site - double d_cs=-0.4, d_cst=+0.34; + double d_cs = ConstantsOxdna::get_d_cs(); + double d_cst = ConstantsOxdna::get_d_cst(); // vectors COM-backbone site, COM-stacking site in lab frame double ra_cs[3],ra_cst[3]; double rb_cs[3],rb_cst[3]; @@ -691,7 +694,7 @@ void PairOxdnaCoaxstk::coeff(int narg, char **arg) { int count; - if (narg != 23) error->all(FLERR,"Incorrect args for pair coefficients in oxdna/coaxstk"); + if (narg != 3 && narg != 23) error->all(FLERR,"Incorrect args for pair coefficients in oxdna/coaxstk"); if (!allocated) allocate(); int ilo,ihi,jlo,jhi; @@ -719,32 +722,109 @@ void PairOxdnaCoaxstk::coeff(int narg, char **arg) double a_cxst3p_one, cosphi_cxst3p_ast_one, b_cxst3p_one, cosphi_cxst3p_c_one; double a_cxst4p_one, cosphi_cxst4p_ast_one, b_cxst4p_one, cosphi_cxst4p_c_one; - k_cxst_one = utils::numeric(FLERR,arg[2],false,lmp); - cut_cxst_0_one = utils::numeric(FLERR,arg[3],false,lmp); - cut_cxst_c_one = utils::numeric(FLERR,arg[4],false,lmp); - cut_cxst_lo_one = utils::numeric(FLERR,arg[5],false,lmp); - cut_cxst_hi_one = utils::numeric(FLERR,arg[6],false,lmp); + if (narg == 23) { + k_cxst_one = utils::numeric(FLERR,arg[2],false,lmp); + cut_cxst_0_one = utils::numeric(FLERR,arg[3],false,lmp); + cut_cxst_c_one = utils::numeric(FLERR,arg[4],false,lmp); + cut_cxst_lo_one = utils::numeric(FLERR,arg[5],false,lmp); + cut_cxst_hi_one = utils::numeric(FLERR,arg[6],false,lmp); - a_cxst1_one = utils::numeric(FLERR,arg[7],false,lmp); - theta_cxst1_0_one = utils::numeric(FLERR,arg[8],false,lmp); - dtheta_cxst1_ast_one = utils::numeric(FLERR,arg[9],false,lmp); + a_cxst1_one = utils::numeric(FLERR,arg[7],false,lmp); + theta_cxst1_0_one = utils::numeric(FLERR,arg[8],false,lmp); + dtheta_cxst1_ast_one = utils::numeric(FLERR,arg[9],false,lmp); - a_cxst4_one = utils::numeric(FLERR,arg[10],false,lmp); - theta_cxst4_0_one = utils::numeric(FLERR,arg[11],false,lmp); - dtheta_cxst4_ast_one = utils::numeric(FLERR,arg[12],false,lmp); + a_cxst4_one = utils::numeric(FLERR,arg[10],false,lmp); + theta_cxst4_0_one = utils::numeric(FLERR,arg[11],false,lmp); + dtheta_cxst4_ast_one = utils::numeric(FLERR,arg[12],false,lmp); - a_cxst5_one = utils::numeric(FLERR,arg[13],false,lmp); - theta_cxst5_0_one = utils::numeric(FLERR,arg[14],false,lmp); - dtheta_cxst5_ast_one = utils::numeric(FLERR,arg[15],false,lmp); + a_cxst5_one = utils::numeric(FLERR,arg[13],false,lmp); + theta_cxst5_0_one = utils::numeric(FLERR,arg[14],false,lmp); + dtheta_cxst5_ast_one = utils::numeric(FLERR,arg[15],false,lmp); - a_cxst6_one = utils::numeric(FLERR,arg[16],false,lmp); - theta_cxst6_0_one = utils::numeric(FLERR,arg[17],false,lmp); - dtheta_cxst6_ast_one = utils::numeric(FLERR,arg[18],false,lmp); + a_cxst6_one = utils::numeric(FLERR,arg[16],false,lmp); + theta_cxst6_0_one = utils::numeric(FLERR,arg[17],false,lmp); + dtheta_cxst6_ast_one = utils::numeric(FLERR,arg[18],false,lmp); - a_cxst3p_one = utils::numeric(FLERR,arg[19],false,lmp); - cosphi_cxst3p_ast_one = utils::numeric(FLERR,arg[20],false,lmp); - a_cxst4p_one = utils::numeric(FLERR,arg[21],false,lmp); - cosphi_cxst4p_ast_one = utils::numeric(FLERR,arg[22],false,lmp); + a_cxst3p_one = utils::numeric(FLERR,arg[19],false,lmp); + cosphi_cxst3p_ast_one = utils::numeric(FLERR,arg[20],false,lmp); + a_cxst4p_one = utils::numeric(FLERR,arg[21],false,lmp); + cosphi_cxst4p_ast_one = utils::numeric(FLERR,arg[22],false,lmp); + } else { + if (comm->me == 0) { // read values from potential file + PotentialFileReader reader(lmp, arg[2], "oxdna potential", " (coaxstk)"); + char * line; + std::string iloc, jloc, potential_name; + + while(line = reader.next_line()) { + try { + ValueTokenizer values(line); + iloc = values.next_string(); + jloc = values.next_string(); + potential_name = values.next_string(); + if (iloc == arg[0] && jloc == arg[1] && potential_name == "coaxstk") { + k_cxst_one = values.next_double(); + cut_cxst_0_one = values.next_double(); + cut_cxst_c_one = values.next_double(); + cut_cxst_lo_one = values.next_double(); + cut_cxst_hi_one = values.next_double(); + + a_cxst1_one = values.next_double(); + theta_cxst1_0_one = values.next_double(); + dtheta_cxst1_ast_one = values.next_double(); + + a_cxst4_one = values.next_double(); + theta_cxst4_0_one = values.next_double(); + dtheta_cxst4_ast_one = values.next_double(); + + a_cxst5_one = values.next_double(); + theta_cxst5_0_one = values.next_double(); + dtheta_cxst5_ast_one = values.next_double(); + + a_cxst6_one = values.next_double(); + theta_cxst6_0_one = values.next_double(); + dtheta_cxst6_ast_one = values.next_double(); + + a_cxst3p_one = values.next_double(); + cosphi_cxst3p_ast_one = values.next_double(); + a_cxst4p_one = values.next_double(); + cosphi_cxst4p_ast_one = values.next_double(); + + break; + } else continue; + } catch (std::exception &e) { + error->one(FLERR, "Problem parsing oxDNA potential file: {}", e.what()); + } + } + if (iloc != arg[0] || jloc != arg[1] || potential_name != "coaxstk") error->one(FLERR, "No corresponding coaxstk potential found in file {} for pair type {} {}", arg[2], arg[0], arg[1]); + } + + MPI_Bcast(&k_cxst_one, 1, MPI_DOUBLE, 0, world); + MPI_Bcast(&cut_cxst_0_one, 1, MPI_DOUBLE, 0, world); + MPI_Bcast(&cut_cxst_c_one, 1, MPI_DOUBLE, 0, world); + MPI_Bcast(&cut_cxst_lo_one, 1, MPI_DOUBLE, 0, world); + MPI_Bcast(&cut_cxst_hi_one, 1, MPI_DOUBLE, 0, world); + + MPI_Bcast(&a_cxst1_one, 1, MPI_DOUBLE, 0, world); + MPI_Bcast(&theta_cxst1_0_one, 1, MPI_DOUBLE, 0, world); + MPI_Bcast(&dtheta_cxst1_ast_one, 1, MPI_DOUBLE, 0, world); + + MPI_Bcast(&a_cxst4_one, 1, MPI_DOUBLE, 0, world); + MPI_Bcast(&theta_cxst4_0_one, 1, MPI_DOUBLE, 0, world); + MPI_Bcast(&dtheta_cxst4_ast_one, 1, MPI_DOUBLE, 0, world); + + MPI_Bcast(&a_cxst5_one, 1, MPI_DOUBLE, 0, world); + MPI_Bcast(&theta_cxst5_0_one, 1, MPI_DOUBLE, 0, world); + MPI_Bcast(&dtheta_cxst5_ast_one, 1, MPI_DOUBLE, 0, world); + + MPI_Bcast(&a_cxst6_one, 1, MPI_DOUBLE, 0, world); + MPI_Bcast(&theta_cxst6_0_one, 1, MPI_DOUBLE, 0, world); + MPI_Bcast(&dtheta_cxst6_ast_one, 1, MPI_DOUBLE, 0, world); + + MPI_Bcast(&a_cxst3p_one, 1, MPI_DOUBLE, 0, world); + MPI_Bcast(&cosphi_cxst3p_ast_one, 1, MPI_DOUBLE, 0, world); + MPI_Bcast(&a_cxst4p_one, 1, MPI_DOUBLE, 0, world); + MPI_Bcast(&cosphi_cxst4p_ast_one, 1, MPI_DOUBLE, 0, world); + } b_cxst_lo_one = 0.25 * (cut_cxst_lo_one - cut_cxst_0_one) * (cut_cxst_lo_one - cut_cxst_0_one)/ (0.5 * (cut_cxst_lo_one - cut_cxst_0_one) * (cut_cxst_lo_one - cut_cxst_0_one) - diff --git a/src/CG-DNA/pair_oxdna_excv.cpp b/src/CG-DNA/pair_oxdna_excv.cpp index 60df42404d..3a98712818 100644 --- a/src/CG-DNA/pair_oxdna_excv.cpp +++ b/src/CG-DNA/pair_oxdna_excv.cpp @@ -20,12 +20,14 @@ #include "atom.h" #include "atom_vec_ellipsoid.h" #include "comm.h" +#include "constants_oxdna.h" #include "error.h" #include "force.h" #include "math_extra.h" #include "memory.h" #include "mf_oxdna.h" #include "neigh_list.h" +#include "potential_file_reader.h" #include #include @@ -97,7 +99,8 @@ PairOxdnaExcv::~PairOxdnaExcv() void PairOxdnaExcv::compute_interaction_sites(double e1[3], double /*e2*/[3], double /*e3*/[3], double rs[3], double rb[3]) { - double d_cs=-0.4, d_cb=+0.4; + double d_cs = ConstantsOxdna::get_d_cs(); + double d_cb = ConstantsOxdna::get_d_cb(); rs[0] = d_cs*e1[0]; rs[1] = d_cs*e1[1]; @@ -500,7 +503,7 @@ void PairOxdnaExcv::coeff(int narg, char **arg) { int count; - if (narg != 11) error->all(FLERR,"Incorrect args for pair coefficients in oxdna/excv"); + if (narg != 3 && narg != 11) error->all(FLERR,"Incorrect args for pair coefficients in oxdna/excv"); if (!allocated) allocate(); int ilo,ihi,jlo,jhi; @@ -518,12 +521,73 @@ void PairOxdnaExcv::coeff(int narg, char **arg) double epsilon_bb_one, sigma_bb_one; double cut_bb_ast_one, cut_bb_c_one, b_bb_one; - // Excluded volume interaction - // LJ parameters - epsilon_ss_one = utils::numeric(FLERR,arg[2],false,lmp); - sigma_ss_one = utils::numeric(FLERR,arg[3],false,lmp); - cut_ss_ast_one = utils::numeric(FLERR,arg[4],false,lmp); + if (narg == 11) { + // Excluded volume interaction + // LJ parameters + epsilon_ss_one = utils::numeric(FLERR,arg[2],false,lmp); + sigma_ss_one = utils::numeric(FLERR,arg[3],false,lmp); + cut_ss_ast_one = utils::numeric(FLERR,arg[4],false,lmp); + // LJ parameters + epsilon_sb_one = utils::numeric(FLERR,arg[5],false,lmp); + sigma_sb_one = utils::numeric(FLERR,arg[6],false,lmp); + cut_sb_ast_one = utils::numeric(FLERR,arg[7],false,lmp); + + // LJ parameters + epsilon_bb_one = utils::numeric(FLERR,arg[8],false,lmp); + sigma_bb_one = utils::numeric(FLERR,arg[9],false,lmp); + cut_bb_ast_one = utils::numeric(FLERR,arg[10],false,lmp); + } else { + if (comm->me == 0) { + PotentialFileReader reader(lmp, arg[2], "oxdna potential", " (excv)"); + char * line; + std::string iloc, jloc, potential_name; + + while(line = reader.next_line()) { + try { + ValueTokenizer values(line); + iloc = values.next_string(); + jloc = values.next_string(); + potential_name = values.next_string(); + if (iloc == arg[0] && jloc == arg[1] && potential_name == "excv") { + // Excluded volume interaction + // LJ parameters + epsilon_ss_one = values.next_double(); + sigma_ss_one = values.next_double(); + cut_ss_ast_one = values.next_double(); + + // LJ parameters + epsilon_sb_one = values.next_double(); + sigma_sb_one = values.next_double(); + cut_sb_ast_one = values.next_double(); + + // LJ parameters + epsilon_bb_one = values.next_double(); + sigma_bb_one = values.next_double(); + cut_bb_ast_one = values.next_double(); + + break; + } else continue; + } catch (std::exception &e) { + error->one(FLERR, "Problem parsing oxDNA potential file: {}", e.what()); + } + } + if (iloc != arg[0] || jloc != arg[1] || potential_name != "excv") error->one(FLERR, "No corresponding excv potential found in file {} for pair type {} {}", arg[2], arg[0], arg[1]); + } + + MPI_Bcast(&epsilon_ss_one, 1, MPI_DOUBLE, 0, world); + MPI_Bcast(&sigma_ss_one, 1, MPI_DOUBLE, 0, world); + MPI_Bcast(&cut_ss_ast_one, 1, MPI_DOUBLE, 0, world); + + MPI_Bcast(&epsilon_sb_one, 1, MPI_DOUBLE, 0, world); + MPI_Bcast(&sigma_sb_one, 1, MPI_DOUBLE, 0, world); + MPI_Bcast(&cut_sb_ast_one, 1, MPI_DOUBLE, 0, world); + + MPI_Bcast(&epsilon_bb_one, 1, MPI_DOUBLE, 0, world); + MPI_Bcast(&sigma_bb_one, 1, MPI_DOUBLE, 0, world); + MPI_Bcast(&cut_bb_ast_one, 1, MPI_DOUBLE, 0, world); + } + // smoothing - determined through continuity and differentiability b_ss_one = 4.0/sigma_ss_one *(6.0*pow(sigma_ss_one/cut_ss_ast_one,7)-12.0*pow(sigma_ss_one/cut_ss_ast_one,13)) @@ -550,11 +614,6 @@ void PairOxdnaExcv::coeff(int narg, char **arg) count = 0; - // LJ parameters - epsilon_sb_one = utils::numeric(FLERR,arg[5],false,lmp); - sigma_sb_one = utils::numeric(FLERR,arg[6],false,lmp); - cut_sb_ast_one = utils::numeric(FLERR,arg[7],false,lmp); - // smoothing - determined through continuity and differentiability b_sb_one = 4.0/sigma_sb_one *(6.0*pow(sigma_sb_one/cut_sb_ast_one,7)-12.0*pow(sigma_sb_one/cut_sb_ast_one,13)) @@ -581,11 +640,6 @@ void PairOxdnaExcv::coeff(int narg, char **arg) count = 0; - // LJ parameters - epsilon_bb_one = utils::numeric(FLERR,arg[8],false,lmp); - sigma_bb_one = utils::numeric(FLERR,arg[9],false,lmp); - cut_bb_ast_one = utils::numeric(FLERR,arg[10],false,lmp); - // smoothing - determined through continuity and differentiability b_bb_one = 4.0/sigma_bb_one *(6.0*pow(sigma_bb_one/cut_bb_ast_one,7)-12.0*pow(sigma_bb_one/cut_bb_ast_one,13)) diff --git a/src/CG-DNA/pair_oxdna_hbond.cpp b/src/CG-DNA/pair_oxdna_hbond.cpp index 2beadc2503..4763f7412d 100644 --- a/src/CG-DNA/pair_oxdna_hbond.cpp +++ b/src/CG-DNA/pair_oxdna_hbond.cpp @@ -19,12 +19,14 @@ #include "atom.h" #include "comm.h" +#include "constants_oxdna.h" #include "error.h" #include "force.h" #include "math_extra.h" #include "memory.h" #include "mf_oxdna.h" #include "neigh_list.h" +#include "potential_file_reader.h" #include #include @@ -145,7 +147,7 @@ void PairOxdnaHbond::compute(int eflag, int vflag) double theta8,t8dir[3],cost8; // distance COM-hbonding site - double d_chb=+0.4; + double d_chb = ConstantsOxdna::get_d_chb(); // vectors COM-h-bonding site in lab frame double ra_chb[3],rb_chb[3]; // Cartesian unit vectors in lab frame @@ -634,7 +636,7 @@ void PairOxdnaHbond::coeff(int narg, char **arg) { int count; - if (narg != 27) error->all(FLERR,"Incorrect args for pair coefficients in oxdna/hbond"); + if (narg != 4 && narg != 27) error->all(FLERR,"Incorrect args for pair coefficients in oxdna/hbond"); if (!allocated) allocate(); int ilo,ihi,jlo,jhi,imod4,jmod4; @@ -671,36 +673,123 @@ void PairOxdnaHbond::coeff(int narg, char **arg) if (strcmp(arg[2],"seqav") == 0) seqdepflag = 0; if (strcmp(arg[2],"seqdep") == 0) seqdepflag = 1; - epsilon_hb_one = utils::numeric(FLERR,arg[3],false,lmp); - a_hb_one = utils::numeric(FLERR,arg[4],false,lmp); - cut_hb_0_one = utils::numeric(FLERR,arg[5],false,lmp); - cut_hb_c_one = utils::numeric(FLERR,arg[6],false,lmp); - cut_hb_lo_one = utils::numeric(FLERR,arg[7],false,lmp); - cut_hb_hi_one = utils::numeric(FLERR,arg[8],false,lmp); + if (narg == 27) { + epsilon_hb_one = utils::numeric(FLERR,arg[3],false,lmp); + a_hb_one = utils::numeric(FLERR,arg[4],false,lmp); + cut_hb_0_one = utils::numeric(FLERR,arg[5],false,lmp); + cut_hb_c_one = utils::numeric(FLERR,arg[6],false,lmp); + cut_hb_lo_one = utils::numeric(FLERR,arg[7],false,lmp); + cut_hb_hi_one = utils::numeric(FLERR,arg[8],false,lmp); - a_hb1_one = utils::numeric(FLERR,arg[9],false,lmp); - theta_hb1_0_one = utils::numeric(FLERR,arg[10],false,lmp); - dtheta_hb1_ast_one = utils::numeric(FLERR,arg[11],false,lmp); + a_hb1_one = utils::numeric(FLERR,arg[9],false,lmp); + theta_hb1_0_one = utils::numeric(FLERR,arg[10],false,lmp); + dtheta_hb1_ast_one = utils::numeric(FLERR,arg[11],false,lmp); - a_hb2_one = utils::numeric(FLERR,arg[12],false,lmp); - theta_hb2_0_one = utils::numeric(FLERR,arg[13],false,lmp); - dtheta_hb2_ast_one = utils::numeric(FLERR,arg[14],false,lmp); + a_hb2_one = utils::numeric(FLERR,arg[12],false,lmp); + theta_hb2_0_one = utils::numeric(FLERR,arg[13],false,lmp); + dtheta_hb2_ast_one = utils::numeric(FLERR,arg[14],false,lmp); - a_hb3_one = utils::numeric(FLERR,arg[15],false,lmp); - theta_hb3_0_one = utils::numeric(FLERR,arg[16],false,lmp); - dtheta_hb3_ast_one = utils::numeric(FLERR,arg[17],false,lmp); + a_hb3_one = utils::numeric(FLERR,arg[15],false,lmp); + theta_hb3_0_one = utils::numeric(FLERR,arg[16],false,lmp); + dtheta_hb3_ast_one = utils::numeric(FLERR,arg[17],false,lmp); - a_hb4_one = utils::numeric(FLERR,arg[18],false,lmp); - theta_hb4_0_one = utils::numeric(FLERR,arg[19],false,lmp); - dtheta_hb4_ast_one = utils::numeric(FLERR,arg[20],false,lmp); + a_hb4_one = utils::numeric(FLERR,arg[18],false,lmp); + theta_hb4_0_one = utils::numeric(FLERR,arg[19],false,lmp); + dtheta_hb4_ast_one = utils::numeric(FLERR,arg[20],false,lmp); - a_hb7_one = utils::numeric(FLERR,arg[21],false,lmp); - theta_hb7_0_one = utils::numeric(FLERR,arg[22],false,lmp); - dtheta_hb7_ast_one = utils::numeric(FLERR,arg[23],false,lmp); + a_hb7_one = utils::numeric(FLERR,arg[21],false,lmp); + theta_hb7_0_one = utils::numeric(FLERR,arg[22],false,lmp); + dtheta_hb7_ast_one = utils::numeric(FLERR,arg[23],false,lmp); + + a_hb8_one = utils::numeric(FLERR,arg[24],false,lmp); + theta_hb8_0_one = utils::numeric(FLERR,arg[25],false,lmp); + dtheta_hb8_ast_one = utils::numeric(FLERR,arg[26],false,lmp); + } else { // read values from potential file + if (comm->me == 0) { + PotentialFileReader reader(lmp, arg[3], "oxdna potential", " (hbond)"); + char * line; + std::string iloc, jloc, potential_name; + + while(line = reader.next_line()) { + try { + ValueTokenizer values(line); + iloc = values.next_string(); + jloc = values.next_string(); + potential_name = values.next_string(); + if (iloc == arg[0] && jloc == arg[1] && potential_name == "hbond") { + + epsilon_hb_one = values.next_double(); + a_hb_one = values.next_double(); + cut_hb_0_one = values.next_double(); + cut_hb_c_one = values.next_double(); + cut_hb_lo_one = values.next_double(); + cut_hb_hi_one = values.next_double(); + + a_hb1_one = values.next_double(); + theta_hb1_0_one = values.next_double(); + dtheta_hb1_ast_one = values.next_double(); + + a_hb2_one = values.next_double(); + theta_hb2_0_one = values.next_double(); + dtheta_hb2_ast_one = values.next_double(); + + a_hb3_one = values.next_double(); + theta_hb3_0_one = values.next_double(); + dtheta_hb3_ast_one = values.next_double(); + + a_hb4_one = values.next_double(); + theta_hb4_0_one = values.next_double(); + dtheta_hb4_ast_one = values.next_double(); + + a_hb7_one = values.next_double(); + theta_hb7_0_one = values.next_double(); + dtheta_hb7_ast_one = values.next_double(); + + a_hb8_one = values.next_double(); + theta_hb8_0_one = values.next_double(); + dtheta_hb8_ast_one = values.next_double(); + + break; + } else continue; + } catch (std::exception &e) { + error->one(FLERR, "Problem parsing oxDNA potential file: {}", e.what()); + } + } + if (iloc != arg[0] || jloc != arg[1] || potential_name != "hbond") error->one(FLERR, "No corresponding hbond potential found in file {} for pair type {} {}", arg[3], arg[0], arg[1]); + } + + MPI_Bcast(&epsilon_hb_one, 1, MPI_DOUBLE, 0, world); + MPI_Bcast(&a_hb_one, 1, MPI_DOUBLE, 0, world); + MPI_Bcast(&cut_hb_0_one, 1, MPI_DOUBLE, 0, world); + MPI_Bcast(&cut_hb_c_one, 1, MPI_DOUBLE, 0, world); + MPI_Bcast(&cut_hb_lo_one, 1, MPI_DOUBLE, 0, world); + MPI_Bcast(&cut_hb_hi_one, 1, MPI_DOUBLE, 0, world); + + MPI_Bcast(&a_hb1_one, 1, MPI_DOUBLE, 0, world); + MPI_Bcast(&theta_hb1_0_one, 1, MPI_DOUBLE, 0, world); + MPI_Bcast(&dtheta_hb1_ast_one, 1, MPI_DOUBLE, 0, world); + + MPI_Bcast(&a_hb2_one, 1, MPI_DOUBLE, 0, world); + MPI_Bcast(&theta_hb2_0_one, 1, MPI_DOUBLE, 0, world); + MPI_Bcast(&dtheta_hb2_ast_one, 1, MPI_DOUBLE, 0, world); + + MPI_Bcast(&a_hb3_one, 1, MPI_DOUBLE, 0, world); + MPI_Bcast(&theta_hb3_0_one, 1, MPI_DOUBLE, 0, world); + MPI_Bcast(&dtheta_hb3_ast_one, 1, MPI_DOUBLE, 0, world); + + MPI_Bcast(&a_hb4_one, 1, MPI_DOUBLE, 0, world); + MPI_Bcast(&theta_hb4_0_one, 1, MPI_DOUBLE, 0, world); + MPI_Bcast(&dtheta_hb4_ast_one, 1, MPI_DOUBLE, 0, world); + + MPI_Bcast(&a_hb7_one, 1, MPI_DOUBLE, 0, world); + MPI_Bcast(&theta_hb7_0_one, 1, MPI_DOUBLE, 0, world); + MPI_Bcast(&dtheta_hb7_ast_one, 1, MPI_DOUBLE, 0, world); + + MPI_Bcast(&a_hb8_one, 1, MPI_DOUBLE, 0, world); + MPI_Bcast(&theta_hb8_0_one, 1, MPI_DOUBLE, 0, world); + MPI_Bcast(&dtheta_hb8_ast_one, 1, MPI_DOUBLE, 0, world); + } - a_hb8_one = utils::numeric(FLERR,arg[24],false,lmp); - theta_hb8_0_one = utils::numeric(FLERR,arg[25],false,lmp); - dtheta_hb8_ast_one = utils::numeric(FLERR,arg[26],false,lmp); b_hb_lo_one = 2*a_hb_one*exp(-a_hb_one*(cut_hb_lo_one-cut_hb_0_one))* 2*a_hb_one*exp(-a_hb_one*(cut_hb_lo_one-cut_hb_0_one))* diff --git a/src/CG-DNA/pair_oxdna_stk.cpp b/src/CG-DNA/pair_oxdna_stk.cpp index cf20d527fe..7991288c1e 100644 --- a/src/CG-DNA/pair_oxdna_stk.cpp +++ b/src/CG-DNA/pair_oxdna_stk.cpp @@ -19,6 +19,7 @@ #include "atom.h" #include "comm.h" +#include "constants_oxdna.h" #include "error.h" #include "force.h" #include "math_extra.h" @@ -26,6 +27,7 @@ #include "mf_oxdna.h" #include "neighbor.h" #include "neigh_list.h" +#include "potential_file_reader.h" #include #include @@ -223,7 +225,8 @@ void PairOxdnaStk::compute(int eflag, int vflag) double cosphi1,cosphi2,cosphi1dir[3],cosphi2dir[3]; // distances COM-backbone site, COM-stacking site - double d_cs=-0.4, d_cst=+0.34; + double d_cs = ConstantsOxdna::get_d_cs(); + double d_cst = ConstantsOxdna::get_d_cst(); // vectors COM-backbone site, COM-stacking site in lab frame double ra_cs[3],ra_cst[3]; double rb_cs[3],rb_cst[3]; @@ -775,7 +778,7 @@ void PairOxdnaStk::coeff(int narg, char **arg) { int count; - if (narg != 24) error->all(FLERR,"Incorrect args for pair coefficients in oxdna/stk"); + if (narg != 7 && narg != 24) error->all(FLERR,"Incorrect args for pair coefficients in oxdna/stk"); if (!allocated) allocate(); int ilo,ihi,jlo,jhi,imod4,jmod4; @@ -812,25 +815,89 @@ void PairOxdnaStk::coeff(int narg, char **arg) kappa_st_one = utils::numeric(FLERR,arg[5],false,lmp); epsilon_st_one = stacking_strength(xi_st_one, kappa_st_one, T); - a_st_one = utils::numeric(FLERR,arg[6],false,lmp); - cut_st_0_one = utils::numeric(FLERR,arg[7],false,lmp); - cut_st_c_one = utils::numeric(FLERR,arg[8],false,lmp); - cut_st_lo_one = utils::numeric(FLERR,arg[9],false,lmp); - cut_st_hi_one = utils::numeric(FLERR,arg[10],false,lmp); + if (narg == 24) { // values are listed in input + a_st_one = utils::numeric(FLERR,arg[6],false,lmp); + cut_st_0_one = utils::numeric(FLERR,arg[7],false,lmp); + cut_st_c_one = utils::numeric(FLERR,arg[8],false,lmp); + cut_st_lo_one = utils::numeric(FLERR,arg[9],false,lmp); + cut_st_hi_one = utils::numeric(FLERR,arg[10],false,lmp); - a_st4_one = utils::numeric(FLERR,arg[11],false,lmp); - theta_st4_0_one = utils::numeric(FLERR,arg[12],false,lmp); - dtheta_st4_ast_one = utils::numeric(FLERR,arg[13],false,lmp); - a_st5_one = utils::numeric(FLERR,arg[14],false,lmp); - theta_st5_0_one = utils::numeric(FLERR,arg[15],false,lmp); - dtheta_st5_ast_one = utils::numeric(FLERR,arg[16],false,lmp); - a_st6_one = utils::numeric(FLERR,arg[17],false,lmp); - theta_st6_0_one = utils::numeric(FLERR,arg[18],false,lmp); - dtheta_st6_ast_one = utils::numeric(FLERR,arg[19],false,lmp); - a_st1_one = utils::numeric(FLERR,arg[20],false,lmp); - cosphi_st1_ast_one = utils::numeric(FLERR,arg[21],false,lmp); - a_st2_one = utils::numeric(FLERR,arg[22],false,lmp); - cosphi_st2_ast_one = utils::numeric(FLERR,arg[23],false,lmp); + a_st4_one = utils::numeric(FLERR,arg[11],false,lmp); + theta_st4_0_one = utils::numeric(FLERR,arg[12],false,lmp); + dtheta_st4_ast_one = utils::numeric(FLERR,arg[13],false,lmp); + a_st5_one = utils::numeric(FLERR,arg[14],false,lmp); + theta_st5_0_one = utils::numeric(FLERR,arg[15],false,lmp); + dtheta_st5_ast_one = utils::numeric(FLERR,arg[16],false,lmp); + a_st6_one = utils::numeric(FLERR,arg[17],false,lmp); + theta_st6_0_one = utils::numeric(FLERR,arg[18],false,lmp); + dtheta_st6_ast_one = utils::numeric(FLERR,arg[19],false,lmp); + a_st1_one = utils::numeric(FLERR,arg[20],false,lmp); + cosphi_st1_ast_one = utils::numeric(FLERR,arg[21],false,lmp); + a_st2_one = utils::numeric(FLERR,arg[22],false,lmp); + cosphi_st2_ast_one = utils::numeric(FLERR,arg[23],false,lmp); + } else { // read values from potential file + if (comm->me == 0) { + PotentialFileReader reader(lmp, arg[6], "oxdna potential", " (stk)"); + char * line; + std::string iloc, jloc, potential_name; + + while(line = reader.next_line()) { + try { + ValueTokenizer values(line); + iloc = values.next_string(); + jloc = values.next_string(); + potential_name = values.next_string(); + if (iloc == arg[0] && jloc == arg[1] && potential_name == "stk") { + + a_st_one = values.next_double(); + cut_st_0_one = values.next_double(); + cut_st_c_one = values.next_double(); + cut_st_lo_one = values.next_double(); + cut_st_hi_one = values.next_double(); + + a_st4_one = values.next_double(); + theta_st4_0_one = values.next_double(); + dtheta_st4_ast_one = values.next_double(); + a_st5_one = values.next_double(); + theta_st5_0_one = values.next_double(); + dtheta_st5_ast_one = values.next_double(); + a_st6_one = values.next_double(); + theta_st6_0_one = values.next_double(); + dtheta_st6_ast_one = values.next_double(); + a_st1_one = values.next_double(); + cosphi_st1_ast_one = values.next_double(); + a_st2_one = values.next_double(); + cosphi_st2_ast_one = values.next_double(); + + break; + } else continue; + } catch (std::exception &e) { + error->one(FLERR, "Problem parsing oxDNA potential file: {}", e.what()); + } + } + if (iloc != arg[0] || jloc != arg[1] || potential_name != "stk") error->one(FLERR, "No corresponding stk potential found in file {} for pair type {} {}", arg[4], arg[0], arg[1]); + } + + MPI_Bcast(&a_st_one, 1, MPI_DOUBLE, 0, world); + MPI_Bcast(&cut_st_0_one, 1, MPI_DOUBLE, 0, world); + MPI_Bcast(&cut_st_c_one, 1, MPI_DOUBLE, 0, world); + MPI_Bcast(&cut_st_lo_one, 1, MPI_DOUBLE, 0, world); + MPI_Bcast(&cut_st_hi_one, 1, MPI_DOUBLE, 0, world); + + MPI_Bcast(&a_st4_one, 1, MPI_DOUBLE, 0, world); + MPI_Bcast(&theta_st4_0_one, 1, MPI_DOUBLE, 0, world); + MPI_Bcast(&dtheta_st4_ast_one, 1, MPI_DOUBLE, 0, world); + MPI_Bcast(&a_st5_one, 1, MPI_DOUBLE, 0, world); + MPI_Bcast(&theta_st5_0_one, 1, MPI_DOUBLE, 0, world); + MPI_Bcast(&dtheta_st5_ast_one, 1, MPI_DOUBLE, 0, world); + MPI_Bcast(&a_st6_one, 1, MPI_DOUBLE, 0, world); + MPI_Bcast(&theta_st6_0_one, 1, MPI_DOUBLE, 0, world); + MPI_Bcast(&dtheta_st6_ast_one, 1, MPI_DOUBLE, 0, world); + MPI_Bcast(&a_st1_one, 1, MPI_DOUBLE, 0, world); + MPI_Bcast(&cosphi_st1_ast_one, 1, MPI_DOUBLE, 0, world); + MPI_Bcast(&a_st2_one, 1, MPI_DOUBLE, 0, world); + MPI_Bcast(&cosphi_st2_ast_one, 1, MPI_DOUBLE, 0, world); + } b_st_lo_one = 2*a_st_one*exp(-a_st_one*(cut_st_lo_one-cut_st_0_one))* 2*a_st_one*exp(-a_st_one*(cut_st_lo_one-cut_st_0_one))* diff --git a/src/CG-DNA/pair_oxdna_xstk.cpp b/src/CG-DNA/pair_oxdna_xstk.cpp index 1eb7caf23a..af80046b50 100644 --- a/src/CG-DNA/pair_oxdna_xstk.cpp +++ b/src/CG-DNA/pair_oxdna_xstk.cpp @@ -19,6 +19,7 @@ #include "atom.h" #include "comm.h" +#include "constants_oxdna.h" #include "error.h" #include "force.h" #include "math_const.h" @@ -26,6 +27,7 @@ #include "memory.h" #include "mf_oxdna.h" #include "neigh_list.h" +#include "potential_file_reader.h" #include #include @@ -122,7 +124,7 @@ void PairOxdnaXstk::compute(int eflag, int vflag) double theta8,theta8p,t8dir[3],cost8; // distance COM-h-bonding site - double d_chb=+0.4; + double d_chb = ConstantsOxdna::get_d_chb(); // vectors COM-h-bonding site in lab frame double ra_chb[3],rb_chb[3]; // Cartesian unit vectors in lab frame @@ -631,7 +633,7 @@ void PairOxdnaXstk::coeff(int narg, char **arg) { int count; - if (narg != 25) error->all(FLERR,"Incorrect args for pair coefficients in oxdna/xstk"); + if (narg != 3 && narg != 25) error->all(FLERR,"Incorrect args for pair coefficients in oxdna/xstk"); if (!allocated) allocate(); int ilo,ihi,jlo,jhi; @@ -662,35 +664,118 @@ void PairOxdnaXstk::coeff(int narg, char **arg) double a_xst8_one, theta_xst8_0_one, dtheta_xst8_ast_one; double b_xst8_one, dtheta_xst8_c_one; - k_xst_one = utils::numeric(FLERR,arg[2],false,lmp); - cut_xst_0_one = utils::numeric(FLERR,arg[3],false,lmp); - cut_xst_c_one = utils::numeric(FLERR,arg[4],false,lmp); - cut_xst_lo_one = utils::numeric(FLERR,arg[5],false,lmp); - cut_xst_hi_one = utils::numeric(FLERR,arg[6],false,lmp); + if (narg == 25) { + k_xst_one = utils::numeric(FLERR,arg[2],false,lmp); + cut_xst_0_one = utils::numeric(FLERR,arg[3],false,lmp); + cut_xst_c_one = utils::numeric(FLERR,arg[4],false,lmp); + cut_xst_lo_one = utils::numeric(FLERR,arg[5],false,lmp); + cut_xst_hi_one = utils::numeric(FLERR,arg[6],false,lmp); - a_xst1_one = utils::numeric(FLERR,arg[7],false,lmp); - theta_xst1_0_one = utils::numeric(FLERR,arg[8],false,lmp); - dtheta_xst1_ast_one = utils::numeric(FLERR,arg[9],false,lmp); + a_xst1_one = utils::numeric(FLERR,arg[7],false,lmp); + theta_xst1_0_one = utils::numeric(FLERR,arg[8],false,lmp); + dtheta_xst1_ast_one = utils::numeric(FLERR,arg[9],false,lmp); - a_xst2_one = utils::numeric(FLERR,arg[10],false,lmp); - theta_xst2_0_one = utils::numeric(FLERR,arg[11],false,lmp); - dtheta_xst2_ast_one = utils::numeric(FLERR,arg[12],false,lmp); + a_xst2_one = utils::numeric(FLERR,arg[10],false,lmp); + theta_xst2_0_one = utils::numeric(FLERR,arg[11],false,lmp); + dtheta_xst2_ast_one = utils::numeric(FLERR,arg[12],false,lmp); - a_xst3_one = utils::numeric(FLERR,arg[13],false,lmp); - theta_xst3_0_one = utils::numeric(FLERR,arg[14],false,lmp); - dtheta_xst3_ast_one = utils::numeric(FLERR,arg[15],false,lmp); + a_xst3_one = utils::numeric(FLERR,arg[13],false,lmp); + theta_xst3_0_one = utils::numeric(FLERR,arg[14],false,lmp); + dtheta_xst3_ast_one = utils::numeric(FLERR,arg[15],false,lmp); - a_xst4_one = utils::numeric(FLERR,arg[16],false,lmp); - theta_xst4_0_one = utils::numeric(FLERR,arg[17],false,lmp); - dtheta_xst4_ast_one = utils::numeric(FLERR,arg[18],false,lmp); + a_xst4_one = utils::numeric(FLERR,arg[16],false,lmp); + theta_xst4_0_one = utils::numeric(FLERR,arg[17],false,lmp); + dtheta_xst4_ast_one = utils::numeric(FLERR,arg[18],false,lmp); - a_xst7_one = utils::numeric(FLERR,arg[19],false,lmp); - theta_xst7_0_one = utils::numeric(FLERR,arg[20],false,lmp); - dtheta_xst7_ast_one = utils::numeric(FLERR,arg[21],false,lmp); + a_xst7_one = utils::numeric(FLERR,arg[19],false,lmp); + theta_xst7_0_one = utils::numeric(FLERR,arg[20],false,lmp); + dtheta_xst7_ast_one = utils::numeric(FLERR,arg[21],false,lmp); - a_xst8_one = utils::numeric(FLERR,arg[22],false,lmp); - theta_xst8_0_one = utils::numeric(FLERR,arg[23],false,lmp); - dtheta_xst8_ast_one = utils::numeric(FLERR,arg[24],false,lmp); + a_xst8_one = utils::numeric(FLERR,arg[22],false,lmp); + theta_xst8_0_one = utils::numeric(FLERR,arg[23],false,lmp); + dtheta_xst8_ast_one = utils::numeric(FLERR,arg[24],false,lmp); + } else { + if (comm->me == 0) { + PotentialFileReader reader(lmp, arg[2], "oxdna potential", " (xstk)"); + char * line; + std::string iloc, jloc, potential_name; + + while(line = reader.next_line()) { + try { + ValueTokenizer values(line); + iloc = values.next_string(); + jloc = values.next_string(); + potential_name = values.next_string(); + if (iloc == arg[0] && jloc == arg[1] && potential_name == "xstk") { + k_xst_one = values.next_double(); + cut_xst_0_one = values.next_double(); + cut_xst_c_one = values.next_double(); + cut_xst_lo_one = values.next_double(); + cut_xst_hi_one = values.next_double(); + + a_xst1_one = values.next_double(); + theta_xst1_0_one = values.next_double(); + dtheta_xst1_ast_one = values.next_double(); + + a_xst2_one = values.next_double(); + theta_xst2_0_one = values.next_double(); + dtheta_xst2_ast_one = values.next_double(); + + a_xst3_one = values.next_double(); + theta_xst3_0_one = values.next_double(); + dtheta_xst3_ast_one = values.next_double(); + + a_xst4_one = values.next_double(); + theta_xst4_0_one = values.next_double(); + dtheta_xst4_ast_one = values.next_double(); + + a_xst7_one = values.next_double(); + theta_xst7_0_one = values.next_double(); + dtheta_xst7_ast_one = values.next_double(); + + a_xst8_one = values.next_double(); + theta_xst8_0_one = values.next_double(); + dtheta_xst8_ast_one = values.next_double(); + + break; + } else continue; + } catch (std::exception &e) { + error->one(FLERR, "Problem parsing oxDNA potential file: {}", e.what()); + } + } + if (iloc != arg[0] || jloc != arg[1] || potential_name != "xstk") error->one(FLERR, "No corresponding xstk potential found in file {} for pair type {} {}", arg[2], arg[0], arg[1]); + } + + MPI_Bcast(&k_xst_one, 1, MPI_DOUBLE, 0, world); + MPI_Bcast(&cut_xst_0_one, 1, MPI_DOUBLE, 0, world); + MPI_Bcast(&cut_xst_c_one, 1, MPI_DOUBLE, 0, world); + MPI_Bcast(&cut_xst_lo_one, 1, MPI_DOUBLE, 0, world); + MPI_Bcast(&cut_xst_hi_one, 1, MPI_DOUBLE, 0, world); + + MPI_Bcast(&a_xst1_one, 1, MPI_DOUBLE, 0, world); + MPI_Bcast(&theta_xst1_0_one, 1, MPI_DOUBLE, 0, world); + MPI_Bcast(&dtheta_xst1_ast_one, 1, MPI_DOUBLE, 0, world); + + MPI_Bcast(&a_xst2_one, 1, MPI_DOUBLE, 0, world); + MPI_Bcast(&theta_xst2_0_one, 1, MPI_DOUBLE, 0, world); + MPI_Bcast(&dtheta_xst2_ast_one, 1, MPI_DOUBLE, 0, world); + + MPI_Bcast(&a_xst3_one, 1, MPI_DOUBLE, 0, world); + MPI_Bcast(&theta_xst3_0_one, 1, MPI_DOUBLE, 0, world); + MPI_Bcast(&dtheta_xst3_ast_one, 1, MPI_DOUBLE, 0, world); + + MPI_Bcast(&a_xst4_one, 1, MPI_DOUBLE, 0, world); + MPI_Bcast(&theta_xst4_0_one, 1, MPI_DOUBLE, 0, world); + MPI_Bcast(&dtheta_xst4_ast_one, 1, MPI_DOUBLE, 0, world); + + MPI_Bcast(&a_xst7_one, 1, MPI_DOUBLE, 0, world); + MPI_Bcast(&theta_xst7_0_one, 1, MPI_DOUBLE, 0, world); + MPI_Bcast(&dtheta_xst7_ast_one, 1, MPI_DOUBLE, 0, world); + + MPI_Bcast(&a_xst8_one, 1, MPI_DOUBLE, 0, world); + MPI_Bcast(&theta_xst8_0_one, 1, MPI_DOUBLE, 0, world); + MPI_Bcast(&dtheta_xst8_ast_one, 1, MPI_DOUBLE, 0, world); + } b_xst_lo_one = 0.25 * (cut_xst_lo_one - cut_xst_0_one) * (cut_xst_lo_one - cut_xst_0_one)/ diff --git a/src/CG-DNA/pair_oxrna2_dh.cpp b/src/CG-DNA/pair_oxrna2_dh.cpp index eceda90768..1495e1f2ee 100644 --- a/src/CG-DNA/pair_oxrna2_dh.cpp +++ b/src/CG-DNA/pair_oxrna2_dh.cpp @@ -16,6 +16,8 @@ #include "pair_oxrna2_dh.h" +#include "constants_oxdna.h" + using namespace LAMMPS_NS; /* ---------------------------------------------------------------------- @@ -24,7 +26,8 @@ using namespace LAMMPS_NS; 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; + double d_cs_x = ConstantsOxdna::get_d_cs(); + double d_cs_z = ConstantsOxdna::get_d_cs_z(); r[0] = d_cs_x * e1[0] + d_cs_z * e3[0]; r[1] = d_cs_x * e1[1] + d_cs_z * e3[1]; diff --git a/src/CG-DNA/pair_oxrna2_excv.cpp b/src/CG-DNA/pair_oxrna2_excv.cpp index a49497ef5a..7dda4f9434 100644 --- a/src/CG-DNA/pair_oxrna2_excv.cpp +++ b/src/CG-DNA/pair_oxrna2_excv.cpp @@ -17,6 +17,8 @@ #include "pair_oxrna2_excv.h" +#include "constants_oxdna.h" + using namespace LAMMPS_NS; /* ---------------------------------------------------------------------- @@ -25,7 +27,9 @@ using namespace LAMMPS_NS; 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; + double d_cs_x = ConstantsOxdna::get_d_cs(); + double d_cs_z = ConstantsOxdna::get_d_cs_z(); + double d_cb = ConstantsOxdna::get_d_cb(); rs[0] = d_cs_x*e1[0] + d_cs_z*e3[0]; rs[1] = d_cs_x*e1[1] + d_cs_z*e3[1]; diff --git a/src/CG-DNA/pair_oxrna2_stk.cpp b/src/CG-DNA/pair_oxrna2_stk.cpp index 8c5c9c22f5..630b8f3c5b 100644 --- a/src/CG-DNA/pair_oxrna2_stk.cpp +++ b/src/CG-DNA/pair_oxrna2_stk.cpp @@ -19,6 +19,7 @@ #include "atom.h" #include "comm.h" +#include "constants_oxdna.h" #include "error.h" #include "force.h" #include "math_const.h" @@ -27,6 +28,7 @@ #include "mf_oxdna.h" #include "neighbor.h" #include "neigh_list.h" +#include "potential_file_reader.h" #include #include @@ -233,9 +235,12 @@ void PairOxrna2Stk::compute(int eflag, int vflag) double cosphi1,cosphi2,cosphi1dir[3],cosphi2dir[3]; // distances COM-backbone site, COM-3' and COM-5' stacking site - double d_cs_x=-0.4, d_cs_z=+0.2; - double d_cst_x_3p=+0.4, d_cst_y_3p=+0.1; - double d_cst_x_5p=+0.124906078525, d_cst_y_5p=-0.00866274917473; + double d_cs_x = ConstantsOxdna::get_d_cs(); + double d_cs_z = ConstantsOxdna::get_d_cs_z(); + double d_cst_x_3p = ConstantsOxdna::get_d_cst_x_3p(); + double d_cst_y_3p = ConstantsOxdna::get_d_cst_y_3p(); + double d_cst_x_5p = ConstantsOxdna::get_d_cst_x_5p(); + double d_cst_y_5p = ConstantsOxdna::get_d_cst_y_5p(); // 3' and p5' auxiliary vectors double d3p_x=-0.462510,d3p_y=-0.528218,d3p_z=+0.712089; @@ -842,7 +847,7 @@ void PairOxrna2Stk::coeff(int narg, char **arg) { int count; - if (narg != 27) error->all(FLERR,"Incorrect args for pair coefficients in oxrna2/stk"); + if (narg != 7 && narg != 27) error->all(FLERR,"Incorrect args for pair coefficients in oxrna2/stk"); if (!allocated) allocate(); int ilo,ihi,jlo,jhi; @@ -882,30 +887,104 @@ void PairOxrna2Stk::coeff(int narg, char **arg) kappa_st_one = utils::numeric(FLERR,arg[5],false,lmp); epsilon_st_one = stacking_strength(xi_st_one, kappa_st_one, T); - a_st_one = utils::numeric(FLERR,arg[6],false,lmp); - cut_st_0_one = utils::numeric(FLERR,arg[7],false,lmp); - cut_st_c_one = utils::numeric(FLERR,arg[8],false,lmp); - cut_st_lo_one = utils::numeric(FLERR,arg[9],false,lmp); - cut_st_hi_one = utils::numeric(FLERR,arg[10],false,lmp); + if (narg == 27) { + a_st_one = utils::numeric(FLERR,arg[6],false,lmp); + cut_st_0_one = utils::numeric(FLERR,arg[7],false,lmp); + cut_st_c_one = utils::numeric(FLERR,arg[8],false,lmp); + cut_st_lo_one = utils::numeric(FLERR,arg[9],false,lmp); + cut_st_hi_one = utils::numeric(FLERR,arg[10],false,lmp); - a_st5_one = utils::numeric(FLERR,arg[11],false,lmp); - theta_st5_0_one = utils::numeric(FLERR,arg[12],false,lmp); - dtheta_st5_ast_one = utils::numeric(FLERR,arg[13],false,lmp); - a_st6_one = utils::numeric(FLERR,arg[14],false,lmp); - theta_st6_0_one = utils::numeric(FLERR,arg[15],false,lmp); - dtheta_st6_ast_one = utils::numeric(FLERR,arg[16],false,lmp); + a_st5_one = utils::numeric(FLERR,arg[11],false,lmp); + theta_st5_0_one = utils::numeric(FLERR,arg[12],false,lmp); + dtheta_st5_ast_one = utils::numeric(FLERR,arg[13],false,lmp); + a_st6_one = utils::numeric(FLERR,arg[14],false,lmp); + theta_st6_0_one = utils::numeric(FLERR,arg[15],false,lmp); + dtheta_st6_ast_one = utils::numeric(FLERR,arg[16],false,lmp); - a_st9_one = utils::numeric(FLERR,arg[17],false,lmp); - theta_st9_0_one = utils::numeric(FLERR,arg[18],false,lmp); - dtheta_st9_ast_one = utils::numeric(FLERR,arg[19],false,lmp); - a_st10_one = utils::numeric(FLERR,arg[20],false,lmp); - theta_st10_0_one = utils::numeric(FLERR,arg[21],false,lmp); - dtheta_st10_ast_one = utils::numeric(FLERR,arg[22],false,lmp); + a_st9_one = utils::numeric(FLERR,arg[17],false,lmp); + theta_st9_0_one = utils::numeric(FLERR,arg[18],false,lmp); + dtheta_st9_ast_one = utils::numeric(FLERR,arg[19],false,lmp); + a_st10_one = utils::numeric(FLERR,arg[20],false,lmp); + theta_st10_0_one = utils::numeric(FLERR,arg[21],false,lmp); + dtheta_st10_ast_one = utils::numeric(FLERR,arg[22],false,lmp); - a_st1_one = utils::numeric(FLERR,arg[23],false,lmp); - cosphi_st1_ast_one = utils::numeric(FLERR,arg[24],false,lmp); - a_st2_one = utils::numeric(FLERR,arg[25],false,lmp); - cosphi_st2_ast_one = utils::numeric(FLERR,arg[26],false,lmp); + a_st1_one = utils::numeric(FLERR,arg[23],false,lmp); + cosphi_st1_ast_one = utils::numeric(FLERR,arg[24],false,lmp); + a_st2_one = utils::numeric(FLERR,arg[25],false,lmp); + cosphi_st2_ast_one = utils::numeric(FLERR,arg[26],false,lmp); + } else { // read values from potential file + if (comm->me == 0) { + PotentialFileReader reader(lmp, arg[6], "oxdna potential", " (stk)"); + char * line; + std::string iloc, jloc, potential_name; + + while(line = reader.next_line()) { + try { + ValueTokenizer values(line); + iloc = values.next_string(); + jloc = values.next_string(); + potential_name = values.next_string(); + if (iloc == arg[0] && jloc == arg[1] && potential_name == "stk") { + + a_st_one = values.next_double(); + cut_st_0_one = values.next_double(); + cut_st_c_one = values.next_double(); + cut_st_lo_one = values.next_double(); + cut_st_hi_one = values.next_double(); + + a_st5_one = values.next_double(); + theta_st5_0_one = values.next_double(); + dtheta_st5_ast_one = values.next_double(); + a_st6_one = values.next_double(); + theta_st6_0_one = values.next_double(); + dtheta_st6_ast_one = values.next_double(); + + a_st9_one = values.next_double(); + theta_st9_0_one = values.next_double(); + dtheta_st9_ast_one = values.next_double(); + a_st10_one = values.next_double(); + theta_st10_0_one = values.next_double(); + dtheta_st10_ast_one = values.next_double(); + + a_st1_one = values.next_double(); + cosphi_st1_ast_one = values.next_double(); + a_st2_one = values.next_double(); + cosphi_st2_ast_one = values.next_double(); + + break; + } else continue; + } catch (std::exception &e) { + error->one(FLERR, "Problem parsing oxDNA potential file: {}", e.what()); + } + } + if (iloc != arg[0] || jloc != arg[1] || potential_name != "stk") error->one(FLERR, "No corresponding stk potential found in file {} for pair type {} {}", arg[4], arg[0], arg[1]); + } + + MPI_Bcast(&a_st_one, 1, MPI_DOUBLE, 0, world); + MPI_Bcast(&cut_st_0_one, 1, MPI_DOUBLE, 0, world); + MPI_Bcast(&cut_st_c_one, 1, MPI_DOUBLE, 0, world); + MPI_Bcast(&cut_st_lo_one, 1, MPI_DOUBLE, 0, world); + MPI_Bcast(&cut_st_hi_one, 1, MPI_DOUBLE, 0, world); + + MPI_Bcast(&a_st5_one, 1, MPI_DOUBLE, 0, world); + MPI_Bcast(&theta_st5_0_one, 1, MPI_DOUBLE, 0, world); + MPI_Bcast(&dtheta_st5_ast_one, 1, MPI_DOUBLE, 0, world); + MPI_Bcast(&a_st6_one, 1, MPI_DOUBLE, 0, world); + MPI_Bcast(&theta_st6_0_one, 1, MPI_DOUBLE, 0, world); + MPI_Bcast(&dtheta_st6_ast_one, 1, MPI_DOUBLE, 0, world); + + MPI_Bcast(&a_st9_one, 1, MPI_DOUBLE, 0, world); + MPI_Bcast(&theta_st9_0_one, 1, MPI_DOUBLE, 0, world); + MPI_Bcast(&dtheta_st9_ast_one, 1, MPI_DOUBLE, 0, world); + MPI_Bcast(&a_st10_one, 1, MPI_DOUBLE, 0, world); + MPI_Bcast(&theta_st10_0_one, 1, MPI_DOUBLE, 0, world); + MPI_Bcast(&dtheta_st10_ast_one, 1, MPI_DOUBLE, 0, world); + + MPI_Bcast(&a_st1_one, 1, MPI_DOUBLE, 0, world); + MPI_Bcast(&cosphi_st1_ast_one, 1, MPI_DOUBLE, 0, world); + MPI_Bcast(&a_st2_one, 1, MPI_DOUBLE, 0, world); + MPI_Bcast(&cosphi_st2_ast_one, 1, MPI_DOUBLE, 0, world); + } b_st_lo_one = 2*a_st_one*exp(-a_st_one*(cut_st_lo_one-cut_st_0_one))* 2*a_st_one*exp(-a_st_one*(cut_st_lo_one-cut_st_0_one))* diff --git a/src/CG-DNA/pair_oxrna2_xstk.cpp b/src/CG-DNA/pair_oxrna2_xstk.cpp index 8a5fa1b1c9..c5f4218d8f 100644 --- a/src/CG-DNA/pair_oxrna2_xstk.cpp +++ b/src/CG-DNA/pair_oxrna2_xstk.cpp @@ -19,6 +19,7 @@ #include "atom.h" #include "comm.h" +#include "constants_oxdna.h" #include "error.h" #include "force.h" #include "math_const.h" @@ -26,6 +27,7 @@ #include "memory.h" #include "mf_oxdna.h" #include "neigh_list.h" +#include "potential_file_reader.h" #include #include @@ -116,7 +118,7 @@ void PairOxrna2Xstk::compute(int eflag, int vflag) double theta8,theta8p,t8dir[3],cost8; // distance COM-h-bonding site - double d_chb=+0.4; + double d_chb = ConstantsOxdna::get_d_chb(); // vectors COM-h-bonding site in lab frame double ra_chb[3],rb_chb[3]; @@ -581,7 +583,7 @@ void PairOxrna2Xstk::coeff(int narg, char **arg) { int count; - if (narg != 22) error->all(FLERR,"Incorrect args for pair coefficients in oxrna2/xstk"); + if (narg != 3 && narg != 22) error->all(FLERR,"Incorrect args for pair coefficients in oxrna2/xstk"); if (!allocated) allocate(); int ilo,ihi,jlo,jhi; @@ -609,32 +611,106 @@ void PairOxrna2Xstk::coeff(int narg, char **arg) double a_xst8_one, theta_xst8_0_one, dtheta_xst8_ast_one; double b_xst8_one, dtheta_xst8_c_one; - k_xst_one = utils::numeric(FLERR,arg[2],false,lmp); - cut_xst_0_one = utils::numeric(FLERR,arg[3],false,lmp); - cut_xst_c_one = utils::numeric(FLERR,arg[4],false,lmp); - cut_xst_lo_one = utils::numeric(FLERR,arg[5],false,lmp); - cut_xst_hi_one = utils::numeric(FLERR,arg[6],false,lmp); + if (narg == 22) { + k_xst_one = utils::numeric(FLERR,arg[2],false,lmp); + cut_xst_0_one = utils::numeric(FLERR,arg[3],false,lmp); + cut_xst_c_one = utils::numeric(FLERR,arg[4],false,lmp); + cut_xst_lo_one = utils::numeric(FLERR,arg[5],false,lmp); + cut_xst_hi_one = utils::numeric(FLERR,arg[6],false,lmp); - a_xst1_one = utils::numeric(FLERR,arg[7],false,lmp); - theta_xst1_0_one = utils::numeric(FLERR,arg[8],false,lmp); - dtheta_xst1_ast_one = utils::numeric(FLERR,arg[9],false,lmp); + a_xst1_one = utils::numeric(FLERR,arg[7],false,lmp); + theta_xst1_0_one = utils::numeric(FLERR,arg[8],false,lmp); + dtheta_xst1_ast_one = utils::numeric(FLERR,arg[9],false,lmp); - a_xst2_one = utils::numeric(FLERR,arg[10],false,lmp); - theta_xst2_0_one = utils::numeric(FLERR,arg[11],false,lmp); - dtheta_xst2_ast_one = utils::numeric(FLERR,arg[12],false,lmp); + a_xst2_one = utils::numeric(FLERR,arg[10],false,lmp); + theta_xst2_0_one = utils::numeric(FLERR,arg[11],false,lmp); + dtheta_xst2_ast_one = utils::numeric(FLERR,arg[12],false,lmp); - a_xst3_one = utils::numeric(FLERR,arg[13],false,lmp); - theta_xst3_0_one = utils::numeric(FLERR,arg[14],false,lmp); - dtheta_xst3_ast_one = utils::numeric(FLERR,arg[15],false,lmp); + a_xst3_one = utils::numeric(FLERR,arg[13],false,lmp); + theta_xst3_0_one = utils::numeric(FLERR,arg[14],false,lmp); + dtheta_xst3_ast_one = utils::numeric(FLERR,arg[15],false,lmp); - a_xst7_one = utils::numeric(FLERR,arg[16],false,lmp); - theta_xst7_0_one = utils::numeric(FLERR,arg[17],false,lmp); - dtheta_xst7_ast_one = utils::numeric(FLERR,arg[18],false,lmp); + a_xst7_one = utils::numeric(FLERR,arg[16],false,lmp); + theta_xst7_0_one = utils::numeric(FLERR,arg[17],false,lmp); + dtheta_xst7_ast_one = utils::numeric(FLERR,arg[18],false,lmp); - a_xst8_one = utils::numeric(FLERR,arg[19],false,lmp); - theta_xst8_0_one = utils::numeric(FLERR,arg[20],false,lmp); - dtheta_xst8_ast_one = utils::numeric(FLERR,arg[21],false,lmp); + a_xst8_one = utils::numeric(FLERR,arg[19],false,lmp); + theta_xst8_0_one = utils::numeric(FLERR,arg[20],false,lmp); + dtheta_xst8_ast_one = utils::numeric(FLERR,arg[21],false,lmp); + } else { + if (comm->me == 0) { + PotentialFileReader reader(lmp, arg[2], "oxdna potential", " (xstk)"); + char * line; + std::string iloc, jloc, potential_name; + while(line = reader.next_line()) { + try { + ValueTokenizer values(line); + iloc = values.next_string(); + jloc = values.next_string(); + potential_name = values.next_string(); + if (iloc == arg[0] && jloc == arg[1] && potential_name == "xstk") { + k_xst_one = values.next_double(); + cut_xst_0_one = values.next_double(); + cut_xst_c_one = values.next_double(); + cut_xst_lo_one = values.next_double(); + cut_xst_hi_one = values.next_double(); + + a_xst1_one = values.next_double(); + theta_xst1_0_one = values.next_double(); + dtheta_xst1_ast_one = values.next_double(); + + a_xst2_one = values.next_double(); + theta_xst2_0_one = values.next_double(); + dtheta_xst2_ast_one = values.next_double(); + + a_xst3_one = values.next_double(); + theta_xst3_0_one = values.next_double(); + dtheta_xst3_ast_one = values.next_double(); + + a_xst7_one = values.next_double(); + theta_xst7_0_one = values.next_double(); + dtheta_xst7_ast_one = values.next_double(); + + a_xst8_one = values.next_double(); + theta_xst8_0_one = values.next_double(); + dtheta_xst8_ast_one = values.next_double(); + + break; + } else continue; + } catch (std::exception &e) { + error->one(FLERR, "Problem parsing oxDNA potential file: {}", e.what()); + } + } + if (iloc != arg[0] || jloc != arg[1] || potential_name != "xstk") error->one(FLERR, "No corresponding xstk potential found in file {} for pair type {} {}", arg[2], arg[0], arg[1]); + } + + MPI_Bcast(&k_xst_one, 1, MPI_DOUBLE, 0, world); + MPI_Bcast(&cut_xst_0_one, 1, MPI_DOUBLE, 0, world); + MPI_Bcast(&cut_xst_c_one, 1, MPI_DOUBLE, 0, world); + MPI_Bcast(&cut_xst_lo_one, 1, MPI_DOUBLE, 0, world); + MPI_Bcast(&cut_xst_hi_one, 1, MPI_DOUBLE, 0, world); + + MPI_Bcast(&a_xst1_one, 1, MPI_DOUBLE, 0, world); + MPI_Bcast(&theta_xst1_0_one, 1, MPI_DOUBLE, 0, world); + MPI_Bcast(&dtheta_xst1_ast_one, 1, MPI_DOUBLE, 0, world); + + MPI_Bcast(&a_xst2_one, 1, MPI_DOUBLE, 0, world); + MPI_Bcast(&theta_xst2_0_one, 1, MPI_DOUBLE, 0, world); + MPI_Bcast(&dtheta_xst2_ast_one, 1, MPI_DOUBLE, 0, world); + + MPI_Bcast(&a_xst3_one, 1, MPI_DOUBLE, 0, world); + MPI_Bcast(&theta_xst3_0_one, 1, MPI_DOUBLE, 0, world); + MPI_Bcast(&dtheta_xst3_ast_one, 1, MPI_DOUBLE, 0, world); + + MPI_Bcast(&a_xst7_one, 1, MPI_DOUBLE, 0, world); + MPI_Bcast(&theta_xst7_0_one, 1, MPI_DOUBLE, 0, world); + MPI_Bcast(&dtheta_xst7_ast_one, 1, MPI_DOUBLE, 0, world); + + MPI_Bcast(&a_xst8_one, 1, MPI_DOUBLE, 0, world); + MPI_Bcast(&theta_xst8_0_one, 1, MPI_DOUBLE, 0, world); + MPI_Bcast(&dtheta_xst8_ast_one, 1, MPI_DOUBLE, 0, world); + } b_xst_lo_one = 0.25 * (cut_xst_lo_one - cut_xst_0_one) * (cut_xst_lo_one - cut_xst_0_one)/ (0.5 * (cut_xst_lo_one - cut_xst_0_one) * (cut_xst_lo_one - cut_xst_0_one) -