Compare commits

...

30 Commits

Author SHA1 Message Date
4b9d0a9566 sync with SVN 2017-03-28 13:26:12 -06:00
0637f23875 patch 28Mar17 2017-03-28 13:12:23 -06:00
9f6e126a2f Merge pull request #437 from ohenrich/user-cgdna
User cgdna
2017-03-28 12:52:26 -06:00
645f56cf70 Merge pull request #436 from Pakketeretet2/better_incorrect_input_handling_nh
Changed the check on initial and final temperature to <= 0 for both.
2017-03-28 12:51:17 -06:00
80e5111dca Merge pull request #434 from akohlmey/imgflags-in-library
improved image flag handling in library interface
2017-03-28 12:50:13 -06:00
7e9f05b617 Merge pull request #433 from akohlmey/fixes-for-stable
More small fixes for stable release
2017-03-28 12:49:09 -06:00
1d8f0c762d Merge branch 'master' into fixes-for-stable 2017-03-28 14:37:30 -04:00
ef6070cbde remove executable permissions for potential files 2017-03-28 14:35:58 -04:00
61f3ff1d2b Merge branch 'master' of github.com:lammps/lammps 2017-03-28 12:35:33 -06:00
111d350a22 fix gcmc units change for chemical potential 2017-03-28 12:34:46 -06:00
1dfd61f532 Merge pull request #432 from Pakketeretet2/user_manifold_fix
Fixed a bug with equal-style variables as manifold params.
2017-03-28 12:33:45 -06:00
5c1f5462e7 Removed contribution line from header files 2017-03-28 19:08:24 +01:00
66a6375405 Resolved merge conflict 2017-03-28 18:58:31 +01:00
604afebf6f Update to oxDNA2 2017-03-28 18:22:02 +01:00
8afed61db1 Upgrade to oxDNA2 2017-03-28 18:16:36 +01:00
ee55a98103 Changed the check on initial and final temperature to <= 0 for both. 2017-03-28 11:22:10 -04:00
f8da9a866a synchronize dump custom/vtk documentation with that of dump custom 2017-03-28 11:00:22 -04:00
28bdebd3c0 avoid segfault when calling PPPM*::memory_usage() before grid communication is initialized 2017-03-28 07:50:48 -04:00
fc51c38abb add some docs for the special treatment of image flags 2017-03-28 02:22:45 -04:00
443ea13eff add image flag packing/unpacking to library/python interface 2017-03-28 02:05:05 -04:00
5feeb79c13 one more line of dead code removed 2017-03-27 15:16:28 -04:00
a241b2d0f7 fix problems with references 2017-03-27 15:01:32 -04:00
61e7595a94 remove references to xmovie, streamline dump and viz descriptions 2017-03-27 14:59:58 -04:00
da9096750e update .gitignore for newly added files 2017-03-27 14:30:21 -04:00
87ea9ba661 bugfix for library interface 2017-03-27 14:29:13 -04:00
c041727e4f remove dead code and reduce trivial compiler warnings (clang++) 2017-03-27 14:28:50 -04:00
3feffbe1de Removed diagnostics. 2017-03-27 13:49:53 -04:00
04fd038d35 Fixed a bug with equal-style variables as manifold params. 2017-03-27 13:46:57 -04:00
af0b5b0e84 Removed dead code 2017-03-22 16:23:29 +00:00
7435084375 Verified oxDNA with modified nucleotide layout 2017-03-22 15:59:10 +00:00
118 changed files with 12921 additions and 5509 deletions

BIN
doc/src/Eqs/fix_gcmc1.jpg Normal file

Binary file not shown.

After

Width:  |  Height:  |  Size: 5.5 KiB

View File

@ -0,0 +1,9 @@
\documentclass[12pt]{article}
\begin{document}
\begin{eqnarray*}
\mu &=&\mu^{id} + \mu^{ex}
\end{eqnarray*}
\end{document}

BIN
doc/src/Eqs/fix_gcmc2.jpg Normal file

Binary file not shown.

After

Width:  |  Height:  |  Size: 10 KiB

10
doc/src/Eqs/fix_gcmc2.tex Normal file
View File

@ -0,0 +1,10 @@
\documentclass[12pt]{article}
\begin{document}
\begin{eqnarray*}
\mu^{id} &=& k T \ln{\rho \Lambda^3} \\
&=& k T \ln{\frac{\phi P \Lambda^3}{k T}}
\end{eqnarray*}
\end{document}

BIN
doc/src/Eqs/fix_gcmc3.jpg Normal file

Binary file not shown.

After

Width:  |  Height:  |  Size: 7.3 KiB

View File

@ -0,0 +1,9 @@
\documentclass[12pt]{article}
\begin{document}
\begin{eqnarray*}
\Lambda &=& \sqrt{ \frac{h^2}{2 \pi m k T}}
\end{eqnarray*}
\end{document}

View File

@ -1,7 +1,7 @@
<!-- HTML_ONLY -->
<HEAD>
<TITLE>LAMMPS Users Manual</TITLE>
<META NAME="docnumber" CONTENT="24 Mar 2017 version">
<META NAME="docnumber" CONTENT="28 Mar 2017 version">
<META NAME="author" CONTENT="http://lammps.sandia.gov - Sandia National Laboratories">
<META NAME="copyright" CONTENT="Copyright (2003) Sandia Corporation. This software and manual is distributed under the GNU General Public License.">
</HEAD>
@ -21,7 +21,7 @@
<H1></H1>
LAMMPS Documentation :c,h3
24 Mar 2017 version :c,h4
28 Mar 2017 version :c,h4
Version info: :h4

Binary file not shown.

View File

@ -1047,6 +1047,10 @@ package"_Section_start.html#start_3.
"oxdna/hbond"_pair_oxdna.html,
"oxdna/stk"_pair_oxdna.html,
"oxdna/xstk"_pair_oxdna.html,
"oxdna2/coaxstk"_pair_oxdna2.html,
"oxdna2/dh"_pair_oxdna2.html,
"oxdna2/excv"_pair_oxdna2.html,
"oxdna2/stk"_pair_oxdna2.html,
"quip"_pair_quip.html,
"reax/c (k)"_pair_reax_c.html,
"smd/hertz"_pair_smd_hertz.html,
@ -1096,7 +1100,8 @@ package"_Section_start.html#start_3.
"harmonic/shift (o)"_bond_harmonic_shift.html,
"harmonic/shift/cut (o)"_bond_harmonic_shift_cut.html,
"oxdna/fene"_bond_oxdna.html :tb(c=4,ea=c)
"oxdna/fene"_bond_oxdna.html,
"oxdna2/fene"_bond_oxdna.html :tb(c=4,ea=c)
:line

View File

