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

Merge remote-tracking branch 'lammps-ro/master' into lammps-icms

Browse Source 
# Resolved Conflicts:
#	doc/Manual.txt
This commit is contained in:
Axel Kohlmeyer
2016-01-15 08:18:31 -05:00
parent 7e1972ef8b ba87515112
commit a066cbb210
26 changed files with 937 additions and 289 deletions
Download Patch File Download Diff File Expand all files Collapse all files

2
doc/Manual.html
View File

@ -135,7 +135,7 @@
<H1></H1><div class="section" id="lammps-documentation">
<h1>LAMMPS-ICMS Documentation<a class="headerlink" href="#lammps-documentation" title="Permalink to this headline"></a></h1>
<div class="section" id="jan-2016-version">
<h2>13 Jan 2016 version<a class="headerlink" href="#jan-2016-version" title="Permalink to this headline"></a></h2>
<h2>14 Jan 2016 version<a class="headerlink" href="#jan-2016-version" title="Permalink to this headline"></a></h2>
</div>
<div class="section" id="version-info">
<h2>Version info:<a class="headerlink" href="#version-info" title="Permalink to this headline"></a></h2>

4
doc/Manual.txt
View File

@ -1,7 +1,7 @@
<!-- HTML_ONLY -->
<HEAD>
<TITLE>LAMMPS-ICMS Users Manual</TITLE>
<META NAME="docnumber" CONTENT="13 Jan 2016 version">
<META NAME="docnumber" CONTENT="14 Jan 2016 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-ICMS Documentation :c,h3
13 Jan 2016 version :c,h4
14 Jan 2016 version :c,h4
Version info: :h4

6
doc/displace_atoms.html
View File

@ -208,8 +208,12 @@ same, regardless of how many processors are being used.</p>
around a rotation axis <em>R</em> = (Rx,Ry,Rz) that goes thru a point <em>P</em> =
(Px,Py,Pz). The direction of rotation for the atoms around the
rotation axis is consistent with the right-hand rule: if your
right-hand&#8217;s thumb points along <em>R</em>, then your fingers wrap around the
right-hand thumb points along <em>R</em>, then your fingers wrap around the
axis in the direction of positive theta.</p>
<p>If the defined <a class="reference internal" href="atom_style.html"><em>atom_style</em></a> assigns an orientation to
each atom (<a class="reference internal" href="atom_style.html"><em>atom styles</em></a> ellipsoid, line, tri, body),
then that property is also updated appropriately to correspond to the
atom&#8217;s rotation.</p>
<p>Distance units for displacements and the origin point of the <em>rotate</em>
style are determined by the setting of <em>box</em> or <em>lattice</em> for the
<em>units</em> keyword. <em>Box</em> means distance units as defined by the

7
doc/displace_atoms.txt
View File

@ -86,9 +86,14 @@ The {rotate} style rotates each atom in the group by the angle {theta}
around a rotation axis {R} = (Rx,Ry,Rz) that goes thru a point {P} =
(Px,Py,Pz). The direction of rotation for the atoms around the
rotation axis is consistent with the right-hand rule: if your
right-hand's thumb points along {R}, then your fingers wrap around the
right-hand thumb points along {R}, then your fingers wrap around the
axis in the direction of positive theta.
If the defined "atom_style"_atom_style.html assigns an orientation to
each atom ("atom styles"_atom_style.html ellipsoid, line, tri, body),
then that property is also updated appropriately to correspond to the
atom's rotation.
Distance units for displacements and the origin point of the {rotate}
style are determined by the setting of {box} or {lattice} for the
{units} keyword. {Box} means distance units as defined by the

97
doc/dump_image.html
View File

@ -143,18 +143,23 @@
<li>color = atom attribute that determines color of each atom</li>
<li>diameter = atom attribute that determines size of each atom</li>
<li>zero or more keyword/value pairs may be appended</li>
<li>keyword = <em>adiam</em> or <em>atom</em> or <em>body</em> or <em>bond</em> or <em>size</em> or <em>view</em> or <em>center</em> or <em>up</em> or <em>zoom</em> or <em>persp</em> or <em>box</em> or <em>axes</em> or <em>subbox</em> or <em>shiny</em> or <em>ssao</em></li>
<li>keyword = <em>atom</em> or <em>adiam</em> or <em>bond</em> or <em>line</em> or <em>tri</em> or <em>body</em> or <em>size</em> or <em>view</em> or <em>center</em> or <em>up</em> or <em>zoom</em> or <em>persp</em> or <em>box</em> or <em>axes</em> or <em>subbox</em> or <em>shiny</em> or <em>ssao</em></li>
</ul>
<pre class="literal-block">
<em>adiam</em> value = number = numeric value for atom diameter (distance units)
<em>atom</em> = yes/no = do or do not draw atoms
<em>body</em> = yes/no bflag1 bflag2
yes/no = do or do not draw atoms as bodies
bflag1,bflag2 = 2 numeric flags to affect how bodies are drawn
<em>adiam</em> size = numeric value for atom diameter (distance units)
<em>bond</em> values = color width = color and width of bonds
color = <em>atom</em> or <em>type</em> or <em>none</em>
width = number or <em>atom</em> or <em>type</em> or <em>none</em>
number = numeric value for bond width (distance units)
<em>line</em> = color width
color = <em>type</em>
width = numeric value for line width (distance units)
<em>tri</em> = color
color = <em>type</em>
<em>body</em> = color bflag1 bflag2
color = <em>type</em>
bflag1,bflag2 = 2 numeric flags to affect how bodies are drawn
<em>size</em> values = width height = size of images
width = width of image in # of pixels
height = height of image in # of pixels
@ -350,26 +355,15 @@ will likely not need to change. The <a class="reference internal" href="dump_mo
also has options specific to the dump image style, particularly for
assigning colors to atoms, bonds, and other image features.</p>
<hr class="docutils" />
<p>The <em>adiam</em> keyword allows you to override the <em>diameter</em> setting to a
per-atom attribute with a specified numeric value. All atoms will be
drawn with that diameter, e.g. 1.5, which is in whatever distance
<a class="reference internal" href="units.html"><em>units</em></a> the input script defines, e.g. Angstroms.</p>
<p>The <em>atom</em> keyword allow you to turn off the drawing of all atoms,
if the specified value is <em>no</em>.</p>
<p>The <em>body</em> keyword can be used when <a class="reference internal" href="atom_style.html"><em>atom_style body</em></a>
is used to define body particles with internal state
(e.g. sub-particles). The <a class="reference internal" href="body.html"><em>body</em></a> doc page descibes the body
styles LAMMPS currently supports, and provides more details as to the
kind of body particles they represent and how they are drawn by this
dump image command. For all the body styles, individual atoms can be
either a body particle or a usual point (non-body) particle. If the
<em>body</em> keyword is set to <em>yes</em>, then atoms which are body particles
are drawn by the method defined by the body style. Non-body particles
the same way they would be if the <em>body</em> keyword is <em>no</em>, i.e. as
spheres. The <em>bflag1</em> and <em>bflag2</em> settings are numerical values
which are passed to the body style to affect how the drawing of a body
particle is done. See the <a class="reference internal" href="body.html"><em>body</em></a> doc page for a description
of what these parameters mean for each body style.</p>
<p>The <em>atom</em> keyword allow you to turn off the drawing of all atoms, if
the specified value is <em>no</em>. Note that this will not turn off the
drawing of particles that are represented as lines, triangles, or
bodies, as discussed below. These particles can be drawn separately
if the <em>line</em>, <em>tri</em>, or <em>body</em> keywords are used.</p>
<p>The <em>adiam</em> keyword allows you to override the <em>diameter</em> setting to
set a single numeric <em>size</em>. All atoms will be drawn with that
diameter, e.g. 1.5, which is in whatever distance <a class="reference internal" href="units.html"><em>units</em></a>
the input script defines, e.g. Angstroms.</p>
<p>The <em>bond</em> keyword allows to you to alter how bonds are drawn. A bond
is only drawn if both atoms in the bond are being drawn due to being
in the specified group and due to other selection criteria
@ -405,6 +399,59 @@ of the 2 atoms in the bond.</p>
<p>If <em>type</em> is specified for the <em>width</em> value then the diameter of each
bond is determined by its bond type. By default all types have
diameter 0.5. This mapping can be changed by the <a class="reference internal" href="dump_modify.html"><em>dump_modify bdiam</em></a> command.</p>
<p>The <em>line</em> keyword can be used when <a class="reference internal" href="atom_style.html"><em>atom_style line</em></a>
is used to define particles as line segments, and will draw them as
lines. If this keyword is not used, such particles will be drawn as
spheres, the same as if they were regular atoms. The only setting
currently allowed for the <em>color</em> value is <em>type</em>, which will color
the lines according to the atom type of the particle. By default the
mapping of types to colors is as follows:</p>
<ul class="simple">
<li>type 1 = red</li>
<li>type 2 = green</li>
<li>type 3 = blue</li>
<li>type 4 = yellow</li>
<li>type 5 = aqua</li>
<li>type 6 = cyan</li>
</ul>
<p>and repeats itself for types &gt; 6. There is not yet an option to
change this via the <a class="reference internal" href="dump_modify.html"><em>dump_modify</em></a> command.</p>
<p>The line <em>width</em> can only be a numeric value, which specifies that all
lines will be drawn as cylinders with that diameter, e.g. 1.0, which
is in whatever distance <a class="reference internal" href="units.html"><em>units</em></a> the input script defines,
e.g. Angstroms.</p>
<p>The <em>tri</em> keyword can be used when <a class="reference internal" href="atom_style.html"><em>atom_style tri</em></a> is
used to define particles as triangles, and will draw them as
triangles. If this keyword is not used, such particles will be drawn
as spheres, the same as if they were regular atoms. The only setting
currently allowed for the <em>color</em> value is <em>type</em>, which will color
the triangles according to the atom type of the particle. By default
the mapping of types to colors is as follows:</p>
<ul class="simple">
<li>type 1 = red</li>
<li>type 2 = green</li>
<li>type 3 = blue</li>
<li>type 4 = yellow</li>
<li>type 5 = aqua</li>
<li>type 6 = cyan</li>
</ul>
<p>and repeats itself for types &gt; 6. There is not yet an option to
change this via the <a class="reference internal" href="dump_modify.html"><em>dump_modify</em></a> command.</p>
<p>The <em>body</em> keyword can be used when <a class="reference internal" href="atom_style.html"><em>atom_style body</em></a>
is used to define body particles with internal state
(e.g. sub-particles), and will drawn them in a manner specific to the
body style. If this keyword is not used, such particles will be drawn
as spheres, the same as if they were regular atoms.</p>
<p>The <a class="reference internal" href="body.html"><em>body</em></a> doc page descibes the body styles LAMMPS
currently supports, and provides more details as to the kind of body
particles they represent and how they are drawn by this dump image
command. For all the body styles, individual atoms can be either a
body particle or a usual point (non-body) particle. Non-body
particles will be drawn the same way they would be as a regular atom.
The <em>bflag1</em> and <em>bflag2</em> settings are numerical values which are
passed to the body style to affect how the drawing of a body particle
is done. See the <a class="reference internal" href="body.html"><em>body</em></a> doc page for a description of what
these parameters mean for each body style.</p>
<hr class="docutils" />
<p>The <em>size</em> keyword sets the width and height of the created images,
i.e. the number of pixels in each direction.</p>

103
doc/dump_image.txt
View File

@ -21,16 +21,21 @@ file = name of file to write image to :l
color = atom attribute that determines color of each atom :l
diameter = atom attribute that determines size of each atom :l
zero or more keyword/value pairs may be appended :l
keyword = {adiam} or {atom} or {body} or {bond} or {size} or {view} or {center} or {up} or {zoom} or {persp} or {box} or {axes} or {subbox} or {shiny} or {ssao} :l
{adiam} value = number = numeric value for atom diameter (distance units)
keyword = {atom} or {adiam} or {bond} or {line} or {tri} or {body} or {size} or {view} or {center} or {up} or {zoom} or {persp} or {box} or {axes} or {subbox} or {shiny} or {ssao} :l
{atom} = yes/no = do or do not draw atoms
{body} = yes/no bflag1 bflag2
yes/no = do or do not draw atoms as bodies
bflag1,bflag2 = 2 numeric flags to affect how bodies are drawn
{adiam} size = numeric value for atom diameter (distance units)
{bond} values = color width = color and width of bonds
color = {atom} or {type} or {none}
width = number or {atom} or {type} or {none}
number = numeric value for bond width (distance units)
{line} = color width
color = {type}
width = numeric value for line width (distance units)
{tri} = color
color = {type}
{body} = color bflag1 bflag2
color = {type}
bflag1,bflag2 = 2 numeric flags to affect how bodies are drawn
{size} values = width height = size of images
width = width of image in # of pixels
height = height of image in # of pixels
@ -234,28 +239,16 @@ assigning colors to atoms, bonds, and other image features.
:line
The {adiam} keyword allows you to override the {diameter} setting to a
per-atom attribute with a specified numeric value. All atoms will be
drawn with that diameter, e.g. 1.5, which is in whatever distance
"units"_units.html the input script defines, e.g. Angstroms.
The {atom} keyword allow you to turn off the drawing of all atoms, if
the specified value is {no}. Note that this will not turn off the
drawing of particles that are represented as lines, triangles, or
bodies, as discussed below. These particles can be drawn separately
if the {line}, {tri}, or {body} keywords are used.
The {atom} keyword allow you to turn off the drawing of all atoms,
if the specified value is {no}.
The {body} keyword can be used when "atom_style body"_atom_style.html
is used to define body particles with internal state
(e.g. sub-particles). The "body"_body.html doc page descibes the body
styles LAMMPS currently supports, and provides more details as to the
kind of body particles they represent and how they are drawn by this
dump image command. For all the body styles, individual atoms can be
either a body particle or a usual point (non-body) particle. If the
{body} keyword is set to {yes}, then atoms which are body particles
are drawn by the method defined by the body style. Non-body particles
the same way they would be if the {body} keyword is {no}, i.e. as
spheres. The {bflag1} and {bflag2} settings are numerical values
which are passed to the body style to affect how the drawing of a body
particle is done. See the "body"_body.html doc page for a description
of what these parameters mean for each body style.
The {adiam} keyword allows you to override the {diameter} setting to
set a single numeric {size}. All atoms will be drawn with that
diameter, e.g. 1.5, which is in whatever distance "units"_units.html
the input script defines, e.g. Angstroms.
The {bond} keyword allows to you to alter how bonds are drawn. A bond
is only drawn if both atoms in the bond are being drawn due to being
@ -300,6 +293,64 @@ bond is determined by its bond type. By default all types have
diameter 0.5. This mapping can be changed by the "dump_modify
bdiam"_dump_modify.html command.
The {line} keyword can be used when "atom_style line"_atom_style.html
is used to define particles as line segments, and will draw them as
lines. If this keyword is not used, such particles will be drawn as
spheres, the same as if they were regular atoms. The only setting
currently allowed for the {color} value is {type}, which will color
the lines according to the atom type of the particle. By default the
mapping of types to colors is as follows:
type 1 = red
type 2 = green
type 3 = blue
type 4 = yellow
type 5 = aqua
type 6 = cyan :ul
and repeats itself for types > 6. There is not yet an option to
change this via the "dump_modify"_dump_modify.html command.
The line {width} can only be a numeric value, which specifies that all
lines will be drawn as cylinders with that diameter, e.g. 1.0, which
is in whatever distance "units"_units.html the input script defines,
e.g. Angstroms.
The {tri} keyword can be used when "atom_style tri"_atom_style.html is
used to define particles as triangles, and will draw them as
triangles. If this keyword is not used, such particles will be drawn
as spheres, the same as if they were regular atoms. The only setting
currently allowed for the {color} value is {type}, which will color
the triangles according to the atom type of the particle. By default
the mapping of types to colors is as follows:
type 1 = red
type 2 = green
type 3 = blue
type 4 = yellow
type 5 = aqua
type 6 = cyan :ul
and repeats itself for types > 6. There is not yet an option to
change this via the "dump_modify"_dump_modify.html command.
The {body} keyword can be used when "atom_style body"_atom_style.html
is used to define body particles with internal state
(e.g. sub-particles), and will drawn them in a manner specific to the
body style. If this keyword is not used, such particles will be drawn
as spheres, the same as if they were regular atoms.
The "body"_body.html doc page descibes the body styles LAMMPS
currently supports, and provides more details as to the kind of body
particles they represent and how they are drawn by this dump image
command. For all the body styles, individual atoms can be either a
body particle or a usual point (non-body) particle. Non-body
particles will be drawn the same way they would be as a regular atom.
The {bflag1} and {bflag2} settings are numerical values which are
passed to the body style to affect how the drawing of a body particle
is done. See the "body"_body.html doc page for a description of what
these parameters mean for each body style.
:line
The {size} keyword sets the width and height of the created images,

21
doc/fix_move.html
View File

@ -249,16 +249,17 @@ fix 1 boundary move variable v_x NULL NULL v_v NULL NULL
</div>
<p>The <em>rotate</em> style rotates atoms around a rotation axis <em>R</em> =
(Rx,Ry,Rz) that goes thru a point <em>P</em> = (Px,Py,Pz). The <em>period</em> of
the rotation is also specified. This style also sets the velocity of
each atom to (omega cross Rperp) where omega is its angular velocity
around the rotation axis and Rperp is a perpendicular vector from the
rotation axis to the atom. If the defined
<a class="reference internal" href="atom_style.html"><em>atom_style</em></a> assigns an angular velocity to each atom,
then each atom&#8217;s angular velocity is also set to omega. Note that the
direction of rotation for the atoms around the rotation axis is
consistent with the right-hand rule: if your right-hand&#8217;s thumb points
along <em>R</em>, then your fingers wrap around the axis in the direction of
rotation.</p>
the rotation is also specified. The direction of rotation for the
atoms around the rotation axis is consistent with the right-hand rule:
if your right-hand thumb points along <em>R</em>, then your fingers wrap
around the axis in the direction of rotation.</p>
<p>This style also sets the velocity of each atom to (omega cross Rperp)
where omega is its angular velocity around the rotation axis and Rperp
is a perpendicular vector from the rotation axis to the atom. If the
defined <a class="reference internal" href="atom_style.html"><em>atom_style</em></a> assigns an angular velocity or
angular moementum or orientation to each atom (<a class="reference internal" href="atom_style.html"><em>atom styles</em></a> sphere, ellipsoid, line, tri, body), then
those properties are also updated appropriately to correspond to the
atom&#8217;s motion and rotation over time.</p>
<p>The <em>variable</em> style allows the position and velocity components of
each atom to be set by formulas specified via the
<a class="reference internal" href="variable.html"><em>variable</em></a> command. Each of the 6 variables is

23
doc/fix_move.txt
View File