@ -759,23 +759,14 @@ LAMMPS itself does not do visualization, but snapshots from LAMMPS
simulations can be visualized (and analyzed) in a variety of ways.
LAMMPS snapshots are created by the "dump"_dump.html command which can
create files in several formats. The native LAMMPS dump format is a
create files in several formats. The native LAMMPS dump format is a
text file (see "dump atom" or "dump custom") which can be visualized
by the "xmovie"_Section_tools.html#xmovie program, included with the
LAMMPS package. This produces simple, fast 2d projections of 3d
systems, and can be useful for rapid debugging of simulation geometry
and atom trajectories.
by several popular visualization tools. The "dump image"_dump_image.html
and "dump movie"_dump_image.html styles can output internally rendered
images and convert a sequence of them to a movie during the MD run.
Several programs included with LAMMPS as auxiliary tools can convert
native LAMMPS dump files to other formats. See the
"Section 9"_Section_tools.html doc page for details. The first is
the "ch2lmp tool"_Section_tools.html#charmm, which contains a
lammps2pdb Perl script which converts LAMMPS dump files into PDB
files. The second is the "lmp2arc tool"_Section_tools.html#arc which
converts LAMMPS dump files into Accelrys' Insight MD program files.
The third is the "lmp2cfg tool"_Section_tools.html#cfg which converts
LAMMPS dump files into CFG files which can be read into the
"AtomEye"_atomeye visualizer.
between LAMMPS format files and other formats.
See the "Section 9"_Section_tools.html doc page for details.
A Python-based toolkit distributed by our group can read native LAMMPS
dump files, including custom dump files with additional columns of
@ -788,22 +779,7 @@ RasMol visualization programs. Pizza.py has tools that do interactive
3d OpenGL visualization and one that creates SVG images of dump file
snapshots.
LAMMPS can create XYZ files directly (via "dump xyz") which is a
simple text-based file format used by many visualization programs
including "VMD"_vmd.
LAMMPS can create DCD files directly (via "dump dcd") which can be
read by "VMD"_vmd in conjunction with a CHARMM PSF file. Using this
form of output avoids the need to convert LAMMPS snapshots to PDB
files. See the "dump"_dump.html command for more information on DCD
files.
LAMMPS can create XTC files directly (via "dump xtc") which is GROMACS
file format which can also be read by "VMD"_vmd for visualization.
See the "dump"_dump.html command for more information on XTC files.
:link(pizza,http://www.sandia.gov/~sjplimp/pizza.html)
:link(vmd,http://www.ks.uiuc.edu/Research/vmd)
:link(ensight,http://www.ensight.com)
:link(atomeye,http://mt.seas.upenn.edu/Archive/Graphics/A)
@ -2013,6 +1989,11 @@ Both methods are thus a means to extract or assign (overwrite) any
peratom quantities within LAMMPS. See the extract() method in the
src/atom.cpp file for a list of valid per-atom properties. New names
could easily be added if the property you want is not listed.
A special treatment is applied for accessing image flags via the
"image" property. Image flags are stored in a packed format with all
three image flags stored in a single integer. When signaling to access
the image flags as 3 individual values per atom instead of 1, the data
is transparently packed or unpacked by the library interface.
The lammps_create_atoms() function takes a list of N atoms as input
with atom types and coords (required), an optionally atom IDs and

View File

@ -338,15 +338,13 @@ dynamics timestepping, particularly if the computations are not
parallel, so it is often better to leave such analysis to
post-processing codes.
A very simple (yet fast) visualizer is provided with the LAMMPS
package - see the "xmovie"_Section_tools.html#xmovie tool in "this
section"_Section_tools.html. It creates xyz projection views of
atomic coordinates and animates them. We find it very useful for
debugging purposes. For high-quality visualization we recommend the
For high-quality visualization we recommend the
following packages:
"VMD"_http://www.ks.uiuc.edu/Research/vmd
"AtomEye"_http://mt.seas.upenn.edu/Archive/Graphics/A
"OVITO"_http://www.ovito.org/
"ParaView"_http://www.paraview.org/
"PyMol"_http://www.pymol.org
"Raster3d"_http://www.bmsc.washington.edu/raster3d/raster3d.html
"RasMol"_http://www.openrasmol.org :ul

View File

@ -1140,7 +1140,7 @@ Package, Description, Author(s), Doc page, Example, Pic/movie, Library
"USER-ATC"_#USER-ATC, atom-to-continuum coupling, Jones & Templeton & Zimmerman (1), "fix atc"_fix_atc.html, USER/atc, "atc"_atc, lib/atc
"USER-AWPMD"_#USER-AWPMD, wave-packet MD, Ilya Valuev (JIHT), "pair_style awpmd/cut"_pair_awpmd.html, USER/awpmd, -, lib/awpmd
"USER-CG-CMM"_#USER-CG-CMM, coarse-graining model, Axel Kohlmeyer (Temple U), "pair_style lj/sdk"_pair_sdk.html, USER/cg-cmm, "cg"_cg, -
"USER-CGDNA"_#USER-CGDNA, coarse-grained DNA force fields, Oliver Henrich (U Edinburgh), src/USER-CGDNA/README, USER/cgdna, -, -
"USER-CGDNA"_#USER-CGDNA, coarse-grained DNA force fields, Oliver Henrich (U Strathclyde Glasgow), src/USER-CGDNA/README, USER/cgdna, -, -
"USER-COLVARS"_#USER-COLVARS, collective variables, Fiorin & Henin & Kohlmeyer (2), "fix colvars"_fix_colvars.html, USER/colvars, "colvars"_colvars, lib/colvars
"USER-DIFFRACTION"_#USER-DIFFRACTION, virutal x-ray and electron diffraction, Shawn Coleman (ARL),"compute xrd"_compute_xrd.html, USER/diffraction, -, -
"USER-DPD"_#USER-DPD, reactive dissipative particle dynamics (DPD), Larentzos & Mattox & Brennan (5), src/USER-DPD/README, USER/dpd, -, -
@ -1288,25 +1288,29 @@ him directly if you have questions.
USER-CGDNA package :link(USER-CGDNA),h5
Contents: The CGDNA package implements coarse-grained force fields for
single- and double-stranded DNA. This is at the moment mainly the
oxDNA model, developed by Doye, Louis and Ouldridge at the University
single- and double-stranded DNA. These are at the moment mainly the
oxDNA and oxDNA2 models, developed by Doye, Louis and Ouldridge at the University
of Oxford. The package also contains Langevin-type rigid-body
integrators with improved stability.
See these doc pages to get started:
"bond_style oxdna/fene"_bond_oxdna.html
"bond_style oxdna2/fene"_bond_oxdna.html
"pair_style oxdna/..."_pair_oxdna.html
"pair_style oxdna2/..."_pair_oxdna2.html
"fix nve/dotc/langevin"_fix_nve_dotc_langevin.html :ul
Supporting info: /src/USER-CGDNA/README, "bond_style
oxdna/fene"_bond_oxdna.html, "pair_style
oxdna/..."_pair_oxdna.html, "fix
oxdna/fene"_bond_oxdna.html, "bond_style
oxdna2/fene"_bond_oxdna.html, "pair_style
oxdna/..."_pair_oxdna.html, "pair_style
oxdna2/..."_pair_oxdna2.html, "fix
nve/dotc/langevin"_fix_nve_dotc_langevin.html
Author: Oliver Henrich at the University of Edinburgh, UK (o.henrich
at epcc.ed.ac.uk or ohenrich at ph.ed.ac.uk). Contact him directly if
you have any questions.
Author: Oliver Henrich at the University of Strathclyde, Glasgow, UK and
University of Edinburgh (ohenrich@ph.ed.ac.uk).
Contact him directly if you have any questions.
:line

View File

@ -698,7 +698,12 @@ method in the src/atom.cpp file for a list of valid names. Again, new
names could easily be added if the property you want is missing. The
vector can be used via normal Python subscripting. If atom IDs are
not consecutively ordered within LAMMPS, a None is returned as
indication of an error.
indication of an error. A special treatment is applied for image flags
stored in the "image" property. All three image flags are stored in
a packed format in a single integer, so count would be 1 to retrieve
that integer, however also a count value of 3 can be used and then
the image flags will be unpacked into 3 individual integers, ordered
in a similar fashion as coordinates.
Note that the data structure gather_atoms("x") returns is different
from the data structure returned by extract_atom("x") in four ways.
@ -727,6 +732,10 @@ corresponding properties for each atom inside LAMMPS. This requires
LAMMPS to have its "map" option enabled; see the
"atom_modify"_atom_modify.html command for details. If it is not, or
if atom IDs are not consecutively ordered, no coordinates are reset.
Similar as for gather_atoms() a special treatment is applied for image
flags, which can be provided in packed (count = 1) or unpacked (count = 3)
format and in the latter case, they will be packed before applied to
atoms.
The array of coordinates passed to scatter_atoms() must be a ctypes
vector of ints or doubles, allocated and initialized something like

View File

@ -7,19 +7,24 @@
:line
bond_style oxdna/fene command :h3
bond_style oxdna2/fene command :h3
[Syntax:]
bond_style oxdna/fene :pre
bond_style oxdna2/fene :pre
[Examples:]
bond_style oxdna/fene
bond_coeff * 2.0 0.25 0.7525 :pre
bond_style oxdna2/fene
bond_coeff * 2.0 0.25 0.7564 :pre
[Description:]
The {oxdna/fene} bond style uses the potential
The {oxdna/fene} and {oxdna2/fene} bond styles use the potential
:c,image(Eqs/bond_oxdna_fene.jpg)
@ -36,13 +41,16 @@ epsilon (energy)
Delta (distance)
r0 (distance) :ul
NOTE: This bond style has to be used together with the corresponding oxDNA pair styles
NOTE: The oxDNA bond style has to be used together with the corresponding oxDNA pair styles
for excluded volume interaction {oxdna/excv}, stacking {oxdna/stk}, cross-stacking {oxdna/xstk}
and coaxial stacking interaction {oxdna/coaxstk} as well as hydrogen-bonding interaction {oxdna/hbond} (see also documentation of
"pair_style oxdna/excv"_pair_oxdna.html). The coefficients
"pair_style oxdna/excv"_pair_oxdna.html). For the oxDNA2 "(Snodin)"_#oxdna2 bond style the analogous pair styles and an additional Debye-Hueckel pair
style {oxdna2/dh} have to be defined.
The coefficients
in the above example have to be kept fixed and cannot be changed without reparametrizing the entire model.
Example input and data files can be found in examples/USER/cgdna/examples/duplex1/ and /duplex2/.
Example input and data files for DNA duplexes can be found in examples/USER/cgdna/examples/oxDNA/ and /oxDNA2/.
A simple python setup tool which creates single straight or helical DNA strands,
DNA duplexes or arrays of DNA duplexes can be found in examples/USER/cgdna/util/.
A technical report with more information on the model, the structure of the input file,
@ -60,7 +68,7 @@ LAMMPS"_Section_start.html#start_3 section for more info on packages.
[Related commands:]
"pair_style oxdna/excv"_pair_oxdna.html, "fix nve/dotc/langevin"_fix_nve_dotc_langevin.html, "bond_coeff"_bond_coeff.html
"pair_style oxdna/excv"_pair_oxdna.html, "pair_style oxdna2/excv"_pair_oxdna2.html, "fix nve/dotc/langevin"_fix_nve_dotc_langevin.html, "bond_coeff"_bond_coeff.html
[Default:] none
@ -68,3 +76,6 @@ LAMMPS"_Section_start.html#start_3 section for more info on packages.
:link(oxdna_fene)
[(Ouldridge)] T.E. Ouldridge, A.A. Louis, J.P.K. Doye, J. Chem. Phys. 134, 085101 (2011).
:link(oxdna2)
[(Snodin)] B.E. Snodin, F. Randisi, M. Mosayebi, et al., J. Chem. Phys. 142, 234901 (2015).

View File

@ -16,6 +16,7 @@ Bond Styles :h1
bond_none
bond_nonlinear
bond_oxdna
bond_oxdna2
bond_quartic
bond_table
bond_zero

View File

@ -54,7 +54,7 @@ adding atoms or molecules to the system (see the "fix
pour"_fix_pour.html, "fix deposit"_fix_deposit.html, and "fix
gcmc"_fix_gcmc.html commands) or expect atoms or molecules to be lost
(e.g. due to exiting the simulation box or via "fix
evaporation"_fix_evaporation.html), then this option should be used to
evaporate"_fix_evaporate.html), then this option should be used to
insure the temperature is correctly normalized.
NOTE: The {extra} and {dynamic} keywords should not be used as they

View File

@ -40,8 +40,11 @@ field.
NOTE: The newer {charmmfsh} style was released in March 2017. We
recommend it be used instead of the older {charmm} style when running
a simulation with the CHARMM force field. See the discussion below
and more details on the "pair_style charmm"_pair_charmm.html doc page.
a simulation with the CHARMM force field and Coulomb cutoffs, via the
"pair_style lj/charmmfsw/coul/charmmfsh"_pair_charmm.html command.
Otherwise the older {charmm} style is fine to use. See the discussion
below and more details on the "pair_style charmm"_pair_charmm.html doc
page.
The following coefficients must be defined for each dihedral type via the
"dihedral_coeff"_dihedral_coeff.html command as in the example above, or in
@ -82,13 +85,18 @@ special_bonds 1-4 scaling factor to 0.0 (which is the
default). Otherwise 1-4 non-bonded interactions in dihedrals will be
computed twice.
For simulations using the CHARMM force field, the difference between
the {charmm} and {charmmfsh} styles is in the computation of the 1-4
non-bond interactions, if the distance between the two atoms is within
the switching distance of the pairwise potential defined by the
corresponding CHARMM pair style, i.e. between the inner and outer
cutoffs specified for the pair style. See the discussion on the
"CHARMM pair_style"_pair_charmm.html doc page for details.
For simulations using the CHARMM force field with a Coulomb cutoff,
the difference between the {charmm} and {charmmfsh} styles is in the
computation of the 1-4 non-bond interactions, though only if the
distance between the two atoms is within the switching region of the
pairwise potential defined by the corresponding CHARMM pair style,
i.e. between the inner and outer cutoffs specified for the pair style.
The {charmmfsh} style should only be used when using the "pair_style
lj/charmmfsw/coul/charmmfsh"_pair_charmm.html to make the Coulombic
pairwise calculations consistent. Use the {charmm} style with
long-range Coulombics or the older "pair_style
lj/charmm/coul/charmm"_pair_charmm.html command. See the discussion
on the "CHARMM pair_style"_pair_charmm.html doc page for details.
Note that for AMBER force fields, which use pair styles with "lj/cut",
the special_bonds 1-4 scaling factor should be set to the AMBER

View File

@ -331,10 +331,7 @@ bonds and colors.
Note that {atom}, {custom}, {dcd}, {xtc}, and {xyz} style dump files
can be read directly by "VMD"_http://www.ks.uiuc.edu/Research/vmd, a
popular molecular viewing program. See
"Section 9"_Section_tools.html#vmd of the manual and the
tools/lmp2vmd/README.txt file for more information about support in
VMD for reading and visualizing LAMMPS dump files.
popular molecular viewing program.
:line

View File

@ -26,7 +26,6 @@ args = list of arguments for a particular style :l
q, mux, muy, muz, mu,
radius, diameter, omegax, omegay, omegaz,
angmomx, angmomy, angmomz, tqx, tqy, tqz,
spin, eradius, ervel, erforce,
c_ID, c_ID\[N\], f_ID, f_ID\[N\], v_name :pre
id = atom ID
@ -51,17 +50,18 @@ args = list of arguments for a particular style :l
angmomx,angmomy,angmomz = angular momentum of aspherical particle
tqx,tqy,tqz = torque on finite-size particles
c_ID = per-atom vector calculated by a compute with ID
c_ID\[N\] = Nth column of per-atom array calculated by a compute with ID
c_ID\[I\] = Ith column of per-atom array calculated by a compute with ID, I can include wildcard (see below)
f_ID = per-atom vector calculated by a fix with ID
f_ID\[N\] = Nth column of per-atom array calculated by a fix with ID
v_name = per-atom vector calculated by an atom-style variable with name :pre
f_ID\[I\] = Ith column of per-atom array calculated by a fix with ID, I can include wildcard (see below)
v_name = per-atom vector calculated by an atom-style variable with name
d_name = per-atom floating point vector with name, managed by fix property/atom
i_name = per-atom integer vector with name, managed by fix property/atom :pre
:ule
[Examples:]
dump dmpvtk all custom/vtk 100 dump*.myforce.vtk id type vx fx
dump dmpvtp flow custom/vtk 100 dump*.%.displace.vtp id type c_myD\[1\] c_myD\[2\] c_myD\[3\] v_ke
dump e_data all custom/vtk 100 dump*.vtu id type spin eradius fx fy fz eforce :pre
dump dmpvtp flow custom/vtk 100 dump*.%.displace.vtp id type c_myD\[1\] c_myD\[2\] c_myD\[3\] v_ke :pre
The style {custom/vtk} is similar to the "custom"_dump.html style but
uses the VTK library to write data to VTK simple legacy or XML format
@ -199,32 +199,38 @@ part of the {custom/vtk} style.
The {id}, {mol}, {proc}, {procp1}, {type}, {element}, {mass}, {vx},
{vy}, {vz}, {fx}, {fy}, {fz}, {q} attributes are self-explanatory.
{id} is the atom ID. {mol} is the molecule ID, included in the data
file for molecular systems. {type} is the atom type. {element} is
typically the chemical name of an element, which you must assign to
each type via the "dump_modify element"_dump_modify.html command.
More generally, it can be any string you wish to associate with an
atom type. {mass} is the atom mass. {vx}, {vy}, {vz}, {fx}, {fy},
{fz}, and {q} are components of atom velocity and force and atomic
charge.
{Id} is the atom ID. {Mol} is the molecule ID, included in the data
file for molecular systems. {Proc} is the ID of the processor (0 to
Nprocs-1) that currently owns the atom. {Procp1} is the proc ID+1,
which can be convenient in place of a {type} attribute (1 to Ntypes)
for coloring atoms in a visualization program. {Type} is the atom
type (1 to Ntypes). {Element} is typically the chemical name of an
element, which you must assign to each type via the "dump_modify
element"_dump_modify.html command. More generally, it can be any
string you wish to associated with an atom type. {Mass} is the atom
mass. {Vx}, {vy}, {vz}, {fx}, {fy}, {fz}, and {q} are components of
atom velocity and force and atomic charge.
There are several options for outputting atom coordinates. The {x},
{y}, {z} attributes are used to write atom coordinates "unscaled", in
the appropriate distance "units"_units.html (Angstroms, sigma, etc).
Additionally, you can use {xs}, {ys}, {zs} if you want to also save the
coordinates "scaled" to the box size, so that each value is 0.0 to
1.0. If the simulation box is triclinic (tilted), then all atom
coords will still be between 0.0 and 1.0. Use {xu}, {yu}, {zu} if you
want the coordinates "unwrapped" by the image flags for each atom.
Unwrapped means that if the atom has passed through a periodic
boundary one or more times, the value is printed for what the
coordinate would be if it had not been wrapped back into the periodic
box. Note that using {xu}, {yu}, {zu} means that the coordinate
values may be far outside the box bounds printed with the snapshot.
Using {xsu}, {ysu}, {zsu} is similar to using {xu}, {yu}, {zu}, except
that the unwrapped coordinates are scaled by the box size. Atoms that
have passed through a periodic boundary will have the corresponding
coordinate increased or decreased by 1.0.
{y}, {z} attributes write atom coordinates "unscaled", in the
appropriate distance "units"_units.html (Angstroms, sigma, etc). Use
{xs}, {ys}, {zs} if you want the coordinates "scaled" to the box size,
so that each value is 0.0 to 1.0. If the simulation box is triclinic
(tilted), then all atom coords will still be between 0.0 and 1.0.
I.e. actual unscaled (x,y,z) = xs*A + ys*B + zs*C, where (A,B,C) are
the non-orthogonal vectors of the simulation box edges, as discussed
in "Section 6.12"_Section_howto.html#howto_12.
Use {xu}, {yu}, {zu} if you want the coordinates "unwrapped" by the
image flags for each atom. Unwrapped means that if the atom has
passed thru a periodic boundary one or more times, the value is
printed for what the coordinate would be if it had not been wrapped
back into the periodic box. Note that using {xu}, {yu}, {zu} means
that the coordinate values may be far outside the box bounds printed
with the snapshot. Using {xsu}, {ysu}, {zsu} is similar to using
{xu}, {yu}, {zu}, except that the unwrapped coordinates are scaled by
the box size. Atoms that have passed through a periodic boundary will
have the corresponding coordinate increased or decreased by 1.0.
The image flags can be printed directly using the {ix}, {iy}, {iz}
attributes. For periodic dimensions, they specify which image of the
@ -255,13 +261,7 @@ The {tqx}, {tqy}, {tqz} attributes are for finite-size particles that
can sustain a rotational torque due to interactions with other
particles.
The {spin}, {eradius}, {ervel}, and {erforce} attributes are for
particles that represent nuclei and electrons modeled with the
electronic force field (EFF). See "atom_style
electron"_atom_style.html and "pair_style eff"_pair_eff.html for more
details.
The {c_ID} and {c_ID\[N\]} attributes allow per-atom vectors or arrays
The {c_ID} and {c_ID\[I\]} attributes allow per-atom vectors or arrays
calculated by a "compute"_compute.html to be output. The ID in the
attribute should be replaced by the actual ID of the compute that has
been defined previously in the input script. See the
@ -275,12 +275,14 @@ command. Instead, global quantities can be output by the
"thermo_style custom"_thermo_style.html command, and local quantities
can be output by the dump local command.
If {c_ID} is used as an attribute, then the per-atom vector calculated
by the compute is printed. If {c_ID\[N\]} is used, then N must be in
the range from 1-M, which will print the Nth column of the M-length
per-atom array calculated by the compute.
If {c_ID} is used as a attribute, then the per-atom vector calculated
by the compute is printed. If {c_ID\[I\]} is used, then I must be in
the range from 1-M, which will print the Ith column of the per-atom
array with M columns calculated by the compute. See the discussion
above for how I can be specified with a wildcard asterisk to
effectively specify multiple values.
The {f_ID} and {f_ID\[N\]} attributes allow vector or array per-atom
The {f_ID} and {f_ID\[I\]} attributes allow vector or array per-atom
quantities calculated by a "fix"_fix.html to be output. The ID in the
attribute should be replaced by the actual ID of the fix that has been
defined previously in the input script. The "fix
@ -291,9 +293,11 @@ any "compute"_compute.html, "fix"_fix.html, or atom-style
be written to a dump file.
If {f_ID} is used as a attribute, then the per-atom vector calculated
by the fix is printed. If {f_ID\[N\]} is used, then N must be in the
range from 1-M, which will print the Nth column of the M-length
per-atom array calculated by the fix.
by the fix is printed. If {f_ID\[I\]} is used, then I must be in the
range from 1-M, which will print the Ith column of the per-atom array
with M columns calculated by the fix. See the discussion above for
how I can be specified with a wildcard asterisk to effectively specify
multiple values.
The {v_name} attribute allows per-atom vectors calculated by a
"variable"_variable.html to be output. The name in the attribute
@ -306,6 +310,10 @@ invoke other computes, fixes, or variables when they are evaluated, so
this is a very general means of creating quantities to output to a
dump file.
The {d_name} and {i_name} attributes allow to output custom per atom
floating point or integer properties that are managed by
"fix property/atom"_fix_property_atom.html.
See "Section 10"_Section_modify.html of the manual for information
on how to add new compute and fix styles to LAMMPS to calculate
per-atom quantities which could then be output into dump files.

View File

@ -27,7 +27,7 @@ fix_modify myCMAP energy yes :pre
This command enables CMAP crossterms to be added to simulations which
use the CHARMM force field. These are relevant for any CHARMM model
of a peptide or protein sequences that is 3 or more amino-acid
residues long; see "(Buck)"_#Buck and "(Brooks)"_#Brooks for details,
residues long; see "(Buck)"_#Buck and "(Brooks)"_#Brooks2 for details,
including the analytic energy expressions for CMAP interactions. The
CMAP crossterms add additional potential energy contributions to pairs
of overlapping phi-psi dihedrals of amino-acids, which are important
@ -128,5 +128,5 @@ LAMMPS"_Section_start.html#start_3 section for more info on packages.
[(Buck)] Buck, Bouguet-Bonnet, Pastor, MacKerell Jr., Biophys J, 90, L36
(2006).
:link(Brooks)
:link(Brooks2)
[(Brooks)] Brooks, Brooks, MacKerell Jr., J Comput Chem, 30, 1545 (2009).

View File

@ -56,26 +56,25 @@ fix 4 my_gas gcmc 1 10 10 1 123456543 300.0 -12.5 1.0 region disk :pre
[Description:]
This fix performs grand canonical Monte Carlo (GCMC) exchanges of
atoms or molecules of the given type with an imaginary ideal gas reservoir at
the specified T and chemical potential (mu) as discussed in
"(Frenkel)"_#Frenkel. If used with the "fix nvt"_fix_nh.html command,
simulations in the grand canonical ensemble (muVT, constant chemical
potential, constant volume, and constant temperature) can be
atoms or molecules of the given type with an imaginary ideal gas
reservoir at the specified T and chemical potential (mu) as discussed
in "(Frenkel)"_#Frenkel. If used with the "fix nvt"_fix_nh.html
command, simulations in the grand canonical ensemble (muVT, constant
chemical potential, constant volume, and constant temperature) can be
performed. Specific uses include computing isotherms in microporous
materials, or computing vapor-liquid coexistence curves.
Every N timesteps the fix attempts a number of GCMC exchanges (insertions
or deletions) of gas atoms or molecules of
the given type between the simulation cell and the imaginary
reservoir. It also attempts a number of Monte Carlo
moves (translations and molecule rotations) of gas of the given type
within the simulation cell or region. The average number of
attempted GCMC exchanges is X. The average number of attempted MC moves is M.
M should typically be chosen to be
approximately equal to the expected number of gas atoms or molecules
of the given type within the simulation cell or region,
which will result in roughly one
MC translation per atom or molecule per MC cycle.
Every N timesteps the fix attempts a number of GCMC exchanges
(insertions or deletions) of gas atoms or molecules of the given type
between the simulation cell and the imaginary reservoir. It also
attempts a number of Monte Carlo moves (translations and molecule
rotations) of gas of the given type within the simulation cell or
region. The average number of attempted GCMC exchanges is X. The
average number of attempted MC moves is M. M should typically be
chosen to be approximately equal to the expected number of gas atoms
or molecules of the given type within the simulation cell or region,
which will result in roughly one MC translation per atom or molecule
per MC cycle.
For MC moves of molecular gasses, rotations and translations are each
attempted with 50% probability. For MC moves of atomic gasses,
@ -83,50 +82,50 @@ translations are attempted 100% of the time. For MC exchanges of
either molecular or atomic gasses, deletions and insertions are each
attempted with 50% probability.
All inserted particles are always assigned to two groups: the default group
"all" and the group specified in the fix gcmc command (which can also
be "all"). In addition, particles are also added to any groups specified
by the {group} and {grouptype} keywords.
If inserted particles are individual atoms, they are
assigned the atom type given by the type argument. If they are molecules,
the type argument has no effect and must be set to zero. Instead,
the type of each atom in the inserted molecule is specified
in the file read by the "molecule"_molecule.html command.
All inserted particles are always assigned to two groups: the default
group "all" and the group specified in the fix gcmc command (which can
also be "all"). In addition, particles are also added to any groups
specified by the {group} and {grouptype} keywords. If inserted
particles are individual atoms, they are assigned the atom type given
by the type argument. If they are molecules, the type argument has no
effect and must be set to zero. Instead, the type of each atom in the
inserted molecule is specified in the file read by the
"molecule"_molecule.html command.
This fix cannot be used to perform MC insertions of gas atoms or
molecules other than the exchanged type, but MC deletions,
translations, and rotations can be performed on any atom/molecule in
the fix group. All atoms in the simulation cell can be moved using
regular time integration translations, e.g. via
"fix nvt"_fix_nh.html, resulting in a hybrid GCMC+MD simulation. A
smaller-than-usual timestep size may be needed when running such a
hybrid simulation, especially if the inserted molecules are not well
equilibrated.
regular time integration translations, e.g. via "fix nvt"_fix_nh.html,
resulting in a hybrid GCMC+MD simulation. A smaller-than-usual
timestep size may be needed when running such a hybrid simulation,
especially if the inserted molecules are not well equilibrated.
This command may optionally use the {region} keyword to define an
exchange and move volume. The specified region must have been
previously defined with a "region"_region.html command. It must be
defined with side = {in}. Insertion attempts occur only within the
specified region. For non-rectangular regions, random trial
points are generated within the rectangular bounding box until a point is found
that lies inside the region. If no valid point is generated after 1000 trials,
no insertion is performed, but it is counted as an attempted insertion.
Move and deletion attempt candidates are selected
from gas atoms or molecules within the region. If there are no candidates,
no move or deletion is performed, but it is counted as an attempt move
or deletion. If an attempted move places the atom or molecule center-of-mass outside
the specified region, a new attempted move is generated. This process is repeated
until the atom or molecule center-of-mass is inside the specified region.
specified region. For non-rectangular regions, random trial points are
generated within the rectangular bounding box until a point is found
that lies inside the region. If no valid point is generated after 1000
trials, no insertion is performed, but it is counted as an attempted
insertion. Move and deletion attempt candidates are selected from gas
atoms or molecules within the region. If there are no candidates, no
move or deletion is performed, but it is counted as an attempt move or
deletion. If an attempted move places the atom or molecule
center-of-mass outside the specified region, a new attempted move is
generated. This process is repeated until the atom or molecule
center-of-mass is inside the specified region.
If used with "fix nvt"_fix_nh.html, the temperature of the imaginary
reservoir, T, should be set to be equivalent to the target temperature
used in fix nvt. Otherwise, the imaginary reservoir
will not be in thermal equilibrium with the simulation cell. Also,
it is important that the temperature used by fix nvt be dynamic,
which can be achieved as follows:
used in fix nvt. Otherwise, the imaginary reservoir will not be in
thermal equilibrium with the simulation cell. Also, it is important
that the temperature used by fix nvt be dynamic/dof, which can be
achieved as follows:
compute mdtemp mdatoms temp
compute_modify mdtemp dynamic yes
compute_modify mdtemp dynamic/dof yes
fix mdnvt mdatoms nvt temp 300.0 300.0 10.0
fix_modify mdnvt temp mdtemp :pre
@ -137,16 +136,16 @@ interactions. Specifically, avoid performing so many MC translations
per timestep that atoms can move beyond the neighbor list skin
distance. See the "neighbor"_neighbor.html command for details.
When an atom or molecule is to be inserted, its
coordinates are chosen at a random position within the current
simulation cell or region, and new atom velocities are randomly chosen from
the specified temperature distribution given by T. The effective
temperature for new atom velocities can be increased or decreased
using the optional keyword {tfac_insert} (see below). Relative
coordinates for atoms in a molecule are taken from the template
molecule provided by the user. The center of mass of the molecule
is placed at the insertion point. The orientation of the molecule
is chosen at random by rotating about this point.
When an atom or molecule is to be inserted, its coordinates are chosen
at a random position within the current simulation cell or region, and
new atom velocities are randomly chosen from the specified temperature
distribution given by T. The effective temperature for new atom
velocities can be increased or decreased using the optional keyword
{tfac_insert} (see below). Relative coordinates for atoms in a
molecule are taken from the template molecule provided by the
user. The center of mass of the molecule is placed at the insertion
point. The orientation of the molecule is chosen at random by rotating
about this point.
Individual atoms are inserted, unless the {mol} keyword is used. It
specifies a {template-ID} previously defined using the
@ -158,15 +157,15 @@ command for details. The only settings required to be in this file
are the coordinates and types of atoms in the molecule.
When not using the {mol} keyword, you should ensure you do not delete
atoms that are bonded to other atoms, or LAMMPS will
soon generate an error when it tries to find bonded neighbors. LAMMPS will
warn you if any of the atoms eligible for deletion have a non-zero
molecule ID, but does not check for this at the time of deletion.
atoms that are bonded to other atoms, or LAMMPS will soon generate an
error when it tries to find bonded neighbors. LAMMPS will warn you if
any of the atoms eligible for deletion have a non-zero molecule ID,
but does not check for this at the time of deletion.
If you wish to insert molecules via the {mol} keyword, that will be
treated as rigid bodies, use the {rigid} keyword, specifying as its
value the ID of a separate "fix rigid/small"_fix_rigid.html
command which also appears in your input script.
value the ID of a separate "fix rigid/small"_fix_rigid.html command
which also appears in your input script.
NOTE: If you wish the new rigid molecules (and other rigid molecules)
to be thermostatted correctly via "fix rigid/small/nvt"_fix_rigid.html
@ -179,43 +178,76 @@ their bonds or angles constrained via SHAKE, use the {shake} keyword,
specifying as its value the ID of a separate "fix
shake"_fix_shake.html command which also appears in your input script.
Optionally, users may specify the maximum rotation angle for
molecular rotations using the {maxangle} keyword and specifying
the angle in degrees. Rotations are performed by generating a random
point on the unit sphere and a random rotation angle on the
range \[0,maxangle). The molecule is then rotated by that angle about an
Optionally, users may specify the maximum rotation angle for molecular
rotations using the {maxangle} keyword and specifying the angle in
degrees. Rotations are performed by generating a random point on the
unit sphere and a random rotation angle on the range
\[0,maxangle). The molecule is then rotated by that angle about an
axis passing through the molecule center of mass. The axis is parallel
to the unit vector defined by the point on the unit sphere.
The same procedure is used for randomly rotating molecules when they
are inserted, except that the maximum angle is 360 degrees.
to the unit vector defined by the point on the unit sphere. The same
procedure is used for randomly rotating molecules when they are
inserted, except that the maximum angle is 360 degrees.
Note that fix GCMC does not use configurational bias
MC or any other kind of sampling of intramolecular degrees of freedom.
Inserted molecules can have different orientations, but they will all
have the same intramolecular configuration,
which was specified in the molecule command input.
Note that fix GCMC does not use configurational bias MC or any other
kind of sampling of intramolecular degrees of freedom. Inserted
molecules can have different orientations, but they will all have the
same intramolecular configuration, which was specified in the molecule
command input.
For atomic gasses, inserted atoms have the specified atom type, but
deleted atoms are any atoms that have been inserted or that belong
to the user-specified fix group. For molecular gasses, exchanged
molecules use the same atom types as in the template molecule
supplied by the user. In both cases, exchanged
atoms/molecules are assigned to two groups: the default group "all"
and the group specified in the fix gcmc command (which can also be
"all").
deleted atoms are any atoms that have been inserted or that belong to
the user-specified fix group. For molecular gasses, exchanged
molecules use the same atom types as in the template molecule supplied
by the user. In both cases, exchanged atoms/molecules are assigned to
two groups: the default group "all" and the group specified in the fix
gcmc command (which can also be "all").
The gas reservoir pressure can be specified using the {pressure}
keyword, in which case the user-specified chemical potential is
ignored. For non-ideal gas reservoirs, the user may also specify the
fugacity coefficient using the {fugacity_coeff} keyword.
The chemical potential is a user-specified input parameter defined
as:
:c,image(Eqs/fix_gcmc1.jpg)
The second term mu_ex is the excess chemical potential due to
energetic interactions and is formally zero for the fictitious gas
reservoir but is non-zero for interacting systems. So, while the
chemical potential of the reservoir and the simulation cell are equal,
mu_ex is not, and as a result, the densities of the two are generally
quite different. The first term mu_id is the ideal gas contribution
to the chemical potential. mu_id can be related to the density or
pressure of the fictitious gas reservoir by:
:c,image(Eqs/fix_gcmc2.jpg)
where k is Boltzman's constant,
T is the user-specified temperature, rho is the number density,
P is the pressure, and phi is the fugacity coefficient.
The constant Lambda is required for dimensional consistency.
For all unit styles except {lj} it is defined as the thermal
de Broglie wavelength
:c,image(Eqs/fix_gcmc3.jpg)
where h is Planck's constant, and m is the mass of the exchanged atom
or molecule. For unit style {lj}, Lambda is simply set to the
unity. Note that prior to March 2017, lambda for unit style {lj} was
calculated using the above formula with h set to the rather specific
value of 0.18292026. Chemical potential under the old definition can
be converted to an equivalent value under the new definition by
subtracting 3kTln(Lambda_old).
As an alternative to specifying mu directly, the ideal gas reservoir
can be defined by its pressure P using the {pressure} keyword, in
which case the user-specified chemical potential is ignored. The user
may also specify the fugacity coefficient phi using the
{fugacity_coeff} keyword, which defaults to unity.
The {full_energy} option means that fix GCMC will compute the total
potential energy of the entire simulated system. The total system
energy before and after the proposed GCMC move is then used in the
Metropolis criterion to determine whether or not to accept the
proposed GCMC move. By default, this option is off, in which case
only partial energies are computed to determine the difference in
energy that would be caused by the proposed GCMC move.
proposed GCMC move. By default, this option is off, in which case only
partial energies are computed to determine the difference in energy
that would be caused by the proposed GCMC move.
The {full_energy} option is needed for systems with complicated
potential energy calculations, including the following:
@ -224,7 +256,7 @@ potential energy calculations, including the following:
many-body pair styles
hybrid pair styles
eam pair styles
triclinic systems
tail corrections
need to include potential energy contributions from other fixes :ul
In these cases, LAMMPS will automatically apply the {full_energy}
@ -233,42 +265,43 @@ keyword and issue a warning message.
When the {mol} keyword is used, the {full_energy} option also includes
the intramolecular energy of inserted and deleted molecules. If this
is not desired, the {intra_energy} keyword can be used to define an
amount of energy that is subtracted from the final energy when a molecule
is inserted, and added to the initial energy when a molecule is
deleted. For molecules that have a non-zero intramolecular energy, this
will ensure roughly the same behavior whether or not the {full_energy}
option is used.
amount of energy that is subtracted from the final energy when a
molecule is inserted, and added to the initial energy when a molecule
is deleted. For molecules that have a non-zero intramolecular energy,
this will ensure roughly the same behavior whether or not the
{full_energy} option is used.
Inserted atoms and molecules are assigned random velocities based on the
specified temperature T. Because the relative velocity of
all atoms in the molecule is zero, this may result in inserted molecules
that are systematically too cold. In addition, the intramolecular potential
energy of the inserted molecule may cause the kinetic energy
of the molecule to quickly increase or decrease after insertion.
The {tfac_insert} keyword allows the user to counteract these effects
by changing the temperature used to assign velocities to
inserted atoms and molecules by a constant factor. For a
particular application, some experimentation may be required
to find a value of {tfac_insert} that results in inserted molecules that
equilibrate quickly to the correct temperature.
Inserted atoms and molecules are assigned random velocities based on
the specified temperature T. Because the relative velocity of all
atoms in the molecule is zero, this may result in inserted molecules
that are systematically too cold. In addition, the intramolecular
potential energy of the inserted molecule may cause the kinetic energy
of the molecule to quickly increase or decrease after insertion. The
{tfac_insert} keyword allows the user to counteract these effects by
changing the temperature used to assign velocities to inserted atoms
and molecules by a constant factor. For a particular application, some
experimentation may be required to find a value of {tfac_insert} that
results in inserted molecules that equilibrate quickly to the correct
temperature.
Some fixes have an associated potential energy. Examples of such fixes
include: "efield"_fix_efield.html, "gravity"_fix_gravity.html,
"addforce"_fix_addforce.html, "langevin"_fix_langevin.html,
"restrain"_fix_restrain.html, "temp/berendsen"_fix_temp_berendsen.html,
"restrain"_fix_restrain.html,
"temp/berendsen"_fix_temp_berendsen.html,
"temp/rescale"_fix_temp_rescale.html, and "wall fixes"_fix_wall.html.
For that energy to be included in the total potential energy of the
system (the quantity used when performing GCMC moves),
you MUST enable the "fix_modify"_fix_modify.html {energy} option for
that fix. The doc pages for individual "fix"_fix.html commands
specify if this should be done.
system (the quantity used when performing GCMC moves), you MUST enable
the "fix_modify"_fix_modify.html {energy} option for that fix. The
doc pages for individual "fix"_fix.html commands specify if this
should be done.
Use the {charge} option to insert atoms with a user-specified point
charge. Note that doing so will cause the system to become non-neutral.
LAMMPS issues a warning when using long-range electrostatics (kspace)
with non-neutral systems. See the
"compute group/group"_compute_group_group.html documentation for more
details about simulating non-neutral systems with kspace on.
charge. Note that doing so will cause the system to become
non-neutral. LAMMPS issues a warning when using long-range
electrostatics (kspace) with non-neutral systems. See the "compute
group/group"_compute_group_group.html documentation for more details
about simulating non-neutral systems with kspace on.
Use of this fix typically will cause the number of atoms to fluctuate,
therefore, you will want to use the
@ -276,16 +309,23 @@ therefore, you will want to use the
current number of atoms is used as a normalizing factor each time
temperature is computed. Here is the necessary command:
NOTE: If the density of the cell is initially very small or zero, and
increases to a much larger density after a period of equilibration,
then certain quantities that are only calculated once at the start
(kspace parameters, tail corrections) may no longer be accurate. The
solution is to start a new simulation after the equilibrium density
has been reached.
With some pair_styles, such as "Buckingham"_pair_buck.html,
"Born-Mayer-Huggins"_pair_born.html and "ReaxFF"_pair_reax_c.html,
two atoms placed close to each other may have an arbitrary large,
negative potential energy due to the functional form of the potential.
While these unphysical configurations are inaccessible
to typical dynamical trajectories,
they can be generated by Monte Carlo moves. The {overlap_cutoff}
keyword suppresses these moves by effectively assigning an
infinite positive energy to all new configurations that place any
pair of atoms closer than the specified overlap cutoff distance.
"Born-Mayer-Huggins"_pair_born.html and "ReaxFF"_pair_reax_c.html, two
atoms placed close to each other may have an arbitrary large, negative
potential energy due to the functional form of the potential. While
these unphysical configurations are inaccessible to typical dynamical
trajectories, they can be generated by Monte Carlo moves. The
{overlap_cutoff} keyword suppresses these moves by effectively
assigning an infinite positive energy to all new configurations that
place any pair of atoms closer than the specified overlap cutoff
distance.
compute_modify thermo_temp dynamic yes :pre
@ -295,10 +335,10 @@ derived from LJ parameters for argon, where h* = h/sqrt(sigma^2 *
epsilon * mass), sigma = 3.429 angstroms, epsilon/k = 121.85 K, and
mass = 39.948 amu.
The {group} keyword assigns all inserted atoms to the "group"_group.html
of the group-ID value. The {grouptype} keyword assigns all
inserted atoms of the specified type to the "group"_group.html
of the group-ID value.
The {group} keyword assigns all inserted atoms to the
"group"_group.html of the group-ID value. The {grouptype} keyword
assigns all inserted atoms of the specified type to the
"group"_group.html of the group-ID value.
[Restart, fix_modify, output, run start/stop, minimize info:]
@ -346,15 +386,15 @@ well in parallel. Only usable for 3D simulations.
Note that very lengthy simulations involving insertions/deletions of
billions of gas molecules may run out of atom or molecule IDs and
trigger an error, so it is better to run multiple shorter-duration
simulations. Likewise, very large molecules have not been tested
and may turn out to be problematic.
simulations. Likewise, very large molecules have not been tested and
may turn out to be problematic.
Use of multiple fix gcmc commands in the same input script can be
problematic if using a template molecule. The issue is that the
user-referenced template molecule in the second fix gcmc command
may no longer exist since it might have been deleted by the first
fix gcmc command. An existing template molecule will need to be
referenced by the user for each subsequent fix gcmc command.
user-referenced template molecule in the second fix gcmc command may
no longer exist since it might have been deleted by the first fix gcmc
command. An existing template molecule will need to be referenced by
the user for each subsequent fix gcmc command.
[Related commands:]
@ -366,7 +406,7 @@ referenced by the user for each subsequent fix gcmc command.
[Default:]
The option defaults are mol = no, maxangle = 10, overlap_cutoff = 0.0,
and full_energy = no,
fugacity_coeff = 1, and full_energy = no,
except for the situations where full_energy is required, as
listed above.

View File

@ -91,7 +91,7 @@ their DOF are assumed to be constant. If you are adding atoms or
molecules to the system (see the "fix pour"_fix_pour.html, "fix
deposit"_fix_deposit.html, and "fix gcmc"_fix_gcmc.html commands) or
expect atoms or molecules to be lost (e.g. due to exiting the
simulation box or via "fix evaporation"_fix_evaporation.html), then
simulation box or via "fix evaporate"_fix_evaporate.html), then
this option should be used to insure the temperature is correctly
normalized.

View File

@ -84,9 +84,9 @@ CHARMM force field.
The styles with {charmm} (not {charmmfsw} or {charmmfsh}) in their
name are the older, original LAMMPS implementations. They compute the
LJ and Coulombic interactions with an energy switching function (esw,
a cubic polynomial, shown in the formula below), which ramps the
energy smoothly to zero between the inner and outer cutoff. This can
cause irregularities in pair-wise forces (due to the discontinuous 2nd
shown in the formula below as S(r)), which ramps the energy smoothly
to zero between the inner and outer cutoff. This can cause
irregularities in pair-wise forces (due to the discontinuous 2nd
derivative of energy at the boundaries of the switching region), which
in some cases can result in detectable artifacts in an MD simulation.
@ -94,7 +94,7 @@ The newer styles with {charmmfsw} or {charmmfsh} in their name replace
the energy switching with force switching (fsw) and force shifting
(fsh) functions, for LJ and Coulombic interactions respectively.
These follow the formulas and description given in
"(Steinbach)"_#Steinbach and "(Brooks)"_#Brooks to minimize these
"(Steinbach)"_#Steinbach and "(Brooks)"_#Brooks1 to minimize these
artifacts.
NOTE: The newer {charmmfsw} or {charmmfsh} styles were released in
@ -248,7 +248,7 @@ the MOLECULE and KSPACE packages are installed by default.
:line
:link(Brooks)
:link(Brooks1)
[(Brooks)] Brooks, et al, J Comput Chem, 30, 1545 (2009).
:link(pair-MacKerell)

View File

@ -14,15 +14,23 @@ pair_style oxdna/coaxstk command :h3
[Syntax:]
pair_style style :pre
pair_style style1 :pre
style = {hybrid/overlay oxdna/excv oxdna/stk oxdna/hbond oxdna/xstk oxdna/coaxstk} :ul
pair_coeff * * style2 args :pre
style1 = {hybrid/overlay oxdna/excv oxdna/stk oxdna/hbond oxdna/xstk oxdna/coaxstk} :ul
style2 = {oxdna/stk}
args = list of arguments for these two particular styles :ul
{oxdna2/stk} args = T 6.0 0.4 0.9 0.32 0.6 1.3 0 0.8 0.9 0 0.95 0.9 0 0.95 2.0 0.65 2.0 0.65
T = temperature (oxDNA units, 0.1 = 300 K) :pre
[Examples:]
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 1.61048 6.0 0.4 0.9 0.32 0.6 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/stk 0.1 6.0 0.4 0.9 0.32 0.6 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 0.0 8.0 0.4 0.75 0.34 0.7 1.5 0 0.7 1.5 0 0.7 1.5 0 0.7 0.46 3.141592653589793 0.7 4.0 1.5707963267948966 0.45 4.0 1.5707963267948966 0.45
pair_coeff 1 4 oxdna/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
pair_coeff 2 3 oxdna/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
@ -42,19 +50,23 @@ The exact functional form of the pair styles is rather complex, which manifests
in the above example. The individual potentials consist of products of modulation factors,
which themselves are constructed from a number of more basic potentials
(Morse, Lennard-Jones, harmonic angle and distance) as well as quadratic smoothing and modulation terms.
We refer to "(Ouldridge-DPhil)"_#Ouldridge-DPhil and "(Ouldridge)"_#Ouldridge
We refer to "(Ouldridge-DPhil)"_#Ouldridge-DPhil1 and "(Ouldridge)"_#Ouldridge1
for a detailed description of the oxDNA force field.
NOTE: These pair styles have to be used together with the related oxDNA bond style
{oxdna/fene} for the connectivity of the phosphate backbone (see also documentation of
"bond_style oxdna/fene"_bond_oxdna.html). The coefficients
"bond_style oxdna/fene"_bond_oxdna.html). With one exception the coefficients
in the above example have to be kept fixed and cannot be changed without reparametrizing the entire model.
The exception is the first coefficient after {oxdna/stk} (T=0.1 in the above example).
When using a Langevin thermostat, e.g. through "fix langevin"_fix_langevin.html
or "fix nve/dotc/langevin"_fix_nve_dotc_langevin.html
the temperature coefficients have to be matched to the one used in the fix.
Example input and data files can be found in examples/USER/cgdna/examples/duplex1/ and /duplex2/.
A simple python setup tool which creates single straight or helical DNA strands,
Example input and data files for DNA duplexes can be found in examples/USER/cgdna/examples/oxDNA/ and /oxDNA2/.
A simple python setup tool which creates single straight or helical DNA strands,
DNA duplexes or arrays of DNA duplexes can be found in examples/USER/cgdna/util/.
A technical report with more information on the model, the structure of the input file,
the setup tool and the performance of the LAMMPS-implementation of oxDNA
the setup tool and the performance of the LAMMPS-implementation of oxDNA
can be found "here"_PDF/USER-CGDNA-overview.pdf.
:line
@ -67,14 +79,14 @@ LAMMPS"_Section_start.html#start_3 section for more info on packages.
[Related commands:]
"bond_style oxdna/fene"_bond_oxdna.html, "fix nve/dotc/langevin"_fix_nve_dotc_langevin.html, "pair_coeff"_pair_coeff.html
"bond_style oxdna/fene"_bond_oxdna.html, "fix nve/dotc/langevin"_fix_nve_dotc_langevin.html, "pair_coeff"_pair_coeff.html,
"bond_style oxdna2/fene"_bond_oxdna.html, "pair_style oxdna2/excv"_pair_oxdna2.html
[Default:] none
:line
:link(Ouldridge-DPhil)
:link(Ouldridge-DPhil1)
[(Ouldrigde-DPhil)] T.E. Ouldridge, Coarse-grained modelling of DNA and DNA self-assembly, DPhil. University of Oxford (2011).
:link(Ouldridge)
:link(Ouldridge1)
[(Ouldridge)] T.E. Ouldridge, A.A. Louis, J.P.K. Doye, J. Chem. Phys. 134, 085101 (2011).

102
doc/src/pair_oxdna2.txt Normal file
View File

@ -0,0 +1,102 @@
"LAMMPS WWW Site"_lws - "LAMMPS Documentation"_ld - "LAMMPS Commands"_lc :c
:link(lws,http://lammps.sandia.gov)
:link(ld,Manual.html)
:link(lc,Section_commands.html#comm)
:line
pair_style oxdna2/excv command :h3
pair_style oxdna2/stk command :h3
pair_style oxdna2/hbond command :h3
pair_style oxdna2/xstk command :h3
pair_style oxdna2/coaxstk command :h3
pair_style oxdna2/dh command :h3
[Syntax:]
pair_style style1 :pre
pair_coeff * * style2 args :pre
style1 = {hybrid/overlay oxdna2/excv oxdna2/stk oxdna2/hbond oxdna2/xstk oxdna2/coaxstk oxdna2/dh} :ul
style2 = {oxdna2/stk} or {oxdna2/dh}
args = list of arguments for these two particular styles :ul
{oxdna2/stk} args = T 6.0 0.4 0.9 0.32 0.6 1.3 0 0.8 0.9 0 0.95 0.9 0 0.95 2.0 0.65 2.0 0.65
T = temperature (oxDNA units, 0.1 = 300 K)
{oxdna2/dh} args = T rhos qeff
T = temperature (oxDNA units, 0.1 = 300 K)
rhos = salt concentration (mole per litre)
qeff = effective charge (elementary charges) :pre
[Examples:]
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 0.1 6.0 0.4 0.9 0.32 0.6 1.3 0 0.8 0.9 0 0.95 0.9 0 0.95 2.0 0.65 2.0 0.65
pair_coeff * * oxdna2/hbond 0.0 8.0 0.4 0.75 0.34 0.7 1.5 0 0.7 1.5 0 0.7 1.5 0 0.7 0.46 3.141592653589793 0.7 4.0 1.5707963267948966 0.45 4.0 1.5707963267948966 0.45
pair_coeff 1 4 oxdna2/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
pair_coeff 2 3 oxdna2/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
pair_coeff * * oxdna2/xstk 47.5 0.575 0.675 0.495 0.655 2.25 0.791592653589793 0.58 1.7 1.0 0.68 1.7 1.0 0.68 1.5 0 0.65 1.7 0.875 0.68 1.7 0.875 0.68
pair_coeff * * oxdna2/coaxstk 58.5 0.4 0.6 0.22 0.58 2.0 2.891592653589793 0.65 1.3 0 0.8 0.9 0 0.95 0.9 0 0.95 40.0 3.116592653589793
pair_coeff * * oxdna2/dh 0.1 1.0 0.815 :pre
[Description:]
The {oxdna2} pair styles compute the pairwise-additive parts of the oxDNA force field
for coarse-grained modelling of DNA. The effective interaction between the nucleotides consists of potentials for the
excluded volume interaction {oxdna2/excv}, the stacking {oxdna2/stk}, cross-stacking {oxdna2/xstk}
and coaxial stacking interaction {oxdna2/coaxstk}, electrostatic Debye-Hueckel interaction {oxdna2/dh}
as well as the hydrogen-bonding interaction {oxdna2/hbond} between complementary pairs of nucleotides on
opposite strands.
The exact functional form of the pair styles is rather complex.
The individual potentials consist of products of modulation factors,
which themselves are constructed from a number of more basic potentials
(Morse, Lennard-Jones, harmonic angle and distance) as well as quadratic smoothing and modulation terms.
We refer to "(Snodin)"_#Snodin and the original oxDNA publications "(Ouldridge-DPhil)"_#Ouldridge-DPhil2
and "(Ouldridge)"_#Ouldridge2 for a detailed description of the oxDNA2 force field.
NOTE: These pair styles have to be used together with the related oxDNA2 bond style
{oxdna2/fene} for the connectivity of the phosphate backbone (see also documentation of
"bond_style oxdna2/fene"_bond_oxdna.html). Almost all coefficients
in the above example have to be kept fixed and cannot be changed without reparametrizing the entire model.
Exceptions are the first coefficient after {oxdna2/stk} (T=0.1 in the above example) and the coefficients
after {oxdna2/dh} (T=0.1, rhos=1.0, qeff=0.815 in the above example). When using a Langevin thermostat
e.g. through "fix langevin"_fix_langevin.html or "fix nve/dotc/langevin"_fix_nve_dotc_langevin.html
the temperature coefficients have to be matched to the one used in the fix.
Example input and data files for DNA duplexes can be found in examples/USER/cgdna/examples/oxDNA/ and /oxDNA2/.
A simple python setup tool which creates single straight or helical DNA strands,
DNA duplexes or arrays of DNA duplexes can be found in examples/USER/cgdna/util/.
A technical report with more information on the model, the structure of the input file,
the setup tool and the performance of the LAMMPS-implementation of oxDNA
can be found "here"_PDF/USER-CGDNA-overview.pdf.
:line
[Restrictions:]
These pair styles can only be used if LAMMPS was built with the
USER-CGDNA package and the MOLECULE and ASPHERE package. See the "Making
LAMMPS"_Section_start.html#start_3 section for more info on packages.
[Related commands:]
"bond_style oxdna2/fene"_bond_oxdna.html, "fix nve/dotc/langevin"_fix_nve_dotc_langevin.html, "pair_coeff"_pair_coeff.html,
"bond_style oxdna/fene"_bond_oxdna.html, "pair_style oxdna/excv"_pair_oxdna.html
[Default:] none
:line
:link(Snodin)
[(Snodin)] B.E. Snodin, F. Randisi, M. Mosayebi, et al., J. Chem. Phys. 142, 234901 (2015).
:link(Ouldridge-DPhil2)
[(Ouldrigde-DPhil)] T.E. Ouldridge, Coarse-grained modelling of DNA and DNA self-assembly, DPhil. University of Oxford (2011).
:link(Ouldridge2)
[(Ouldridge)] T.E. Ouldridge, A.A. Louis, J.P.K. Doye, J. Chem. Phys. 134, 085101 (2011).

View File

@ -68,6 +68,7 @@ Pair Styles :h1
pair_nm
pair_none
pair_oxdna
pair_oxdna2
pair_peri
pair_polymorphic
pair_quip

View File

@ -7,9 +7,9 @@ Input, data and log files for a DNA duplex (double-stranded DNA)
consisiting of 5 base pairs. The duplex contains two strands with
complementary base pairs. The topology is
A - A - A - A - A
A - C - G - T - A
| | | | |
T - T - T - T - T
T - G - C - A - T
/examples/duplex2:
Input, data and log files for a nicked DNA duplex (double-stranded DNA)
@ -18,9 +18,9 @@ complementary base pairs, but the backbone on one side is not continuous:
two individual strands on one side form a duplex with a longer single
strand on the other side. The topology is
A - A - A - A - A - A - A - A
A - C - G - T - A - C - G - T
| | | | | | | |
T - T - T T - T - T - T - T
T - G - C - A T - G - C - A
/util:
This directory contains a simple python setup tool which creates

View File

@ -1,74 +0,0 @@
# LAMMPS data file
10 atoms
10 ellipsoids
8 bonds
4 atom types
1 bond types
# System size
-20.000000 20.000000 xlo xhi
-20.000000 20.000000 ylo yhi
-20.000000 20.000000 zlo zhi
# Atom masses for each atom type
Masses
1 3.1575
2 3.1575
3 3.1575
4 3.1575
# Atom-ID, type, position, molecule-ID, ellipsoid flag, density
Atoms
1 1 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 1 1 1
2 1 1.3274493266864451e-01 -4.2912827978022683e-01 3.7506163469402809e-01 1 1 1
3 1 4.8460810659772807e-01 -7.0834970533509178e-01 7.5012326938805618e-01 1 1 1
4 1 9.3267359196674593e-01 -7.4012419946742802e-01 1.1251849040820843e+00 1 1 1
5 1 1.3204192238113461e+00 -5.1335201721887447e-01 1.5002465387761124e+00 1 1 1
6 4 1.9958077618865377e-01 5.1335201721887447e-01 1.5002465387761124e+00 1 1 1
7 4 5.8732640803325409e-01 7.4012419946742802e-01 1.1251849040820843e+00 1 1 1
8 4 1.0353918934022719e+00 7.0834970533509178e-01 7.5012326938805618e-01 1 1 1
9 4 1.3872550673313555e+00 4.2912827978022683e-01 3.7506163469402809e-01 1 1 1
10 4 1.5200000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 1 1 1
# Atom-ID, translational, rotational velocity
Velocities
1 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00
2 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00
3 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00
4 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00
5 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00
6 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00
7 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00
8 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00
9 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00
10 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00
# Atom-ID, shape, quaternion
Ellipsoids
1 1.1739845031423408e+00 1.1739845031423408e+00 1.1739845031423408e+00 1.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00
2 1.1739845031423408e+00 1.1739845031423408e+00 1.1739845031423408e+00 9.5533648912560598e-01 0.0000000000000000e+00 0.0000000000000000e+00 2.9552020666133955e-01
3 1.1739845031423408e+00 1.1739845031423408e+00 1.1739845031423408e+00 8.2533561490967822e-01 0.0000000000000000e+00 0.0000000000000000e+00 5.6464247339503526e-01
4 1.1739845031423408e+00 1.1739845031423408e+00 1.1739845031423408e+00 6.2160996827066439e-01 0.0000000000000000e+00 0.0000000000000000e+00 7.8332690962748319e-01
5 1.1739845031423408e+00 1.1739845031423408e+00 1.1739845031423408e+00 3.6235775447667351e-01 0.0000000000000000e+00 0.0000000000000000e+00 9.3203908596722607e-01
6 1.1739845031423408e+00 1.1739845031423408e+00 1.1739845031423408e+00 0.0000000000000000e+00 9.3203908596722607e-01 -3.6235775447667351e-01 0.0000000000000000e+00
7 1.1739845031423408e+00 1.1739845031423408e+00 1.1739845031423408e+00 0.0000000000000000e+00 7.8332690962748319e-01 -6.2160996827066439e-01 0.0000000000000000e+00
8 1.1739845031423408e+00 1.1739845031423408e+00 1.1739845031423408e+00 0.0000000000000000e+00 5.6464247339503526e-01 -8.2533561490967822e-01 0.0000000000000000e+00
9 1.1739845031423408e+00 1.1739845031423408e+00 1.1739845031423408e+00 0.0000000000000000e+00 2.9552020666133955e-01 -9.5533648912560598e-01 0.0000000000000000e+00
10 1.1739845031423408e+00 1.1739845031423408e+00 1.1739845031423408e+00 0.0000000000000000e+00 0.0000000000000000e+00 -1.0000000000000000e+00 0.0000000000000000e+00
# Bond topology
Bonds
1 1 1 2
2 1 2 3
3 1 3 4
4 1 4 5
5 1 6 7
6 1 7 8
7 1 8 9
8 1 9 10

View File

@ -1,75 +0,0 @@
variable number equal 1
variable ofreq equal 1000
variable efreq equal 1000
units lj
dimension 3
newton off
boundary p p p
atom_style hybrid bond ellipsoid
atom_modify sort 0 1.0
# Pair interactions require lists of neighbours to be calculated
neighbor 1.0 bin
neigh_modify every 1 delay 0 check yes
read_data data.duplex1
set atom * mass 3.1575
group all type 1 4
# oxDNA bond interactions - FENE backbone
bond_style oxdna_fene
bond_coeff * 2.0 0.25 0.7525
# oxDNA pair interactions
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 1.61048 6.0 0.4 0.9 0.32 0.6 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 0.0 8.0 0.4 0.75 0.34 0.7 1.5 0 0.7 1.5 0 0.7 1.5 0 0.7 0.46 3.141592653589793 0.7 4.0 1.5707963267948966 0.45 4.0 1.5707963267948966 0.45
pair_coeff 1 4 oxdna_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
pair_coeff 2 3 oxdna_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
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
# NVE ensemble
#fix 1 all nve/dotc/langevin 0.1 0.1 0.03 457145 angmom 10
fix 1 all nve/dot
timestep 1e-5
#comm_style tiled
#fix 3 all balance 10000 1.1 rcb
#compute mol all chunk/atom molecule
#compute mychunk all vcm/chunk mol
#fix 4 all ave/time 10000 1 10000 c_mychunk[1] c_mychunk[2] c_mychunk[3] file vcm.txt mode vector
dump pos all xyz ${ofreq} traj.${number}.xyz
compute quat all property/atom quatw quati quatj quatk
dump quat all custom ${ofreq} quat.${number}.txt id c_quat[1] c_quat[2] c_quat[3] c_quat[4]
dump_modify quat sort id
dump_modify quat format line "%d %13.6le %13.6le %13.6le %13.6le"
compute erot all erotate/asphere
compute ekin all ke
compute epot all pe
variable erot equal c_erot
variable ekin equal c_ekin
variable epot equal c_epot
variable etot equal c_erot+c_ekin+c_epot
fix 5 all print ${efreq} "$(step) ekin = ${ekin} | erot = ${erot} | epot = ${epot} | etot = ${etot}" screen yes
dump out all custom ${ofreq} out.${number}.txt id x y z vx vy vz fx fy fz tqx tqy tqz
dump_modify out sort id
dump_modify out format line "%d %13.6le %13.6le %13.6le %13.6le %13.6le %13.6le %13.6le %13.6le %13.6le %13.6le %13.6le %13.6le"
run 1000000
#write_restart config.${number}.*

File diff suppressed because it is too large Load Diff

File diff suppressed because it is too large Load Diff

View File

@ -1,97 +0,0 @@
# LAMMPS data file
16 atoms
16 ellipsoids
13 bonds
4 atom types
1 bond types
# System size
-20.0 20.0 xlo xhi
-20.0 20.0 ylo yhi
-20.0 20.0 zlo zhi
# Atom masses for each atom type
Masses
1 3.1575
2 3.1575
3 3.1575
4 3.1575
# Atom-ID, type, position, molecule-ID, ellipsoid flag, density
Atoms
1 1 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 1 1 1
2 1 1.327449326686445e-01 -4.291282797802268e-01 3.750616346940281e-01 1 1 1
3 1 4.846081065977281e-01 -7.083497053350921e-01 7.501232693880562e-01 1 1 1
4 1 9.326735919667459e-01 -7.401241994674285e-01 1.125184904082084e+00 1 1 1
5 1 1.320419223811347e+00 -5.133520172188747e-01 1.500246538776112e+00 1 1 1
6 1 1.512394297416339e+00 -1.072512061254991e-01 1.875308173470140e+00 1 1 1
7 1 1.441536396413952e+00 3.363155369040876e-01 2.250369808164169e+00 1 1 1
8 1 1.132598224218932e+00 6.623975870343269e-01 2.625431442858197e+00 1 1 1
9 4 5.873264080332541e-01 7.401241994674285e-01 1.125184904082084e+00 1 1 1
10 4 1.035391893402272e+00 7.083497053350921e-01 7.501232693880562e-01 1 1 1
11 4 1.387255067331356e+00 4.291282797802267e-01 3.750616346940281e-01 1 1 1
12 4 1.520000000000000e+00 1.260981291332700e-33 0.000000000000000e+00 1 1 1
13 4 3.874017757810680e-01 -6.623975870343268e-01 2.625431442858197e+00 1 1 1
14 4 7.846360358604798e-02 -3.363155369040874e-01 2.250369808164169e+00 1 1 1
15 4 7.605702583661333e-03 1.072512061254995e-01 1.875308173470140e+00 1 1 1
16 4 1.995807761886533e-01 5.133520172188748e-01 1.500246538776112e+00 1 1 1
# Atom-ID, translational, rotational velocity
Velocities
1 0.0 0.0 0.0 0.0 0.0 0.0
2 0.0 0.0 0.0 0.0 0.0 0.0
3 0.0 0.0 0.0 0.0 0.0 0.0
4 0.0 0.0 0.0 0.0 0.0 0.0
5 0.0 0.0 0.0 0.0 0.0 0.0
6 0.0 0.0 0.0 0.0 0.0 0.0
7 0.0 0.0 0.0 0.0 0.0 0.0
8 0.0 0.0 0.0 0.0 0.0 0.0
9 0.0 0.0 0.0 0.0 0.0 0.0
10 0.0 0.0 0.0 0.0 0.0 0.0
11 0.0 0.0 0.0 0.0 0.0 0.0
12 0.0 0.0 0.0 0.0 0.0 0.0
13 0.0 0.0 0.0 0.0 0.0 0.0
14 0.0 0.0 0.0 0.0 0.0 0.0
15 0.0 0.0 0.0 0.0 0.0 0.0
16 0.0 0.0 0.0 0.0 0.0 0.0
# Atom-ID, shape, quaternion
Ellipsoids
1 1.1739845031423408 1.1739845031423408 1.1739845031423408 1.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00
2 1.1739845031423408 1.1739845031423408 1.1739845031423408 9.553364891256060e-01 0.000000000000000e+00 0.000000000000000e+00 2.955202066613395e-01
3 1.1739845031423408 1.1739845031423408 1.1739845031423408 8.253356149096783e-01 0.000000000000000e+00 0.000000000000000e+00 5.646424733950354e-01
4 1.1739845031423408 1.1739845031423408 1.1739845031423408 6.216099682706646e-01 0.000000000000000e+00 0.000000000000000e+00 7.833269096274833e-01
5 1.1739845031423408 1.1739845031423408 1.1739845031423408 3.623577544766736e-01 0.000000000000000e+00 0.000000000000000e+00 9.320390859672263e-01
6 1.1739845031423408 1.1739845031423408 1.1739845031423408 7.073720166770291e-02 0.000000000000000e+00 0.000000000000000e+00 9.974949866040544e-01
7 1.1739845031423408 1.1739845031423408 1.1739845031423408 -2.272020946930869e-01 -0.000000000000000e+00 0.000000000000000e+00 9.738476308781953e-01
8 1.1739845031423408 1.1739845031423408 1.1739845031423408 -5.048461045998575e-01 -0.000000000000000e+00 0.000000000000000e+00 8.632093666488738e-01
9 1.1739845031423408 1.1739845031423408 1.1739845031423408 4.796493962806427e-17 7.833269096274833e-01 -6.216099682706646e-01 3.806263289803786e-17
10 1.1739845031423408 1.1739845031423408 1.1739845031423408 5.707093416549944e-17 5.646424733950354e-01 -8.253356149096784e-01 2.218801320830406e-17
11 1.1739845031423408 1.1739845031423408 1.1739845031423408 6.107895212550935e-17 2.955202066613394e-01 -9.553364891256061e-01 4.331404380149668e-18
12 1.1739845031423408 1.1739845031423408 1.1739845031423408 5.963096920061075e-17 0.000000000000000e+00 -1.000000000000000e+00 -1.391211590127312e-17
13 1.1739845031423408 1.1739845031423408 1.1739845031423408 5.285632939302787e-17 8.632093666488739e-01 5.048461045998572e-01 -3.091290830301125e-17
14 1.1739845031423408 1.1739845031423408 1.1739845031423408 4.136019110019290e-17 9.738476308781953e-01 2.272020946930868e-01 -4.515234267244800e-17
15 1.1739845031423408 1.1739845031423408 1.1739845031423408 2.616947011741696e-17 9.974949866040544e-01 -7.073720166770313e-02 -5.535845274597425e-17
16 1.1739845031423408 1.1739845031423408 1.1739845031423408 8.641108308308281e-18 9.320390859672264e-01 -3.623577544766736e-01 -6.061955710708163e-17
# Bond-ID, type, atom pairs
Bonds
1 1 1 2
2 1 2 3
3 1 3 4
4 1 4 5
5 1 5 6
6 1 6 7
7 1 7 8
8 1 13 14
9 1 14 15
10 1 15 16
11 1 9 10
12 1 10 11
13 1 11 12

View File

@ -1,75 +0,0 @@
variable number equal 2
variable ofreq equal 1000
variable efreq equal 1000
units lj
dimension 3
newton off
boundary p p p
atom_style hybrid bond ellipsoid
atom_modify sort 0 1.0
# Pair interactions require lists of neighbours to be calculated
neighbor 1.0 bin
neigh_modify every 1 delay 0 check yes
read_data data.duplex2
set atom * mass 3.1575
group all type 1 4
# oxDNA bond interactions - FENE backbone
bond_style oxdna_fene
bond_coeff * 2.0 0.25 0.7525
# oxDNA pair interactions
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 1.61048 6.0 0.4 0.9 0.32 0.6 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 0.0 8.0 0.4 0.75 0.34 0.7 1.5 0 0.7 1.5 0 0.7 1.5 0 0.7 0.46 3.141592653589793 0.7 4.0 1.5707963267948966 0.45 4.0 1.5707963267948966 0.45
pair_coeff 1 4 oxdna_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
pair_coeff 2 3 oxdna_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
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
# NVE ensemble
fix 1 all nve/dotc/langevin 0.1 0.1 0.03 457145 angmom 10
#fix 1 all nve/dot
timestep 1e-5
#comm_style tiled
#fix 3 all balance 10000 1.1 rcb
#compute mol all chunk/atom molecule
#compute mychunk all vcm/chunk mol
#fix 4 all ave/time 10000 1 10000 c_mychunk[1] c_mychunk[2] c_mychunk[3] file vcm.txt mode vector
dump pos all xyz ${ofreq} traj.${number}.xyz
compute quat all property/atom quatw quati quatj quatk
dump quat all custom ${ofreq} quat.${number}.txt id c_quat[1] c_quat[2] c_quat[3] c_quat[4]
dump_modify quat sort id
dump_modify quat format line "%d %13.6le %13.6le %13.6le %13.6le"
compute erot all erotate/asphere
compute ekin all ke
compute epot all pe
variable erot equal c_erot
variable ekin equal c_ekin
variable epot equal c_epot
variable etot equal c_erot+c_ekin+c_epot
fix 5 all print ${efreq} "$(step) ekin = ${ekin} | erot = ${erot} | epot = ${epot} | etot = ${etot}" screen yes
dump out all custom ${ofreq} out.${number}.txt id x y z vx vy vz fx fy fz tqx tqy tqz
dump_modify out sort id
dump_modify out format line "%d %13.6le %13.6le %13.6le %13.6le %13.6le %13.6le %13.6le %13.6le %13.6le %13.6le %13.6le %13.6le"
run 1000000
#write_restart config.${number}.*

File diff suppressed because it is too large Load Diff

File diff suppressed because it is too large Load Diff

View File

@ -0,0 +1,73 @@
# LAMMPS data file
10 atoms
10 ellipsoids
8 bonds
4 atom types
1 bond types
# System size
-20.000000 20.000000 xlo xhi
-20.000000 20.000000 ylo yhi
-20.000000 20.000000 zlo zhi
Masses
1 3.1575
2 3.1575
3 3.1575
4 3.1575
# Atom-ID, type, position, molecule-ID, ellipsoid flag, density
Atoms
1 1 -6.000000000000001e-01 0.000000000000000e+00 0.000000000000000e+00 1 1 1
2 2 -4.860249842674776e-01 -3.518234140414736e-01 3.897628551303122e-01 1 1 1
3 3 -1.874009511073395e-01 -5.699832309147915e-01 7.795257102606244e-01 1 1 1
4 4 1.824198365552941e-01 -5.715968887521518e-01 1.169288565390937e+00 1 1 1
5 1 4.829362784135484e-01 -3.560513319622209e-01 1.559051420521249e+00 1 1 1
6 4 -4.829362784135484e-01 3.560513319622209e-01 1.559051420521249e+00 2 1 1
7 1 -1.824198365552941e-01 5.715968887521516e-01 1.169288565390937e+00 2 1 1
8 2 1.874009511073395e-01 5.699832309147913e-01 7.795257102606243e-01 2 1 1
9 3 4.860249842674775e-01 3.518234140414733e-01 3.897628551303121e-01 2 1 1
10 4 5.999999999999996e-01 -1.332267629550188e-16 -1.110223024625157e-16 2 1 1
# Atom-ID, translational, rotational velocity
Velocities
1 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00
2 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00
3 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00
4 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00
5 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00
6 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00
7 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00
8 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00
9 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00
10 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00
# Atom-ID, shape, quaternion
Ellipsoids
1 1.173984503142341e+00 1.173984503142341e+00 1.173984503142341e+00 1.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00
2 1.173984503142341e+00 1.173984503142341e+00 1.173984503142341e+00 9.513258223252946e-01 0.000000000000000e+00 0.000000000000000e+00 3.081869234362515e-01
3 1.173984503142341e+00 1.173984503142341e+00 1.173984503142341e+00 8.100416404457962e-01 0.000000000000000e+00 0.000000000000000e+00 5.863723567357894e-01
4 1.173984503142341e+00 1.173984503142341e+00 1.173984503142341e+00 5.899012371043606e-01 0.000000000000000e+00 0.000000000000000e+00 8.074754054847398e-01
5 1.173984503142341e+00 1.173984503142341e+00 1.173984503142341e+00 3.123349185122326e-01 0.000000000000000e+00 0.000000000000000e+00 9.499720515246527e-01
6 1.173984503142341e+00 1.173984503142341e+00 1.173984503142341e+00 0.000000000000000e+00 9.499720515246527e-01 -3.123349185122326e-01 -0.000000000000000e+00
7 1.173984503142341e+00 1.173984503142341e+00 1.173984503142341e+00 0.000000000000000e+00 8.074754054847401e-01 -5.899012371043604e-01 0.000000000000000e+00
8 1.173984503142341e+00 1.173984503142341e+00 1.173984503142341e+00 0.000000000000000e+00 5.863723567357896e-01 -8.100416404457959e-01 0.000000000000000e+00
9 1.173984503142341e+00 1.173984503142341e+00 1.173984503142341e+00 -0.000000000000000e+00 -3.081869234362514e-01 9.513258223252947e-01 0.000000000000000e+00
10 1.173984503142341e+00 1.173984503142341e+00 1.173984503142341e+00 -0.000000000000000e+00 1.110223024625157e-16 1.000000000000000e+00 -0.000000000000000e+00
# Bond topology
Bonds
1 1 1 2
2 1 2 3
3 1 3 4
4 1 4 5
5 1 6 7
6 1 7 8
7 1 8 9
8 1 9 10

View File

@ -0,0 +1,77 @@
variable number equal 1
variable ofreq equal 1000
variable efreq equal 1000
units lj
dimension 3
newton off
boundary p p p
atom_style hybrid bond ellipsoid
atom_modify sort 0 1.0
# Pair interactions require lists of neighbours to be calculated
neighbor 1.0 bin
neigh_modify every 1 delay 0 check yes
read_data data.duplex1
set atom * mass 3.1575
group all type 1 4
# oxDNA bond interactions - FENE backbone
bond_style oxdna/fene
bond_coeff * 2.0 0.25 0.7525
# oxDNA pair interactions
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 0.1 6.0 0.4 0.9 0.32 0.6 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 0.0 8.0 0.4 0.75 0.34 0.7 1.5 0 0.7 1.5 0 0.7 1.5 0 0.7 0.46 3.141592653589793 0.7 4.0 1.5707963267948966 0.45 4.0 1.5707963267948966 0.45
pair_coeff 1 4 oxdna/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
pair_coeff 2 3 oxdna/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
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
# NVE ensemble
fix 1 all nve/dot
#fix 1 all nve/dotc/langevin 0.1 0.1 0.03 457145 angmom 10
#fix 1 all nve/asphere
#fix 2 all langevin 0.1 0.1 0.03 457145 angmom 10
timestep 1e-5
#comm_style tiled
#fix 3 all balance 10000 1.1 rcb
#compute mol all chunk/atom molecule
#compute mychunk all vcm/chunk mol
#fix 4 all ave/time 10000 1 10000 c_mychunk[1] c_mychunk[2] c_mychunk[3] file vcm.txt mode vector
#dump pos all xyz ${ofreq} traj.${number}.xyz
#compute quat all property/atom quatw quati quatj quatk
#dump quat all custom ${ofreq} quat.${number}.txt id c_quat[1] c_quat[2] c_quat[3] c_quat[4]
#dump_modify quat sort id
#dump_modify quat format line "%d %13.6le %13.6le %13.6le %13.6le"
compute erot all erotate/asphere
compute ekin all ke
compute epot all pe
variable erot equal c_erot
variable ekin equal c_ekin
variable epot equal c_epot
variable etot equal c_erot+c_ekin+c_epot
fix 5 all print ${efreq} "$(step) ekin = ${ekin} | erot = ${erot} | epot = ${epot} | etot = ${etot}" screen yes
#dump out all custom ${ofreq} out.${number}.txt id x y z vx vy vz fx fy fz tqx tqy tqz
#dump_modify out sort id
#dump_modify out format line "%d %13.6le %13.6le %13.6le %13.6le %13.6le %13.6le %13.6le %13.6le %13.6le %13.6le %13.6le %13.6le"
run 1000000
#write_restart config.${number}.*

File diff suppressed because it is too large Load Diff

File diff suppressed because it is too large Load Diff

View File

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

View File

@ -0,0 +1,77 @@
variable number equal 2
variable ofreq equal 1000
variable efreq equal 1000
units lj
dimension 3
newton off
boundary p p p
atom_style hybrid bond ellipsoid
atom_modify sort 0 1.0
# Pair interactions require lists of neighbours to be calculated
neighbor 1.0 bin
neigh_modify every 1 delay 0 check yes
read_data data.duplex2
set atom * mass 3.1575
group all type 1 4
# oxDNA bond interactions - FENE backbone
bond_style oxdna/fene
bond_coeff * 2.0 0.25 0.7525
# oxDNA pair interactions
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 0.1 6.0 0.4 0.9 0.32 0.6 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 0.0 8.0 0.4 0.75 0.34 0.7 1.5 0 0.7 1.5 0 0.7 1.5 0 0.7 0.46 3.141592653589793 0.7 4.0 1.5707963267948966 0.45 4.0 1.5707963267948966 0.45
pair_coeff 1 4 oxdna/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
pair_coeff 2 3 oxdna/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
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
# NVE ensemble
#fix 1 all nve/dot
fix 1 all nve/dotc/langevin 0.1 0.1 0.03 457145 angmom 10
#fix 1 all nve/asphere
#fix 2 all langevin 0.1 0.1 0.03 457145 angmom 10
timestep 1e-5
#comm_style tiled
#fix 3 all balance 10000 1.1 rcb
#compute mol all chunk/atom molecule
#compute mychunk all vcm/chunk mol
#fix 4 all ave/time 10000 1 10000 c_mychunk[1] c_mychunk[2] c_mychunk[3] file vcm.txt mode vector
#dump pos all xyz ${ofreq} traj.${number}.xyz
#compute quat all property/atom quatw quati quatj quatk
#dump quat all custom ${ofreq} quat.${number}.txt id c_quat[1] c_quat[2] c_quat[3] c_quat[4]
#dump_modify quat sort id
#dump_modify quat format line "%d %13.6le %13.6le %13.6le %13.6le"
compute erot all erotate/asphere
compute ekin all ke
compute epot all pe
variable erot equal c_erot
variable ekin equal c_ekin
variable epot equal c_epot
variable etot equal c_erot+c_ekin+c_epot
fix 5 all print ${efreq} "$(step) ekin = ${ekin} | erot = ${erot} | epot = ${epot} | etot = ${etot}" screen yes
#dump out all custom ${ofreq} out.${number}.txt id x y z vx vy vz fx fy fz tqx tqy tqz
#dump_modify out sort id
#dump_modify out format line "%d %13.6le %13.6le %13.6le %13.6le %13.6le %13.6le %13.6le %13.6le %13.6le %13.6le %13.6le %13.6le"
run 1000000
#write_restart config.${number}.*

File diff suppressed because it is too large Load Diff

File diff suppressed because it is too large Load Diff

View File

@ -0,0 +1,73 @@
# LAMMPS data file
10 atoms
10 ellipsoids
8 bonds
4 atom types
1 bond types
# System size
-20.000000 20.000000 xlo xhi
-20.000000 20.000000 ylo yhi
-20.000000 20.000000 zlo zhi
Masses
1 3.1575
2 3.1575
3 3.1575
4 3.1575
# Atom-ID, type, position, molecule-ID, ellipsoid flag, density
Atoms
1 1 -6.000000000000001e-01 0.000000000000000e+00 0.000000000000000e+00 1 1 1
2 2 -4.860249842674776e-01 -3.518234140414736e-01 3.897628551303122e-01 1 1 1
3 3 -1.874009511073395e-01 -5.699832309147915e-01 7.795257102606244e-01 1 1 1
4 4 1.824198365552941e-01 -5.715968887521518e-01 1.169288565390937e+00 1 1 1
5 1 4.829362784135484e-01 -3.560513319622209e-01 1.559051420521249e+00 1 1 1
6 4 -4.829362784135484e-01 3.560513319622209e-01 1.559051420521249e+00 2 1 1
7 1 -1.824198365552941e-01 5.715968887521516e-01 1.169288565390937e+00 2 1 1
8 2 1.874009511073395e-01 5.699832309147913e-01 7.795257102606243e-01 2 1 1
9 3 4.860249842674775e-01 3.518234140414733e-01 3.897628551303121e-01 2 1 1
10 4 5.999999999999996e-01 -1.332267629550188e-16 -1.110223024625157e-16 2 1 1
# Atom-ID, translational, rotational velocity
Velocities
1 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00
2 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00
3 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00
4 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00
5 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00
6 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00
7 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00
8 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00
9 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00
10 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00
# Atom-ID, shape, quaternion
Ellipsoids
1 1.173984503142341e+00 1.173984503142341e+00 1.173984503142341e+00 1.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00
2 1.173984503142341e+00 1.173984503142341e+00 1.173984503142341e+00 9.513258223252946e-01 0.000000000000000e+00 0.000000000000000e+00 3.081869234362515e-01
3 1.173984503142341e+00 1.173984503142341e+00 1.173984503142341e+00 8.100416404457962e-01 0.000000000000000e+00 0.000000000000000e+00 5.863723567357894e-01
4 1.173984503142341e+00 1.173984503142341e+00 1.173984503142341e+00 5.899012371043606e-01 0.000000000000000e+00 0.000000000000000e+00 8.074754054847398e-01
5 1.173984503142341e+00 1.173984503142341e+00 1.173984503142341e+00 3.123349185122326e-01 0.000000000000000e+00 0.000000000000000e+00 9.499720515246527e-01
6 1.173984503142341e+00 1.173984503142341e+00 1.173984503142341e+00 0.000000000000000e+00 9.499720515246527e-01 -3.123349185122326e-01 -0.000000000000000e+00
7 1.173984503142341e+00 1.173984503142341e+00 1.173984503142341e+00 0.000000000000000e+00 8.074754054847401e-01 -5.899012371043604e-01 0.000000000000000e+00
8 1.173984503142341e+00 1.173984503142341e+00 1.173984503142341e+00 0.000000000000000e+00 5.863723567357896e-01 -8.100416404457959e-01 0.000000000000000e+00
9 1.173984503142341e+00 1.173984503142341e+00 1.173984503142341e+00 -0.000000000000000e+00 -3.081869234362514e-01 9.513258223252947e-01 0.000000000000000e+00
10 1.173984503142341e+00 1.173984503142341e+00 1.173984503142341e+00 -0.000000000000000e+00 1.110223024625157e-16 1.000000000000000e+00 -0.000000000000000e+00
# Bond topology
Bonds
1 1 1 2
2 1 2 3
3 1 3 4
4 1 4 5
5 1 6 7
6 1 7 8
7 1 8 9
8 1 9 10

View File

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

File diff suppressed because it is too large Load Diff

File diff suppressed because it is too large Load Diff

View File

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

View File

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

File diff suppressed because it is too large Load Diff

File diff suppressed because it is too large Load Diff

View File

@ -1,4 +1,4 @@
single 0,0,0:AAAAA
single_helix 0,0,0:AAAAA
duplex 0,0,0:AAAAA
duplex_array 10,10:-112.0:AAAAA
single 0,0,0:ACGTA
single_helix 0,0,0:ACGTA
duplex 0,0,0:ACGTA
duplex_array 10,10:-112.0:ACGTA

View File

@ -1,18 +1,23 @@
# GCMC for LJ simple fluid, no dynamics
# T = 2.0
# rho ~ 0.5
# p ~ 1.5
# mu_ex ~ 0.0
# comparable to Frenkel and Smit GCMC Case Study, Figure 5.8
# variables available on command line
# variables modifiable using -var command line switch
variable mu index -21.0
variable disp index 1.0
variable mu index -1.25
variable temp index 2.0
variable lbox index 10.0
variable disp index 1.0
variable lbox index 5.0
# global model settings
units lj
atom_style atomic
pair_style lj/cut 3.0
pair_modify tail yes
pair_style lj/cut 3.0
pair_modify tail no # turn of to avoid triggering full_energy
# box
@ -28,15 +33,27 @@ mass * 1.0
fix mygcmc all gcmc 1 100 100 1 29494 ${temp} ${mu} ${disp}
# averaging
variable rho equal density
variable p equal press
variable nugget equal 1.0e-8
variable lambda equal 1.0
variable muex equal ${mu}-${temp}*ln(density*${lambda}+${nugget})
fix ave all ave/time 10 100 1000 v_rho v_p v_muex ave one file rho_vs_p.dat
variable rhoav equal f_ave[1]
variable pav equal f_ave[2]
variable muexav equal f_ave[3]
# output
variable tacc equal f_mygcmc[2]/(f_mygcmc[1]+0.1)
variable iacc equal f_mygcmc[4]/(f_mygcmc[3]+0.1)
variable dacc equal f_mygcmc[6]/(f_mygcmc[5]+0.1)
variable tacc equal f_mygcmc[2]/(f_mygcmc[1]+${nugget})
variable iacc equal f_mygcmc[4]/(f_mygcmc[3]+${nugget})
variable dacc equal f_mygcmc[6]/(f_mygcmc[5]+${nugget})
compute_modify thermo_temp dynamic yes
thermo_style custom step temp press pe ke density atoms v_iacc v_dacc v_tacc
thermo 100
thermo_style custom step temp press pe ke density atoms v_iacc v_dacc v_tacc v_rhoav v_pav v_muexav
thermo 1000
# run
run 1000
run 10000

View File

@ -1,28 +1,35 @@
LAMMPS (17 Mar 2017)
OMP_NUM_THREADS environment is not set. Defaulting to 1 thread. (../comm.cpp:90)
using 1 OpenMP thread(s) per MPI task
# GCMC for LJ simple fluid, no dynamics
# T = 2.0
# rho ~ 0.5
# p ~ 1.5
# mu_ex ~ 0.0
# comparable to Frenkel and Smit GCMC Case Study, Figure 5.8
# variables available on command line
# variables modifiable using -var command line switch
variable mu index -21.0
variable disp index 1.0
variable mu index -1.25
variable temp index 2.0
variable lbox index 10.0
variable disp index 1.0
variable lbox index 5.0
# global model settings
units lj
atom_style atomic
pair_style lj/cut 3.0
pair_modify tail yes
pair_modify tail no # turn of to avoid triggering full_energy
# box
region box block 0 ${lbox} 0 ${lbox} 0 ${lbox}
region box block 0 10.0 0 ${lbox} 0 ${lbox}
region box block 0 10.0 0 10.0 0 ${lbox}
region box block 0 10.0 0 10.0 0 10.0
region box block 0 5.0 0 ${lbox} 0 ${lbox}
region box block 0 5.0 0 5.0 0 ${lbox}
region box block 0 5.0 0 5.0 0 5.0
create_box 1 box
Created orthogonal box = (0 0 0) to (10 10 10)
Created orthogonal box = (0 0 0) to (5 5 5)
1 by 1 by 1 MPI processor grid
# lj parameters
@ -34,70 +41,89 @@ mass * 1.0
fix mygcmc all gcmc 1 100 100 1 29494 ${temp} ${mu} ${disp}
fix mygcmc all gcmc 1 100 100 1 29494 2.0 ${mu} ${disp}
fix mygcmc all gcmc 1 100 100 1 29494 2.0 -21.0 ${disp}
fix mygcmc all gcmc 1 100 100 1 29494 2.0 -21.0 1.0
fix mygcmc all gcmc 1 100 100 1 29494 2.0 -1.25 ${disp}
fix mygcmc all gcmc 1 100 100 1 29494 2.0 -1.25 1.0
# averaging
variable rho equal density
variable p equal press
variable nugget equal 1.0e-8
variable lambda equal 1.0
variable muex equal ${mu}-${temp}*ln(density*${lambda}+${nugget})
variable muex equal -1.25-${temp}*ln(density*${lambda}+${nugget})
variable muex equal -1.25-2.0*ln(density*${lambda}+${nugget})
variable muex equal -1.25-2.0*ln(density*1+${nugget})
variable muex equal -1.25-2.0*ln(density*1+1e-08)
fix ave all ave/time 10 100 1000 v_rho v_p v_muex ave one file rho_vs_p.dat
variable rhoav equal f_ave[1]
variable pav equal f_ave[2]
variable muexav equal f_ave[3]
# output
variable tacc equal f_mygcmc[2]/(f_mygcmc[1]+0.1)
variable iacc equal f_mygcmc[4]/(f_mygcmc[3]+0.1)
variable dacc equal f_mygcmc[6]/(f_mygcmc[5]+0.1)
variable tacc equal f_mygcmc[2]/(f_mygcmc[1]+${nugget})
variable tacc equal f_mygcmc[2]/(f_mygcmc[1]+1e-08)
variable iacc equal f_mygcmc[4]/(f_mygcmc[3]+${nugget})
variable iacc equal f_mygcmc[4]/(f_mygcmc[3]+1e-08)
variable dacc equal f_mygcmc[6]/(f_mygcmc[5]+${nugget})
variable dacc equal f_mygcmc[6]/(f_mygcmc[5]+1e-08)
compute_modify thermo_temp dynamic yes
thermo_style custom step temp press pe ke density atoms v_iacc v_dacc v_tacc
thermo 100
thermo_style custom step temp press pe ke density atoms v_iacc v_dacc v_tacc v_rhoav v_pav v_muexav
thermo 1000
# run
run 1000
run 10000
Neighbor list info ...
update every 1 steps, delay 10 steps, check yes
max neighbors/atom: 2000, page size: 100000
master list distance cutoff = 3.3
ghost atom cutoff = 3.3
binsize = 1.65, bins = 7 7 7
binsize = 1.65, bins = 4 4 4
1 neighbor lists, perpetual/occasional/extra = 1 0 0
(1) pair lj/cut, perpetual
attributes: half, newton on
pair build: half/bin/atomonly/newton
stencil: half/bin/3d/newton
bin: standard
Per MPI rank memory allocation (min/avg/max) = 0.4369 | 0.4369 | 0.4369 Mbytes
Step Temp Press PotEng KinEng Density Atoms v_iacc v_dacc v_tacc
0 0 0 0 -0 0 0 0 0 0
100 1.9042848 0.39026453 -1.7692765 2.8466449 0.292 292 0.3619855 0.30247792 0.40278761
200 1.8651924 0.47815517 -1.8494955 2.7886155 0.305 305 0.34021109 0.30357196 0.37759189
300 2.0626994 0.52068504 -1.8197295 3.0834166 0.291 291 0.32055605 0.3003043 0.36103862
400 2.0394818 0.53751435 -1.7636699 3.0482184 0.278 278 0.31698808 0.29995864 0.35441275
500 1.9628066 0.54594742 -1.7145336 2.9339513 0.287 287 0.31211861 0.29724228 0.35161407
600 1.9845913 0.40846162 -1.8199325 2.9669308 0.299 299 0.30976643 0.29612711 0.34933559
700 1.8582606 0.53445462 -1.7869306 2.777974 0.296 296 0.30642103 0.29446478 0.34633665
800 2.0340641 0.66057698 -1.7075279 3.0403148 0.283 283 0.30730979 0.29746793 0.34768045
900 2.0830765 0.63731971 -1.894775 3.114911 0.322 322 0.30636338 0.29737705 0.34737644
1000 1.9688933 0.50024802 -1.7013944 2.9428299 0.281 281 0.3053174 0.29772245 0.34788254
Loop time of 3.98286 on 1 procs for 1000 steps with 281 atoms
Per MPI rank memory allocation (min/avg/max) = 0.433 | 0.433 | 0.433 Mbytes
Step Temp Press PotEng KinEng Density Atoms v_iacc v_dacc v_tacc v_rhoav v_pav v_muexav
0 0 0 0 -0 0 0 0 0 0 0 0 0
1000 2.4038954 2.1735585 -2.7041368 3.5476844 0.496 62 0.064790036 0.06313096 0.1081294 0.54304 1.4513524 -0.025479219
2000 2.0461168 1.1913842 -2.9880181 3.0212194 0.512 64 0.067416408 0.066335853 0.11306166 0.52736 1.3274665 0.034690004
3000 1.7930436 1.3788681 -3.2212667 2.6505861 0.552 69 0.067733191 0.066877836 0.1133516 0.5344 1.3834744 0.0070582537
4000 1.981449 1.2541054 -2.8222868 2.9217977 0.472 59 0.068546991 0.067856412 0.11442807 0.52504 1.3815629 0.043309657
5000 2.0946818 1.0701629 -3.5213291 3.0977688 0.568 71 0.06813743 0.067567891 0.11342906 0.53824 1.4049567 -0.0054539777
6000 1.9793484 0.68224187 -3.410211 2.9247088 0.536 67 0.067797628 0.067420108 0.11295333 0.5384 1.401683 -0.0066894359
7000 2.1885798 1.6745012 -3.185499 3.2345922 0.544 68 0.068630201 0.068261832 0.11403705 0.5244 1.449239 0.045987399
8000 2.2175324 1.5897263 -3.078898 3.2759002 0.528 66 0.068180395 0.067899629 0.11332691 0.53928 1.5488388 -0.01075766
9000 1.8610779 1.0396231 -2.923262 2.7465908 0.496 62 0.068346453 0.068028117 0.1134132 0.52912 1.4352871 0.027082544
10000 2.1079271 1.1746643 -2.9112062 3.1091925 0.48 60 0.068352878 0.068054948 0.11335434 0.5316 1.4462327 0.018503094
Loop time of 13.05 on 1 procs for 10000 steps with 60 atoms
Performance: 108464.750 tau/day, 251.076 timesteps/s
99.9% CPU use with 1 MPI tasks x no OpenMP threads
Performance: 331035.016 tau/day, 766.285 timesteps/s
100.0% CPU use with 1 MPI tasks x 1 OpenMP threads
MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total
---------------------------------------------------------------
Pair | 0.10563 | 0.10563 | 0.10563 | 0.0 | 2.65
Neigh | 0.33428 | 0.33428 | 0.33428 | 0.0 | 8.39
Comm | 0.027969 | 0.027969 | 0.027969 | 0.0 | 0.70
Output | 0.00017285 | 0.00017285 | 0.00017285 | 0.0 | 0.00
Modify | 3.5096 | 3.5096 | 3.5096 | 0.0 | 88.12
Other | | 0.005197 | | | 0.13
Pair | 0.37239 | 0.37239 | 0.37239 | 0.0 | 2.85
Neigh | 0.94764 | 0.94764 | 0.94764 | 0.0 | 7.26
Comm | 0.092473 | 0.092473 | 0.092473 | 0.0 | 0.71
Output | 0.00023365 | 0.00023365 | 0.00023365 | 0.0 | 0.00
Modify | 11.627 | 11.627 | 11.627 | 0.0 | 89.09
Other | | 0.01054 | | | 0.08
Nlocal: 281 ave 281 max 281 min
Nlocal: 60 ave 60 max 60 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Nghost: 977 ave 977 max 977 min
Nghost: 663 ave 663 max 663 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Neighs: 5902 ave 5902 max 5902 min
Neighs: 2133 ave 2133 max 2133 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Total # of neighbors = 5902
Ave neighs/atom = 21.0036
Neighbor list builds = 1000
Total # of neighbors = 2133
Ave neighs/atom = 35.55
Neighbor list builds = 10000
Dangerous builds = 0
Total wall time: 0:00:03
Total wall time: 0:00:13

View File

@ -1,28 +1,35 @@
LAMMPS (17 Mar 2017)
OMP_NUM_THREADS environment is not set. Defaulting to 1 thread. (../comm.cpp:90)
using 1 OpenMP thread(s) per MPI task
# GCMC for LJ simple fluid, no dynamics
# T = 2.0
# rho ~ 0.5
# p ~ 1.5
# mu_ex ~ 0.0
# comparable to Frenkel and Smit GCMC Case Study, Figure 5.8
# variables available on command line
# variables modifiable using -var command line switch
variable mu index -21.0
variable disp index 1.0
variable mu index -1.25
variable temp index 2.0
variable lbox index 10.0
variable disp index 1.0
variable lbox index 5.0
# global model settings
units lj
atom_style atomic
pair_style lj/cut 3.0
pair_modify tail yes
pair_modify tail no # turn of to avoid triggering full_energy
# box
region box block 0 ${lbox} 0 ${lbox} 0 ${lbox}
region box block 0 10.0 0 ${lbox} 0 ${lbox}
region box block 0 10.0 0 10.0 0 ${lbox}
region box block 0 10.0 0 10.0 0 10.0
region box block 0 5.0 0 ${lbox} 0 ${lbox}
region box block 0 5.0 0 5.0 0 ${lbox}
region box block 0 5.0 0 5.0 0 5.0
create_box 1 box
Created orthogonal box = (0 0 0) to (10 10 10)
Created orthogonal box = (0 0 0) to (5 5 5)
1 by 2 by 2 MPI processor grid
# lj parameters
@ -34,70 +41,89 @@ mass * 1.0
fix mygcmc all gcmc 1 100 100 1 29494 ${temp} ${mu} ${disp}
fix mygcmc all gcmc 1 100 100 1 29494 2.0 ${mu} ${disp}
fix mygcmc all gcmc 1 100 100 1 29494 2.0 -21.0 ${disp}
fix mygcmc all gcmc 1 100 100 1 29494 2.0 -21.0 1.0
fix mygcmc all gcmc 1 100 100 1 29494 2.0 -1.25 ${disp}
fix mygcmc all gcmc 1 100 100 1 29494 2.0 -1.25 1.0
# averaging
variable rho equal density
variable p equal press
variable nugget equal 1.0e-8
variable lambda equal 1.0
variable muex equal ${mu}-${temp}*ln(density*${lambda}+${nugget})
variable muex equal -1.25-${temp}*ln(density*${lambda}+${nugget})
variable muex equal -1.25-2.0*ln(density*${lambda}+${nugget})
variable muex equal -1.25-2.0*ln(density*1+${nugget})
variable muex equal -1.25-2.0*ln(density*1+1e-08)
fix ave all ave/time 10 100 1000 v_rho v_p v_muex ave one file rho_vs_p.dat
variable rhoav equal f_ave[1]
variable pav equal f_ave[2]
variable muexav equal f_ave[3]
# output
variable tacc equal f_mygcmc[2]/(f_mygcmc[1]+0.1)
variable iacc equal f_mygcmc[4]/(f_mygcmc[3]+0.1)
variable dacc equal f_mygcmc[6]/(f_mygcmc[5]+0.1)
variable tacc equal f_mygcmc[2]/(f_mygcmc[1]+${nugget})
variable tacc equal f_mygcmc[2]/(f_mygcmc[1]+1e-08)
variable iacc equal f_mygcmc[4]/(f_mygcmc[3]+${nugget})
variable iacc equal f_mygcmc[4]/(f_mygcmc[3]+1e-08)
variable dacc equal f_mygcmc[6]/(f_mygcmc[5]+${nugget})
variable dacc equal f_mygcmc[6]/(f_mygcmc[5]+1e-08)
compute_modify thermo_temp dynamic yes
thermo_style custom step temp press pe ke density atoms v_iacc v_dacc v_tacc
thermo 100
thermo_style custom step temp press pe ke density atoms v_iacc v_dacc v_tacc v_rhoav v_pav v_muexav
thermo 1000
# run
run 1000
run 10000
Neighbor list info ...
update every 1 steps, delay 10 steps, check yes
max neighbors/atom: 2000, page size: 100000
master list distance cutoff = 3.3
ghost atom cutoff = 3.3
binsize = 1.65, bins = 7 7 7
binsize = 1.65, bins = 4 4 4
1 neighbor lists, perpetual/occasional/extra = 1 0 0
(1) pair lj/cut, perpetual
attributes: half, newton on
pair build: half/bin/atomonly/newton
stencil: half/bin/3d/newton
bin: standard
Per MPI rank memory allocation (min/avg/max) = 0.434 | 0.434 | 0.434 Mbytes
Step Temp Press PotEng KinEng Density Atoms v_iacc v_dacc v_tacc
0 0 0 0 -0 0 0 0 0 0
100 2.0328045 0.58661762 -1.6812724 3.0385824 0.287 287 0.35917318 0.30067507 0.38663622
200 1.9594279 0.50682399 -1.7308396 2.9287927 0.284 284 0.33788365 0.30337335 0.37300293
300 2.0602937 0.7028247 -1.9278541 3.0806296 0.315 315 0.31882007 0.29697498 0.36167185
400 1.995183 0.4328246 -1.8715454 2.983026 0.307 307 0.31527654 0.29681901 0.35673374
500 2.1390101 0.48232215 -1.554138 3.1960306 0.257 257 0.31372975 0.30003067 0.35558858
600 2.0584244 0.4929049 -1.6995569 3.0767263 0.283 283 0.31114213 0.29801665 0.35160109
700 1.9155066 0.49654243 -1.5770611 2.8624174 0.265 265 0.31056419 0.29944173 0.35157337
800 2.0883562 0.52731947 -1.8261112 3.1220925 0.3 300 0.30730979 0.29704354 0.34898892
900 2.0470677 0.5605993 -2.0130053 3.0610656 0.322 322 0.30484441 0.29586719 0.34678883
1000 2.004135 0.50642204 -1.6956257 2.9955798 0.283 283 0.30396929 0.29634309 0.34770304
Loop time of 3.688 on 4 procs for 1000 steps with 283 atoms
Per MPI rank memory allocation (min/avg/max) = 0.4477 | 0.4477 | 0.4477 Mbytes
Step Temp Press PotEng KinEng Density Atoms v_iacc v_dacc v_tacc v_rhoav v_pav v_muexav
0 0 0 0 -0 0 0 0 0 0 0 0 0
1000 1.956397 1.7699101 -2.7889468 2.8864874 0.488 61 0.068894746 0.067229075 0.1141726 0.53288 1.3832798 0.013392866
2000 2.040943 0.56060899 -2.8001647 3.0077055 0.456 57 0.069858594 0.068831934 0.11629114 0.5232 1.3587174 0.049995794
3000 2.0004866 1.5736515 -3.3098044 2.9572411 0.552 69 0.069594029 0.068727791 0.11592543 0.53096 1.4129434 0.020022578
4000 2.1127942 2.642809 -2.8865084 3.1211733 0.528 66 0.070268697 0.069533235 0.11693806 0.52424 1.3444615 0.046884078
5000 2.3663648 1.354269 -3.1917346 3.4957662 0.528 66 0.070519633 0.069960064 0.11710321 0.52688 1.3595814 0.036270867
6000 1.9224136 0.82756699 -3.1965 2.839257 0.52 65 0.06985018 0.069474645 0.11628632 0.536 1.47062 0.00141549
7000 2.0266192 1.5593811 -2.9972341 2.9931606 0.52 65 0.070244693 0.069880791 0.11666541 0.52528 1.3246332 0.040754793
8000 1.7790467 1.8680568 -2.8028819 2.6275151 0.52 65 0.070454494 0.070172368 0.11736806 0.524 1.4213649 0.047985191
9000 1.7968847 1.3195587 -3.261001 2.6550983 0.536 67 0.069952011 0.069618327 0.11650087 0.53904 1.4624201 -0.01069837
10000 2.1566109 1.1015729 -3.4999837 3.1880335 0.552 69 0.069603309 0.069284134 0.11625548 0.53128 1.3587249 0.02075238
Loop time of 13.0611 on 4 procs for 10000 steps with 69 atoms
Performance: 117136.751 tau/day, 271.150 timesteps/s
99.2% CPU use with 4 MPI tasks x no OpenMP threads
Performance: 330753.007 tau/day, 765.632 timesteps/s
99.7% CPU use with 4 MPI tasks x 1 OpenMP threads
MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total
---------------------------------------------------------------
Pair | 0.024644 | 0.026027 | 0.027483 | 0.6 | 0.71
Neigh | 0.085449 | 0.088998 | 0.092893 | 0.9 | 2.41
Comm | 0.045756 | 0.051296 | 0.056578 | 1.7 | 1.39
Output | 0.00028491 | 0.00030857 | 0.00035262 | 0.0 | 0.01
Modify | 3.5189 | 3.5191 | 3.5194 | 0.0 | 95.42
Other | | 0.002221 | | | 0.06
Pair | 0.08888 | 0.09443 | 0.099889 | 1.4 | 0.72
Neigh | 0.27721 | 0.28532 | 0.29177 | 1.1 | 2.18
Comm | 0.27648 | 0.28875 | 0.30268 | 1.9 | 2.21
Output | 0.00029635 | 0.00043058 | 0.00048113 | 0.0 | 0.00
Modify | 12.384 | 12.384 | 12.384 | 0.0 | 94.82
Other | | 0.008055 | | | 0.06
Nlocal: 70.75 ave 77 max 68 min
Histogram: 1 2 0 0 0 0 0 0 0 1
Nghost: 514.25 ave 520 max 507 min
Histogram: 1 0 0 0 1 0 0 1 0 1
Neighs: 1483.5 ave 1715 max 1359 min
Histogram: 2 0 0 1 0 0 0 0 0 1
Nlocal: 17.25 ave 23 max 10 min
Histogram: 1 0 0 0 0 0 2 0 0 1
Nghost: 506.5 ave 519 max 490 min
Histogram: 1 0 1 0 0 0 0 0 0 2
Neighs: 705.75 ave 998 max 369 min
Histogram: 1 0 0 0 0 1 1 0 0 1
Total # of neighbors = 5934
Ave neighs/atom = 20.9682
Neighbor list builds = 1000
Total # of neighbors = 2823
Ave neighs/atom = 40.913
Neighbor list builds = 10000
Dangerous builds = 0
Total wall time: 0:00:03
Total wall time: 0:00:13

0
potentials/InP.vashishta Executable file → Normal file
View File

0
potentials/SiC.vashishta Executable file → Normal file
View File

0
potentials/SiO.1990.vashishta Executable file → Normal file
View File

0
potentials/SiO.1994.vashishta Executable file → Normal file
View File

0
potentials/SiO.1997.vashishta Executable file → Normal file
View File

0
potentials/ffield.smtbq.Al Executable file → Normal file
View File

0
potentials/ffield.smtbq.Al2O3 Executable file → Normal file
View File

0
potentials/ffield.smtbq.TiO2 Executable file → Normal file
View File

7
src/.gitignore vendored
View File

@ -76,6 +76,13 @@
/pair_awpmd_cut.cpp
/pair_awpmd_cut.h
/dihedral_charmmfsh.cpp
/dihedral_charmmfsh.h
/pair_lj_charmmfsw_coul_charmmfsh.cpp
/pair_lj_charmmfsw_coul_charmmfsh.h
/pair_lj_charmmfsw_coul_long.cpp
/pair_lj_charmmfsw_coul_long.h
/angle_cg_cmm.cpp
/angle_cg_cmm.h
/angle_charmm.cpp

View File

@ -131,9 +131,8 @@ void FixWallGranRegion::init()
void FixWallGranRegion::post_force(int vflag)
{
int i,m,nc,iwall;
double rinv,fx,fy,fz,tooclose;
double dx,dy,dz,rsq,meff;
double xc[3],vwall[3];
double vwall[3];
// do not update shear history during setup

View File

@ -3158,7 +3158,7 @@ double PPPMKokkos<DeviceType>::memory_usage()
if (peratom_allocate_flag)
bytes += 6 * nbrick * sizeof(FFT_SCALAR);
bytes += cg->memory_usage();
if (cg) bytes += cg->memory_usage();
return bytes;
}

View File

@ -3085,7 +3085,7 @@ double PPPM::memory_usage()
bytes += 2 * nfft_both * sizeof(FFT_SCALAR);;
}
bytes += cg->memory_usage();
if (cg) bytes += cg->memory_usage();
return bytes;
}

View File

@ -8240,7 +8240,7 @@ double PPPMDisp::memory_usage()
bytes += 6 * nfft_both * sizeof(double); // vg
bytes += nfft_both * sizeof(double); // greensfn
bytes += nfft_both * 3 * sizeof(FFT_SCALAR); // density_FFT, work1, work2
bytes += cg->memory_usage();
if (cg) bytes += cg->memory_usage();
}
if (function[1] + function[2] + function[3]) {
@ -8250,7 +8250,7 @@ double PPPMDisp::memory_usage()
bytes += 6 * nfft_both_6 * sizeof(double); // vg
bytes += nfft_both_6 * sizeof(double); // greensfn
bytes += nfft_both_6 * (mixing + 2) * sizeof(FFT_SCALAR); // density_FFT, work1, work2
bytes += cg_6->memory_usage();
if (cg_6) bytes += cg_6->memory_usage();
}
return bytes;
}

View File

@ -72,7 +72,7 @@ FixGCMC::FixGCMC(LAMMPS *lmp, int narg, char **arg) :
Fix(lmp, narg, arg),
idregion(NULL), full_flag(0), ngroups(0), groupstrings(NULL), ngrouptypes(0), grouptypestrings(NULL),
grouptypebits(NULL), grouptypes(NULL), local_gas_list(NULL), atom_coord(NULL), random_equal(NULL), random_unequal(NULL),
coords(NULL), imageflags(NULL), idrigid(NULL), idshake(NULL), fixrigid(NULL), fixshake(NULL)
coords(NULL), imageflags(NULL), fixrigid(NULL), fixshake(NULL), idrigid(NULL), idshake(NULL)
{
if (narg < 11) error->all(FLERR,"Illegal fix gcmc command");
@ -432,7 +432,8 @@ void FixGCMC::init()
(force->pair == NULL) ||
(force->pair->single_enable == 0) ||
(force->pair_match("hybrid",0)) ||
(force->pair_match("eam",0))
(force->pair_match("eam",0)) ||
(force->pair->tail_flag)
) {
full_flag = true;
if (comm->me == 0)
@ -618,13 +619,18 @@ void FixGCMC::init()
}
// compute beta, lambda, sigma, and the zz factor
// For LJ units, lambda=1
beta = 1.0/(force->boltz*reservoir_temperature);
double lambda = sqrt(force->hplanck*force->hplanck/
(2.0*MY_PI*gas_mass*force->mvv2e*
if (strcmp(update->unit_style,"lj") == 0)
zz = exp(beta*chemical_potential);
else {
double lambda = sqrt(force->hplanck*force->hplanck/
(2.0*MY_PI*gas_mass*force->mvv2e*
force->boltz*reservoir_temperature));
zz = exp(beta*chemical_potential)/(pow(lambda,3.0));
}
sigma = sqrt(force->boltz*reservoir_temperature*tfac_insert/gas_mass/force->mvv2e);
zz = exp(beta*chemical_potential)/(pow(lambda,3.0));
if (pressure_flag) zz = pressure*fugacity_coeff*beta/force->nktv2p;
imagezero = ((imageint) IMGMAX << IMG2BITS) |

View File

@ -964,22 +964,22 @@ void FixCMAP::bc_coeff(double *gs, double *d1gs, double *d2gs, double *d12gs)
// calculate the bicubic interpolation coefficients c_ij
static int wt[16][16] =
{ 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0,
-3, 0, 0, 3, 0, 0, 0, 0,-2, 0, 0,-1, 0, 0, 0, 0,
2, 0, 0,-2, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0,
0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0,
0, 0, 0, 0,-3, 0, 0, 3, 0, 0, 0, 0,-2, 0, 0,-1,
0, 0, 0, 0, 2, 0, 0,-2, 0, 0, 0, 0, 1, 0, 0, 1,
-3, 3, 0, 0,-2,-1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0,-3, 3, 0, 0,-2,-1, 0, 0,
9,-9, 9,-9, 6, 3,-3,-6, 6,-6,-3, 3, 4, 2, 1, 2,
-6, 6,-6, 6,-4,-2, 2, 4,-3, 3, 3,-3,-2,-1,-1,-2,
2,-2, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 2,-2, 0, 0, 1, 1, 0, 0,
-6, 6,-6, 6,-3,-3, 3, 3,-4, 4, 2,-2,-2,-2,-1,-1,
4,-4, 4,-4, 2, 2,-2,-2, 2,-2,-2, 2, 1, 1, 1, 1
{ {1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
{0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0},
{-3, 0, 0, 3, 0, 0, 0, 0,-2, 0, 0,-1, 0, 0, 0, 0},
{2, 0, 0,-2, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0},
{0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0},
{0, 0, 0, 0,-3, 0, 0, 3, 0, 0, 0, 0,-2, 0, 0,-1},
{0, 0, 0, 0, 2, 0, 0,-2, 0, 0, 0, 0, 1, 0, 0, 1},
{-3, 3, 0, 0,-2,-1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
{0, 0, 0, 0, 0, 0, 0, 0,-3, 3, 0, 0,-2,-1, 0, 0},
{9,-9, 9,-9, 6, 3,-3,-6, 6,-6,-3, 3, 4, 2, 1, 2},
{-6, 6,-6, 6,-4,-2, 2, 4,-3, 3, 3,-3,-2,-1,-1,-2},
{2,-2, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
{0, 0, 0, 0, 0, 0, 0, 0, 2,-2, 0, 0, 1, 1, 0, 0},
{-6, 6,-6, 6,-3,-3, 3, 3,-4, 4, 2,-2,-2,-2,-1,-1},
{4,-4, 4,-4, 2, 2,-2,-2, 2,-2,-2, 2, 1, 1, 1, 1}
};
int i, j, k, in;

View File

@ -67,11 +67,6 @@ void PairEAMOpt::eval()
double rhor0j,rhor1j,rhor2j,rhor3j;
} fast_alpha_t;
typedef struct {
double frho0,frho1,frho2,frho3,frho4,frho5,frho6;
double _pad[1];
} fast_beta_t;
typedef struct {
double rhor4i,rhor5i,rhor6i;
double rhor4j,rhor5j,rhor6j;

View File

@ -249,7 +249,7 @@ double FixQEqDynamic::compute_eneg()
int FixQEqDynamic::pack_forward_comm(int n, int *list, double *buf,
int pbc_flag, int *pbc)
{
int m;
int m=0;
if( pack_flag == 1 )
for(m = 0; m < n; m++) buf[m] = atom->q[list[m]];

View File

@ -0,0 +1,49 @@
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
Contributing author: Oliver Henrich (University of Strathclyde, Glasgow)
------------------------------------------------------------------------- */
#include <math.h>
#include <stdlib.h>
#include "bond_oxdna2_fene.h"
using namespace LAMMPS_NS;
/* ---------------------------------------------------------------------- */
BondOxdna2Fene::BondOxdna2Fene(LAMMPS *lmp) : BondOxdnaFene(lmp)
{
}
/* ---------------------------------------------------------------------- */
BondOxdna2Fene::~BondOxdna2Fene()
{
}
/* ----------------------------------------------------------------------
compute vector COM-sugar-phosphate backbone interaction site in oxDNA2
------------------------------------------------------------------------- */
void BondOxdna2Fene::compute_interaction_sites(double e1[3],
double e2[3], double r[3])
{
double d_cs_x=-0.34, d_cs_y=+0.3408;
r[0] = d_cs_x*e1[0] + d_cs_y*e2[0];
r[1] = d_cs_x*e1[1] + d_cs_y*e2[1];
r[2] = d_cs_x*e1[2] + d_cs_y*e2[2];
}

View File

@ -0,0 +1,64 @@
/* -*- c++ -*- ----------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#ifdef BOND_CLASS
BondStyle(oxdna2/fene,BondOxdna2Fene)
#else
#ifndef LMP_BOND_OXDNA2_FENE_H
#define LMP_BOND_OXDNA2_FENE_H
#include "bond_oxdna_fene.h"
namespace LAMMPS_NS {
class BondOxdna2Fene : public BondOxdnaFene {
public:
BondOxdna2Fene(class LAMMPS *);
virtual ~BondOxdna2Fene();
virtual void compute_interaction_sites(double *, double *, double *);
};
}
#endif
#endif
/* ERROR/WARNING messages:
W: FENE bond too long: %ld %d %d %g
A FENE bond has stretched dangerously far. It's interaction strength
will be truncated to attempt to prevent the bond from blowing up.
E: Bad FENE bond
Two atoms in a FENE bond have become so far apart that the bond cannot
be computed.
E: Incorrect args for bond coefficients
Self-explanatory. Check the input script or data file.
W: Use special bonds = 0,1,1 with bond style oxdna
Most FENE models need this setting for the special_bonds command.
W: FENE bond too long: %ld %g
A FENE bond has stretched dangerously far. It's interaction strength
will be truncated to attempt to prevent the bond from blowing up.
*/

View File

@ -11,7 +11,7 @@
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
Contributing author: Oliver Henrich (EPCC, University of Edinburgh)
Contributing author: Oliver Henrich (University of Strathclyde, Glasgow)
------------------------------------------------------------------------- */
#include <math.h>
@ -51,6 +51,21 @@ BondOxdnaFene::~BondOxdnaFene()
}
}
/* ----------------------------------------------------------------------
compute vector COM-sugar-phosphate backbone interaction site in oxDNA
------------------------------------------------------------------------- */
void BondOxdnaFene::compute_interaction_sites(double e1[3],
double e2[3], double r[3])
{
double d_cs=-0.4;
r[0] = d_cs*e1[0];
r[1] = d_cs*e1[1];
r[2] = d_cs*e1[2];
}
/* ----------------------------------------------------------------------
compute function for oxDNA FENE-bond interaction
s=sugar-phosphate backbone site, b=base site, st=stacking site
@ -62,8 +77,6 @@ void BondOxdnaFene::compute(int eflag, int vflag)
double delr[3],ebond,fbond;
double rsq,Deltasq,rlogarg;
double r,rr0,rr0sq;
// distances COM-backbone site
double d_cs=-0.24;
// vectors COM-backbone site in lab frame
double ra_cs[3],rb_cs[3];
@ -100,12 +113,8 @@ void BondOxdnaFene::compute(int eflag, int vflag)
MathExtra::q_to_exyz(qb,bx,by,bz);
// vector COM-backbone site a and b
ra_cs[0] = d_cs*ax[0];
ra_cs[1] = d_cs*ax[1];
ra_cs[2] = d_cs*ax[2];
rb_cs[0] = d_cs*bx[0];
rb_cs[1] = d_cs*bx[1];
rb_cs[2] = d_cs*bx[2];
compute_interaction_sites(ax,ay,ra_cs);
compute_interaction_sites(bx,by,rb_cs);
// vector backbone site b to a
delr[0] = x[a][0] + ra_cs[0] - x[b][0] - rb_cs[0];
@ -202,7 +211,7 @@ 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 != 4) error->all(FLERR,"Incorrect args for bond coefficients in oxdna/fene");
if (!allocated) allocate();
int ilo,ihi;
@ -222,7 +231,7 @@ void BondOxdnaFene::coeff(int narg, char **arg)
count++;
}
if (count == 0) error->all(FLERR,"Incorrect args for bond coefficients in oxdna_fene");
if (count == 0) error->all(FLERR,"Incorrect args for bond coefficients in oxdna/fene");
}
@ -252,7 +261,7 @@ void BondOxdnaFene::init_style()
force->special_coul[1] != 1.0 || force->special_coul[2] != 1.0 || force->special_coul[3] != 1.0)
{
if (comm->me == 0)
error->warning(FLERR,"Use special bonds lj = 0,1,1 and coul = 1,1,1 with bond style oxdna_fene");
error->warning(FLERR,"Use special bonds lj = 0,1,1 and coul = 1,1,1 with bond style oxdna/fene");
}
}

View File

@ -10,13 +10,10 @@
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
Contributing author: Oliver Henrich (EPCC, University of Edinburgh)
------------------------------------------------------------------------- */
#ifdef BOND_CLASS
BondStyle(oxdna_fene,BondOxdnaFene)
BondStyle(oxdna/fene,BondOxdnaFene)
#else
@ -31,6 +28,7 @@ class BondOxdnaFene : public Bond {
public:
BondOxdnaFene(class LAMMPS *);
virtual ~BondOxdnaFene();
virtual void compute_interaction_sites(double *, double *, double *);
virtual void compute(int, int);
void coeff(int, char **);
void init_style();

View File

@ -10,9 +10,8 @@
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
Contributing author: Oliver Henrich (EPCC, University of Edinburgh)
Contributing author: Oliver Henrich (University of Strathclyde, Glasgow)
------------------------------------------------------------------------- */
#include <math.h>

View File

@ -12,7 +12,7 @@
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
Contributing author: Oliver Henrich (EPCC, University of Edinburgh)
Contributing author: Oliver Henrich (University of Strathclyde, Glasgow)
------------------------------------------------------------------------- */
#include <math.h>

View File

@ -10,9 +10,6 @@
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
Contributing author: Oliver Henrich (EPCC, University of Edinburgh)
------------------------------------------------------------------------- */
#ifndef MF_OXDNA_H
#define MF_OXDNA_H
@ -31,6 +28,8 @@ namespace MFOxdna {
inline double DF4(double, double, double, double, double, double);
inline double F5(double, double, double, double, double);
inline double DF5(double, double, double, double, double);
inline double F6(double, double, double);
inline double DF6(double, double, double);
inline double is_3pto5p(const double *, const double *);
}
@ -252,6 +251,32 @@ inline double MFOxdna::DF5(double x, double a, double x_ast,
return 0;
}
/* ----------------------------------------------------------------------
f6 modulation factor
------------------------------------------------------------------------- */
inline double MFOxdna::F6(double theta, double a, double b)
{
if (theta < b) {
return 0.0;
}
else {
return 0.5 * a * (theta-b)*(theta-b);
}
}
/* ----------------------------------------------------------------------
derivative of f6 modulation factor
------------------------------------------------------------------------- */
inline double MFOxdna::DF6(double theta, double a, double b)
{
if (theta < b) {
return 0.0;
}
else {
return a * (theta-b);
}
}
/* ----------------------------------------------------------------------
test for directionality by projecting base normal n onto delr = a - b,
returns 1 if nucleotide b to nucleotide a is 3' to 5', otherwise -1

File diff suppressed because it is too large Load Diff

View File

@ -0,0 +1,85 @@
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#ifdef PAIR_CLASS
PairStyle(oxdna2/coaxstk,PairOxdna2Coaxstk)
#else
#ifndef LMP_PAIR_OXDNA2_COAXSTK_H
#define LMP_PAIR_OXDNA2_COAXSTK_H
#include "pair.h"
namespace LAMMPS_NS {
class PairOxdna2Coaxstk : public Pair {
public:
PairOxdna2Coaxstk(class LAMMPS *);
virtual ~PairOxdna2Coaxstk();
virtual void compute(int, int);
void settings(int, char **);
void coeff(int, char **);
void init_style();
void init_list(int, class NeighList *);
double init_one(int, int);
void write_restart(FILE *);
void read_restart(FILE *);
void write_restart_settings(FILE *);
void read_restart_settings(FILE *);
void write_data(FILE *);
void write_data_all(FILE *);
void *extract(const char *, int &);
protected:
// coaxial stacking interaction
double **k_cxst, **cut_cxst_0, **cut_cxst_c, **cut_cxst_lo, **cut_cxst_hi;
double **cut_cxst_lc, **cut_cxst_hc, **b_cxst_lo, **b_cxst_hi;
double **cutsq_cxst_hc;
double **a_cxst1, **theta_cxst1_0, **dtheta_cxst1_ast;
double **b_cxst1, **dtheta_cxst1_c;
double **a_cxst4, **theta_cxst4_0, **dtheta_cxst4_ast;
double **b_cxst4, **dtheta_cxst4_c;
double **a_cxst5, **theta_cxst5_0, **dtheta_cxst5_ast;
double **b_cxst5, **dtheta_cxst5_c;
double **a_cxst6, **theta_cxst6_0, **dtheta_cxst6_ast;
double **b_cxst6, **dtheta_cxst6_c;
double **AA_cxst1, **BB_cxst1;
virtual void allocate();
};
}
#endif
#endif
/* ERROR/WARNING messages:
E: Illegal ... command
Self-explanatory. Check the input script syntax and compare to the
documentation for the command. You can use -echo screen as a
command-line option when running LAMMPS to see the offending line.
E: Incorrect args for pair coefficients
Self-explanatory. Check the input script or data file.
*/

View File

@ -0,0 +1,553 @@
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
Contributing author: Oliver Henrich (University of Strathclyde, Glasgow)
------------------------------------------------------------------------- */
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include "pair_oxdna2_dh.h"
#include "mf_oxdna.h"
#include "atom.h"
#include "comm.h"
#include "force.h"
#include "neighbor.h"
#include "neigh_list.h"
#include "neigh_request.h"
#include "update.h"
#include "integrate.h"
#include "math_const.h"
#include "memory.h"
#include "error.h"
#include "atom_vec_ellipsoid.h"
#include "math_extra.h"
using namespace LAMMPS_NS;
using namespace MathConst;
using namespace MFOxdna;
/* ---------------------------------------------------------------------- */
PairOxdna2Dh::PairOxdna2Dh(LAMMPS *lmp) : Pair(lmp)
{
single_enable = 0;
writedata = 1;
}
/* ---------------------------------------------------------------------- */
PairOxdna2Dh::~PairOxdna2Dh()
{
if (allocated) {
memory->destroy(setflag);
memory->destroy(cutsq);
memory->destroy(qeff_dh_pf);
memory->destroy(kappa_dh);
memory->destroy(b_dh);
memory->destroy(cut_dh_ast);
memory->destroy(cutsq_dh_ast);
memory->destroy(cut_dh_c);
memory->destroy(cutsq_dh_c);
}
}
/* ----------------------------------------------------------------------
compute vector COM-sugar-phosphate backbone interaction site in oxDNA2
------------------------------------------------------------------------- */
void PairOxdna2Dh::compute_interaction_sites(double e1[3],
double e2[3], double r[3])
{
double d_cs_x=-0.34, d_cs_y=+0.3408;
r[0] = d_cs_x*e1[0] + d_cs_y*e2[0];
r[1] = d_cs_x*e1[1] + d_cs_y*e2[1];
r[2] = d_cs_x*e1[2] + d_cs_y*e2[2];
}
/* ----------------------------------------------------------------------
compute function for oxDNA pair interactions
s=sugar-phosphate backbone site, b=base site, st=stacking site
------------------------------------------------------------------------- */
void PairOxdna2Dh::compute(int eflag, int vflag)
{
double delf[3],delta[3],deltb[3]; // force, torque increment;;
double rtmp_s[3],delr[3];
double evdwl,fpair,factor_lj;
double r,rsq,rinv;
// vectors COM-backbone sites in lab frame
double ra_cs[3],rb_cs[3];
// quaternions and Cartesian unit vectors in lab frame
double *qa,ax[3],ay[3],az[3];
double *qb,bx[3],by[3],bz[3];
double *special_lj = force->special_lj;
double **x = atom->x;
double **f = atom->f;
double **torque = atom->torque;
int *type = atom->type;
int nlocal = atom->nlocal;
int newton_pair = force->newton_pair;
int *alist,*blist,*numneigh,**firstneigh;
AtomVecEllipsoid *avec = (AtomVecEllipsoid *) atom->style_match("ellipsoid");
AtomVecEllipsoid::Bonus *bonus = avec->bonus;
int a,b,ia,ib,anum,bnum,atype,btype;
evdwl = 0.0;
if (eflag || vflag) ev_setup(eflag,vflag);
else evflag = vflag_fdotr = 0;
anum = list->inum;
alist = list->ilist;
numneigh = list->numneigh;
firstneigh = list->firstneigh;
// loop over pair interaction neighbours of my atoms
for (ia = 0; ia < anum; ia++) {
a = alist[ia];
atype = type[a];
qa=bonus[a].quat;
MathExtra::q_to_exyz(qa,ax,ay,az);
// vector COM-backbone site a
compute_interaction_sites(ax,ay,ra_cs);
rtmp_s[0] = x[a][0] + ra_cs[0];
rtmp_s[1] = x[a][1] + ra_cs[1];
rtmp_s[2] = x[a][2] + ra_cs[2];
blist = firstneigh[a];
bnum = numneigh[a];
for (ib = 0; ib < bnum; ib++) {
b = blist[ib];
factor_lj = special_lj[sbmask(b)]; // = 0 for nearest neighbours
b &= NEIGHMASK;
btype = type[b];
qb=bonus[b].quat;
MathExtra::q_to_exyz(qb,bx,by,bz);
// vector COM-backbone site b
compute_interaction_sites(bx,by,rb_cs);
// vector backbone site b to a
delr[0] = rtmp_s[0] - x[b][0] - rb_cs[0];
delr[1] = rtmp_s[1] - x[b][1] - rb_cs[1];
delr[2] = rtmp_s[2] - x[b][2] - rb_cs[2];
rsq = delr[0]*delr[0] + delr[1]*delr[1] + delr[2]*delr[2];
if (rsq <= cutsq_dh_c[atype][btype]) {
r = sqrt(rsq);
rinv = 1.0/r;
if (r <= cut_dh_ast[atype][btype]) {
fpair = qeff_dh_pf[atype][btype] * exp(-kappa_dh[atype][btype] * r) *
(kappa_dh[atype][btype] + rinv) * rinv * rinv;
if (eflag) {
evdwl = qeff_dh_pf[atype][btype] * exp(-kappa_dh[atype][btype]*r) * rinv;
}
}
else {
fpair = 2.0 * b_dh[atype][btype] * (cut_dh_c[atype][btype] - r) * rinv;
if (eflag) {
evdwl = b_dh[atype][btype] * (r - cut_dh_c[atype][btype]) * (r - cut_dh_c[atype][btype]);
}
}
// knock out nearest-neighbour interaction between adjacent backbone sites
fpair *= factor_lj;
evdwl *= factor_lj;
delf[0] = delr[0] * fpair;
delf[1] = delr[1] * fpair;
delf[2] = delr[2] * fpair;
// apply force and torque to each of 2 atoms
if (newton_pair || a < nlocal) {
f[a][0] += delf[0];
f[a][1] += delf[1];
f[a][2] += delf[2];
MathExtra::cross3(ra_cs,delf,delta);
torque[a][0] += delta[0];
torque[a][1] += delta[1];
torque[a][2] += delta[2];
}
if (newton_pair || b < nlocal) {
f[b][0] -= delf[0];
f[b][1] -= delf[1];
f[b][2] -= delf[2];
MathExtra::cross3(rb_cs,delf,deltb);
torque[b][0] -= deltb[0];
torque[b][1] -= deltb[1];
torque[b][2] -= deltb[2];
}
// increment energy and virial
if (evflag) ev_tally(a,b,nlocal,newton_pair,
evdwl,0.0,fpair,delr[0],delr[1],delr[2]);
}
}
}
}
/* ----------------------------------------------------------------------
allocate all arrays
------------------------------------------------------------------------- */
void PairOxdna2Dh::allocate()
{
allocated = 1;
int n = atom->ntypes;
memory->create(setflag,n+1,n+1,"pair:setflag");
for (int i = 1; i <= n; i++)
for (int j = i; j <= n; j++)
setflag[i][j] = 0;
memory->create(cutsq,n+1,n+1,"pair:cutsq");
memory->create(kappa_dh,n+1,n+1,"pair:kappa_dh");
memory->create(qeff_dh_pf,n+1,n+1,"pair:qeff_dh_pf");
memory->create(b_dh,n+1,n+1,"pair:b_dh");
memory->create(cut_dh_ast,n+1,n+1,"pair:cut_dh_ast");
memory->create(cutsq_dh_ast,n+1,n+1,"pair:cutsq_dh_ast");
memory->create(cut_dh_c,n+1,n+1,"pair:cut_dh_c");
memory->create(cutsq_dh_c,n+1,n+1,"pair:cutsq_dh_c");
}
/* ----------------------------------------------------------------------
global settings
------------------------------------------------------------------------- */
void PairOxdna2Dh::settings(int narg, char **arg)
{
if (narg != 0) error->all(FLERR,"Illegal pair_style command");
}
/* ----------------------------------------------------------------------
set coeffs for one or more type pairs
------------------------------------------------------------------------- */
void PairOxdna2Dh::coeff(int narg, char **arg)
{
int count;
if (narg != 5) error->all(FLERR,"Incorrect args for pair coefficients in oxdna/dh");
if (!allocated) allocate();
int ilo,ihi,jlo,jhi;
force->bounds(FLERR,arg[0],atom->ntypes,ilo,ihi);
force->bounds(FLERR,arg[1],atom->ntypes,jlo,jhi);
count = 0;
double T, rhos_dh_one, qeff_dh_one;
T = force->numeric(FLERR,arg[2]);
rhos_dh_one = force->numeric(FLERR,arg[3]);
qeff_dh_one = force->numeric(FLERR,arg[4]);
double lambda_dh_one, kappa_dh_one, qeff_dh_pf_one;
double b_dh_one, cut_dh_ast_one, cut_dh_c_one;
// Debye length and inverse Debye length
/*
NOTE:
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
(see B. Snodin et al., J. Chem. Phys. 142, 234901 (2015).)
We use
eps_0 = vacuum permittivity = 8.854187817e-12 F/m
eps_r = relative permittivity of water = 80
k_B = Boltzmann constant = 1.3806485279e-23 J/K
T = absolute temperature = 300 K
N_A = Avogadro constant = 6.02214085774e23 / mol
e = elementary charge = 1.6021766208e-19 C
oxDNA_length_unit = 8.518e-10 m
*/
lambda_dh_one = 0.3616455075438555*sqrt(T/0.1/rhos_dh_one);
kappa_dh_one = 1.0/lambda_dh_one;
// prefactor in DH interaction containing qeff^2
/*
NOTE:
The numerical factor is
qeff_dh_pf = e^2/(4 * pi * eps_0 * eps_r)
* 1/(oxDNA_energy_unit * oxDNA_length_unit)
(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;
// smoothing parameters - determined through continuity and differentiability
cut_dh_ast_one = 3.0*lambda_dh_one;
b_dh_one = -(exp(-cut_dh_ast_one/lambda_dh_one) * qeff_dh_pf_one * qeff_dh_pf_one *
(cut_dh_ast_one + lambda_dh_one) * (cut_dh_ast_one + lambda_dh_one))/
(-4.0 * cut_dh_ast_one * cut_dh_ast_one * cut_dh_ast_one *
lambda_dh_one * lambda_dh_one * qeff_dh_pf_one);
cut_dh_c_one = cut_dh_ast_one * (qeff_dh_pf_one*cut_dh_ast_one +
3.0*qeff_dh_pf_one * lambda_dh_one)/
(qeff_dh_pf_one * (cut_dh_ast_one+lambda_dh_one));
for (int i = ilo; i <= ihi; i++) {
for (int j = MAX(jlo,i); j <= jhi; j++) {
kappa_dh[i][j] = kappa_dh_one;
qeff_dh_pf[i][j] = qeff_dh_pf_one;
b_dh[i][j] = b_dh_one;
cut_dh_ast[i][j] = cut_dh_ast_one;
cut_dh_c[i][j] = cut_dh_c_one;
setflag[i][j] = 1;
count++;
}
}
if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients in oxdna/dh");
}
/* ----------------------------------------------------------------------
init specific to this pair style
------------------------------------------------------------------------- */
void PairOxdna2Dh::init_style()
{
int irequest;
// request regular neighbor lists
irequest = neighbor->request(this,instance_me);
}
/* ----------------------------------------------------------------------
neighbor callback to inform pair style of neighbor list to use regular
------------------------------------------------------------------------- */
void PairOxdna2Dh::init_list(int id, NeighList *ptr)
{
if (id == 0) list = ptr;
if (id > 0) error->all(FLERR,"Respa not supported");
}
/* ----------------------------------------------------------------------
init for one type pair i,j and corresponding j,i
------------------------------------------------------------------------- */
double PairOxdna2Dh::init_one(int i, int j)
{
if (setflag[i][j] == 0) {
error->all(FLERR,"Coefficient mixing not defined in oxDNA");
}
if (offset_flag) {
error->all(FLERR,"Offset not supported in oxDNA");
}
kappa_dh[j][i] = kappa_dh[i][j];
qeff_dh_pf[j][i] = qeff_dh_pf[i][j];
b_dh[j][i] = b_dh[i][j];
cut_dh_ast[j][i] = cut_dh_ast[i][j];
cut_dh_c[j][i] = cut_dh_c[i][j];
cutsq_dh_ast[i][j] = cut_dh_ast[i][j]*cut_dh_ast[i][j];
cutsq_dh_ast[j][i] = cutsq_dh_ast[i][j];
cutsq_dh_c[i][j] = cut_dh_c[i][j]*cut_dh_c[i][j];
cutsq_dh_c[j][i] = cutsq_dh_c[i][j];
// set the master list distance cutoff
return cut_dh_c[i][j];
}
/* ----------------------------------------------------------------------
proc 0 writes to restart file
------------------------------------------------------------------------- */
void PairOxdna2Dh::write_restart(FILE *fp)
{
write_restart_settings(fp);
int i,j;
for (i = 1; i <= atom->ntypes; i++)
for (j = i; j <= atom->ntypes; j++) {
fwrite(&setflag[i][j],sizeof(int),1,fp);
if (setflag[i][j]) {
fwrite(&kappa_dh[i][j],sizeof(double),1,fp);
fwrite(&qeff_dh_pf[i][j],sizeof(double),1,fp);
fwrite(&b_dh[i][j],sizeof(double),1,fp);
fwrite(&cut_dh_ast[i][j],sizeof(double),1,fp);
fwrite(&cut_dh_c[i][j],sizeof(double),1,fp);
}
}
}
/* ----------------------------------------------------------------------
proc 0 reads from restart file, bcasts
------------------------------------------------------------------------- */
void PairOxdna2Dh::read_restart(FILE *fp)
{
read_restart_settings(fp);
allocate();
int i,j;
int me = comm->me;
for (i = 1; i <= atom->ntypes; i++)
for (j = i; j <= atom->ntypes; j++) {
if (me == 0) fread(&setflag[i][j],sizeof(int),1,fp);
MPI_Bcast(&setflag[i][j],1,MPI_INT,0,world);
if (setflag[i][j]) {
if (me == 0) {
fread(&kappa_dh[i][j],sizeof(double),1,fp);
fread(&qeff_dh_pf[i][j],sizeof(double),1,fp);
fread(&b_dh[i][j],sizeof(double),1,fp);
fread(&cut_dh_ast[i][j],sizeof(double),1,fp);
fread(&cut_dh_c[i][j],sizeof(double),1,fp);
}
MPI_Bcast(&kappa_dh[i][j],1,MPI_DOUBLE,0,world);
MPI_Bcast(&qeff_dh_pf[i][j],1,MPI_DOUBLE,0,world);
MPI_Bcast(&b_dh[i][j],1,MPI_DOUBLE,0,world);
MPI_Bcast(&cut_dh_ast[i][j],1,MPI_DOUBLE,0,world);
MPI_Bcast(&cut_dh_c[i][j],1,MPI_DOUBLE,0,world);
}
}
}
/* ----------------------------------------------------------------------
proc 0 writes to restart file
------------------------------------------------------------------------- */
void PairOxdna2Dh::write_restart_settings(FILE *fp)
{
fwrite(&offset_flag,sizeof(int),1,fp);
fwrite(&mix_flag,sizeof(int),1,fp);
fwrite(&tail_flag,sizeof(int),1,fp);
}
/* ----------------------------------------------------------------------
proc 0 reads from restart file, bcasts
------------------------------------------------------------------------- */
void PairOxdna2Dh::read_restart_settings(FILE *fp)
{
int me = comm->me;
if (me == 0) {
fread(&offset_flag,sizeof(int),1,fp);
fread(&mix_flag,sizeof(int),1,fp);
fread(&tail_flag,sizeof(int),1,fp);
}
MPI_Bcast(&offset_flag,1,MPI_INT,0,world);
MPI_Bcast(&mix_flag,1,MPI_INT,0,world);
MPI_Bcast(&tail_flag,1,MPI_INT,0,world);
}
/* ----------------------------------------------------------------------
proc 0 writes to data file
------------------------------------------------------------------------- */
void PairOxdna2Dh::write_data(FILE *fp)
{
for (int i = 1; i <= atom->ntypes; i++)
fprintf(fp,"%d\
%g %g\
%g %g %g\
\n",i,
kappa_dh[i][i],qeff_dh_pf[i][i],
b_dh[i][i],cut_dh_ast[i][i],cut_dh_c[i][i]);
}
/* ----------------------------------------------------------------------
proc 0 writes all pairs to data file
------------------------------------------------------------------------- */
void PairOxdna2Dh::write_data_all(FILE *fp)
{
for (int i = 1; i <= atom->ntypes; i++)
for (int j = i; j <= atom->ntypes; j++)
fprintf(fp,"%d %d\
%g %g\
%g %g %g\
\n",i,j,
kappa_dh[i][j],qeff_dh_pf[i][j],
b_dh[i][j],cut_dh_ast[i][j],cut_dh_c[i][j]);
}
/* ---------------------------------------------------------------------- */
void *PairOxdna2Dh::extract(const char *str, int &dim)
{
dim = 2;
if (strcmp(str,"kappa_dh") == 0) return (void *) kappa_dh;
if (strcmp(str,"qeff_dh_pf") == 0) return (void *) qeff_dh_pf;
if (strcmp(str,"b_dh") == 0) return (void *) b_dh;
if (strcmp(str,"cut_dh_ast") == 0) return (void *) cut_dh_ast;
if (strcmp(str,"cut_dh_c") == 0) return (void *) cut_dh_c;
return NULL;
}

View File

@ -0,0 +1,72 @@
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#ifdef PAIR_CLASS
PairStyle(oxdna2/dh,PairOxdna2Dh)
#else
#ifndef LMP_PAIR_OXDNA2_DH_H
#define LMP_PAIR_OXDNA2_DH_H
#include "pair.h"
namespace LAMMPS_NS {
class PairOxdna2Dh : public Pair {
public:
PairOxdna2Dh(class LAMMPS *);
virtual ~PairOxdna2Dh();
virtual void compute_interaction_sites(double *, double *, double *);
virtual void compute(int, int);
void settings(int, char **);
void coeff(int, char **);
void init_style();
void init_list(int, class NeighList *);
double init_one(int, int);
void write_restart(FILE *);
void read_restart(FILE *);
void write_restart_settings(FILE *);
void read_restart_settings(FILE *);
void write_data(FILE *);
void write_data_all(FILE *);
void *extract(const char *, int &);
protected:
double **qeff_dh_pf,**kappa_dh;
double **b_dh,**cut_dh_ast,**cutsq_dh_ast,**cut_dh_c,**cutsq_dh_c;
virtual void allocate();
};
}
#endif
#endif
/* ERROR/WARNING messages:
E: Illegal ... command
Self-explanatory. Check the input script syntax and compare to the
documentation for the command. You can use -echo screen as a
command-line option when running LAMMPS to see the offending line.
E: Incorrect args for pair coefficients
Self-explanatory. Check the input script or data file.
*/

View File

@ -0,0 +1,55 @@
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
Contributing author: Oliver Henrich (University of Strathclyde, Glasgow)
------------------------------------------------------------------------- */
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include "pair_oxdna2_excv.h"
using namespace LAMMPS_NS;
/* ---------------------------------------------------------------------- */
PairOxdna2Excv::PairOxdna2Excv(LAMMPS *lmp) : PairOxdnaExcv(lmp)
{
}
/* ---------------------------------------------------------------------- */
PairOxdna2Excv::~PairOxdna2Excv()
{
}
/* ----------------------------------------------------------------------
compute vector COM-excluded volume interaction sites in oxDNA2
------------------------------------------------------------------------- */
void PairOxdna2Excv::compute_interaction_sites(double e1[3],
double e2[3], double rs[3], double rb[3])
{
double d_cs_x=-0.34, d_cs_y=+0.3408, d_cb=+0.4;
rs[0] = d_cs_x*e1[0] + d_cs_y*e2[0];
rs[1] = d_cs_x*e1[1] + d_cs_y*e2[1];
rs[2] = d_cs_x*e1[2] + d_cs_y*e2[2];
rb[0] = d_cb*e1[0];
rb[1] = d_cb*e1[1];
rb[2] = d_cb*e1[2];
}

View File

@ -0,0 +1,52 @@
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#ifdef PAIR_CLASS
PairStyle(oxdna2/excv,PairOxdna2Excv)
#else
#ifndef LMP_PAIR_OXDNA2_EXCV_H
#define LMP_PAIR_OXDNA2_EXCV_H
#include "pair_oxdna_excv.h"
namespace LAMMPS_NS {
class PairOxdna2Excv : public PairOxdnaExcv {
public:
PairOxdna2Excv(class LAMMPS *);
virtual ~PairOxdna2Excv();
virtual void compute_interaction_sites(double *,
double *, double *, double *);
};
}
#endif
#endif
/* ERROR/WARNING messages:
E: Illegal ... command
Self-explanatory. Check the input script syntax and compare to the
documentation for the command. You can use -echo screen as a
command-line option when running LAMMPS to see the offending line.
E: Incorrect args for pair coefficients
Self-explanatory. Check the input script or data file.
*/

View File

@ -0,0 +1,50 @@
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
Contributing author: Oliver Henrich (University of Strathclyde, Glasgow)
------------------------------------------------------------------------- */
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include "pair_oxdna2_stk.h"
using namespace LAMMPS_NS;
/* ---------------------------------------------------------------------- */
PairOxdna2Stk::PairOxdna2Stk(LAMMPS *lmp) : PairOxdnaStk(lmp)
{
}
/* ---------------------------------------------------------------------- */
PairOxdna2Stk::~PairOxdna2Stk()
{
}
/* ----------------------------------------------------------------------
return temperature dependent oxDNA2 stacking strength
------------------------------------------------------------------------- */
double PairOxdna2Stk::stacking_strength(double T)
{
double eps;
eps = 1.3523 + 2.6717 * T;
return eps;
}

View File

@ -0,0 +1,53 @@
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#ifdef PAIR_CLASS
PairStyle(oxdna2/stk,PairOxdna2Stk)
#else
#ifndef LMP_PAIR_OXDNA2_STK_H
#define LMP_PAIR_OXDNA2_STK_H
#include "pair_oxdna_stk.h"
namespace LAMMPS_NS {
class PairOxdna2Stk : public PairOxdnaStk {
public:
PairOxdna2Stk(class LAMMPS *);
virtual ~PairOxdna2Stk();
protected:
virtual double stacking_strength(double);
};
}
#endif
#endif
/* ERROR/WARNING messages:
E: Illegal ... command
Self-explanatory. Check the input script syntax and compare to the
documentation for the command. You can use -echo screen as a
command-line option when running LAMMPS to see the offending line.
E: Incorrect args for pair coefficients
Self-explanatory. Check the input script or data file.
*/

View File

@ -11,7 +11,7 @@
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
Contributing author: Oliver Henrich (EPCC, University of Edinburgh)
Contributing author: Oliver Henrich (University of Strathclyde, Glasgow)
------------------------------------------------------------------------- */
#include <math.h>
@ -127,7 +127,7 @@ void PairOxdnaCoaxstk::compute(int eflag, int vflag)
double dcdrax,dcdray,dcdraz;
// distances COM-backbone site, COM-stacking site
double d_cs=-0.24, d_cst=0.5;
double d_cs=-0.4, d_cst=+0.34;
// vectors COM-backbone site, COM-stacking site in lab frame
double ra_cs[3],ra_cst[3];
double rb_cs[3],rb_cst[3];
@ -680,7 +680,7 @@ void PairOxdnaCoaxstk::coeff(int narg, char **arg)
{
int count;
if (narg != 23) error->all(FLERR,"Incorrect args for pair coefficients in oxdna_xstack");
if (narg != 23) error->all(FLERR,"Incorrect args for pair coefficients in oxdna/coaxstk");
if (!allocated) allocate();
int ilo,ihi,jlo,jhi;
@ -818,7 +818,7 @@ void PairOxdnaCoaxstk::coeff(int narg, char **arg)
}
}
if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients in oxdna_xstack");
if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients in oxdna/coaxstk");
}

View File

@ -10,13 +10,10 @@
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
Contributing author: Oliver Henrich (EPCC, University of Edinburgh)
------------------------------------------------------------------------- */
#ifdef PAIR_CLASS
PairStyle(oxdna_coaxstk,PairOxdnaCoaxstk)
PairStyle(oxdna/coaxstk,PairOxdnaCoaxstk)
#else

View File

@ -11,7 +11,7 @@
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
Contributing author: Oliver Henrich (EPCC, University of Edinburgh)
Contributing author: Oliver Henrich (University of Strathclyde, Glasgow)
------------------------------------------------------------------------- */
#include <math.h>
@ -88,6 +88,24 @@ PairOxdnaExcv::~PairOxdnaExcv()
}
}
/* ----------------------------------------------------------------------
compute vector COM-excluded volume interaction sites in oxDNA
------------------------------------------------------------------------- */
void PairOxdnaExcv::compute_interaction_sites(double e1[3],
double e2[3], double rs[3], double rb[3])
{
double d_cs=-0.4, d_cb=+0.4;
rs[0] = d_cs*e1[0];
rs[1] = d_cs*e1[1];
rs[2] = d_cs*e1[2];
rb[0] = d_cb*e1[0];
rb[1] = d_cb*e1[1];
rb[2] = d_cb*e1[2];
}
/* ----------------------------------------------------------------------
compute function for oxDNA pair interactions
s=sugar-phosphate backbone site, b=base site, st=stacking site
@ -102,8 +120,6 @@ void PairOxdnaExcv::compute(int eflag, int vflag)
double delr_ss[3],rsq_ss,delr_sb[3],rsq_sb;
double delr_bs[3],rsq_bs,delr_bb[3],rsq_bb;
// distances COM-backbone site, COM-base site
double d_cs=-0.24, d_cb=0.56;
// vectors COM-backbone site, COM-base site in lab frame
double ra_cs[3],ra_cb[3];
double rb_cs[3],rb_cb[3];
@ -146,18 +162,13 @@ void PairOxdnaExcv::compute(int eflag, int vflag)
qa=bonus[a].quat;
MathExtra::q_to_exyz(qa,ax,ay,az);
// position of backbone site a
ra_cs[0] = d_cs*ax[0];
ra_cs[1] = d_cs*ax[1];
ra_cs[2] = d_cs*ax[2];
// vector COM - backbone and base site a
compute_interaction_sites(ax,ay,ra_cs,ra_cb);
rtmp_s[0] = x[a][0] + ra_cs[0];
rtmp_s[1] = x[a][1] + ra_cs[1];
rtmp_s[2] = x[a][2] + ra_cs[2];
// position of base site a
ra_cb[0] = d_cb*ax[0];
ra_cb[1] = d_cb*ax[1];
ra_cb[2] = d_cb*ax[2];
rtmp_b[0] = x[a][0] + ra_cb[0];
rtmp_b[1] = x[a][1] + ra_cb[1];
rtmp_b[2] = x[a][2] + ra_cb[2];
@ -176,12 +187,8 @@ void PairOxdnaExcv::compute(int eflag, int vflag)
qb=bonus[b].quat;
MathExtra::q_to_exyz(qb,bx,by,bz);
rb_cs[0] = d_cs*bx[0];
rb_cs[1] = d_cs*bx[1];
rb_cs[2] = d_cs*bx[2];
rb_cb[0] = d_cb*bx[0];
rb_cb[1] = d_cb*bx[1];
rb_cb[2] = d_cb*bx[2];
// vector COM - backbone and base site b
compute_interaction_sites(bx,by,rb_cs,rb_cb);
// vector backbone site b to a
delr_ss[0] = rtmp_s[0] - (x[b][0] + rb_cs[0]);
@ -448,7 +455,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 != 11) error->all(FLERR,"Incorrect args for pair coefficients in oxdna/excv");
if (!allocated) allocate();
int ilo,ihi,jlo,jhi;
@ -494,7 +501,7 @@ void PairOxdnaExcv::coeff(int narg, char **arg)
}
}
if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients in oxdna_excv");
if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients in oxdna/excv");
count = 0;
@ -525,7 +532,7 @@ void PairOxdnaExcv::coeff(int narg, char **arg)
}
}
if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients in oxdna_excv");
if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients in oxdna/excv");
count = 0;
@ -556,7 +563,7 @@ void PairOxdnaExcv::coeff(int narg, char **arg)
}
}
if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients in oxdna_excv");
if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients in oxdna/excv");
}
@ -657,7 +664,7 @@ double PairOxdnaExcv::init_one(int i, int j)
cutsq_bb_c[j][i] = cutsq_bb_c[i][j];
// set the master list distance cutoff
return cut_ss_ast[i][j];
return cut_ss_c[i][j];
}

View File

@ -10,13 +10,10 @@
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
Contributing author: Oliver Henrich (EPCC, University of Edinburgh)
------------------------------------------------------------------------- */
#ifdef PAIR_CLASS
PairStyle(oxdna_excv,PairOxdnaExcv)
PairStyle(oxdna/excv,PairOxdnaExcv)
#else
@ -31,6 +28,8 @@ class PairOxdnaExcv : public Pair {
public:
PairOxdnaExcv(class LAMMPS *);
virtual ~PairOxdnaExcv();
virtual void compute_interaction_sites(double *, double *,
double *, double *);
virtual void compute(int, int);
void settings(int, char **);
void coeff(int, char **);

View File

@ -11,7 +11,7 @@
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
Contributing author: Oliver Henrich (EPCC, University of Edinburgh)
Contributing author: Oliver Henrich (University of Strathclyde, Glasgow)
------------------------------------------------------------------------- */
#include <math.h>
@ -125,7 +125,7 @@ void PairOxdnaHbond::compute(int eflag, int vflag)
double theta8,t8dir[3],cost8;
// distance COM-hbonding site
double d_chb=0.56;
double d_chb=+0.4;
// vectors COM-h-bonding site in lab frame
double ra_chb[3],rb_chb[3];
@ -607,7 +607,7 @@ void PairOxdnaHbond::coeff(int narg, char **arg)
{
int count;
if (narg != 26) error->all(FLERR,"Incorrect args for pair coefficients in oxdna_hbond");
if (narg != 26) error->all(FLERR,"Incorrect args for pair coefficients in oxdna/hbond");
if (!allocated) allocate();
int ilo,ihi,jlo,jhi;
@ -770,7 +770,7 @@ void PairOxdnaHbond::coeff(int narg, char **arg)
}
}
if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients in oxdna_hbond");
if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients in oxdna/hbond");
}

View File

@ -10,13 +10,11 @@
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
Contributing author: Oliver Henrich (EPCC, University of Edinburgh)
------------------------------------------------------------------------- */
#ifdef PAIR_CLASS
PairStyle(oxdna_hbond,PairOxdnaHbond)
PairStyle(oxdna/hbond,PairOxdnaHbond)
PairStyle(oxdna2/hbond,PairOxdnaHbond)
#else

View File

@ -11,7 +11,7 @@
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
Contributing author: Oliver Henrich (EPCC, University of Edinburgh)
Contributing author: Oliver Henrich (University of Strathclyde, Glasgow)
------------------------------------------------------------------------- */
#include <math.h>
@ -116,7 +116,7 @@ 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.24, d_cst=0.5;
double d_cs=-0.4, d_cst=+0.34;
// vectors COM-backbone site, COM-stacking site in lab frame
double ra_cs[3],ra_cst[3];
double rb_cs[3],rb_cst[3];
@ -644,6 +644,19 @@ void PairOxdnaStk::settings(int narg, char **arg)
}
/* ----------------------------------------------------------------------
return temperature dependent oxDNA stacking strength
------------------------------------------------------------------------- */
double PairOxdnaStk::stacking_strength(double T)
{
double eps;
eps = 1.3448 + 2.6568 * T;
return eps;
}
/* ----------------------------------------------------------------------
set coeffs for one or more type pairs
------------------------------------------------------------------------- */
@ -652,7 +665,7 @@ void PairOxdnaStk::coeff(int narg, char **arg)
{
int count;
if (narg != 21) error->all(FLERR,"Incorrect args for pair coefficients in oxdna_stk");
if (narg != 21) error->all(FLERR,"Incorrect args for pair coefficients in oxdna/stk");
if (!allocated) allocate();
int ilo,ihi,jlo,jhi;
@ -662,7 +675,7 @@ void PairOxdnaStk::coeff(int narg, char **arg)
// stacking interaction
count = 0;
double epsilon_st_one, a_st_one, b_st_lo_one, b_st_hi_one;
double T, epsilon_st_one, a_st_one, b_st_lo_one, b_st_hi_one;
double cut_st_0_one, cut_st_c_one, cut_st_lo_one, cut_st_hi_one;
double cut_st_lc_one, cut_st_hc_one, tmp, shift_st_one;
@ -678,7 +691,9 @@ void PairOxdnaStk::coeff(int narg, char **arg)
double a_st1_one, cosphi_st1_ast_one, b_st1_one, cosphi_st1_c_one;
double a_st2_one, cosphi_st2_ast_one, b_st2_one, cosphi_st2_c_one;
epsilon_st_one = force->numeric(FLERR,arg[2]);
T = force->numeric(FLERR,arg[2]);
epsilon_st_one = stacking_strength(T);
a_st_one = force->numeric(FLERR,arg[3]);
cut_st_0_one = force->numeric(FLERR,arg[4]);
cut_st_c_one = force->numeric(FLERR,arg[5]);
@ -790,7 +805,7 @@ void PairOxdnaStk::coeff(int narg, char **arg)
}
}
if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients in oxdna_stk");
if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients in oxdna/stk");
}

View File

@ -10,13 +10,10 @@
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
Contributing author: Oliver Henrich (EPCC, University of Edinburgh)
------------------------------------------------------------------------- */
#ifdef PAIR_CLASS
PairStyle(oxdna_stk,PairOxdnaStk)
PairStyle(oxdna/stk,PairOxdnaStk)
#else
@ -47,6 +44,7 @@ class PairOxdnaStk : public Pair {
protected:
// stacking interaction
virtual double stacking_strength(double);
double **epsilon_st, **a_st, **cut_st_0, **cut_st_c;
double **cut_st_lo, **cut_st_hi;
double **cut_st_lc, **cut_st_hc, **b_st_lo, **b_st_hi, **shift_st;

View File

@ -11,7 +11,7 @@
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
Contributing author: Oliver Henrich (EPCC, University of Edinburgh)
Contributing author: Oliver Henrich (University of Strathclyde, Glasgow)
------------------------------------------------------------------------- */
#include <math.h>
@ -125,7 +125,7 @@ void PairOxdnaXstk::compute(int eflag, int vflag)
double theta8,theta8p,t8dir[3],cost8;
// distance COM-h-bonding site
double d_chb=0.56;
double d_chb=+0.4;
// vectors COM-h-bonding site in lab frame
double ra_chb[3],rb_chb[3];
@ -625,7 +625,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 != 25) error->all(FLERR,"Incorrect args for pair coefficients in oxdna/xstk");
if (!allocated) allocate();
int ilo,ihi,jlo,jhi;
@ -772,7 +772,7 @@ void PairOxdnaXstk::coeff(int narg, char **arg)
}
}
if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients in oxdna_xstk");
if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients in oxdna/xstk");
}

View File

@ -10,13 +10,11 @@
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
Contributing author: Oliver Henrich (EPCC, University of Edinburgh)
------------------------------------------------------------------------- */
#ifdef PAIR_CLASS
PairStyle(oxdna_xstk,PairOxdnaXstk)
PairStyle(oxdna/xstk,PairOxdnaXstk)
PairStyle(oxdna2/xstk,PairOxdnaXstk)
#else

View File

@ -666,7 +666,6 @@ void FixRX::setup_pre_force(int vflag)
int *mask = atom->mask;
int newton_pair = force->newton_pair;
double tmp;
int ii;
if(restartFlag){
restartFlag = 0;

View File

@ -64,7 +64,6 @@ void NPairHalfBinNewtonSSA::build(NeighList *list)
tagint **special = atom->special;
int **nspecial = atom->nspecial;
int nlocal = atom->nlocal;
int nall = nlocal + atom->nghost;
if (includegroup) nlocal = atom->nfirst;
int *ssaAIR = atom->ssaAIR;

View File

@ -640,10 +640,12 @@ double PairLJCutTholeLong::single(int i, int j, int itype, int jtype,
if (j != di_closest){
if (drudetype[i] == CORE_TYPE) dqi = -atom->q[di];
else if (drudetype[i] == DRUDE_TYPE) dqi = atom->q[i];
else dqi = 0.0;
if (drudetype[j] == CORE_TYPE) {
dj = atom->map(drudeid[j]);
dqj = -atom->q[dj];
} else if (drudetype[j] == DRUDE_TYPE) dqj = atom->q[j];
else dqj = 0.0;
asr = ascreen[itype][jtype] * r;
exp_asr = exp(-asr);
dcoul = force->qqrd2e * dqi * dqj / r;

Some files were not shown because too many files have changed in this diff Show More