@ -122,16 +122,19 @@ fix 1 boundary move variable v_x NULL NULL v_v NULL NULL :pre
The {rotate} style rotates atoms around a rotation axis {R} =
(Rx,Ry,Rz) that goes thru a point {P} = (Px,Py,Pz). The {period} of
the rotation is also specified. This style also sets the velocity of
each atom to (omega cross Rperp) where omega is its angular velocity
around the rotation axis and Rperp is a perpendicular vector from the
rotation axis to the atom. If the defined
"atom_style"_atom_style.html assigns an angular velocity to each atom,
then each atom's angular velocity is also set to omega. Note that the
direction of rotation for the atoms around the rotation axis is
consistent with the right-hand rule: if your right-hand's thumb points
along {R}, then your fingers wrap around the axis in the direction of
rotation.
the rotation is also specified. The direction of rotation for the
atoms around the rotation axis is consistent with the right-hand rule:
if your right-hand thumb points along {R}, then your fingers wrap
around the axis in the direction of rotation.
This style also sets the velocity of each atom to (omega cross Rperp)
where omega is its angular velocity around the rotation axis and Rperp
is a perpendicular vector from the rotation axis to the atom. If the
defined "atom_style"_atom_style.html assigns an angular velocity or
angular moementum or orientation to each atom ("atom
styles"_atom_style.html sphere, ellipsoid, line, tri, body), then
those properties are also updated appropriately to correspond to the
atom's motion and rotation over time.
The {variable} style allows the position and velocity components of
each atom to be set by formulas specified via the

18
doc/fix_store_state.html
View File

@ -200,15 +200,15 @@ time for future use in a calculation or output. See the discussion of
inputs.</p>
<p>If <em>N</em> is not zero, then the attributes will be updated every <em>N</em>
steps.</p>
<div class="admonition warning">
<p class="first admonition-title">Warning</p>
<p class="last">Actually, only atom attributes specified by keywords
like <em>xu</em> or <em>vy</em> or <em>radius</em> are initially stored immediately at the
point in your input script when the fix is defined. Attributes
specified by a compute, fix, or variable are not initially stored
until the first run following the fix definition begins. This is
because calculating those attributes may require quantities that are
not defined in between runs.</p>
<div class="admonition note">
<p class="first admonition-title">Note</p>
<p class="last">Actually, only atom attributes specified by keywords like <em>xu</em>
or <em>vy</em> or <em>radius</em> are initially stored immediately at the point in
your input script when the fix is defined. Attributes specified by a
compute, fix, or variable are not initially stored until the first run
following the fix definition begins. This is because calculating
those attributes may require quantities that are not defined in
between runs.</p>
</div>
<p>The list of possible attributes is the same as that used by the <a class="reference internal" href="dump.html"><em>dump custom</em></a> command, which describes their meaning.</p>
<p>If the <em>com</em> keyword is set to <em>yes</em> then the <em>xu</em>, <em>yu</em>, and <em>zu</em>

14
doc/fix_store_state.txt
View File

@ -74,13 +74,13 @@ inputs.
If {N} is not zero, then the attributes will be updated every {N}
steps.
IMPORTANT NOTE: Actually, only atom attributes specified by keywords
like {xu} or {vy} or {radius} are initially stored immediately at the
point in your input script when the fix is defined. Attributes
specified by a compute, fix, or variable are not initially stored
until the first run following the fix definition begins. This is
because calculating those attributes may require quantities that are
not defined in between runs.
NOTE: Actually, only atom attributes specified by keywords like {xu}
or {vy} or {radius} are initially stored immediately at the point in
your input script when the fix is defined. Attributes specified by a
compute, fix, or variable are not initially stored until the first run
following the fix definition begins. This is because calculating
those attributes may require quantities that are not defined in
between runs.
The list of possible attributes is the same as that used by the "dump
custom"_dump.html command, which describes their meaning.

2
doc/pair_beck.html
View File

@ -189,7 +189,7 @@ more instructions on how to use the accelerated styles effectively.</p>
<hr class="docutils" />
<p><strong>Mixing, shift, table, tail correction, restart, rRESPA info</strong>:</p>
<p>For atom type pairs I,J and I != J, coeffiecients must be specified.
No default miture rules are used.</p>
No default mixing rules are used.</p>
<p>This pair style does not support the <a class="reference internal" href="pair_modify.html"><em>pair_modify</em></a> shift
option for the energy of the pair interaction.</p>
<p>The <a class="reference internal" href="pair_modify.html"><em>pair_modify</em></a> table option is not relevant

2
doc/pair_beck.txt
View File

@ -74,7 +74,7 @@ more instructions on how to use the accelerated styles effectively.
[Mixing, shift, table, tail correction, restart, rRESPA info]:
For atom type pairs I,J and I != J, coeffiecients must be specified.
No default miture rules are used.
No default mixing rules are used.
This pair style does not support the "pair_modify"_pair_modify.html shift
option for the energy of the pair interaction.

102
doc/pair_line_lj.html
View File

@ -139,8 +139,9 @@
<div class="section" id="examples">
<h2>Examples<a class="headerlink" href="#examples" title="Permalink to this headline"></a></h2>
<div class="highlight-python"><div class="highlight"><pre>pair_style line/lj 3.0
pair_coeff * * 1.0 1.0
pair_coeff 1 1 1.0 1.5 2.5
pair_coeff * * 1.0 1.0 1.0 0.8 1.12
pair_coeff 1 2 1.0 2.0 1.0 1.5 1.12 5.0
pair_coeff 1 2 1.0 0.0 1.0 1.0 2.5
</pre></div>
</div>
</div>
@ -154,44 +155,77 @@ N1*N2 Lennard-Jones interactions. Interactions between a line segment
with N spherical particles and a point particle are treated as the
pairwise sum of N Lennard-Jones interactions. See the <a class="reference internal" href="pair_lj.html"><em>pair_style lj/cut</em></a> doc page for the definition of Lennard-Jones
interactions.</p>
<p>The cutoff distance for an interaction between 2 line segments, or
between a line segment and a point particle, is calculated from the
position of the line segment (its center), not between pairs of
individual spheres comprising the line segment. Thus an interaction
is either calculated in its entirety or not at all.</p>
<p>The set of non-overlapping spherical particles that represent a line
segment, for purposes of this pair style, are generated in the
following manner. Their size is a function of the line segment length
and the specified sigma for that particle type. If a line segment has
a length L and is of type I, then the number of spheres N that
represent the segment is calculated as N = L/sigma_II, rounded up to
an integer value. Thus if L is not evenly divisibly by sigam_II, N is
incremented to include one extra sphere. In this case, the spheres
must be slightly smaller than sigma_II so as not to overlap, so a new
sigma-prime is chosen as the sphere diameter, such that L/N =
sigma-prime. Thus the line segment interacts with other segments or
point particles as a collection of N spheres of diameter sigma-prime,
evenly spaced along the line segment, so as to exactly cover its
length.</p>
<p>The LJ interaction between 2 spheres on different line segments of
types I,J is computed with an arithmetic mixing of the sigma values of
the 2 spheres and using the specified epsilon value for I,J atom
types. Note that because the sigma values for line segment spheres is
computed using only sigma_II values, specific to the line segment&#8217;s
type, this means that any specified sigma_IJ values (for I != J) are
effectively ignored.</p>
<p>The set of non-overlapping spherical sub-particles that represent a
line segment are generated in the following manner. Their size is a
function of the line segment length and the specified sub-particle
size for that particle type. If a line segment has a length L and is
of type I, then the number of spheres N that represent the segment is
calculated as N = L/sizeI, rounded up to an integer value. Thus if L
is not evenly divisibly by sizeI, N is incremented to include one
extra sphere. The centers of the spheres are spaced equally along the
line segment. Imagine N+1 equally-space points, which include the 2
end points of the segment. The sphere centers are halfway between
each pair of points.</p>
<p>The LJ interaction between 2 spheres on different line segments (or a
sphere on a line segment and a point particles) is computed with
sub-particle epsilon, sigma, and cutoff values that are set by the
pair_coeff command, as described below. If the distance bewteen the 2
spheres is greater than the sub-particle cutoff, there is no
interaction. This means that some pairs of sub-particles on 2 line
segments may interact, but others may not.</p>
<p>For purposes of creating the neighbor list for pairs of interacting
line segments or lines/point particles, a regular particle-particle
cutoff is used, as defined by the <em>cutoff</em> setting above in the
pair_style command or overridden with an optional argument in the
pair_coeff command for a type pair as discussed below. The distance
between the centers of 2 line segments, or the center of a line
segment and a point particle, must be less than this distance (plus
the neighbor skin; see the <a class="reference external" href="neighbor">neighbor</a> command), for the pair
of particles to be included in the neighbor list.</p>
<div class="admonition note">
<p class="first admonition-title">Note</p>
<p class="last">This means that a too-short value for the <em>cutoff</em> setting can
exclude a pair of particles from the neighbor list even if pairs of
their sub-particle spheres would interact, based on the sub-particle
cutoff specified in the pair_coeff command. E.g. sub-particles at the
ends of the line segments that are close to each other. Which may not
be what you want, since it means the ends of 2 line segments could
pass through each other. It is up to you to specify a <em>cutoff</em>
setting that is consistent with the length of the line segments you
are using and the sub-particle cutoff settings.</p>
</div>
<p>For style <em>line/lj</em>, the following coefficients must be defined for
each pair of atoms types via the <a class="reference internal" href="pair_coeff.html"><em>pair_coeff</em></a> command
each pair of atom types via the <a class="reference internal" href="pair_coeff.html"><em>pair_coeff</em></a> command
as in the examples above, or in the data file or restart files read by
the <a class="reference internal" href="read_data.html"><em>read_data</em></a> or <a class="reference internal" href="read_restart.html"><em>read_restart</em></a>
commands:</p>
<ul class="simple">
<li>sizeI (distance units)</li>
<li>sizeJ (distance units)</li>
<li>epsilon (energy units)</li>
<li>sigma (distance units)</li>
<li>subcutoff (distance units)</li>
<li>cutoff (distance units)</li>
</ul>
<p>The last coefficient is optional. If not specified, the global cutoff
is used.</p>
<p>The <em>sizeI</em> and <em>sizeJ</em> coefficients are the sub-particle sizes for
line particles of type I and type J. They are used to define the N
sub-particles per segment as described above. These coefficients are
actually stored on a per-type basis. Thus if there are multiple
pair_coeff commmands that involve type I, as either the first or
second atom type, you should use consistent values for sizeI or sizeJ
in all of them. If you do not do this, the last value specified for
sizeI will apply to all segments of type I. If typeI or typeJ refers
to point particles, the corresponding sizeI or sizeJ is ignored; it
can be set to 0.0.</p>
<p>The <em>epsilon</em>, <em>sigma</em>, and <em>subcutoff</em> coefficients are used to
compute an LJ interactions between a pair of sub-particles on 2 line
segments (of type I and J), or between a sub particle/point particle
pair. As discussed above, the <em>subcutoff</em> and <em>cutoff</em> params are
different. The latter is only used for building the neighbor list
when the distance between centers of two line segments or one segment
and a point particle is calculated.</p>
<p>The <em>cutoff</em> coefficient is optional. If not specified, the global
cutoff is used.</p>
<hr class="docutils" />
<p>Styles with a <em>cuda</em>, <em>gpu</em>, <em>intel</em>, <em>kk</em>, <em>omp</em>, or <em>opt</em> suffix are
functionally the same as the corresponding style without the suffix.
@ -210,10 +244,8 @@ use the <a class="reference internal" href="suffix.html"><em>suffix</em></a> com
more instructions on how to use the accelerated styles effectively.</p>
<hr class="docutils" />
<p><strong>Mixing, shift, table, tail correction, restart, rRESPA info</strong>:</p>
<p>For atom type pairs I,J and I != J, the epsilon and sigma coefficients
and cutoff distance for all of this pair style can be mixed. The
default mix value is <em>geometric</em>. See the &#8220;pair_modify&#8221; command for
details.</p>
<p>For atom type pairs I,J and I != J, coeffiecients must be specified.
No default mixing rules are used.</p>
<p>This pair style does not support the <a class="reference internal" href="pair_modify.html"><em>pair_modify</em></a>
shift, table, and tail options.</p>
<p>This pair style does not write its information to <a class="reference internal" href="restart.html"><em>binary restart files</em></a>.</p>

102
doc/pair_line_lj.txt
View File

@ -18,8 +18,9 @@ cutoff = global cutoff for interactions (distance units)
[Examples:]
pair_style line/lj 3.0
pair_coeff * * 1.0 1.0
pair_coeff 1 1 1.0 1.5 2.5 :pre
pair_coeff * * 1.0 1.0 1.0 0.8 1.12
pair_coeff 1 2 1.0 2.0 1.0 1.5 1.12 5.0
pair_coeff 1 2 1.0 0.0 1.0 1.0 2.5 :pre
[Description:]
@ -33,47 +34,80 @@ pairwise sum of N Lennard-Jones interactions. See the "pair_style
lj/cut"_pair_lj.html doc page for the definition of Lennard-Jones
interactions.
The cutoff distance for an interaction between 2 line segments, or
between a line segment and a point particle, is calculated from the
position of the line segment (its center), not between pairs of
individual spheres comprising the line segment. Thus an interaction
is either calculated in its entirety or not at all.
The set of non-overlapping spherical sub-particles that represent a
line segment are generated in the following manner. Their size is a
function of the line segment length and the specified sub-particle
size for that particle type. If a line segment has a length L and is
of type I, then the number of spheres N that represent the segment is
calculated as N = L/sizeI, rounded up to an integer value. Thus if L
is not evenly divisibly by sizeI, N is incremented to include one
extra sphere. The centers of the spheres are spaced equally along the
line segment. Imagine N+1 equally-space points, which include the 2
end points of the segment. The sphere centers are halfway between
each pair of points.
The set of non-overlapping spherical particles that represent a line
segment, for purposes of this pair style, are generated in the
following manner. Their size is a function of the line segment length
and the specified sigma for that particle type. If a line segment has
a length L and is of type I, then the number of spheres N that
represent the segment is calculated as N = L/sigma_II, rounded up to
an integer value. Thus if L is not evenly divisibly by sigam_II, N is
incremented to include one extra sphere. In this case, the spheres
must be slightly smaller than sigma_II so as not to overlap, so a new
sigma-prime is chosen as the sphere diameter, such that L/N =
sigma-prime. Thus the line segment interacts with other segments or
point particles as a collection of N spheres of diameter sigma-prime,
evenly spaced along the line segment, so as to exactly cover its
length.
The LJ interaction between 2 spheres on different line segments (or a
sphere on a line segment and a point particles) is computed with
sub-particle epsilon, sigma, and cutoff values that are set by the
pair_coeff command, as described below. If the distance bewteen the 2
spheres is greater than the sub-particle cutoff, there is no
interaction. This means that some pairs of sub-particles on 2 line
segments may interact, but others may not.
The LJ interaction between 2 spheres on different line segments of
types I,J is computed with an arithmetic mixing of the sigma values of
the 2 spheres and using the specified epsilon value for I,J atom
types. Note that because the sigma values for line segment spheres is
computed using only sigma_II values, specific to the line segment's
type, this means that any specified sigma_IJ values (for I != J) are
effectively ignored.
For purposes of creating the neighbor list for pairs of interacting
line segments or lines/point particles, a regular particle-particle
cutoff is used, as defined by the {cutoff} setting above in the
pair_style command or overridden with an optional argument in the
pair_coeff command for a type pair as discussed below. The distance
between the centers of 2 line segments, or the center of a line
segment and a point particle, must be less than this distance (plus
the neighbor skin; see the "neighbor"_neighbor command), for the pair
of particles to be included in the neighbor list.
NOTE: This means that a too-short value for the {cutoff} setting can
exclude a pair of particles from the neighbor list even if pairs of
their sub-particle spheres would interact, based on the sub-particle
cutoff specified in the pair_coeff command. E.g. sub-particles at the
ends of the line segments that are close to each other. Which may not
be what you want, since it means the ends of 2 line segments could
pass through each other. It is up to you to specify a {cutoff}
setting that is consistent with the length of the line segments you
are using and the sub-particle cutoff settings.
For style {line/lj}, the following coefficients must be defined for
each pair of atoms types via the "pair_coeff"_pair_coeff.html command
each pair of atom types via the "pair_coeff"_pair_coeff.html command
as in the examples above, or in the data file or restart files read by
the "read_data"_read_data.html or "read_restart"_read_restart.html
commands:
sizeI (distance units)
sizeJ (distance units)
epsilon (energy units)
sigma (distance units)
subcutoff (distance units)
cutoff (distance units) :ul
The last coefficient is optional. If not specified, the global cutoff
is used.
The {sizeI} and {sizeJ} coefficients are the sub-particle sizes for
line particles of type I and type J. They are used to define the N
sub-particles per segment as described above. These coefficients are
actually stored on a per-type basis. Thus if there are multiple
pair_coeff commmands that involve type I, as either the first or
second atom type, you should use consistent values for sizeI or sizeJ
in all of them. If you do not do this, the last value specified for
sizeI will apply to all segments of type I. If typeI or typeJ refers
to point particles, the corresponding sizeI or sizeJ is ignored; it
can be set to 0.0.
The {epsilon}, {sigma}, and {subcutoff} coefficients are used to
compute an LJ interactions between a pair of sub-particles on 2 line
segments (of type I and J), or between a sub particle/point particle
pair. As discussed above, the {subcutoff} and {cutoff} params are
different. The latter is only used for building the neighbor list
when the distance between centers of two line segments or one segment
and a point particle is calculated.
The {cutoff} coefficient is optional. If not specified, the global
cutoff is used.
:line
@ -102,10 +136,8 @@ more instructions on how to use the accelerated styles effectively.
[Mixing, shift, table, tail correction, restart, rRESPA info]:
For atom type pairs I,J and I != J, the epsilon and sigma coefficients
and cutoff distance for all of this pair style can be mixed. The
default mix value is {geometric}. See the "pair_modify" command for
details.
For atom type pairs I,J and I != J, coeffiecients must be specified.
No default mixing rules are used.
This pair style does not support the "pair_modify"_pair_modify.html
shift, table, and tail options.

2
doc/searchindex.js
View File

File diff suppressed because one or more lines are too long

93
src/ASPHERE/pair_line_lj.cpp
View File

@ -52,7 +52,10 @@ PairLineLJ::~PairLineLJ()
memory->destroy(setflag);
memory->destroy(cutsq);
memory->destroy(subsize);
memory->destroy(cut);
memory->destroy(cutsub);
memory->destroy(cutsubsq);
memory->destroy(epsilon);
memory->destroy(sigma);
memory->destroy(lj1);
@ -66,7 +69,7 @@ PairLineLJ::~PairLineLJ()
void PairLineLJ::compute(int eflag, int vflag)
{
int i,j,ii,jj,inum,jnum,itype,jtype;
int i,j,ii,jj,inum,jnum,itype,jtype,tmp;
int ni,nj,npi,npj,ifirst,jfirst;
double xtmp,ytmp,ztmp,delx,dely,delz,evdwl,fpair;
double rsq,r2inv,r6inv,term1,term2,sig,sig3,forcelj;
@ -130,10 +133,10 @@ void PairLineLJ::compute(int eflag, int vflag)
evdwl = 0.0;
if (line[i] >= 0 && line[j] >= 0) {
if (dnum[i] == 0) discretize(i,sigma[itype][itype]);
if (dnum[i] == 0) discretize(i,subsize[itype]);
npi = dnum[i];
ifirst = dfirst[i];
if (dnum[j] == 0) discretize(j,sigma[jtype][jtype]);
if (dnum[j] == 0) discretize(j,subsize[jtype]);
npj = dnum[j];
jfirst = dfirst[j];
@ -154,7 +157,11 @@ void PairLineLJ::compute(int eflag, int vflag)
dely = xi[1] - xj[1];
rsq = delx*delx + dely*dely;
sig = 0.5 * (discrete[ifirst+ni].sigma+discrete[jfirst+nj].sigma);
// skip this pair of sub-particles if outside sub cutoff
if (rsq >= cutsubsq[itype][jtype]) continue;
sig = sigma[itype][jtype];
sig3 = sig*sig*sig;
term2 = 24.0*epsilon[itype][jtype] * sig3*sig3;
term1 = 2.0 * term2 * sig3*sig3;
@ -167,16 +174,15 @@ void PairLineLJ::compute(int eflag, int vflag)
fi[0] = delx*fpair;
fi[1] = dely*fpair;
f[i][0] += fi[0];
f[i][1] += fi[1];
torque[i][2] += dxi*fi[1] - dyi*fi[0];
if (newton_pair || j < nlocal) {
fj[0] = -delx*fpair;
fj[1] = -dely*fpair;
f[j][0] += fj[0];
f[j][1] += fj[1];
torque[j][2] += dxj*fj[1] - dyj*fj[0];
f[j][0] -= fi[0];
f[j][1] -= fi[1];
torque[j][2] -= dxj*fi[1] - dyj*fi[0];
}
}
}
@ -185,7 +191,7 @@ void PairLineLJ::compute(int eflag, int vflag)
// convert line into Np particles based on sigma and line length
} else if (line[i] >= 0) {
if (dnum[i] == 0) discretize(i,sigma[itype][itype]);
if (dnum[i] == 0) discretize(i,subsize[itype]);
npi = dnum[i];
ifirst = dfirst[i];
@ -202,7 +208,11 @@ void PairLineLJ::compute(int eflag, int vflag)
dely = xi[1] - xj[1];
rsq = delx*delx + dely*dely;
sig = 0.5 * (discrete[ifirst+ni].sigma+sigma[jtype][jtype]);
// skip this pair of sub-particles if outside sub cutoff
if (rsq >= cutsubsq[itype][jtype]) continue;
sig = sigma[itype][jtype];
sig3 = sig*sig*sig;
term2 = 24.0*epsilon[itype][jtype] * sig3*sig3;
term1 = 2.0 * term2 * sig3*sig3;
@ -220,10 +230,8 @@ void PairLineLJ::compute(int eflag, int vflag)
torque[i][2] += dxi*fi[1] - dyi*fi[0];
if (newton_pair || j < nlocal) {
fj[0] = -delx*fpair;
fj[1] = -dely*fpair;
f[j][0] += fj[0];
f[j][1] += fj[1];
f[j][0] -= fi[0];
f[j][1] -= fi[1];
}
}
@ -231,7 +239,7 @@ void PairLineLJ::compute(int eflag, int vflag)
// convert line into Np particles based on sigma and line length
} else if (line[j] >= 0) {
if (dnum[j] == 0) discretize(j,sigma[jtype][jtype]);
if (dnum[j] == 0) discretize(j,subsize[jtype]);
npj = dnum[j];
jfirst = dfirst[j];
@ -248,7 +256,11 @@ void PairLineLJ::compute(int eflag, int vflag)
dely = xi[1] - xj[1];
rsq = delx*delx + dely*dely;
sig = 0.5 * (sigma[itype][itype]+discrete[jfirst+nj].sigma);
// skip this pair of sub-particles if outside sub cutoff
if (rsq >= cutsubsq[itype][jtype]) continue;
sig = sigma[itype][jtype];
sig3 = sig*sig*sig;
term2 = 24.0*epsilon[itype][jtype] * sig3*sig3;
term1 = 2.0 * term2 * sig3*sig3;
@ -265,11 +277,9 @@ void PairLineLJ::compute(int eflag, int vflag)
f[i][1] += fi[1];
if (newton_pair || j < nlocal) {
f[j][0] += fj[0];
f[j][1] += fj[1];
fj[0] = -delx*fpair;
fj[1] = -dely*fpair;
torque[j][2] += dxj*fj[1] - dyj*fj[0];
f[j][0] -= fi[0];
f[j][1] -= fi[1];
torque[j][2] -= dxj*fi[1] - dyj*fi[0];
}
}
@ -318,7 +328,10 @@ void PairLineLJ::allocate()
memory->create(cutsq,n+1,n+1,"pair:cutsq");
memory->create(subsize,n+1,"pair:subsize");
memory->create(cut,n+1,n+1,"pair:cut");
memory->create(cutsub,n+1,n+1,"pair:cutsub");
memory->create(cutsubsq,n+1,n+1,"pair:cutsubsq");
memory->create(epsilon,n+1,n+1,"pair:epsilon");
memory->create(sigma,n+1,n+1,"pair:sigma");
memory->create(lj1,n+1,n+1,"pair:lj1");
@ -353,7 +366,7 @@ void PairLineLJ::settings(int narg, char **arg)
void PairLineLJ::coeff(int narg, char **arg)
{
if (narg < 4 || narg > 5)
if (narg < 7 || narg > 8)
error->all(FLERR,"Incorrect args for pair coefficients");
if (!allocated) allocate();
@ -361,17 +374,23 @@ void PairLineLJ::coeff(int narg, char **arg)
force->bounds(arg[0],atom->ntypes,ilo,ihi);
force->bounds(arg[1],atom->ntypes,jlo,jhi);
double epsilon_one = force->numeric(FLERR,arg[2]);
double sigma_one = force->numeric(FLERR,arg[3]);
double size_itype = force->numeric(FLERR,arg[2]);
double size_jtype = force->numeric(FLERR,arg[3]);
double epsilon_one = force->numeric(FLERR,arg[4]);
double sigma_one = force->numeric(FLERR,arg[5]);
double cutsub_one = force->numeric(FLERR,arg[6]);
double cut_one = cut_global;
if (narg == 5) cut_one = force->numeric(FLERR,arg[4]);
if (narg == 8) cut_one = force->numeric(FLERR,arg[7]);
int count = 0;
for (int i = ilo; i <= ihi; i++) {
for (int j = MAX(jlo,i); j <= jhi; j++) {
subsize[i] = size_itype;
subsize[j] = size_jtype;
epsilon[i][j] = epsilon_one;
sigma[i][j] = sigma_one;
cutsub[i][j] = cutsub_one;
cut[i][j] = cut_one;
setflag[i][j] = 1;
count++;
@ -399,12 +418,9 @@ void PairLineLJ::init_style()
double PairLineLJ::init_one(int i, int j)
{
if (setflag[i][j] == 0) {
epsilon[i][j] = mix_energy(epsilon[i][i],epsilon[j][j],
sigma[i][i],sigma[j][j]);
sigma[i][j] = mix_distance(sigma[i][i],sigma[j][j]);
cut[i][j] = mix_distance(cut[i][i],cut[j][j]);
}
if (setflag[i][j] == 0) error->all(FLERR,"All pair coeffs are not set");
cutsubsq[i][j] = cutsub[i][j] * cutsub[i][j];
lj1[i][j] = 48.0 * epsilon[i][j] * pow(sigma[i][j],12.0);
lj2[i][j] = 24.0 * epsilon[i][j] * pow(sigma[i][j],6.0);
@ -413,6 +429,8 @@ double PairLineLJ::init_one(int i, int j)
epsilon[j][i] = epsilon[i][j];
sigma[j][i] = sigma[i][j];
cutsubsq[j][i] = cutsubsq[i][j];
lj1[j][i] = lj1[i][j];
lj2[j][i] = lj2[i][j];
lj3[j][i] = lj3[i][j];
@ -422,16 +440,17 @@ double PairLineLJ::init_one(int i, int j)
}
/* ----------------------------------------------------------------------
discretize line segment I into N sub-segments no more than sigma in length
store new discrete particles in Discrete list
discretize line segment I into N sub-particles with <= size separation
store displacement dx,dy of discrete particles in Discrete list
------------------------------------------------------------------------- */
void PairLineLJ::discretize(int i, double sigma)
void PairLineLJ::discretize(int i, double size)
{
AtomVecLine::Bonus *bonus = avec->bonus;
double length = bonus[atom->line[i]].length;
double theta = bonus[atom->line[i]].theta;
int n = static_cast<int> (length/sigma) + 1;
int n = static_cast<int> (length/size) + 1;
dnum[i] = n;
dfirst[i] = ndiscrete;
@ -441,14 +460,12 @@ void PairLineLJ::discretize(int i, double sigma)
memory->srealloc(discrete,dmax*sizeof(Discrete),"pair:discrete");
}
sigma = length/n;
double delta;
for (int m = 0; m < n; m++) {
delta = -0.5 + (2*m+1)/(2.0*n);
discrete[ndiscrete].dx = delta*length*cos(theta);
discrete[ndiscrete].dy = delta*length*sin(theta);
discrete[ndiscrete].sigma = sigma;
ndiscrete++;
}
}

8
src/ASPHERE/pair_line_lj.h
View File

@ -36,14 +36,16 @@ class PairLineLJ : public Pair {
protected:
double cut_global;
double *subsize;
double **epsilon,**sigma,**cutsub,**cutsubsq;
double **cut;
double **epsilon,**sigma;
double **lj1,**lj2,**lj3,**lj4;
double **lj1,**lj2,**lj3,**lj4; // for sphere/sphere interactions
class AtomVecLine *avec;
double *size; // per-type size of sub-particles to tile line segment
struct Discrete {
double dx,dy;
double sigma;
};
Discrete *discrete; // list of all discretes for all lines
int ndiscrete; // number of discretes in list

2
src/ASPHERE/pair_tri_lj.cpp
View File

@ -513,7 +513,7 @@ double PairTriLJ::init_one(int i, int j)
------------------------------------------------------------------------- */
void PairTriLJ::discretize(int i, double sigma,
double *c1, double *c2, double *c3)
double *c1, double *c2, double *c3)
{
double centroid[3],dc1[3],dc2[3],dc3[3];

70
src/displace_atoms.cpp
View File

@ -24,11 +24,16 @@
#include "group.h"
#include "math_const.h"
#include "random_park.h"
#include "error.h"
#include "force.h"
#include "input.h"
#include "variable.h"
#include "atom_vec_ellipsoid.h"
#include "atom_vec_line.h"
#include "atom_vec_tri.h"
#include "atom_vec_body.h"
#include "math_extra.h"
#include "memory.h"
#include "error.h"
using namespace LAMMPS_NS;
using namespace MathConst;
@ -203,8 +208,10 @@ void DisplaceAtoms::command(int narg, char **arg)
// X = P + C + A cos(theta) + B sin(theta)
if (style == ROTATE) {
double axis[3],point[3];
double theta_new;
double axis[3],point[3],qrotate[4],qnew[4];
double a[3],b[3],c[3],d[3],disp[3],runit[3];
double *quat;
int dim = domain->dimension;
point[0] = xscale*force->numeric(FLERR,arg[2]);
@ -224,11 +231,36 @@ void DisplaceAtoms::command(int narg, char **arg)
runit[1] = axis[1]/len;
runit[2] = axis[2]/len;
double sine = sin(MY_PI*theta/180.0);
double cosine = cos(MY_PI*theta/180.0);
double angle = MY_PI*theta/180.0;
double sine = sin(angle);
double cosine = cos(angle);
double ddotr;
// flags for additional orientation info stored by some atom styles
int ellipsoid_flag = atom->ellipsoid_flag;
int line_flag = atom->line_flag;
int tri_flag = atom->tri_flag;
int body_flag = atom->body_flag;
int theta_flag = 0;
int quat_flag = 0;
if (line_flag) theta_flag = 1;
if (ellipsoid_flag || tri_flag || body_flag) quat_flag = 1;
// AtomVec pointers to retrieve per-atom storage of extra quantities
AtomVecEllipsoid *avec_ellipsoid =
(AtomVecEllipsoid *) atom->style_match("ellipsoid");
AtomVecLine *avec_line = (AtomVecLine *) atom->style_match("line");
AtomVecTri *avec_tri = (AtomVecTri *) atom->style_match("tri");
AtomVecBody *avec_body = (AtomVecBody *) atom->style_match("body");
double **x = atom->x;
int *ellipsoid = atom->ellipsoid;
int *line = atom->line;
int *tri = atom->tri;
int *body = atom->body;
int *mask = atom->mask;
int nlocal = atom->nlocal;
@ -253,6 +285,36 @@ void DisplaceAtoms::command(int narg, char **arg)
x[i][0] = point[0] + c[0] + disp[0];
x[i][1] = point[1] + c[1] + disp[1];
if (dim == 3) x[i][2] = point[2] + c[2] + disp[2];
// theta for lines
if (theta_flag && line[i] >= 0.0) {
theta_new = fmod(avec_line->bonus[line[i]].theta+angle,MY_2PI);
avec_line->bonus[atom->line[i]].theta = theta_new;
}
// quats for ellipsoids, tris, and bodies
if (quat_flag) {
quat = NULL;
if (ellipsoid_flag && ellipsoid[i] >= 0)
quat = avec_ellipsoid->bonus[ellipsoid[i]].quat;
else if (tri_flag && tri[i] >= 0)
quat = avec_tri->bonus[tri[i]].quat;
else if (body_flag && body[i] >= 0)
quat = avec_body->bonus[body[i]].quat;
if (quat) {
qrotate[0] = cosine;
qrotate[1] = runit[0]*sine;
qrotate[2] = runit[1]*sine;
qrotate[3] = runit[2]*sine;
MathExtra::quatquat(qrotate,quat,qnew);
quat[0] = qnew[0];
quat[1] = qnew[1];
quat[2] = qnew[2];
quat[3] = qnew[3];
}
}
}
}
}

200
src/dump_image.cpp
View File

@ -18,7 +18,8 @@
#include "dump_image.h"
#include "image.h"
#include "atom.h"
#include "atom_vec.h"
#include "atom_vec_line.h"
#include "atom_vec_tri.h"
#include "atom_vec_body.h"
#include "body.h"
#include "molecule.h"
@ -29,6 +30,7 @@
#include "input.h"
#include "variable.h"
#include "math_const.h"
#include "math_extra.h"
#include "error.h"
#include "memory.h"
@ -106,7 +108,7 @@ DumpImage::DumpImage(LAMMPS *lmp, int narg, char **arg) :
// set defaults for optional args
atomflag = YES;
bodyflag = NO;
lineflag = triflag = bodyflag = NO;
if (atom->nbondtypes == 0) bondflag = NO;
else {
bondflag = YES;
@ -132,28 +134,19 @@ DumpImage::DumpImage(LAMMPS *lmp, int narg, char **arg) :
int iarg = ioptional;
while (iarg < narg) {
if (strcmp(arg[iarg],"adiam") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal dump image command");
adiam = NUMERIC;
adiamvalue = force->numeric(FLERR,arg[iarg+1]);
if (adiamvalue <= 0.0) error->all(FLERR,"Illegal dump image command");
iarg += 2;
} else if (strcmp(arg[iarg],"atom") == 0) {
if (strcmp(arg[iarg],"atom") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal dump image command");
if (strcmp(arg[iarg+1],"yes") == 0) atomflag = YES;
else if (strcmp(arg[iarg+1],"no") == 0) atomflag = NO;
else error->all(FLERR,"Illegal dump image command");
iarg += 2;
} else if (strcmp(arg[iarg],"body") == 0) {
if (iarg+4 > narg) error->all(FLERR,"Illegal dump image command");
if (strcmp(arg[iarg+1],"yes") == 0) bodyflag = YES;
else if (strcmp(arg[iarg+1],"no") == 0) bodyflag = NO;
else error->all(FLERR,"Illegal dump image command");
bodyflag1 = force->numeric(FLERR,arg[iarg+2]);
bodyflag2 = force->numeric(FLERR,arg[iarg+3]);
iarg += 4;
} else if (strcmp(arg[iarg],"adiam") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal dump image command");
adiam = NUMERIC;
adiamvalue = force->numeric(FLERR,arg[iarg+1]);
if (adiamvalue <= 0.0) error->all(FLERR,"Illegal dump image command");
iarg += 2;
} else if (strcmp(arg[iarg],"bond") == 0) {
if (iarg+3 > narg) error->all(FLERR,"Illegal dump image command");
@ -174,6 +167,32 @@ DumpImage::DumpImage(LAMMPS *lmp, int narg, char **arg) :
else error->all(FLERR,"Illegal dump image command");
iarg += 3;
} else if (strcmp(arg[iarg],"line") == 0) {
if (iarg+3 > narg) error->all(FLERR,"Illegal dump image command");
lineflag = YES;
if (strcmp(arg[iarg+1],"type") == 0) lcolor = TYPE;
else error->all(FLERR,"Illegal dump image command");
ldiam = NUMERIC;
ldiamvalue = force->numeric(FLERR,arg[iarg+2]);
iarg += 3;
} else if (strcmp(arg[iarg],"tri") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal dump image command");
triflag = YES;
if (strcmp(arg[iarg+1],"type") == 0) tcolor = TYPE;
else error->all(FLERR,"Illegal dump image command");
iarg += 2;
} else if (strcmp(arg[iarg],"body") == 0) {
if (iarg+4 > narg) error->all(FLERR,"Illegal dump image command");
bodyflag = YES;
if (strcmp(arg[iarg+1],"type") == 0) bodycolor = TYPE;
else error->all(FLERR,"Illegal dump image command");
bodyflag1 = force->numeric(FLERR,arg[iarg+2]);
bodyflag2 = force->numeric(FLERR,arg[iarg+3]);
iarg += 4;
} else if (strcmp(arg[iarg],"size") == 0) {
if (iarg+3 > narg) error->all(FLERR,"Illegal dump image command");
int width = force->inumeric(FLERR,arg[iarg+1]);
@ -333,13 +352,26 @@ DumpImage::DumpImage(LAMMPS *lmp, int narg, char **arg) :
} else error->all(FLERR,"Illegal dump image command");
}
// error check for bodyflag
// error checks and setup for lineflag, triflag, bodyflag
if (bodyflag) {
AtomVecBody *avec = (AtomVecBody *) atom->style_match("body");
if (!avec) error->all(FLERR,"Dump image body yes requires atom style body");
bptr = avec->bptr;
if (lineflag) {
avec_line = (AtomVecLine *) atom->style_match("line");
if (!avec_line)
error->all(FLERR,"Dump image line requires atom style line");
}
if (triflag) {
avec_tri = (AtomVecTri *) atom->style_match("tri");
if (!avec_tri)
error->all(FLERR,"Dump image tri requires atom style tri");
}
if (bodyflag) {
avec_body = (AtomVecBody *) atom->style_match("body");
if (!avec_body)
error->all(FLERR,"Dump image body yes requires atom style body");
}
extraflag = 0;
if (lineflag || triflag || bodyflag) extraflag = 1;
// allocate image buffer now that image size is known
@ -679,18 +711,21 @@ void DumpImage::view_params()
void DumpImage::create_image()
{
int i,j,k,m,n,itype,atom1,atom2,imol,iatom,btype,ibonus;
int i,j,k,m,n,itype,atom1,atom2,imol,iatom,btype,ibonus,drawflag;
tagint tagprev;
double diameter,delx,dely,delz;
int *bodyvec;
double **bodyarray;
double *color,*color1,*color2;
double xmid[3];
double xmid[3],pt1[3],pt2[3],pt3[3];
double mat[3][3];
// render my atoms
if (atomflag) {
double **x = atom->x;
int *line = atom->line;
int *tri = atom->tri;
int *body = atom->body;
m = 0;
@ -719,18 +754,109 @@ void DumpImage::create_image()
diameter = buf[m+1];
}
if (!body || !bodyflag || body[j] < 0)
image->draw_sphere(x[j],color,diameter);
else {
ibonus = body[i];
n = bptr->image(ibonus,bodyflag1,bodyflag2,bodyvec,bodyarray);
for (k = 0; k < n; k++) {
if (bodyvec[k] == SPHERE)
image->draw_sphere(bodyarray[k],color,bodyarray[k][3]);
else if (bodyvec[k] == LINE)
image->draw_cylinder(&bodyarray[k][0],&bodyarray[k][3],
color,bodyarray[k][6],3);
}
// do not draw if line,tri,body keywords enabled and atom is one of those
drawflag = 1;
if (extraflag) {
if (lineflag && line[j] >= 0) drawflag = 0;
if (triflag && tri[j] >= 0) drawflag = 0;
if (bodyflag && body[j] >= 0) drawflag = 0;
}
if (drawflag) image->draw_sphere(x[j],color,diameter);
m += size_one;
}
}
// render atoms that are lines
if (lineflag) {
double length,theta,dx,dy;
double **x = atom->x;
int *line = atom->line;
int *type = atom->type;
for (i = 0; i < nchoose; i++) {
j = clist[i];
if (line[j] < 0) continue;
if (lcolor == TYPE) {
color = colortype[type[j]];
}
if (ldiam == NUMERIC) {
diameter = ldiamvalue;
}
length = avec_line->bonus[line[j]].length;
theta = avec_line->bonus[line[j]].theta;
dx = 0.5*length*cos(theta);
dy = 0.5*length*sin(theta);
pt1[0] = x[j][0] + dx;
pt1[1] = x[j][1] + dy;
pt1[2] = 0.0;
pt2[0] = x[j][0] - dx;
pt2[1] = x[j][1] - dy;
pt2[2] = 0.0;
image->draw_cylinder(pt1,pt2,color,ldiamvalue,3);
}
}
// render atoms that are triangles
if (triflag) {
double **x = atom->x;
int *tri = atom->tri;
int *type = atom->type;
for (i = 0; i < nchoose; i++) {
j = clist[i];
if (tri[j] < 0) continue;
if (tcolor == TYPE) {
color = colortype[type[j]];
}
MathExtra::quat_to_mat(avec_tri->bonus[tri[i]].quat,mat);
MathExtra::matvec(mat,avec_tri->bonus[tri[i]].c1,pt1);
MathExtra::matvec(mat,avec_tri->bonus[tri[i]].c2,pt2);
MathExtra::matvec(mat,avec_tri->bonus[tri[i]].c3,pt3);
MathExtra::add3(pt1,x[i],pt1);
MathExtra::add3(pt2,x[i],pt2);
MathExtra::add3(pt3,x[i],pt3);
image->draw_triangle(pt1,pt2,pt3,color);
}
}
// render atoms that are bodies
if (bodyflag) {
Body *bptr = avec_body->bptr;
double **x = atom->x;
int *body = atom->body;
m = 0;
for (i = 0; i < nchoose; i++) {
j = clist[i];
if (body[j] < 0) continue;
if (bodycolor == TYPE) {
itype = static_cast<int> (buf[m]);
color = colortype[itype];
}
ibonus = body[i];
n = bptr->image(ibonus,bodyflag1,bodyflag2,bodyvec,bodyarray);
for (k = 0; k < n; k++) {
if (bodyvec[k] == SPHERE)
image->draw_sphere(bodyarray[k],color,bodyarray[k][3]);
else if (bodyvec[k] == LINE)
image->draw_cylinder(&bodyarray[k][0],&bodyarray[k][3],
color,bodyarray[k][6],3);
}
m += size_one;

20
src/dump_image.h
View File

@ -37,13 +37,24 @@ class DumpImage : public DumpCustom {
int filetype;
enum{PPM,JPG,PNG};
int atomflag; // 0/1 for draw atoms
int acolor,adiam; // what determines color/diam of atoms
double adiamvalue; // atom diameter value
int atomflag,bondflag; // 0/1 for draw atoms,bonds
int lineflag; // 0/1 for draw atoms as lines
int lcolor,ldiam; // what determines color/diam of lines
double ldiamvalue; // line diameter value
int triflag; // 0/1 for draw atoms as triangles
int tcolor; // what determines color of tris
int bodyflag; // 0/1 for draw atoms as bodies
double bodyflag1,bodyflag2; // user params for drawing bodies
int bodycolor; // what determines color of bodies
double bodyflag1,bodyflag2; // user-specified params for drawing bodies
int bondflag; // 0/1 for draw bonds
int bcolor,bdiam; // what determines color/diam of bonds
double bdiamvalue; // bond diameter value
int extraflag; // 0/1 for any of line/tri/body flag set
char *thetastr,*phistr; // variables for view theta,phi
int thetavar,phivar; // index to theta,phi vars
int cflag; // static/dynamic box center
@ -64,8 +75,11 @@ class DumpImage : public DumpCustom {
double *diamtype,*diamelement,*bdiamtype; // per-type diameters
double **colortype,**colorelement,**bcolortype; // per-type colors
class AtomVecLine *avec_line; // ptrs to atom style (sub)classes
class AtomVecTri *avec_tri;
class AtomVecBody *avec_body;
class Image *image; // class that renders each image
class Body *bptr; // class for Body particles
int *chooseghost; // extended choose array for comm
double **bufcopy; // buffer for communicating bond/atom info

284
src/fix_move.cpp
View File

@ -26,7 +26,12 @@
#include "respa.h"
#include "input.h"
#include "variable.h"
#include "atom_vec_ellipsoid.h"
#include "atom_vec_line.h"
#include "atom_vec_tri.h"
#include "atom_vec_body.h"
#include "math_const.h"
#include "math_extra.h"
#include "memory.h"
#include "error.h"
@ -37,6 +42,8 @@ using namespace MathConst;
enum{LINEAR,WIGGLE,ROTATE,VARIABLE};
enum{EQUAL,ATOM};
#define INERTIA 0.2 // moment of inertia prefactor for ellipsoid
/* ---------------------------------------------------------------------- */
FixMove::FixMove(LAMMPS *lmp, int narg, char **arg) :
@ -216,7 +223,7 @@ FixMove::FixMove(LAMMPS *lmp, int narg, char **arg) :
// set omega_rotate from period
if (mstyle == WIGGLE || mstyle == ROTATE) omega_rotate = 2.0*MY_PI / period;
if (mstyle == WIGGLE || mstyle == ROTATE) omega_rotate = MY_2PI / period;
// runit = unit vector along rotation axis
@ -229,24 +236,54 @@ FixMove::FixMove(LAMMPS *lmp, int narg, char **arg) :
runit[2] = axis[2]/len;
}
// set omega_flag if particles store omega
// set flags for extra attributes particles may store
// relevant extra attributes = omega, angmom, theta, quat
omega_flag = atom->omega_flag;
angmom_flag = atom->angmom_flag;
radius_flag = atom->radius_flag;
ellipsoid_flag = atom->ellipsoid_flag;
line_flag = atom->line_flag;
tri_flag = atom->tri_flag;
body_flag = atom->body_flag;
theta_flag = quat_flag = 0;
if (line_flag) theta_flag = 1;
if (ellipsoid_flag || tri_flag || body_flag) quat_flag = 1;
extra_flag = 0;
if (omega_flag || angmom_flag || theta_flag || quat_flag) extra_flag = 1;
// perform initial allocation of atom-based array
// register with Atom class
xoriginal = NULL;
toriginal = NULL;
qoriginal = NULL;
grow_arrays(atom->nmax);
atom->add_callback(0);
atom->add_callback(1);
displace = velocity = NULL;
// AtomVec pointers to retrieve per-atom storage of extra quantities
avec_ellipsoid = (AtomVecEllipsoid *) atom->style_match("ellipsoid");
avec_line = (AtomVecLine *) atom->style_match("line");
avec_tri = (AtomVecTri *) atom->style_match("tri");
avec_body = (AtomVecBody *) atom->style_match("body");
// xoriginal = initial unwrapped positions of atoms
// toriginal = initial theta of lines
// qoriginal = initial quat of extended particles
double **x = atom->x;
imageint *image = atom->image;
int *ellipsoid = atom->ellipsoid;
int *line = atom->line;
int *tri = atom->tri;
int *body = atom->body;
int *mask = atom->mask;
int nlocal = atom->nlocal;
@ -255,6 +292,45 @@ FixMove::FixMove(LAMMPS *lmp, int narg, char **arg) :
else xoriginal[i][0] = xoriginal[i][1] = xoriginal[i][2] = 0.0;
}
if (theta_flag) {
for (int i = 0; i < nlocal; i++) {
if ((mask[i] & groupbit) && line[i] >= 0)
toriginal[i] = avec_line->bonus[line[i]].theta;
else toriginal[i] = 0.0;
}
}
if (quat_flag) {
double *quat;
for (int i = 0; i < nlocal; i++) {
quat = NULL;
if (mask[i] & groupbit) {
if (ellipsoid_flag && ellipsoid[i] >= 0)
quat = avec_ellipsoid->bonus[ellipsoid[i]].quat;
else if (tri_flag && tri[i] >= 0)
quat = avec_tri->bonus[tri[i]].quat;
else if (body_flag && body[i] >= 0)
quat = avec_body->bonus[body[i]].quat;
}
if (quat) {
qoriginal[i][0] = quat[0];
qoriginal[i][1] = quat[1];
qoriginal[i][2] = quat[2];
qoriginal[i][3] = quat[3];
} else qoriginal[i][0] = qoriginal[i][1] =
qoriginal[i][2] = qoriginal[i][3] = 0.0;
}
}
// nrestart = size of per-atom restart data
// nrestart = 1 + xorig + torig + qorig
nrestart = 4;
if (theta_flag) nrestart++;
if (quat_flag) nrestart += 4;
// time origin for movement = current timestep
time_origin = update->ntimestep;
}
@ -270,6 +346,8 @@ FixMove::~FixMove()
// delete locally stored arrays
memory->destroy(xoriginal);
memory->destroy(toriginal);
memory->destroy(qoriginal);
memory->destroy(displace);
memory->destroy(velocity);
@ -381,9 +459,12 @@ void FixMove::init()
void FixMove::initial_integrate(int vflag)
{
double dtfm;
double xold[3],a[3],b[3],c[3],d[3],disp[3];
int flag;
double ddotr,dx,dy,dz;
double dtfm,theta_new;
double xold[3],a[3],b[3],c[3],d[3],disp[3],w[3],ex[3],ey[3],ez[3];
double inertia_ellipsoid[3],qrotate[4];
double *quat,*inertia,*shape;
double delta = (update->ntimestep - time_origin) * dt;
@ -391,10 +472,17 @@ void FixMove::initial_integrate(int vflag)
double **v = atom->v;
double **f = atom->f;
double **omega = atom->omega;
double **angmom = atom->angmom;
double *radius = atom->radius;
double *rmass = atom->rmass;
double *mass = atom->mass;
int *type = atom->type;
int *ellipsoid = atom->ellipsoid;
int *line = atom->line;
int *tri = atom->tri;
int *body = atom->body;
int *mask = atom->mask;
int nlocal = atom->nlocal;
// for linear: X = X0 + V*dt
@ -553,10 +641,79 @@ void FixMove::initial_integrate(int vflag)
v[i][0] = omega_rotate * (runit[1]*disp[2] - runit[2]*disp[1]);
v[i][1] = omega_rotate * (runit[2]*disp[0] - runit[0]*disp[2]);
v[i][2] = omega_rotate * (runit[0]*disp[1] - runit[1]*disp[0]);
if (omega_flag) {
omega[i][0] = omega_rotate*runit[0];
omega[i][1] = omega_rotate*runit[1];
omega[i][2] = omega_rotate*runit[2];
// set any extra attributes affected by rotation
if (extra_flag) {
// omega for spheres and lines
if (omega_flag) {
flag = 0;
if (radius_flag && radius[i] > 0.0) flag = 1;
if (line_flag && line[i] >= 0.0) flag = 1;
if (flag) {
omega[i][0] = omega_rotate*runit[0];
omega[i][1] = omega_rotate*runit[1];
omega[i][2] = omega_rotate*runit[2];
}
}
// angmom for ellipsoids, tris, and bodies
if (angmom_flag) {
quat = inertia = NULL;
if (ellipsoid_flag && ellipsoid[i] >= 0) {
quat = avec_ellipsoid->bonus[ellipsoid[i]].quat;
shape = avec_ellipsoid->bonus[ellipsoid[i]].shape;
inertia_ellipsoid[0] =
INERTIA*rmass[i] * (shape[1]*shape[1]+shape[2]*shape[2]);
inertia_ellipsoid[1] =
INERTIA*rmass[i] * (shape[0]*shape[0]+shape[2]*shape[2]);
inertia_ellipsoid[2] =
INERTIA*rmass[i] * (shape[0]*shape[0]+shape[1]*shape[1]);
inertia = inertia_ellipsoid;
} else if (tri_flag && tri[i] >= 0) {
quat = avec_tri->bonus[tri[i]].quat;
inertia = avec_tri->bonus[tri[i]].inertia;
} else if (body_flag && body[i] >= 0) {
quat = avec_body->bonus[body[i]].quat;
inertia = avec_body->bonus[body[i]].inertia;
}
if (quat) {
w[0] = omega_rotate*runit[0];
w[1] = omega_rotate*runit[1];
w[2] = omega_rotate*runit[2];
MathExtra::q_to_exyz(quat,ex,ey,ez);
MathExtra::omega_to_angmom(w,ex,ey,ez,inertia,angmom[i]);
}
}
// theta for lines
if (theta_flag && line[i] >= 0.0) {
theta_new = fmod(toriginal[i]+arg,MY_2PI);
avec_line->bonus[atom->line[i]].theta = theta_new;
}
// quats for ellipsoids, tris, and bodies
if (quat_flag) {
quat = NULL;
if (ellipsoid_flag && ellipsoid[i] >= 0)
quat = avec_ellipsoid->bonus[ellipsoid[i]].quat;
else if (tri_flag && tri[i] >= 0)
quat = avec_tri->bonus[tri[i]].quat;
else if (body_flag && body[i] >= 0)
quat = avec_body->bonus[body[i]].quat;
if (quat) {
qrotate[0] = cosine;
qrotate[1] = runit[0]*sine;
qrotate[2] = runit[1]*sine;
qrotate[3] = runit[2]*sine;
MathExtra::quatquat(qrotate,qoriginal[i],quat);
}
}
}
domain->remap_near(x[i],xold);
@ -564,6 +721,9 @@ void FixMove::initial_integrate(int vflag)
}
// for variable: compute x,v from variables
// NOTE: also allow for changes to extra attributes?
// omega, angmom, theta, quat
// only necessary if prescribed motion involves rotation
} else if (mstyle == VARIABLE) {
@ -809,6 +969,8 @@ void FixMove::final_integrate_respa(int ilevel, int iloop)
double FixMove::memory_usage()
{
double bytes = atom->nmax*3 * sizeof(double);
if (theta_flag) bytes += atom->nmax * sizeof(double);
if (quat_flag) bytes += atom->nmax*4 * sizeof(double);
if (displaceflag) bytes += atom->nmax*3 * sizeof(double);
if (velocityflag) bytes += atom->nmax*3 * sizeof(double);
return bytes;
@ -850,6 +1012,8 @@ void FixMove::restart(char *buf)
void FixMove::grow_arrays(int nmax)
{
memory->grow(xoriginal,nmax,3,"move:xoriginal");
if (theta_flag) memory->grow(toriginal,nmax,"move:toriginal");
if (quat_flag) memory->grow(qoriginal,nmax,4,"move:qoriginal");
array_atom = xoriginal;
}
@ -862,6 +1026,13 @@ void FixMove::copy_arrays(int i, int j, int delflag)
xoriginal[j][0] = xoriginal[i][0];
xoriginal[j][1] = xoriginal[i][1];
xoriginal[j][2] = xoriginal[i][2];
if (theta_flag) toriginal[j] = toriginal[i];
if (quat_flag) {
qoriginal[j][0] = qoriginal[i][0];
qoriginal[j][1] = qoriginal[i][1];
qoriginal[j][2] = qoriginal[i][2];
qoriginal[j][3] = qoriginal[i][3];
}
}
/* ----------------------------------------------------------------------
@ -870,8 +1041,15 @@ void FixMove::copy_arrays(int i, int j, int delflag)
void FixMove::set_arrays(int i)
{
double theta;
double *quat;
double **x = atom->x;
imageint *image = atom->image;
int *ellipsoid = atom->ellipsoid;
int *line = atom->line;
int *tri = atom->tri;
int *body = atom->body;
int *mask = atom->mask;
// particle not in group
@ -932,6 +1110,33 @@ void FixMove::set_arrays(int i)
xoriginal[i][0] = point[0] + c[0] + disp[0];
xoriginal[i][1] = point[1] + c[1] + disp[1];
xoriginal[i][2] = point[2] + c[2] + disp[2];
// set theta and quat extra attributes affected by rotation
if (extra_flag) {
// theta for lines
if (theta_flag && line[i] >= 0.0) {
theta = avec_line->bonus[atom->line[i]].theta;
toriginal[i] = theta - 0.0; // NOTE: edit this line
}
// quats for ellipsoids, tris, and bodies
if (quat_flag) {
quat = NULL;
if (ellipsoid_flag && ellipsoid[i] >= 0)
quat = avec_ellipsoid->bonus[ellipsoid[i]].quat;
else if (tri_flag && tri[i] >= 0)
quat = avec_tri->bonus[tri[i]].quat;
else if (body_flag && body[i] >= 0)
quat = avec_body->bonus[body[i]].quat;
if (quat) {
// qoriginal = f(quat,-delta); // NOTE: edit this line
}
}
}
}
}
@ -941,10 +1146,18 @@ void FixMove::set_arrays(int i)
int FixMove::pack_exchange(int i, double *buf)
{
buf[0] = xoriginal[i][0];
buf[1] = xoriginal[i][1];
buf[2] = xoriginal[i][2];
return 3;
int n = 0;
buf[n++] = xoriginal[i][0];
buf[n++] = xoriginal[i][1];
buf[n++] = xoriginal[i][2];
if (theta_flag) buf[n++] = toriginal[i];
if (quat_flag) {
buf[n++] = qoriginal[i][0];
buf[n++] = qoriginal[i][1];
buf[n++] = qoriginal[i][2];
buf[n++] = qoriginal[i][3];
}
return n;
}
/* ----------------------------------------------------------------------
@ -953,10 +1166,18 @@ int FixMove::pack_exchange(int i, double *buf)
int FixMove::unpack_exchange(int nlocal, double *buf)
{
xoriginal[nlocal][0] = buf[0];
xoriginal[nlocal][1] = buf[1];
xoriginal[nlocal][2] = buf[2];
return 3;
int n = 0;
xoriginal[nlocal][0] = buf[n++];
xoriginal[nlocal][1] = buf[n++];
xoriginal[nlocal][2] = buf[n++];
if (theta_flag) toriginal[nlocal] = buf[n++];
if (quat_flag) {
qoriginal[nlocal][0] = buf[n++];
qoriginal[nlocal][1] = buf[n++];
qoriginal[nlocal][2] = buf[n++];
qoriginal[nlocal][3] = buf[n++];
}
return n;
}
/* ----------------------------------------------------------------------
@ -965,11 +1186,19 @@ int FixMove::unpack_exchange(int nlocal, double *buf)
int FixMove::pack_restart(int i, double *buf)
{
buf[0] = 4;
buf[1] = xoriginal[i][0];
buf[2] = xoriginal[i][1];
buf[3] = xoriginal[i][2];
return 4;
int n = 1;
buf[n++] = xoriginal[i][0];
buf[n++] = xoriginal[i][1];
buf[n++] = xoriginal[i][2];
if (theta_flag) buf[n++] = toriginal[i];
if (quat_flag) {
buf[n++] = qoriginal[i][0];
buf[n++] = qoriginal[i][1];
buf[n++] = qoriginal[i][2];
buf[n++] = qoriginal[i][3];
}
buf[0] = n;
return n;
}
/* ----------------------------------------------------------------------
@ -989,6 +1218,13 @@ void FixMove::unpack_restart(int nlocal, int nth)
xoriginal[nlocal][0] = extra[nlocal][m++];
xoriginal[nlocal][1] = extra[nlocal][m++];
xoriginal[nlocal][2] = extra[nlocal][m++];
if (theta_flag) toriginal[nlocal] = extra[nlocal][m++];
if (quat_flag) {
qoriginal[nlocal][0] = extra[nlocal][m++];
qoriginal[nlocal][1] = extra[nlocal][m++];
qoriginal[nlocal][2] = extra[nlocal][m++];
qoriginal[nlocal][3] = extra[nlocal][m++];
}
}
/* ----------------------------------------------------------------------
@ -997,7 +1233,7 @@ void FixMove::unpack_restart(int nlocal, int nth)
int FixMove::maxsize_restart()
{
return 4;
return nrestart;
}
/* ----------------------------------------------------------------------
@ -1006,7 +1242,7 @@ int FixMove::maxsize_restart()
int FixMove::size_restart(int nlocal)
{
return 4;
return nrestart;
}
/* ---------------------------------------------------------------------- */

12
src/fix_move.h
View File

@ -61,13 +61,23 @@ class FixMove : public Fix {
double dt,dtv,dtf;
int xvar,yvar,zvar,vxvar,vyvar,vzvar;
int xvarstyle,yvarstyle,zvarstyle,vxvarstyle,vyvarstyle,vzvarstyle;
int omega_flag,nlevels_respa;
int extra_flag,omega_flag,angmom_flag;
int radius_flag,ellipsoid_flag,line_flag,tri_flag,body_flag;
int theta_flag,quat_flag;
int nlevels_respa,nrestart;
int time_origin;
double **xoriginal; // original coords of atoms
double *toriginal; // original theta of atoms
double **qoriginal; // original quat of atoms
int displaceflag,velocityflag;
int maxatom;
double **displace,**velocity;
class AtomVecEllipsoid *avec_ellipsoid;
class AtomVecLine *avec_line;
class AtomVecTri *avec_tri;
class AtomVecBody *avec_body;
};
}

1
src/math_extra.cpp
View File

@ -230,6 +230,7 @@ void richardson(double *q, double *m, double *w, double *moments, double dtq)
apply evolution operators to quat, quat momentum
Miller et al., J Chem Phys. 116, 8649-8659 (2002)
------------------------------------------------------------------------- */
void no_squish_rotate(int k, double *p, double *q, double *inertia,
double dt)
{

29
src/math_extra.h
View File

@ -76,7 +76,8 @@ namespace MathExtra {
void rotate(double matrix[3][3], int i, int j, int k, int l,
double s, double tau);
void richardson(double *q, double *m, double *w, double *moments, double dtq);
void no_squish_rotate(int k, double *p, double *q, double *inertia, double dt);
void no_squish_rotate(int k, double *p, double *q, double *inertia,
double dt);
// shape matrix operations
// upper-triangular 3x3 matrix stored in Voigt notation as 6-vector
@ -151,7 +152,8 @@ inline void MathExtra::normalize3(const double *v, double *ans)
scale a vector to length
------------------------------------------------------------------------- */
inline void MathExtra::snormalize3(const double length, const double *v, double *ans)
inline void MathExtra::snormalize3(const double length, const double *v,
double *ans)
{
double scale = length/sqrt(v[0]*v[0]+v[1]*v[1]+v[2]*v[2]);
ans[0] = v[0]*scale;
@ -276,7 +278,7 @@ inline double MathExtra::det3(const double m[3][3])
------------------------------------------------------------------------- */
inline void MathExtra::diag_times3(const double *d, const double m[3][3],
double ans[3][3])
double ans[3][3])
{
ans[0][0] = d[0]*m[0][0];
ans[0][1] = d[0]*m[0][1];
@ -312,7 +314,7 @@ void MathExtra::times3_diag(const double m[3][3], const double *d,
------------------------------------------------------------------------- */
inline void MathExtra::plus3(const double m[3][3], const double m2[3][3],
double ans[3][3])
double ans[3][3])
{
ans[0][0] = m[0][0]+m2[0][0];
ans[0][1] = m[0][1]+m2[0][1];
@ -330,7 +332,7 @@ inline void MathExtra::plus3(const double m[3][3], const double m2[3][3],
------------------------------------------------------------------------- */
inline void MathExtra::times3(const double m[3][3], const double m2[3][3],
double ans[3][3])
double ans[3][3])
{
ans[0][0] = m[0][0]*m2[0][0] + m[0][1]*m2[1][0] + m[0][2]*m2[2][0];
ans[0][1] = m[0][0]*m2[0][1] + m[0][1]*m2[1][1] + m[0][2]*m2[2][1];
@ -405,7 +407,8 @@ inline void MathExtra::invert3(const double m[3][3], double ans[3][3])
matrix times vector
------------------------------------------------------------------------- */
inline void MathExtra::matvec(const double m[3][3], const double *v, double *ans)
inline void MathExtra::matvec(const double m[3][3], const double *v,
double *ans)
{
ans[0] = m[0][0]*v[0] + m[0][1]*v[1] + m[0][2]*v[2];
ans[1] = m[1][0]*v[0] + m[1][1]*v[1] + m[1][2]*v[2];
@ -416,8 +419,8 @@ inline void MathExtra::matvec(const double m[3][3], const double *v, double *ans
matrix times vector
------------------------------------------------------------------------- */
inline void MathExtra::matvec(const double *ex, const double *ey, const double *ez,
const double *v, double *ans)
inline void MathExtra::matvec(const double *ex, const double *ey,
const double *ez, const double *v, double *ans)
{
ans[0] = ex[0]*v[0] + ey[0]*v[1] + ez[0]*v[2];
ans[1] = ex[1]*v[0] + ey[1]*v[1] + ez[1]*v[2];
@ -471,7 +474,8 @@ inline void MathExtra::transpose_diag3(const double m[3][3], const double *d,
row vector times matrix
------------------------------------------------------------------------- */
inline void MathExtra::vecmat(const double *v, const double m[3][3], double *ans)
inline void MathExtra::vecmat(const double *v, const double m[3][3],
double *ans)
{
ans[0] = v[0]*m[0][0] + v[1]*m[1][0] + v[2]*m[2][0];
ans[1] = v[0]*m[0][1] + v[1]*m[1][1] + v[2]*m[2][1];
@ -494,8 +498,8 @@ inline void MathExtra::scalar_times3(const double f, double m[3][3])
upper-triangular 3x3, stored as 6-vector in Voigt notation
------------------------------------------------------------------------- */
inline void MathExtra::multiply_shape_shape(const double *one, const double *two,
double *ans)
inline void MathExtra::multiply_shape_shape(const double *one,
const double *two, double *ans)
{
ans[0] = one[0]*two[0];
ans[1] = one[1]*two[1];
@ -601,7 +605,8 @@ inline void MathExtra::axisangle_to_quat(const double *v, const double angle,
Apply principal rotation generator about x to rotation matrix m
------------------------------------------------------------------------- */
inline void MathExtra::rotation_generator_x(const double m[3][3], double ans[3][3])
inline void MathExtra::rotation_generator_x(const double m[3][3],
double ans[3][3])
{
ans[0][0] = 0;
ans[0][1] = -m[0][2];

2
src/version.h
View File

@ -1 +1 @@
#define LAMMPS_VERSION "13 Jan 2016"
#define LAMMPS_VERSION "14 Jan 2016"