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

# Resolved Conflicts:
#	doc/Manual.html
#	doc/Manual.txt
#	src/USER-SMTBQ/pair_smtbq.cpp
This commit is contained in:
Axel Kohlmeyer
2016-01-12 07:02:29 -05:00
48 changed files with 1769 additions and 321 deletions

View File

@ -134,8 +134,8 @@
<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="dec-2015-version">
<h2>24 Dec 2015 version<a class="headerlink" href="#dec-2015-version" title="Permalink to this headline"></a></h2>
<div class="section" id="jan-2016-version">
<h2>11 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>

View File

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

View File

@ -171,7 +171,7 @@ style = <em>bin/1d</em> or <em>bin/2d</em> or <em>bin/3d</em> or <em>bin/sphere<
</pre>
<ul class="simple">
<li>zero or more keyword/values pairs may be appended</li>
<li>keyword = <em>region</em> or <em>nchunk</em> or <em>static</em> or <em>compress</em> or <em>bound</em> or <em>discard</em> or <em>units</em></li>
<li>keyword = <em>region</em> or <em>nchunk</em> or <em>static</em> or <em>compress</em> or <em>bound</em> or <em>discard</em> or <em>pbc</em> or <em>units</em></li>
</ul>
<pre class="literal-block">
<em>region</em> value = region-ID
@ -198,6 +198,8 @@ style = <em>bin/1d</em> or <em>bin/2d</em> or <em>bin/3d</em> or <em>bin/sphere<
x/y/z = <em>x</em> or <em>y</em> or <em>z</em> to bound sptial bins in this dimension
lo = <em>lower</em> or coordinate value (distance units)
hi = <em>upper</em> or coordinate value (distance units)
<em>pbc</em> value = <em>no</em> or <em>yes</em>
yes = use periodic distance for bin/sphere and bin/cylinder styles
<em>units</em> value = <em>box</em> or <em>lattice</em> or <em>reduced</em>
</pre>
</div>
@ -616,6 +618,19 @@ value, which is assumed to be inside (or at least near) the simulation
box boundaries, though LAMMPS does not check for this. Note that
using the <em>bound</em> keyword typically reduces the total number of bins
and thus the number of chunks <em>Nchunk</em>.</p>
<p>The <em>pbc</em> keyword only applies to the <em>bin/sphere</em> and <em>bin/cylinder</em>
styles. If set to <em>yes</em>, the distance an atom is from the sphere
origin or cylinder axis is calculated in a minimum image sense with
respect to periodic dimensions, when determining which bin the atom is
in. I.e. if x is a periodic dimension and the distance between the
atom and the sphere center in the x dimension is greater than 0.5 *
simulation box length in x, then a box length is subtracted to give a
distance &lt; 0.5 * simulation box length. This allosws the sphere or
cylinder center to be near a box edge, and atoms on the other side of
the periodic box will still be close to the center point/axis. Note
that with a setting of <em>yes</em>, the outer sphere or cylinder radius must
also be &lt;= 0.5 * simulation box length in any periodic dimension
except for the cylinder axis dimension, or an error is generated.</p>
<p>The <em>units</em> keyword only applies to the <em>binning</em> styles; otherwise it
is ignored. For the <em>bin/1d</em>, <em>bin/2d</em>, <em>bin/3d</em> styles, it
determines the meaning of the distance units used for the bin sizes
@ -683,6 +698,7 @@ the restarted simulation begins.</p>
<li>discard = yes, for all styles except binning</li>
<li>discard = mixed, for binning styles</li>
<li>bound = lower and upper in all dimensions</li>
<li>pbc = no</li>
<li>units = lattice</li>
</ul>
</div>

View File

@ -48,7 +48,7 @@ style = {bin/1d} or {bin/2d} or {bin/3d} or {bin/sphere} or {type} or {molecule}
v_name = per-atom vector calculated by an atom-style variable with name :pre
zero or more keyword/values pairs may be appended :l
keyword = {region} or {nchunk} or {static} or {compress} or {bound} or {discard} or {units} :l
keyword = {region} or {nchunk} or {static} or {compress} or {bound} or {discard} or {pbc} or {units} :l
{region} value = region-ID
region-ID = ID of region atoms must be in to be part of a chunk
{nchunk} value = {once} or {every}
@ -73,6 +73,8 @@ keyword = {region} or {nchunk} or {static} or {compress} or {bound} or {discard}
x/y/z = {x} or {y} or {z} to bound sptial bins in this dimension
lo = {lower} or coordinate value (distance units)
hi = {upper} or coordinate value (distance units)
{pbc} value = {no} or {yes}
yes = use periodic distance for bin/sphere and bin/cylinder styles
{units} value = {box} or {lattice} or {reduced} :pre
:ule
@ -570,6 +572,20 @@ box boundaries, though LAMMPS does not check for this. Note that
using the {bound} keyword typically reduces the total number of bins
and thus the number of chunks {Nchunk}.
The {pbc} keyword only applies to the {bin/sphere} and {bin/cylinder}
styles. If set to {yes}, the distance an atom is from the sphere
origin or cylinder axis is calculated in a minimum image sense with
respect to periodic dimensions, when determining which bin the atom is
in. I.e. if x is a periodic dimension and the distance between the
atom and the sphere center in the x dimension is greater than 0.5 *
simulation box length in x, then a box length is subtracted to give a
distance < 0.5 * simulation box length. This allosws the sphere or
cylinder center to be near a box edge, and atoms on the other side of
the periodic box will still be close to the center point/axis. Note
that with a setting of {yes}, the outer sphere or cylinder radius must
also be <= 0.5 * simulation box length in any periodic dimension
except for the cylinder axis dimension, or an error is generated.
The {units} keyword only applies to the {binning} styles; otherwise it
is ignored. For the {bin/1d}, {bin/2d}, {bin/3d} styles, it
determines the meaning of the distance units used for the bin sizes
@ -645,4 +661,5 @@ compress = no
discard = yes, for all styles except binning
discard = mixed, for binning styles
bound = lower and upper in all dimensions
pbc = no
units = lattice :ul

View File

@ -169,7 +169,7 @@ fix saed/vtk 1 1 1 c_2 file Ni_000.saed
<div class="section" id="description">
<h2>Description<a class="headerlink" href="#description" title="Permalink to this headline"></a></h2>
<p>Define a computation that calculates electron diffraction intensity as
described in <a class="reference internal" href="fix_saed_vtk.html#coleman"><span>(Coleman)</span></a> on a mesh of reciprocal lattice nodes
described in <a class="reference internal" href="compute_xrd.html#coleman"><span>(Coleman)</span></a> on a mesh of reciprocal lattice nodes
defined by the entire simulation domain (or manually) using simulated
radiation of wavelength lambda.</p>
<p>The electron diffraction intensity I at each reciprocal lattice point
@ -212,13 +212,14 @@ the <em>Kmax</em>, <em>Zone</em>, and <em>dR_Ewald</em> parameters. The rectili
created about the origin of reciprocal space is terminated at the
boundary of a sphere of radius <em>Kmax</em> centered at the origin. If
<em>Zone</em> parameters z1=z2=z3=0 are used, diffraction intensities are
computed throughout the entire spherical volume - note this can greatly
increase the cost of computation. Otherwise, <em>Zone</em> parameters will
denote the z1=h, z2=k, and z3=l (in a global since) zone axis of an
intersecting Ewald sphere. Diffraction intensities will only be
computed at the intersection of the reciprocal lattice mesh and a
<em>dR_Ewald</em> thick surface of the Ewald sphere. See the example 3D
intestiety data and the intersection of a [010] zone axis in the below image.</p>
computed throughout the entire spherical volume - note this can
greatly increase the cost of computation. Otherwise, <em>Zone</em>
parameters will denote the z1=h, z2=k, and z3=l (in a global since)
zone axis of an intersecting Ewald sphere. Diffraction intensities
will only be computed at the intersection of the reciprocal lattice
mesh and a <em>dR_Ewald</em> thick surface of the Ewald sphere. See the
example 3D intestiety data and the intersection of a [010] zone axis
in the below image.</p>
<a data-lightbox="group-default"
href="_images/saed_ewald_intersect.jpg"
class=""
@ -279,6 +280,8 @@ options.</p>
</div>
<div class="section" id="restrictions">
<h2>Restrictions<a class="headerlink" href="#restrictions" title="Permalink to this headline"></a></h2>
<p>This compute is part of the USER-DIFFRACTION package. It is only
enabled if LAMMPS was built with that package. See the <a class="reference internal" href="Section_start.html#start-3"><span>Making LAMMPS</span></a> section for more info.</p>
<p>The compute_saed command does not work for triclinic cells.</p>
</div>
<div class="section" id="related-commands">

View File

@ -87,13 +87,14 @@ the {Kmax}, {Zone}, and {dR_Ewald} parameters. The rectilinear mesh
created about the origin of reciprocal space is terminated at the
boundary of a sphere of radius {Kmax} centered at the origin. If
{Zone} parameters z1=z2=z3=0 are used, diffraction intensities are
computed throughout the entire spherical volume - note this can greatly
increase the cost of computation. Otherwise, {Zone} parameters will
denote the z1=h, z2=k, and z3=l (in a global since) zone axis of an
intersecting Ewald sphere. Diffraction intensities will only be
computed at the intersection of the reciprocal lattice mesh and a
{dR_Ewald} thick surface of the Ewald sphere. See the example 3D
intestiety data and the intersection of a \[010\] zone axis in the below image.
computed throughout the entire spherical volume - note this can
greatly increase the cost of computation. Otherwise, {Zone}
parameters will denote the z1=h, z2=k, and z3=l (in a global since)
zone axis of an intersecting Ewald sphere. Diffraction intensities
will only be computed at the intersection of the reciprocal lattice
mesh and a {dR_Ewald} thick surface of the Ewald sphere. See the
example 3D intestiety data and the intersection of a \[010\] zone axis
in the below image.
:c,image(JPG/saed_ewald_intersect_small.jpg,JPG/saed_ewald_intersect.jpg)
@ -151,6 +152,10 @@ All array values calculated by this compute are "intensive".
[Restrictions:]
This compute is part of the USER-DIFFRACTION package. It is only
enabled if LAMMPS was built with that package. See the "Making
LAMMPS"_Section_start.html#start_3 section for more info.
The compute_saed command does not work for triclinic cells.
[Related commands:]

View File

@ -167,7 +167,7 @@ fix 2 all ave/histo/weight 1 1 1 10 100 250 c_2[1] c_2[2] mode vector file Deg2T
<div class="section" id="description">
<h2>Description<a class="headerlink" href="#description" title="Permalink to this headline"></a></h2>
<p>Define a computation that calculates x-ray diffraction intensity as described
in <a class="reference internal" href="fix_saed_vtk.html#coleman"><span>(Coleman)</span></a> on a mesh of reciprocal lattice nodes defined
in <a class="reference internal" href="#coleman"><span>(Coleman)</span></a> on a mesh of reciprocal lattice nodes defined
by the entire simulation domain (or manually) using a simulated radiation
of wavelength lambda.</p>
<p>The x-ray diffraction intensity, I, at each reciprocal lattice point, k,
@ -297,13 +297,15 @@ which by the mesh. The global array has 2 columns.</p>
<p>The first column contains the diffraction angle in the units (radians
or degrees) provided with the <em>2Theta</em> values. The second column contains
the computed diffraction intensities as described above.</p>
<p>The array can be accessed by any command that uses global values
from a compute as input. See <a class="reference internal" href="Section_howto.html#howto-15"><span>this section</span></a> for an overview of LAMMPS output
options.</p>
<p>The array can be accessed by any command that uses global values from
a compute as input. See <a class="reference internal" href="Section_howto.html#howto-15"><span>this section</span></a>
for an overview of LAMMPS output options.</p>
<p>All array values calculated by this compute are &#8220;intensive&#8221;.</p>
</div>
<div class="section" id="restrictions">
<h2>Restrictions<a class="headerlink" href="#restrictions" title="Permalink to this headline"></a></h2>
<p>This compute is part of the USER-DIFFRACTION package. It is only
enabled if LAMMPS was built with that package. See the <a class="reference internal" href="Section_start.html#start-3"><span>Making LAMMPS</span></a> section for more info.</p>
<p>The compute_xrd command does not work for triclinic cells.</p>
</div>
<div class="section" id="related-commands">

View File

@ -163,14 +163,18 @@ The first column contains the diffraction angle in the units (radians
or degrees) provided with the {2Theta} values. The second column contains
the computed diffraction intensities as described above.
The array can be accessed by any command that uses global values
from a compute as input. See "this section"_Section_howto.html#howto_15 for an overview of LAMMPS output
options.
The array can be accessed by any command that uses global values from
a compute as input. See "this section"_Section_howto.html#howto_15
for an overview of LAMMPS output options.
All array values calculated by this compute are "intensive".
[Restrictions:]
This compute is part of the USER-DIFFRACTION package. It is only
enabled if LAMMPS was built with that package. See the "Making
LAMMPS"_Section_start.html#start_3 section for more info.
The compute_xrd command does not work for triclinic cells.
[Related commands:]

View File

@ -138,10 +138,10 @@
<li>input = one or more atom attributes</li>
</ul>
<div class="highlight-python"><div class="highlight"><pre>possible attributes = id, mol, type, mass,
x, y, z, xs, ys, zs, xu, yu, zu, ix, iy, iz,
x, y, z, xs, ys, zs, xu, yu, zu, xsu, ysu, zsu, ix, iy, iz,
vx, vy, vz, fx, fy, fz,
q, mux, muy, muz,
radius, omegax, omegay, omegaz,
q, mux, muy, muz, mu,
radius, diameter, omegax, omegay, omegaz,
angmomx, angmomy, angmomz, tqx, tqy, tqz,
c_ID, c_ID[N], f_ID, f_ID[N], v_name,
d_name, i_name
@ -154,12 +154,14 @@ mass = atom mass
x,y,z = unscaled atom coordinates
xs,ys,zs = scaled atom coordinates
xu,yu,zu = unwrapped atom coordinates
xsu,ysu,zsu = scaled unwrapped atom coordinates
ix,iy,iz = box image that the atom is in
vx,vy,vz = atom velocities
fx,fy,fz = forces on atoms
q = atom charge
mux,muy,muz = orientation of dipolar atom
radius = radius of spherical particle
mu = magnitued of dipole moment of atom
radius,diameter = radius.diameter of spherical particle
omegax,omegay,omegaz = angular velocity of spherical particle
angmomx,angmomy,angmomz = angular momentum of aspherical particle
tqx,tqy,tqz = torque on finite-size particles
@ -195,25 +197,23 @@ time the fix is defined. If <em>N</em> is 0, then the values are never
updated, so this is a way of archiving an atom attribute at a given
time for future use in a calculation or output. See the discussion of
<a class="reference internal" href="Section_howto.html#howto-15"><span>output commands</span></a> that take fixes as
inputs. And see for example, the <a class="reference internal" href="compute_reduce.html"><em>compute reduce</em></a>, <a class="reference internal" href="fix_ave_atom.html"><em>fix ave/atom</em></a>, <a class="reference internal" href="fix_ave_histo.html"><em>fix ave/histo</em></a>, <a class="reference internal" href="fix_ave_spatial.html"><em>fix ave/spatial</em></a>,
and <a class="reference internal" href="variable.html"><em>atom-style variable</em></a> commands.</p>
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 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> 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 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>
<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>
inputs store the position of each atom relative to the center-of-mass
of the group of atoms, instead of storing the absolute position. This
option is used by the <a class="reference internal" href="compute_msd.html"><em>compute msd</em></a> command.</p>
of the group of atoms, instead of storing the absolute position.</p>
<p>The requested values are stored in a per-atom vector or array as
discussed below. Zeroes are stored for atoms not in the specified
group.</p>

View File

@ -17,10 +17,10 @@ store/state = style name of this fix command :l
N = store atom attributes every N steps, N = 0 for initial store only :l
input = one or more atom attributes :l
possible attributes = id, mol, type, mass,
x, y, z, xs, ys, zs, xu, yu, zu, ix, iy, iz,
x, y, z, xs, ys, zs, xu, yu, zu, xsu, ysu, zsu, ix, iy, iz,
vx, vy, vz, fx, fy, fz,
q, mux, muy, muz,
radius, omegax, omegay, omegaz,
q, mux, muy, muz, mu,
radius, diameter, omegax, omegay, omegaz,
angmomx, angmomy, angmomz, tqx, tqy, tqz,
c_ID, c_ID\[N\], f_ID, f_ID\[N\], v_name,
d_name, i_name :pre
@ -32,12 +32,14 @@ input = one or more atom attributes :l
x,y,z = unscaled atom coordinates
xs,ys,zs = scaled atom coordinates
xu,yu,zu = unwrapped atom coordinates
xsu,ysu,zsu = scaled unwrapped atom coordinates
ix,iy,iz = box image that the atom is in
vx,vy,vz = atom velocities
fx,fy,fz = forces on atoms
q = atom charge
mux,muy,muz = orientation of dipolar atom
radius = radius of spherical particle
mu = magnitued of dipole moment of atom
radius,diameter = radius.diameter of spherical particle
omegax,omegay,omegaz = angular velocity of spherical particle
angmomx,angmomy,angmomz = angular momentum of aspherical particle
tqx,tqy,tqz = torque on finite-size particles
@ -67,29 +69,25 @@ time the fix is defined. If {N} is 0, then the values are never
updated, so this is a way of archiving an atom attribute at a given
time for future use in a calculation or output. See the discussion of
"output commands"_Section_howto.html#howto_15 that take fixes as
inputs. And see for example, the "compute
reduce"_compute_reduce.html, "fix ave/atom"_fix_ave_atom.html, "fix
ave/histo"_fix_ave_histo.html, "fix ave/spatial"_fix_ave_spatial.html,
and "atom-style variable"_variable.html commands.
inputs.
If {N} is not zero, then the attributes will be updated every {N}
steps.
NOTE: Actually, only atom attributes specified by keywords like {xu}
or {vy} 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.
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.
The list of possible attributes is the same as that used by the "dump
custom"_dump.html command, which describes their meaning.
If the {com} keyword is set to {yes} then the {xu}, {yu}, and {zu}
inputs store the position of each atom relative to the center-of-mass
of the group of atoms, instead of storing the absolute position. This
option is used by the "compute msd"_compute_msd.html command.
of the group of atoms, instead of storing the absolute position.
The requested values are stored in a per-atom vector or array as
discussed below. Zeroes are stored for atoms not in the specified

File diff suppressed because one or more lines are too long

View File

@ -139,7 +139,7 @@
<li><p class="first">one or more keyword/value pairs may be appended</p>
</li>
<li><dl class="first docutils">
<dt>keyword = <em>type</em> or <em>type/fraction</em> or <em>mol</em> or <em>x</em> or <em>y</em> or <em>z</em> or <em>charge</em> or <em>dipole</em> or <em>dipole/random</em> or <em>quat</em> or <em>quat/random</em> or <em>diameter</em> or <em>shape</em> or <em>length</em> or <em>tri</em> or <em>theta</em> or <em>angmom</em> or <em>omega</em> or <em>mass</em> or <em>density</em> or <em>volume</em> or <em>image</em> or</dt>
<dt>keyword = <em>type</em> or <em>type/fraction</em> or <em>mol</em> or <em>x</em> or <em>y</em> or <em>z</em> or <em>charge</em> or <em>dipole</em> or <em>dipole/random</em> or <em>quat</em> or <em>quat/random</em> or <em>diameter</em> or <em>shape</em> or <em>length</em> or <em>tri</em> or <em>theta</em> or <em>theta/random</em> or <em>angmom</em> or <em>omega</em> or <em>mass</em> or <em>density</em> or <em>volume</em> or <em>image</em> or</dt>
<dd><p class="first last"><em>bond</em> or <em>angle</em> or <em>dihedral</em> or <em>improper</em> or
<em>meso_e</em> or <em>meso_cv</em> or <em>meso_rho</em> or <em>smd_contact_radius</em> or <em>smd_mass_density</em> or <em>i_name</em> or <em>d_name</em></p>
</dd>
@ -184,6 +184,8 @@
<em>theta</em> value = angle (degrees)
angle = orientation of line segment with respect to x-axis
angle can be an atom-style variable (see below)
<em>theta/random</em> value = seed
seed = random # seed (positive integer) for line segment orienations
<em>angmom</em> values = Lx Ly Lz
Lx,Ly,Lz = components of angular momentum vector (distance-mass-velocity units)
any of Lx,Ly,Lz can be an atom-style variable (see below)
@ -318,19 +320,20 @@ vector to set as the orientation of the dipole moment vectors of the
selected atoms. The magnitude of the dipole moment is set
by the length of this orientation vector.</p>
<p>Keyword <em>dipole/random</em> randomizes the orientation of the dipole
moment vectors of the selected atoms and sets the magnitude of each to
the specified <em>Dlen</em> value. For 2d systems, the z component of the
moment vectors for the selected atoms and sets the magnitude of each
to the specified <em>Dlen</em> value. For 2d systems, the z component of the
orientation is set to 0.0. Random numbers are used in such a way that
the orientation of a particular atom is the same, regardless of how
many processors are being used. This keyword does not allow use of an
atom-style variable.</p>
<p>Keyword <em>quat</em> uses the specified values to create a quaternion
(4-vector) that represents the orientation of the selected atoms. The
particles must be ellipsoids as defined by the <a class="reference internal" href="atom_style.html"><em>atom_style ellipsoid</em></a> command or triangles as defined by the
<a class="reference internal" href="atom_style.html"><em>atom_style tri</em></a> command. Note that particles defined
by <a class="reference internal" href="atom_style.html"><em>atom_style ellipsoid</em></a> have 3 shape parameters.
The 3 values must be non-zero for each particle set by this command.
They are used to specify the aspect ratios of an ellipsoidal particle,
particles must define a quaternion for their orientation
(e.g. ellipsoids, triangles, body particles) as defined by the
<a class="reference internal" href="atom_style.html"><em>atom_style</em></a> command. Note that particles defined by
<a class="reference internal" href="atom_style.html"><em>atom_style ellipsoid</em></a> have 3 shape parameters. The 3
values must be non-zero for each particle set by this command. They
are used to specify the aspect ratios of an ellipsoidal particle,
which is oriented by default with its x-axis along the simulation
box&#8217;s x-axis, and similarly for y and z. If this body is rotated (via
the right-hand rule) by an angle theta around a unit rotation vector
@ -340,16 +343,16 @@ c*sin(theta/2)). The theta and a,b,c values are the arguments to the
<em>quat</em> keyword. LAMMPS normalizes the quaternion in case (a,b,c) was
not specified as a unit vector. For 2d systems, the a,b,c values are
ignored, since a rotation vector of (0,0,1) is the only valid choice.</p>
<p>Keyword <em>quat/random</em> randomizes the orientation of the quaternion of
the selected atoms. The particles must be ellipsoids as defined by
the <a class="reference internal" href="atom_style.html"><em>atom_style ellipsoid</em></a> command or triangles as
defined by the <a class="reference internal" href="atom_style.html"><em>atom_style tri</em></a> command. Random
numbers are used in such a way that the orientation of a particular
atom is the same, regardless of how many processors are being used.
For 2d systems, only orientations in the xy plane are generated. As
with keyword <em>quat</em>, for ellipsoidal particles, the 3 shape values
must be non-zero for each particle set by this command. This keyword
does not allow use of an atom-style variable.</p>
<p>Keyword <em>quat/random</em> randomizes the orientation of the quaternion for
the selected atoms. The particles must define a quaternion for their
orientation (e.g. ellipsoids, triangles, body particles) as defined by
the <a class="reference internal" href="atom_style.html"><em>atom_style</em></a> command. Random numbers are used in
such a way that the orientation of a particular atom is the same,
regardless of how many processors are being used. For 2d systems,
only orientations in the xy plane are generated. As with keyword
<em>quat</em>, for ellipsoidal particles, the 3 shape values must be non-zero
for each particle set by this command. This keyword does not allow
use of an atom-style variable.</p>
<p>Keyword <em>diameter</em> sets the size of the selected atoms. The particles
must be finite-size spheres as defined by the <a class="reference internal" href="atom_style.html"><em>atom_style sphere</em></a> command. The diameter of a particle can be
set to 0.0, which means they will be treated as point particles. Note
@ -385,6 +388,12 @@ density, e.g. via the <a class="reference internal" href="read_data.html"><em>re
<p>Keyword <em>theta</em> sets the orientation of selected atoms. The particles
must be line segments as defined by the <a class="reference internal" href="atom_style.html"><em>atom_style line</em></a> command. The specified value is used to set the
orientation angle of the line segments with respect to the x axis.</p>
<p>Keyword <em>theta/random</em> randomizes the orientation of theta for the
selected atoms. The particles must be line segments as defined by the
<a class="reference internal" href="atom_style.html"><em>atom_style line</em></a> command. Random numbers are used in
such a way that the orientation of a particular atom is the same,
regardless of how many processors are being used. This keyword does
not allow use of an atom-style variable.</p>
<p>Keyword <em>angmom</em> sets the angular momentum of selected atoms. The
particles must be ellipsoids as defined by the <a class="reference internal" href="atom_style.html"><em>atom_style ellipsoid</em></a> command or triangles as defined by the
<a class="reference internal" href="atom_style.html"><em>atom_style tri</em></a> command. The angular momentum vector

View File

@ -18,7 +18,7 @@ one or more keyword/value pairs may be appended :l
keyword = {type} or {type/fraction} or {mol} or {x} or {y} or {z} or \
{charge} or {dipole} or {dipole/random} or {quat} or \
{quat/random} or {diameter} or {shape} or \
{length} or {tri} or {theta} or {angmom} or {omega} or \
{length} or {tri} or {theta} or {theta/random} or {angmom} or {omega} or \
{mass} or {density} or {volume} or {image} or
{bond} or {angle} or {dihedral} or {improper} or
{meso_e} or {meso_cv} or {meso_rho} or \
@ -61,6 +61,8 @@ keyword = {type} or {type/fraction} or {mol} or {x} or {y} or {z} or \
{theta} value = angle (degrees)
angle = orientation of line segment with respect to x-axis
angle can be an atom-style variable (see below)
{theta/random} value = seed
seed = random # seed (positive integer) for line segment orienations
{angmom} values = Lx Ly Lz
Lx,Ly,Lz = components of angular momentum vector (distance-mass-velocity units)
any of Lx,Ly,Lz can be an atom-style variable (see below)
@ -209,8 +211,8 @@ selected atoms. The magnitude of the dipole moment is set
by the length of this orientation vector.
Keyword {dipole/random} randomizes the orientation of the dipole
moment vectors of the selected atoms and sets the magnitude of each to
the specified {Dlen} value. For 2d systems, the z component of the
moment vectors for the selected atoms and sets the magnitude of each
to the specified {Dlen} value. For 2d systems, the z component of the
orientation is set to 0.0. Random numbers are used in such a way that
the orientation of a particular atom is the same, regardless of how
many processors are being used. This keyword does not allow use of an
@ -218,12 +220,12 @@ atom-style variable.
Keyword {quat} uses the specified values to create a quaternion
(4-vector) that represents the orientation of the selected atoms. The
particles must be ellipsoids as defined by the "atom_style
ellipsoid"_atom_style.html command or triangles as defined by the
"atom_style tri"_atom_style.html command. Note that particles defined
by "atom_style ellipsoid"_atom_style.html have 3 shape parameters.
The 3 values must be non-zero for each particle set by this command.
They are used to specify the aspect ratios of an ellipsoidal particle,
particles must define a quaternion for their orientation
(e.g. ellipsoids, triangles, body particles) as defined by the
"atom_style"_atom_style.html command. Note that particles defined by
"atom_style ellipsoid"_atom_style.html have 3 shape parameters. The 3
values must be non-zero for each particle set by this command. They
are used to specify the aspect ratios of an ellipsoidal particle,
which is oriented by default with its x-axis along the simulation
box's x-axis, and similarly for y and z. If this body is rotated (via
the right-hand rule) by an angle theta around a unit rotation vector
@ -234,16 +236,16 @@ c*sin(theta/2)). The theta and a,b,c values are the arguments to the
not specified as a unit vector. For 2d systems, the a,b,c values are
ignored, since a rotation vector of (0,0,1) is the only valid choice.
Keyword {quat/random} randomizes the orientation of the quaternion of
the selected atoms. The particles must be ellipsoids as defined by
the "atom_style ellipsoid"_atom_style.html command or triangles as
defined by the "atom_style tri"_atom_style.html command. Random
numbers are used in such a way that the orientation of a particular
atom is the same, regardless of how many processors are being used.
For 2d systems, only orientations in the xy plane are generated. As
with keyword {quat}, for ellipsoidal particles, the 3 shape values
must be non-zero for each particle set by this command. This keyword
does not allow use of an atom-style variable.
Keyword {quat/random} randomizes the orientation of the quaternion for
the selected atoms. The particles must define a quaternion for their
orientation (e.g. ellipsoids, triangles, body particles) as defined by
the "atom_style"_atom_style.html command. Random numbers are used in
such a way that the orientation of a particular atom is the same,
regardless of how many processors are being used. For 2d systems,
only orientations in the xy plane are generated. As with keyword
{quat}, for ellipsoidal particles, the 3 shape values must be non-zero
for each particle set by this command. This keyword does not allow
use of an atom-style variable.
Keyword {diameter} sets the size of the selected atoms. The particles
must be finite-size spheres as defined by the "atom_style
@ -289,6 +291,13 @@ must be line segments as defined by the "atom_style
line"_atom_style.html command. The specified value is used to set the
orientation angle of the line segments with respect to the x axis.
Keyword {theta/random} randomizes the orientation of theta for the
selected atoms. The particles must be line segments as defined by the
"atom_style line"_atom_style.html command. Random numbers are used in
such a way that the orientation of a particular atom is the same,
regardless of how many processors are being used. This keyword does
not allow use of an atom-style variable.
Keyword {angmom} sets the angular momentum of selected atoms. The
particles must be ellipsoids as defined by the "atom_style
ellipsoid"_atom_style.html command or triangles as defined by the

View File

@ -1,4 +1,4 @@
LAMMPS (23 Oct 2015)
LAMMPS (24 Dec 2015)
# Al2O3 crystal, qeq on, minimizes, then calculates elastic constants
variable T_depart equal 300
@ -64,6 +64,7 @@ velocity all create 300 277387
pair_style smtbq
pair_coeff * * ffield.smtbq.Al Al
Reading potential file ffield.smtbq.Al with DATE: 2015-10-22
neighbor 0.5 bin
neigh_modify every 20 delay 0 check yes
@ -100,20 +101,20 @@ Step Temp Press PotEng KinEng TotEng Lx Ly Lz Volume
8 298.92698 770.098 -5600.6212 64.875459 -5535.7458 28.637825 34.721517 28.059223 27900.653
9 298.64244 780.92261 -5600.5595 64.813707 -5535.7458 28.637825 34.721517 28.059223 27900.653
10 298.32468 793.00943 -5600.4905 64.744743 -5535.7458 28.637825 34.721517 28.059223 27900.653
Loop time of 5.10336 on 1 procs for 10 steps with 1680 atoms
Loop time of 5.10475 on 1 procs for 10 steps with 1680 atoms
Performance: 0.034 ns/day, 708.800 hours/ns, 1.959 timesteps/s
Performance: 0.034 ns/day, 708.993 hours/ns, 1.959 timesteps/s
99.9% CPU use with 1 MPI tasks x no OpenMP threads
MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total
---------------------------------------------------------------
Pair | 5.1023 | 5.1023 | 5.1023 | 0.0 | 99.98
Pair | 5.1037 | 5.1037 | 5.1037 | 0.0 | 99.98
Neigh | 0 | 0 | 0 | 0.0 | 0.00
Comm | 0.00045538 | 0.00045538 | 0.00045538 | 0.0 | 0.01
Output | 0.00024509 | 0.00024509 | 0.00024509 | 0.0 | 0.00
Modify | 0.00026131 | 0.00026131 | 0.00026131 | 0.0 | 0.01
Other | | 0.0001056 | | | 0.00
Comm | 0.00044179 | 0.00044179 | 0.00044179 | 0.0 | 0.01
Output | 0.00024772 | 0.00024772 | 0.00024772 | 0.0 | 0.00
Modify | 0.00026011 | 0.00026011 | 0.00026011 | 0.0 | 0.01
Other | | 0.0001063 | | | 0.00
Nlocal: 1680 ave 1680 max 1680 min
Histogram: 1 0 0 0 0 0 0 0 0 0
@ -186,7 +187,7 @@ Step Temp Press PotEng KinEng TotEng Lx Ly Lz Volume
50 298.32468 2486.1481 -5600.8719 64.744743 -5536.1271 28.615167 34.699435 28.038181 27839.955
51 298.32468 2498.3384 -5600.8721 64.744743 -5536.1274 28.616069 34.696071 28.039596 27839.537
52 298.32468 2443.1247 -5600.8722 64.744743 -5536.1274 28.617616 34.696737 28.039439 27841.421
Loop time of 36.9389 on 1 procs for 42 steps with 1680 atoms
Loop time of 36.8807 on 1 procs for 42 steps with 1680 atoms
99.9% CPU use with 1 MPI tasks x no OpenMP threads
@ -202,12 +203,12 @@ Minimization stats:
MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total
---------------------------------------------------------------
Pair | 36.918 | 36.918 | 36.918 | 0.0 | 99.94
Pair | 36.86 | 36.86 | 36.86 | 0.0 | 99.95
Neigh | 0 | 0 | 0 | 0.0 | 0.00
Comm | 0.0033035 | 0.0033035 | 0.0033035 | 0.0 | 0.01
Output | 0.0011785 | 0.0011785 | 0.0011785 | 0.0 | 0.00
Comm | 0.0032628 | 0.0032628 | 0.0032628 | 0.0 | 0.01
Output | 0.001102 | 0.001102 | 0.001102 | 0.0 | 0.00
Modify | 0 | 0 | 0 | 0.0 | 0.00
Other | | 0.01603 | | | 0.04
Other | | 0.0159 | | | 0.04
Nlocal: 1680 ave 1680 max 1680 min
Histogram: 1 0 0 0 0 0 0 0 0 0
@ -247,20 +248,20 @@ Step Temp Press PotEng KinEng TotEng Lx Ly Lz Volume
60 297.05416 2489.245 -5600.5964 64.469006 -5536.1274 28.617616 34.696737 28.039439 27841.421
61 296.74469 2500.7159 -5600.5293 64.401842 -5536.1274 28.617616 34.696737 28.039439 27841.421
62 296.40194 2513.4457 -5600.4549 64.327455 -5536.1274 28.617616 34.696737 28.039439 27841.421
Loop time of 5.13028 on 1 procs for 10 steps with 1680 atoms
Loop time of 5.13584 on 1 procs for 10 steps with 1680 atoms
Performance: 0.034 ns/day, 712.539 hours/ns, 1.949 timesteps/s
Performance: 0.034 ns/day, 713.311 hours/ns, 1.947 timesteps/s
99.9% CPU use with 1 MPI tasks x no OpenMP threads
MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total
---------------------------------------------------------------
Pair | 5.1292 | 5.1292 | 5.1292 | 0.0 | 99.98
Pair | 5.1346 | 5.1346 | 5.1346 | 0.0 | 99.98
Neigh | 0 | 0 | 0 | 0.0 | 0.00
Comm | 0.00045729 | 0.00045729 | 0.00045729 | 0.0 | 0.01
Output | 0.00024509 | 0.00024509 | 0.00024509 | 0.0 | 0.00
Modify | 0.00026488 | 0.00026488 | 0.00026488 | 0.0 | 0.01
Other | | 0.0001016 | | | 0.00
Comm | 0.00047112 | 0.00047112 | 0.00047112 | 0.0 | 0.01
Output | 0.00036144 | 0.00036144 | 0.00036144 | 0.0 | 0.01
Modify | 0.00029016 | 0.00029016 | 0.00029016 | 0.0 | 0.01
Other | | 0.0001101 | | | 0.00
Nlocal: 1680 ave 1680 max 1680 min
Histogram: 1 0 0 0 0 0 0 0 0 0

View File

@ -0,0 +1,278 @@
LAMMPS (24 Dec 2015)
# Al2O3 crystal, qeq on, minimizes, then calculates elastic constants
variable T_depart equal 300
variable dt equal 0.0002
#Constante
variable rac3 equal sqrt(3.0)
variable rac1_2 equal sqrt(0.5)
variable rac3_2 equal sqrt(1.5)
#Structure
variable a equal 4.05
variable nx equal 10
variable ny equal 7
variable nz equal 4
variable bx equal ${a}*${nx}*${rac1_2}
variable bx equal 4.05*${nx}*${rac1_2}
variable bx equal 4.05*10*${rac1_2}
variable bx equal 4.05*10*0.707106781186548
variable by equal ${a}*${ny}*${rac3_2}
variable by equal 4.05*${ny}*${rac3_2}
variable by equal 4.05*7*${rac3_2}
variable by equal 4.05*7*1.22474487139159
variable bz equal ${a}*${nz}*${rac3}
variable bz equal 4.05*${nz}*${rac3}
variable bz equal 4.05*4*${rac3}
variable bz equal 4.05*4*1.73205080756888
# =======================================================================
units metal
atom_style charge
dimension 3
boundary p p p
lattice sc 1.0
Lattice spacing in x,y,z = 1 1 1
region box_vide prism 0 ${bx} 0 ${by} 0 ${bz} 0.0 0.0 0.0
region box_vide prism 0 28.6378246380552 0 ${by} 0 ${bz} 0.0 0.0 0.0
region box_vide prism 0 28.6378246380552 0 34.7215171039516 0 ${bz} 0.0 0.0 0.0
region box_vide prism 0 28.6378246380552 0 34.7215171039516 0 28.0592230826159 0.0 0.0 0.0
create_box 1 box_vide
Created triclinic box = (0 0 0) to (28.6378 34.7215 28.0592) with tilt (0 0 0)
2 by 2 by 1 MPI processor grid
# Aluminium atoms z = [111]
lattice custom ${a} a1 ${rac1_2} 0.0 0.0 a2 0.0 ${rac3_2} 0.0 a3 0.0 0.0 ${rac3} basis 0.0 0.0 0.0 basis 0.5 0.5 0.0 basis 0.5 0.166666667 0.33333 basis 0.0 0.666666667 0.33333 basis 0.0 0.333333333 0.66667 basis 0.5 0.833333333 0.66667
lattice custom 4.05 a1 ${rac1_2} 0.0 0.0 a2 0.0 ${rac3_2} 0.0 a3 0.0 0.0 ${rac3} basis 0.0 0.0 0.0 basis 0.5 0.5 0.0 basis 0.5 0.166666667 0.33333 basis 0.0 0.666666667 0.33333 basis 0.0 0.333333333 0.66667 basis 0.5 0.833333333 0.66667
lattice custom 4.05 a1 0.707106781186548 0.0 0.0 a2 0.0 ${rac3_2} 0.0 a3 0.0 0.0 ${rac3} basis 0.0 0.0 0.0 basis 0.5 0.5 0.0 basis 0.5 0.166666667 0.33333 basis 0.0 0.666666667 0.33333 basis 0.0 0.333333333 0.66667 basis 0.5 0.833333333 0.66667
lattice custom 4.05 a1 0.707106781186548 0.0 0.0 a2 0.0 1.22474487139159 0.0 a3 0.0 0.0 ${rac3} basis 0.0 0.0 0.0 basis 0.5 0.5 0.0 basis 0.5 0.166666667 0.33333 basis 0.0 0.666666667 0.33333 basis 0.0 0.333333333 0.66667 basis 0.5 0.833333333 0.66667
lattice custom 4.05 a1 0.707106781186548 0.0 0.0 a2 0.0 1.22474487139159 0.0 a3 0.0 0.0 1.73205080756888 basis 0.0 0.0 0.0 basis 0.5 0.5 0.0 basis 0.5 0.166666667 0.33333 basis 0.0 0.666666667 0.33333 basis 0.0 0.333333333 0.66667 basis 0.5 0.833333333 0.66667
Lattice spacing in x,y,z = 2.86378 4.96022 7.01481
create_atoms 1 region box_vide
Created 1680 atoms
mass 1 26.98
velocity all create ${T_depart} 277387
velocity all create 300 277387
pair_style smtbq
pair_coeff * * ffield.smtbq.Al Al
Reading potential file ffield.smtbq.Al with DATE: 2015-10-22
neighbor 0.5 bin
neigh_modify every 20 delay 0 check yes
timestep ${dt}
timestep 0.0002
thermo_style custom step temp press pe ke etotal lx ly lz vol
thermo_modify flush yes
thermo 1
#dump 5 all custom 1 box_Al.lammpstrj id type q x y z
fix 3 all nve
run 10
Neighbor list info ...
1 neighbor list requests
update every 20 steps, delay 0 steps, check yes
max neighbors/atom: 2000, page size: 100000
master list distance cutoff = 11.6714
ghost atom cutoff = 11.6714
binsize = 5.8357 -> bins = 5 6 5
Memory usage per processor = 4.2858 Mbytes
Step Temp Press PotEng KinEng TotEng Lx Ly Lz Volume
0 300 729.26605 -5600.8541 65.108335 -5535.7458 28.637825 34.721517 28.059223 27900.653
1 299.98339 729.897 -5600.8505 65.104731 -5535.7458 28.637825 34.721517 28.059223 27900.653
2 299.93356 731.7909 -5600.8397 65.093915 -5535.7458 28.637825 34.721517 28.059223 27900.653
3 299.8505 734.94708 -5600.8217 65.07589 -5535.7458 28.637825 34.721517 28.059223 27900.653
4 299.73426 739.36449 -5600.7964 65.050662 -5535.7458 28.637825 34.721517 28.059223 27900.653
5 299.58486 745.04171 -5600.764 65.018238 -5535.7458 28.637825 34.721517 28.059223 27900.653
6 299.40234 751.97694 -5600.7244 64.978627 -5535.7458 28.637825 34.721517 28.059223 27900.653
7 299.18677 760.16798 -5600.6776 64.931841 -5535.7458 28.637825 34.721517 28.059223 27900.653
8 298.9382 769.6122 -5600.6237 64.877894 -5535.7458 28.637825 34.721517 28.059223 27900.653
9 298.65669 780.3065 -5600.5626 64.8168 -5535.7458 28.637825 34.721517 28.059223 27900.653
10 298.34234 792.24735 -5600.4943 64.748578 -5535.7458 28.637825 34.721517 28.059223 27900.653
Loop time of 1.49703 on 4 procs for 10 steps with 1680 atoms
Performance: 0.115 ns/day, 207.920 hours/ns, 6.680 timesteps/s
98.4% CPU use with 4 MPI tasks x no OpenMP threads
MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total
---------------------------------------------------------------
Pair | 1.4956 | 1.4957 | 1.4957 | 0.0 | 99.91
Neigh | 0 | 0 | 0 | 0.0 | 0.00
Comm | 0.00078988 | 0.00090927 | 0.0010064 | 0.3 | 0.06
Output | 0.00022459 | 0.00027347 | 0.00034189 | 0.3 | 0.02
Modify | 6.4373e-05 | 7.8619e-05 | 9.3222e-05 | 0.1 | 0.01
Other | | 8.106e-05 | | | 0.01
Nlocal: 420 ave 484 max 360 min
Histogram: 1 0 1 0 0 0 1 0 0 1
Nghost: 4305 ave 4365 max 4241 min
Histogram: 1 0 0 1 0 0 0 1 0 1
Neighs: 0 ave 0 max 0 min
Histogram: 4 0 0 0 0 0 0 0 0 0
FullNghs: 159600 ave 183920 max 136800 min
Histogram: 1 0 1 0 0 0 1 0 0 1
Total # of neighbors = 638400
Ave neighs/atom = 380
Neighbor list builds = 0
Dangerous builds = 0
unfix 3
#thermo 15
fix 1 all box/relax tri 0.0 vmax 0.001
minimize 1.0e-8 1.0e-10 1000 10000
WARNING: Resetting reneighboring criteria during minimization (../min.cpp:168)
Neighbor list info ...
1 neighbor list requests
update every 1 steps, delay 0 steps, check yes
max neighbors/atom: 2000, page size: 100000
master list distance cutoff = 11.6714
ghost atom cutoff = 11.6714
binsize = 5.8357 -> bins = 5 6 5
Memory usage per processor = 5.4127 Mbytes
Step Temp Press PotEng KinEng TotEng Lx Ly Lz Volume
10 298.34234 792.24735 -5600.4943 64.748578 -5535.7458 28.637825 34.721517 28.059223 27900.653
11 298.34234 2483.41 -5600.5251 64.748578 -5535.7765 28.617939 34.697415 28.039749 27842.588
12 298.34234 3212.2068 -5600.6294 64.748578 -5535.8809 28.605195 34.685635 28.035928 27816.95
13 298.34234 1756.312 -5600.6529 64.748578 -5535.9044 28.62806 34.701751 28.05036 27866.456
14 298.34234 845.30512 -5600.6668 64.748578 -5535.9182 28.652683 34.704776 28.055059 27897.528
15 298.34234 1950.7901 -5600.7433 64.748578 -5535.9948 28.66644 34.670055 28.031028 27859.115
16 298.34234 2860.9521 -5600.7579 64.748578 -5536.0093 28.637802 34.669634 28.028024 27827.963
17 298.34234 2406.0293 -5600.7642 64.748578 -5536.0156 28.624953 34.689277 28.040313 27843.439
18 298.34234 2531.6658 -5600.7645 64.748578 -5536.0159 28.619984 34.690929 28.039526 27839.151
19 298.34234 2410.8263 -5600.7651 64.748578 -5536.0166 28.615355 34.70117 28.039934 27843.27
20 298.34234 2538.4705 -5600.7654 64.748578 -5536.0169 28.613142 34.700755 28.038049 27838.913
21 298.34234 2455.2909 -5600.7657 64.748578 -5536.0171 28.615273 34.700334 28.039158 27841.749
22 298.34234 2578.5637 -5600.766 64.748578 -5536.0174 28.614997 34.695849 28.038814 27837.541
23 298.34234 2413.5858 -5600.7667 64.748578 -5536.0181 28.617591 34.695982 28.041829 27843.164
24 298.34234 2575.34 -5600.7677 64.748578 -5536.0191 28.614711 34.696953 28.038299 27837.636
25 298.34234 2313.5501 -5600.7689 64.748578 -5536.0203 28.618801 34.702208 28.039031 27846.559
26 298.34234 2589.6554 -5600.7724 64.748578 -5536.0239 28.622807 34.692685 28.033282 27837.106
27 298.34234 2258.4841 -5600.7755 64.748578 -5536.0269 28.633037 34.683012 28.042438 27848.382
28 298.34234 2668.7002 -5600.7789 64.748578 -5536.0303 28.632179 34.665065 28.043666 27834.358
29 298.34234 2339.309 -5600.8141 64.748578 -5536.0655 28.626885 34.675842 28.051214 27845.356
30 298.34234 2846.4227 -5600.8223 64.748578 -5536.0737 28.612492 34.69364 28.033449 27828.006
31 298.34234 2534.6985 -5600.8246 64.748578 -5536.076 28.611965 34.709562 28.031793 27838.619
32 298.34234 2278.6053 -5600.8262 64.748578 -5536.0776 28.617042 34.703828 28.040241 27847.349
33 298.34234 2477.381 -5600.8269 64.748578 -5536.0783 28.616395 34.695005 28.041167 27840.56
34 298.34234 2561.8128 -5600.8271 64.748578 -5536.0785 28.615917 34.69647 28.037549 27837.678
35 298.34234 2470.4145 -5600.8273 64.748578 -5536.0787 28.617208 34.699332 28.037111 27840.795
36 298.34234 2415.0452 -5600.8276 64.748578 -5536.079 28.617651 34.69562 28.041578 27842.683
37 298.34234 2514.0703 -5600.8277 64.748578 -5536.0791 28.616929 34.693138 28.040886 27839.303
38 298.34234 2515.7338 -5600.828 64.748578 -5536.0794 28.619947 34.694884 28.03646 27839.244
39 298.34234 2402.2948 -5600.8281 64.748578 -5536.0796 28.620623 34.69709 28.037911 27843.113
40 298.34234 2457.1081 -5600.8289 64.748578 -5536.0804 28.612617 34.697892 28.043218 27841.237
41 298.34234 2729.4589 -5600.8299 64.748578 -5536.0814 28.610746 34.694899 28.038106 27831.941
42 298.34234 2487.8313 -5600.8331 64.748578 -5536.0845 28.620239 34.701333 28.031882 27840.157
43 298.34234 1959.1329 -5600.8393 64.748578 -5536.0908 28.609631 34.719775 28.045509 27858.164
44 298.34234 2535.5621 -5600.8461 64.748578 -5536.0975 28.587496 34.720535 28.046731 27838.433
45 298.34234 2978.3628 -5600.858 64.748578 -5536.1095 28.609197 34.699148 28.027451 27823.264
46 298.34234 2351.5532 -5600.8632 64.748578 -5536.1146 28.624986 34.701661 28.031447 27844.605
47 298.34234 2294.1629 -5600.8704 64.748578 -5536.1218 28.614261 34.700247 28.045027 27846.522
48 298.34234 2528.0548 -5600.871 64.748578 -5536.1225 28.612265 34.696612 28.041878 27838.536
49 298.34234 2516.2649 -5600.8715 64.748578 -5536.1229 28.618047 34.695729 28.037328 27838.936
50 298.34234 2461.5711 -5600.8715 64.748578 -5536.1229 28.61865 34.696255 28.038191 27840.802
Loop time of 10.0107 on 4 procs for 40 steps with 1680 atoms
98.5% CPU use with 4 MPI tasks x no OpenMP threads
Minimization stats:
Stopping criterion = energy tolerance
Energy initial, next-to-last, final =
-5600.49433365 -5600.87146764 -5600.87150075
Force two-norm initial, final = 50.9055 1.01476
Force max component initial, final = 29.378 0.811083
Final line search alpha, max atom move = 0.00128157 0.00103946
Iterations, force evaluations = 40 75
MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total
---------------------------------------------------------------
Pair | 9.9932 | 9.9933 | 9.9933 | 0.0 | 99.83
Neigh | 0 | 0 | 0 | 0.0 | 0.00
Comm | 0.0059471 | 0.0059789 | 0.0059991 | 0.0 | 0.06
Output | 0.00088 | 0.0010666 | 0.0015972 | 0.9 | 0.01
Modify | 0 | 0 | 0 | 0.0 | 0.00
Other | | 0.01039 | | | 0.10
Nlocal: 420 ave 426 max 416 min
Histogram: 1 1 0 0 0 1 0 0 0 1
Nghost: 4305 ave 4309 max 4299 min
Histogram: 1 0 0 0 0 1 0 0 0 2
Neighs: 0 ave 0 max 0 min
Histogram: 4 0 0 0 0 0 0 0 0 0
FullNghs: 159600 ave 161880 max 158080 min
Histogram: 1 1 0 0 0 1 0 0 0 1
Total # of neighbors = 638400
Ave neighs/atom = 380
Neighbor list builds = 0
Dangerous builds = 0
unfix 1
thermo 1
fix 3 all nve
run 10
Neighbor list info ...
1 neighbor list requests
update every 20 steps, delay 0 steps, check yes
max neighbors/atom: 2000, page size: 100000
master list distance cutoff = 11.6714
ghost atom cutoff = 11.6714
binsize = 5.8357 -> bins = 5 6 5
Memory usage per processor = 4.2877 Mbytes
Step Temp Press PotEng KinEng TotEng Lx Ly Lz Volume
50 298.34234 2461.5711 -5600.8715 64.748578 -5536.1229 28.61865 34.696255 28.038191 27840.802
51 298.30489 2462.7206 -5600.8634 64.740449 -5536.1229 28.61865 34.696255 28.038191 27840.802
52 298.23415 2465.1307 -5600.848 64.725096 -5536.1229 28.61865 34.696255 28.038191 27840.802
53 298.13012 2468.8007 -5600.8254 64.702519 -5536.1229 28.61865 34.696255 28.038191 27840.802
54 297.99284 2473.7295 -5600.7956 64.672726 -5536.1229 28.61865 34.696255 28.038191 27840.802
55 297.82234 2479.9155 -5600.7586 64.635722 -5536.1229 28.61865 34.696255 28.038191 27840.802
56 297.61866 2487.3567 -5600.7144 64.591518 -5536.1229 28.61865 34.696255 28.038191 27840.802
57 297.38185 2496.0508 -5600.663 64.540124 -5536.1229 28.61865 34.696255 28.038191 27840.802
58 297.11198 2505.9952 -5600.6045 64.481554 -5536.1229 28.61865 34.696255 28.038191 27840.802
59 296.80911 2517.1867 -5600.5387 64.415823 -5536.1229 28.61865 34.696255 28.038191 27840.802
60 296.47332 2529.6221 -5600.4659 64.342947 -5536.1229 28.61865 34.696255 28.038191 27840.802
Loop time of 1.2969 on 4 procs for 10 steps with 1680 atoms
Performance: 0.133 ns/day, 180.125 hours/ns, 7.711 timesteps/s
98.4% CPU use with 4 MPI tasks x no OpenMP threads
MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total
---------------------------------------------------------------
Pair | 1.2955 | 1.2955 | 1.2955 | 0.0 | 99.89
Neigh | 0 | 0 | 0 | 0.0 | 0.00
Comm | 0.00080395 | 0.00095397 | 0.0010161 | 0.3 | 0.07
Output | 0.00023341 | 0.00028217 | 0.00041556 | 0.5 | 0.02
Modify | 7.2241e-05 | 7.993e-05 | 8.8692e-05 | 0.1 | 0.01
Other | | 8.166e-05 | | | 0.01
Nlocal: 420 ave 424 max 416 min
Histogram: 1 0 0 0 0 2 0 0 0 1
Nghost: 4305 ave 4309 max 4301 min
Histogram: 1 0 0 0 0 2 0 0 0 1
Neighs: 0 ave 0 max 0 min
Histogram: 4 0 0 0 0 0 0 0 0 0
FullNghs: 159600 ave 161120 max 158080 min
Histogram: 1 0 0 0 0 2 0 0 0 1
Total # of neighbors = 638400
Ave neighs/atom = 380
Neighbor list builds = 0
Dangerous builds = 0
Total wall time: 0:00:13

View File

@ -0,0 +1,210 @@
LAMMPS (24 Dec 2015)
# Al2O3 crystal, qeq on, minimizes, then calculates elastic constants
variable T_depart equal 300
variable dt equal 0.0002
# =======================================================================
units metal
atom_style charge
dimension 3
boundary p p p
read_data data.Alpha
triclinic box = (0 0 0) to (23.769 24.7015 25.9564) with tilt (0 0 0)
1 by 2 by 2 MPI processor grid
reading atoms ...
1800 atoms
# ^ Orthorombic box of corundum strcture
mass 1 16.00
group Oxy type 1
1080 atoms in group Oxy
compute chargeOxy Oxy property/atom q
compute q_Oxy Oxy reduce ave c_chargeOxy
mass 2 26.98
group Al type 2
720 atoms in group Al
compute chargeAl Al property/atom q
compute q_Al Al reduce ave c_chargeAl
velocity all create ${T_depart} 277387
velocity all create 300 277387
pair_style smtbq
pair_coeff * * ffield.smtbq.Al2O3 O Al
Reading potential file ffield.smtbq.Al2O3 with DATE: 2015-10-22
neighbor 0.5 bin
neigh_modify every 20 delay 0 check yes
timestep ${dt}
timestep 0.0002
thermo_style custom step temp press pe ke etotal c_q_Al c_q_Oxy lx ly lz vol
thermo_modify flush yes
thermo 1
#dump 5 all custom 500 boxAlpha_alumina.lammpstrj id type q x y z
fix 3 all nve
run 10
Neighbor list info ...
1 neighbor list requests
update every 20 steps, delay 0 steps, check yes
max neighbors/atom: 2000, page size: 100000
master list distance cutoff = 11.6714
ghost atom cutoff = 11.6714
binsize = 5.8357 -> bins = 5 5 5
Memory usage per processor = 3.94008 Mbytes
Step Temp Press PotEng KinEng TotEng q_Al q_Oxy Lx Ly Lz Volume
0 300 91921.482 -11494.543 69.7617 -11424.781 2.6095997 -1.7397331 23.769 24.7015 25.9564 15239.78
1 299.96467 91922.303 -11494.535 69.753485 -11424.781 2.6095996 -1.739733 23.769 24.7015 25.9564 15239.78
2 299.75126 91933.246 -11494.485 69.703859 -11424.781 2.6095978 -1.7397318 23.769 24.7015 25.9564 15239.78
3 299.36045 91954.835 -11494.394 69.61298 -11424.781 2.6095941 -1.7397294 23.769 24.7015 25.9564 15239.78
4 298.79335 91986.343 -11494.262 69.481107 -11424.781 2.6095886 -1.7397257 23.769 24.7015 25.9564 15239.78
5 298.05151 92027.62 -11494.09 69.3086 -11424.781 2.6095812 -1.7397208 23.769 24.7015 25.9564 15239.78
6 297.13689 92078.615 -11493.877 69.095915 -11424.781 2.6095721 -1.7397147 23.769 24.7015 25.9564 15239.78
7 296.05187 92139.141 -11493.625 68.843606 -11424.781 2.6095613 -1.7397075 23.769 24.7015 25.9564 15239.78
8 294.79923 92209.15 -11493.334 68.552319 -11424.781 2.6095488 -1.7396992 23.769 24.7015 25.9564 15239.78
9 293.38215 92288.12 -11493.004 68.222793 -11424.781 2.6095347 -1.7396898 23.769 24.7015 25.9564 15239.78
10 291.80421 92376.81 -11492.637 67.855859 -11424.781 2.6095191 -1.7396794 23.769 24.7015 25.9564 15239.78
Loop time of 43.8203 on 4 procs for 10 steps with 1800 atoms
Performance: 0.004 ns/day, 6086.154 hours/ns, 0.228 timesteps/s
99.9% CPU use with 4 MPI tasks x no OpenMP threads
MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total
---------------------------------------------------------------
Pair | 43.817 | 43.817 | 43.817 | 0.0 | 99.99
Neigh | 0 | 0 | 0 | 0.0 | 0.00
Comm | 0.0020506 | 0.0024325 | 0.0026755 | 0.5 | 0.01
Output | 0.00047874 | 0.00057358 | 0.00084901 | 0.7 | 0.00
Modify | 9.8944e-05 | 0.00010616 | 0.00011468 | 0.1 | 0.00
Other | | 0.0001459 | | | 0.00
Nlocal: 450 ave 450 max 450 min
Histogram: 4 0 0 0 0 0 0 0 0 0
Nghost: 6820 ave 6820 max 6820 min
Histogram: 4 0 0 0 0 0 0 0 0 0
Neighs: 0 ave 0 max 0 min
Histogram: 4 0 0 0 0 0 0 0 0 0
FullNghs: 361800 ave 361800 max 361800 min
Histogram: 4 0 0 0 0 0 0 0 0 0
Total # of neighbors = 1447200
Ave neighs/atom = 804
Neighbor list builds = 0
Dangerous builds = 0
unfix 3
thermo 1
fix 1 all box/relax tri 0.0 vmax 0.001
minimize 1.0e-3 1.0e-5 1000 10000
WARNING: Resetting reneighboring criteria during minimization (../min.cpp:168)
Neighbor list info ...
1 neighbor list requests
update every 1 steps, delay 0 steps, check yes
max neighbors/atom: 2000, page size: 100000
master list distance cutoff = 11.6714
ghost atom cutoff = 11.6714
binsize = 5.8357 -> bins = 5 5 5
Memory usage per processor = 5.31508 Mbytes
Step Temp Press PotEng KinEng TotEng q_Al q_Oxy Lx Ly Lz Volume
10 291.80421 92376.81 -11492.637 67.855859 -11424.781 2.6095191 -1.7396794 23.769 24.7015 25.9564 15239.78
11 291.80421 84416.246 -11494.722 67.855859 -11426.866 2.6087748 -1.7391832 23.787835 24.721015 25.982356 15279.17
Loop time of 6.48552 on 4 procs for 1 steps with 1800 atoms
99.9% CPU use with 4 MPI tasks x no OpenMP threads
Minimization stats:
Stopping criterion = energy tolerance
Energy initial, next-to-last, final =
-11492.6369832 -11492.6369832 -11494.7221261
Force two-norm initial, final = 1453.27 1325.26
Force max component initial, final = 968.201 892.249
Final line search alpha, max atom move = 1.03284e-06 0.000921553
Iterations, force evaluations = 1 1
MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total
---------------------------------------------------------------
Pair | 6.4848 | 6.4848 | 6.4848 | 0.0 | 99.99
Neigh | 0 | 0 | 0 | 0.0 | 0.00
Comm | 0.00037098 | 0.00037396 | 0.00037694 | 0.0 | 0.01
Output | 0 | 0 | 0 | 0.0 | 0.00
Modify | 0 | 0 | 0 | 0.0 | 0.00
Other | | 0.0003674 | | | 0.01
Nlocal: 450 ave 450 max 450 min
Histogram: 4 0 0 0 0 0 0 0 0 0
Nghost: 6759 ave 6772 max 6750 min
Histogram: 1 1 0 0 1 0 0 0 0 1
Neighs: 0 ave 0 max 0 min
Histogram: 4 0 0 0 0 0 0 0 0 0
FullNghs: 361140 ave 361150 max 361130 min
Histogram: 1 0 0 1 0 0 0 1 0 1
Total # of neighbors = 1444562
Ave neighs/atom = 802.534
Neighbor list builds = 0
Dangerous builds = 0
unfix 1
thermo 1
fix 3 all nve
run 10
Neighbor list info ...
1 neighbor list requests
update every 20 steps, delay 0 steps, check yes
max neighbors/atom: 2000, page size: 100000
master list distance cutoff = 11.6714
ghost atom cutoff = 11.6714
binsize = 5.8357 -> bins = 5 5 5
Memory usage per processor = 4.19008 Mbytes
Step Temp Press PotEng KinEng TotEng q_Al q_Oxy Lx Ly Lz Volume
11 291.80421 84416.246 -11494.722 67.855859 -11426.866 2.6087748 -1.7391832 23.787835 24.721015 25.982356 15279.17
12 290.08293 84514.767 -11494.322 67.455594 -11426.866 2.6087578 -1.7391718 23.787835 24.721015 25.982356 15279.17
13 288.21041 84622.406 -11493.886 67.020161 -11426.866 2.6087394 -1.7391596 23.787835 24.721015 25.982356 15279.17
14 286.19128 84738.689 -11493.417 66.550634 -11426.866 2.6087199 -1.7391466 23.787835 24.721015 25.982356 15279.17
15 284.03049 84864.242 -11492.914 66.048166 -11426.866 2.6086993 -1.7391329 23.787835 24.721015 25.982356 15279.17
16 281.73331 84998.125 -11492.38 65.513983 -11426.866 2.6086776 -1.7391184 23.787835 24.721015 25.982356 15279.17
17 279.30534 85140.233 -11491.815 64.949384 -11426.866 2.6086551 -1.7391034 23.787835 24.721015 25.982356 15279.17
18 276.75244 85290.405 -11491.221 64.355737 -11426.866 2.6086319 -1.7390879 23.787835 24.721015 25.982356 15279.17
19 274.08079 85448.449 -11490.6 63.734472 -11426.866 2.608608 -1.739072 23.787835 24.721015 25.982356 15279.17
20 271.29678 85614.064 -11489.953 63.087082 -11426.866 2.6085837 -1.7390558 23.787835 24.721015 25.982356 15279.17
21 268.40708 85786.72 -11489.281 62.415114 -11426.865 2.608559 -1.7390393 23.787835 24.721015 25.982356 15279.17
Loop time of 47.1154 on 4 procs for 10 steps with 1800 atoms
Performance: 0.004 ns/day, 6543.802 hours/ns, 0.212 timesteps/s
99.9% CPU use with 4 MPI tasks x no OpenMP threads
MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total
---------------------------------------------------------------
Pair | 47.112 | 47.112 | 47.112 | 0.0 | 99.99
Neigh | 0 | 0 | 0 | 0.0 | 0.00
Comm | 0.0019891 | 0.0023952 | 0.0026395 | 0.5 | 0.01
Output | 0.0004952 | 0.00062233 | 0.00096345 | 0.8 | 0.00
Modify | 0.00010109 | 0.00010723 | 0.00012183 | 0.1 | 0.00
Other | | 0.0001532 | | | 0.00
Nlocal: 450 ave 450 max 450 min
Histogram: 4 0 0 0 0 0 0 0 0 0
Nghost: 6611 ave 6620 max 6600 min
Histogram: 1 0 0 0 0 0 2 0 0 1
Neighs: 0 ave 0 max 0 min
Histogram: 4 0 0 0 0 0 0 0 0 0
FullNghs: 360316 ave 360333 max 360300 min
Histogram: 1 1 0 0 0 0 1 0 0 1
Total # of neighbors = 1441262
Ave neighs/atom = 800.701
Neighbor list builds = 0
Dangerous builds = 0
Total wall time: 0:01:48

View File

@ -0,0 +1,256 @@
LAMMPS (24 Dec 2015)
# Al2O3 crystal, qeq on, minimizes, then calculates elastic constants
variable T_depart equal 300
variable dt equal 0.0002
variable a equal 4.5937
variable c equal 2.9587
variable ca equal ${c}/${a}
variable ca equal 2.9587/${a}
variable ca equal 2.9587/4.5937
variable nx equal 6
variable ny equal 6
variable nz equal 11
variable bx equal ${a}*${nx}
variable bx equal 4.5937*${nx}
variable bx equal 4.5937*6
variable by equal ${a}*${ny}
variable by equal 4.5937*${ny}
variable by equal 4.5937*6
variable bz equal ${c}*${nz}
variable bz equal 2.9587*${nz}
variable bz equal 2.9587*11
# =======================================================================
units metal
atom_style charge
dimension 3
boundary p p p
lattice sc 1.0
Lattice spacing in x,y,z = 1 1 1
region box_vide prism 0 ${bx} 0 ${by} 0 ${bz} 0.0 0.0 0.0
region box_vide prism 0 27.5622 0 ${by} 0 ${bz} 0.0 0.0 0.0
region box_vide prism 0 27.5622 0 27.5622 0 ${bz} 0.0 0.0 0.0
region box_vide prism 0 27.5622 0 27.5622 0 32.5457 0.0 0.0 0.0
create_box 2 box_vide
Created triclinic box = (0 0 0) to (27.5622 27.5622 32.5457) with tilt (0 0 0)
1 by 2 by 2 MPI processor grid
#lattice sc 1.0
#region box_TiO2 block 0 ${bx} 0 ${by} 0 ${bz}
# titanium atoms
lattice custom ${a} origin 0.0 0.0 0.0 orient x 1 0 0 orient y 0 1 0 orient z 0 0 1 a1 1.0 0.0 0.0 a2 0.0 1.0 0.0 a3 0.0 0.0 ${ca} basis 0.0 0.0 0.0 basis 0.5 0.5 0.5
lattice custom 4.5937 origin 0.0 0.0 0.0 orient x 1 0 0 orient y 0 1 0 orient z 0 0 1 a1 1.0 0.0 0.0 a2 0.0 1.0 0.0 a3 0.0 0.0 ${ca} basis 0.0 0.0 0.0 basis 0.5 0.5 0.5
lattice custom 4.5937 origin 0.0 0.0 0.0 orient x 1 0 0 orient y 0 1 0 orient z 0 0 1 a1 1.0 0.0 0.0 a2 0.0 1.0 0.0 a3 0.0 0.0 0.644077758669482 basis 0.0 0.0 0.0 basis 0.5 0.5 0.5
Lattice spacing in x,y,z = 4.5937 4.5937 2.9587
create_atoms 2 region box_vide
Created 792 atoms
# Oxygen atoms
lattice custom ${a} origin 0.0 0.0 0.0 orient x 1 0 0 orient y 0 1 0 orient z 0 0 1 a1 1.0 0.0 0.0 a2 0.0 1.0 0.0 a3 0.0 0.0 ${ca} basis 0.30478 0.30478 0.0 basis 0.69522 0.69522 0.0 basis 0.19522 0.80478 0.5 basis 0.80478 0.19522 0.5
lattice custom 4.5937 origin 0.0 0.0 0.0 orient x 1 0 0 orient y 0 1 0 orient z 0 0 1 a1 1.0 0.0 0.0 a2 0.0 1.0 0.0 a3 0.0 0.0 ${ca} basis 0.30478 0.30478 0.0 basis 0.69522 0.69522 0.0 basis 0.19522 0.80478 0.5 basis 0.80478 0.19522 0.5
lattice custom 4.5937 origin 0.0 0.0 0.0 orient x 1 0 0 orient y 0 1 0 orient z 0 0 1 a1 1.0 0.0 0.0 a2 0.0 1.0 0.0 a3 0.0 0.0 0.644077758669482 basis 0.30478 0.30478 0.0 basis 0.69522 0.69522 0.0 basis 0.19522 0.80478 0.5 basis 0.80478 0.19522 0.5
Lattice spacing in x,y,z = 4.5937 4.5937 2.9587
create_atoms 1 region box_vide
Created 1584 atoms
mass 1 16.00
group Oxy type 1
1584 atoms in group Oxy
compute chargeOxy Oxy property/atom q
compute q_Oxy Oxy reduce ave c_chargeOxy
mass 2 47.867
group Ti type 2
792 atoms in group Ti
compute chargeTi Ti property/atom q
compute q_Ti Ti reduce ave c_chargeTi
velocity all create ${T_depart} 277387
velocity all create 300 277387
pair_style smtbq
pair_coeff * * ffield.smtbq.TiO2 O Ti
Reading potential file ffield.smtbq.TiO2 with DATE: 2015-10-22
neighbor 0.5 bin
neigh_modify every 20 delay 0 check yes
timestep ${dt}
timestep 0.0002
thermo_style custom step temp press pe ke etotal c_q_Ti c_q_Oxy lx ly lz vol
thermo_modify flush yes
thermo 1
#dump 5 all custom 500 boxAlpha_alumina.lammpstrj id type q x y z
fix 3 all nve
run 10
Neighbor list info ...
1 neighbor list requests
update every 20 steps, delay 0 steps, check yes
max neighbors/atom: 2000, page size: 100000
master list distance cutoff = 12.6744
ghost atom cutoff = 12.6744
binsize = 6.3372 -> bins = 5 5 6
Memory usage per processor = 4.60415 Mbytes
Step Temp Press PotEng KinEng TotEng q_Ti q_Oxy Lx Ly Lz Volume
0 300 44365.066 -15815.239 92.097853 -15723.142 2.5521775 -1.2760888 27.5622 27.5622 32.5457 24724.15
1 299.90799 44372.657 -15815.211 92.069608 -15723.142 2.5521771 -1.2760885 27.5622 27.5622 32.5457 24724.15
2 299.64406 44386.645 -15815.13 91.988581 -15723.142 2.5521705 -1.2760852 27.5622 27.5622 32.5457 24724.15
3 299.20863 44406.97 -15814.997 91.854908 -15723.142 2.5521584 -1.2760792 27.5622 27.5622 32.5457 24724.15
4 298.60246 44433.509 -15814.812 91.668818 -15723.143 2.552141 -1.2760705 27.5622 27.5622 32.5457 24724.15
5 297.82659 44466.206 -15814.574 91.430631 -15723.144 2.5521181 -1.276059 27.5622 27.5622 32.5457 24724.15
6 296.88235 44505.016 -15814.285 91.140758 -15723.144 2.5520898 -1.2760449 27.5622 27.5622 32.5457 24724.15
7 295.77139 44549.856 -15813.945 90.7997 -15723.145 2.5520562 -1.2760281 27.5622 27.5622 32.5457 24724.15
8 294.49562 44600.663 -15813.555 90.408048 -15723.147 2.5520173 -1.2760087 27.5622 27.5622 32.5457 24724.15
9 293.05725 44657.353 -15813.114 89.966478 -15723.148 2.5519732 -1.2759866 27.5622 27.5622 32.5457 24724.15
10 291.45876 44719.821 -15812.625 89.475755 -15723.149 2.551924 -1.275962 27.5622 27.5622 32.5457 24724.15
Loop time of 161.071 on 4 procs for 10 steps with 2376 atoms
Performance: 0.001 ns/day, 22370.960 hours/ns, 0.062 timesteps/s
99.8% CPU use with 4 MPI tasks x no OpenMP threads
MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total
---------------------------------------------------------------
Pair | 161.07 | 161.07 | 161.07 | 0.0 |100.00
Neigh | 0 | 0 | 0 | 0.0 | 0.00
Comm | 0.0021119 | 0.0024694 | 0.0026705 | 0.4 | 0.00
Output | 0.00060225 | 0.00070065 | 0.00098777 | 0.6 | 0.00
Modify | 0.00012326 | 0.00013238 | 0.00014782 | 0.1 | 0.00
Other | | 0.0001593 | | | 0.00
Nlocal: 594 ave 630 max 558 min
Histogram: 2 0 0 0 0 0 0 0 0 2
Nghost: 7638 ave 7674 max 7602 min
Histogram: 2 0 0 0 0 0 0 0 0 2
Neighs: 0 ave 0 max 0 min
Histogram: 4 0 0 0 0 0 0 0 0 0
FullNghs: 492624 ave 522720 max 462528 min
Histogram: 2 0 0 0 0 0 0 0 0 2
Total # of neighbors = 1970496
Ave neighs/atom = 829.333
Neighbor list builds = 0
Dangerous builds = 0
unfix 3
#thermo 15
fix 1 all box/relax tri 0.0 vmax 0.001
minimize 1.0e-3 1.0e-5 1000 10000
WARNING: Resetting reneighboring criteria during minimization (../min.cpp:168)
Neighbor list info ...
1 neighbor list requests
update every 1 steps, delay 0 steps, check yes
max neighbors/atom: 2000, page size: 100000
master list distance cutoff = 12.6744
ghost atom cutoff = 12.6744
binsize = 6.3372 -> bins = 5 5 6
Memory usage per processor = 5.97915 Mbytes
Step Temp Press PotEng KinEng TotEng q_Ti q_Oxy Lx Ly Lz Volume
10 291.45876 44719.821 -15812.625 89.475755 -15723.149 2.551924 -1.275962 27.5622 27.5622 32.5457 24724.15
11 291.45876 38966.735 -15814.114 89.475755 -15724.638 2.5514114 -1.2757057 27.582769 27.582778 32.578246 24785.835
Loop time of 22.3212 on 4 procs for 1 steps with 2376 atoms
99.9% CPU use with 4 MPI tasks x no OpenMP threads
Minimization stats:
Stopping criterion = energy tolerance
Energy initial, next-to-last, final =
-15812.6250198 -15812.6250198 -15814.1138993
Force two-norm initial, final = 1103.27 950.466
Force max component initial, final = 758.704 657.184
Final line search alpha, max atom move = 1.31804e-06 0.000866193
Iterations, force evaluations = 1 1
MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total
---------------------------------------------------------------
Pair | 22.32 | 22.32 | 22.32 | 0.0 |100.00
Neigh | 0 | 0 | 0 | 0.0 | 0.00
Comm | 0.00042605 | 0.00043178 | 0.00043893 | 0.0 | 0.00
Output | 0 | 0 | 0 | 0.0 | 0.00
Modify | 0 | 0 | 0 | 0.0 | 0.00
Other | | 0.0004204 | | | 0.00
Nlocal: 594 ave 600 max 587 min
Histogram: 1 0 0 0 1 0 1 0 0 1
Nghost: 7638 ave 7645 max 7632 min
Histogram: 1 0 0 1 0 1 0 0 0 1
Neighs: 0 ave 0 max 0 min
Histogram: 4 0 0 0 0 0 0 0 0 0
FullNghs: 492126 ave 497133 max 486366 min
Histogram: 1 0 0 0 1 0 1 0 0 1
Total # of neighbors = 1968506
Ave neighs/atom = 828.496
Neighbor list builds = 0
Dangerous builds = 0
unfix 1
thermo 1
fix 3 all nve
run 10
Neighbor list info ...
1 neighbor list requests
update every 20 steps, delay 0 steps, check yes
max neighbors/atom: 2000, page size: 100000
master list distance cutoff = 12.6744
ghost atom cutoff = 12.6744
binsize = 6.3372 -> bins = 5 5 6
Memory usage per processor = 4.85415 Mbytes
Step Temp Press PotEng KinEng TotEng q_Ti q_Oxy Lx Ly Lz Volume
11 291.45876 38966.735 -15814.114 89.475755 -15724.638 2.5514114 -1.2757057 27.582769 27.582778 32.578246 24785.835
12 289.71602 39034.428 -15813.58 88.940745 -15724.639 2.5513568 -1.2756784 27.582769 27.582778 32.578246 24785.835
13 287.81994 39107.705 -15813 88.358662 -15724.641 2.5512976 -1.2756488 27.582769 27.582778 32.578246 24785.835
14 285.77378 39186.429 -15812.373 87.730505 -15724.642 2.5512335 -1.2756168 27.582769 27.582778 32.578246 24785.835
15 283.58105 39270.434 -15811.702 87.057353 -15724.644 2.5511647 -1.2755824 27.582769 27.582778 32.578246 24785.835
16 281.24552 39359.544 -15810.986 86.340362 -15724.646 2.5510913 -1.2755457 27.582769 27.582778 32.578246 24785.835
17 278.77119 39453.574 -15810.229 85.580761 -15724.648 2.5510134 -1.2755067 27.582769 27.582778 32.578246 24785.835
18 276.16232 39552.323 -15809.43 84.779855 -15724.65 2.5509311 -1.2754655 27.582769 27.582778 32.578246 24785.835
19 273.42337 39655.586 -15808.592 83.939017 -15724.652 2.5508445 -1.2754223 27.582769 27.582778 32.578246 24785.835
20 270.55904 39763.139 -15807.715 83.05969 -15724.655 2.5507538 -1.2753769 27.582769 27.582778 32.578246 24785.835
21 267.57425 39874.754 -15806.801 82.14338 -15724.658 2.5506591 -1.2753296 27.582769 27.582778 32.578246 24785.835
Loop time of 170.172 on 4 procs for 10 steps with 2376 atoms
Performance: 0.001 ns/day, 23634.979 hours/ns, 0.059 timesteps/s
99.9% CPU use with 4 MPI tasks x no OpenMP threads
MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total
---------------------------------------------------------------
Pair | 170.17 | 170.17 | 170.17 | 0.0 |100.00
Neigh | 0 | 0 | 0 | 0.0 | 0.00
Comm | 0.0025029 | 0.0028628 | 0.0030224 | 0.4 | 0.00
Output | 0.00063372 | 0.00081372 | 0.0011685 | 0.7 | 0.00
Modify | 0.00012541 | 0.0001539 | 0.00021005 | 0.3 | 0.00
Other | | 0.0001655 | | | 0.00
Nlocal: 594 ave 600 max 587 min
Histogram: 1 0 0 0 1 0 1 0 0 1
Nghost: 7638 ave 7645 max 7632 min
Histogram: 1 0 0 1 0 1 0 0 0 1
Neighs: 0 ave 0 max 0 min
Histogram: 4 0 0 0 0 0 0 0 0 0
FullNghs: 490100 ave 495079 max 484284 min
Histogram: 1 0 0 0 1 0 1 0 0 1
Total # of neighbors = 1960398
Ave neighs/atom = 825.083
Neighbor list builds = 0
Dangerous builds = 0
Total wall time: 0:06:39

View File

@ -98,25 +98,16 @@ int BodyNparticle::unpack_border_body(AtomVecBody::Bonus *bonus, double *buf)
------------------------------------------------------------------------- */
void BodyNparticle::data_body(int ibonus, int ninteger, int ndouble,
char **ifile, char **dfile)
int *ifile, double *dfile)
{
AtomVecBody::Bonus *bonus = &avec->bonus[ibonus];
// error in data file if any values are NULL
for (int i = 0; i < ninteger; i++)
if (ifile[i] == NULL)
error->one(FLERR,"Invalid format in Bodies section of data file");
for (int i = 0; i < ndouble; i++)
if (dfile[i] == NULL)
error->one(FLERR,"Invalid format in Bodies section of data file");
// set ninteger, ndouble in bonus and allocate 2 vectors of ints, doubles
if (ninteger != 1)
error->one(FLERR,"Incorrect # of integer values in "
"Bodies section of data file");
int nsub = atoi(ifile[0]);
int nsub = ifile[0];
if (nsub < 1)
error->one(FLERR,"Incorrect integer value in "
"Bodies section of data file");
@ -133,12 +124,12 @@ void BodyNparticle::data_body(int ibonus, int ninteger, int ndouble,
// diagonalize inertia tensor
double tensor[3][3];
tensor[0][0] = atof(dfile[0]);
tensor[1][1] = atof(dfile[1]);
tensor[2][2] = atof(dfile[2]);
tensor[0][1] = tensor[1][0] = atof(dfile[3]);
tensor[0][2] = tensor[2][0] = atof(dfile[4]);
tensor[1][2] = tensor[2][1] = atof(dfile[5]);
tensor[0][0] = dfile[0];
tensor[1][1] = dfile[1];
tensor[2][2] = dfile[2];
tensor[0][1] = tensor[1][0] = dfile[3];
tensor[0][2] = tensor[2][0] = dfile[4];
tensor[1][2] = tensor[2][1] = dfile[5];
double *inertia = bonus->inertia;
double evectors[3][3];
@ -188,9 +179,9 @@ void BodyNparticle::data_body(int ibonus, int ninteger, int ndouble,
int j = 6;
int k = 0;
for (int i = 0; i < nsub; i++) {
delta[0] = atof(dfile[j]);
delta[1] = atof(dfile[j+1]);
delta[2] = atof(dfile[j+2]);
delta[0] = dfile[j];
delta[1] = dfile[j+1];
delta[2] = dfile[j+2];
MathExtra::transpose_matvec(ex_space,ey_space,ez_space,
delta,&bonus->dvalue[k]);
j += 3;
@ -198,6 +189,43 @@ void BodyNparticle::data_body(int ibonus, int ninteger, int ndouble,
}
}
/* ----------------------------------------------------------------------
return radius of body particle defined by ifile/dfile params
params are ordered as in data file
called by Molecule class which needs single body size
------------------------------------------------------------------------- */
double BodyNparticle::radius_body(int ninteger, int ndouble,
int *ifile, double *dfile)
{
int nsub = ifile[0];
if (nsub < 1)
error->one(FLERR,"Incorrect integer value in "
"Bodies section of data file");
if (ndouble != 6 + 3*nsub)
error->one(FLERR,"Incorrect # of floating-point values in "
"Bodies section of data file");
// sub-particle coords are relative to body center at (0,0,0)
// offset = 6 for sub-particle coords
double onerad;
double maxrad = 0.0;
double delta[3];
int offset = 6;
for (int i = 0; i < nsub; i++) {
delta[0] = dfile[offset];
delta[1] = dfile[offset+1];
delta[2] = dfile[offset+2];
offset += 3;
onerad = MathExtra::len3(delta);
maxrad = MAX(maxrad,onerad);
}
return maxrad;
}
/* ---------------------------------------------------------------------- */
int BodyNparticle::noutcol()

View File

@ -34,7 +34,8 @@ class BodyNparticle : public Body {
int pack_border_body(struct AtomVecBody::Bonus *, double *);
int unpack_border_body(struct AtomVecBody::Bonus *, double *);
void data_body(int, int, int, char **, char **);
void data_body(int, int, int, int *, double *);
double radius_body(int, int, int *, double *);
int noutrow(int);
int noutcol();

View File

@ -74,13 +74,14 @@ ComputeTempCS::ComputeTempCS(LAMMPS *lmp, int narg, char **arg) :
strcpy(id_fix,id);
strcat(id_fix,"_COMPUTE_STORE");
char **newarg = new char*[5];
char **newarg = new char*[6];
newarg[0] = id_fix;
newarg[1] = group->names[igroup];
newarg[2] = (char *) "STORE";
newarg[3] = (char *) "0";
newarg[4] = (char *) "1";
modify->add_fix(5,newarg);
newarg[3] = (char *) "peratom";
newarg[4] = (char *) "0";
newarg[5] = (char *) "1";
modify->add_fix(6,newarg);
fix = (FixStore *) modify->fix[modify->nfix-1];
delete [] newarg;

View File

@ -378,6 +378,13 @@ void FixPour::pre_exchange()
if (next_reneighbor != update->ntimestep) return;
// clear ghost count and any ghost bonus data internal to AtomVec
// same logic as beginning of Comm::exchange()
// do it now b/c inserting atoms will overwrite ghost atoms
atom->nghost = 0;
atom->avec->clear_bonus();
// find current max atom and molecule IDs if necessary
if (!idnext) find_maxid();
@ -638,7 +645,11 @@ void FixPour::pre_exchange()
radtmp = newcoord[3];
atom->radius[n] = radtmp;
atom->rmass[n] = 4.0*MY_PI/3.0 * radtmp*radtmp*radtmp * denstmp;
} else atom->add_molecule_atom(onemols[imol],m,n,maxtag_all);
} else {
onemols[imol]->quat_external = quat;
atom->add_molecule_atom(onemols[imol],m,n,maxtag_all);
}
modify->create_attribute(n);
}
}
@ -667,7 +678,8 @@ void FixPour::pre_exchange()
// reset global natoms,nbonds,etc
// increment maxtag_all and maxmol_all if necessary
// if global map exists, reset it now instead of waiting for comm
// since adding atoms messes up ghosts
// since other pre-exchange fixes may use it
// invoke map_init() b/c atom count has grown
if (ninserted_atoms) {
atom->natoms += ninserted_atoms;
@ -682,7 +694,6 @@ void FixPour::pre_exchange()
if (maxtag_all >= MAXTAGINT)
error->all(FLERR,"New atom IDs exceed maximum allowed ID");
if (atom->map_style) {
atom->nghost = 0;
atom->map_init();
atom->map_set();
}
@ -1008,23 +1019,24 @@ void *FixPour::extract(const char *str, int &itype)
} else {
// find a molecule in template with matching type
// loop over onemols molecules
// skip a molecule with no atoms as large as itype
oneradius = 0.0;
for (int i = 0; i < nmol; i++) {
if (itype-ntype > onemols[i]->ntypes) continue;
if (itype > ntype+onemols[i]->ntypes) continue;
double *radius = onemols[i]->radius;
int *type = onemols[i]->type;
int natoms = onemols[i]->natoms;
// check radii of matching types in Molecule
// check radii of atoms in Molecule with matching types
// default to 0.5, if radii not defined in Molecule
// same as atom->avec->create_atom(), invoked in pre_exchange()
oneradius = 0.0;
for (int i = 0; i < natoms; i++)
if (type[i] == itype-ntype) {
if (type[i]+ntype == itype) {
if (radius) oneradius = MAX(oneradius,radius[i]);
else oneradius = 0.5;
else oneradius = MAX(oneradius,0.5);
}
}
}

View File

@ -301,6 +301,13 @@ void FixDeposit::pre_exchange()
if (next_reneighbor != update->ntimestep) return;
// clear ghost count and any ghost bonus data internal to AtomVec
// same logic as beginning of Comm::exchange()
// do it now b/c inserting atoms will overwrite ghost atoms
atom->nghost = 0;
atom->avec->clear_bonus();
// compute current offset = bottom of insertion volume
double offset = 0.0;
@ -524,8 +531,10 @@ void FixDeposit::pre_exchange()
atom->v[n][0] = vnew[0];
atom->v[n][1] = vnew[1];
atom->v[n][2] = vnew[2];
if (mode == MOLECULE)
if (mode == MOLECULE) {
onemols[imol]->quat_external = quat;
atom->add_molecule_atom(onemols[imol],m,n,maxtag_all);
}
modify->create_attribute(n);
}
}
@ -559,7 +568,8 @@ void FixDeposit::pre_exchange()
// reset global natoms,nbonds,etc
// increment maxtag_all and maxmol_all if necessary
// if global map exists, reset it now instead of waiting for comm
// since adding atoms messes up ghosts
// since other pre-exchange fixes may use it
// invoke map_init() b/c atom count has grown
if (success) {
atom->natoms += natom;
@ -576,7 +586,6 @@ void FixDeposit::pre_exchange()
error->all(FLERR,"New atom IDs exceed maximum allowed ID");
if (mode == MOLECULE && atom->molecule_flag) maxmol_all++;
if (atom->map_style) {
atom->nghost = 0;
atom->map_init();
atom->map_set();
}
@ -811,23 +820,24 @@ void *FixDeposit::extract(const char *str, int &itype)
} else {
// find a molecule in template with matching type
// loop over onemols molecules
// skip a molecule with no atoms as large as itype
oneradius = 0.0;
for (int i = 0; i < nmol; i++) {
if (itype-ntype > onemols[i]->ntypes) continue;
if (itype > ntype+onemols[i]->ntypes) continue;
double *radius = onemols[i]->radius;
int *type = onemols[i]->type;
int natoms = onemols[i]->natoms;
// check radii of matching types in Molecule
// check radii of atoms in Molecule with matching types
// default to 0.5, if radii not defined in Molecule
// same as atom->avec->create_atom(), invoked in pre_exchange()
oneradius = 0.0;
for (int i = 0; i < natoms; i++)
if (type[i] == itype-ntype) {
if (type[i]+ntype == itype) {
if (radius) oneradius = MAX(oneradius,radius[i]);
else oneradius = 0.5;
else oneradius = MAX(oneradius,0.5);
}
}
}

View File

@ -154,13 +154,14 @@ void TAD::command(int narg, char **arg)
// create FixStore object to store revert state
narg2 = 5;
narg2 = 6;
args = new char*[narg2];
args[0] = (char *) "tad_revert";
args[1] = (char *) "all";
args[2] = (char *) "STORE";
args[3] = (char *) "0";
args[4] = (char *) "7";
args[3] = (char *) "peratom";
args[4] = (char *) "0";
args[5] = (char *) "7";
modify->add_fix(narg2,args);
fix_revert = (FixStore *) modify->fix[modify->nfix-1];
delete [] args;

View File

@ -224,11 +224,12 @@ void FixAdaptFEP::post_constructor()
id_fix_diam = NULL;
id_fix_chg = NULL;
char **newarg = new char*[5];
char **newarg = new char*[6];
newarg[1] = group->names[igroup];
newarg[2] = (char *) "STORE";
newarg[3] = (char *) "1";
newarg[3] = (char *) "peratom";
newarg[4] = (char *) "1";
newarg[5] = (char *) "1";
if (diamflag) {
int n = strlen(id) + strlen("_FIX_STORE_DIAM") + 1;
@ -236,7 +237,7 @@ void FixAdaptFEP::post_constructor()
strcpy(id_fix_diam,id);
strcat(id_fix_diam,"_FIX_STORE_DIAM");
newarg[0] = id_fix_diam;
modify->add_fix(5,newarg);
modify->add_fix(6,newarg);
fix_diam = (FixStore *) modify->fix[modify->nfix-1];
if (fix_diam->restart_reset) fix_diam->restart_reset = 0;
@ -259,7 +260,7 @@ void FixAdaptFEP::post_constructor()
strcpy(id_fix_chg,id);
strcat(id_fix_chg,"_FIX_STORE_CHG");
newarg[0] = id_fix_chg;
modify->add_fix(5,newarg);
modify->add_fix(6,newarg);
fix_chg = (FixStore *) modify->fix[modify->nfix-1];
if (fix_chg->restart_reset) fix_chg->restart_reset = 0;

View File

@ -828,9 +828,7 @@ void PairSMTBQ::read_file(char *file)
}
//A adapter au STO
ncov = params[0].sto * params[0].n0;
for (i=1; i<num_atom_types; i++)
ncov = min(ncov,params[i].sto * params[i].n0);
ncov = min((params[0].sto)*(params[0].n0),(params[1].sto)*(params[1].n0));
if (verbose) printf (" Parametre ncov = %f\n",ncov);
if (verbose) printf (" ********************************************* \n");
@ -914,7 +912,7 @@ void PairSMTBQ::compute(int eflag, int vflag)
// Charges Communication
// ----------------------
forward(q) ; reverse(q);
forward(q) ; // reverse(q);
memory->create(nmol,nteam+1,"pair:nmol");
memory->create(tmp,nteam+1,7,"pair:tmp");
@ -1271,7 +1269,9 @@ void PairSMTBQ::tabqeq()
double aCoeff,bCoeff,rcoupe,nang;
int n = atom->ntypes;
int nmax = atom->nmax;
int nlocal = atom->nlocal;
int nghost = atom->nghost;
nmax = atom->nmax;
verbose = 1;
verbose = 0;
@ -1284,6 +1284,7 @@ void PairSMTBQ::tabqeq()
if (verbose) printf ("kmax %d, ds %f, nmax %d\n",kmax,ds,nmax);
if (verbose) printf ("nlocal = %d, nghost = %d\n",nlocal,nghost);
if (verbose) printf ("nntypes %d, kmax %d, rc %f, n %d\n",nntype,kmax,rc,n);
// allocate arrays
@ -2396,9 +2397,9 @@ void PairSMTBQ::QForce_charge(int loop)
// to calculate the N-body forces
// ================================================
forward (sbcov) ; reverse (sbcov);
forward (coord) ; reverse (coord);
forward (sbmet) ; reverse (sbmet);
forward (sbcov) ; // reverse (sbcov);
forward (coord) ; // reverse (coord);
forward (sbmet) ; // reverse (sbmet);
if (nteam == 0) return; //no oxide
@ -2655,7 +2656,7 @@ void PairSMTBQ::Charge()
// Communication des charges
// ---------------------------
forward(q) ; reverse(q);
forward(q) ; // reverse(q);
// Calcul des potentiel
@ -2735,7 +2736,7 @@ void PairSMTBQ::Charge()
// =======================================
// Charge Communication.
// =======================================
forward(q); reverse(q);
forward(q); // reverse(q);
//:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
// ==========================================
@ -3124,7 +3125,7 @@ void PairSMTBQ::groupQEqAllParallel_QEq()
// Groups communication
// OOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOO
for (i = 0; i < nproc; i++) {
forward_int(flag_gp[i]); reverse_int(flag_gp[i]);
forward_int(flag_gp[i]); //reverse_int(flag_gp[i]);
}
// ---
@ -3214,7 +3215,7 @@ void PairSMTBQ::groupQEqAllParallel_QEq()
Group Communication between proc for unification
================================================== */
for (i = 0; i < nproc; i++) {
forward_int(flag_gp[i]); reverse_int(flag_gp[i]);
forward_int(flag_gp[i]);// reverse_int(flag_gp[i]);
}
// =============== End of COMM =================
@ -3404,7 +3405,7 @@ int PairSMTBQ::pack_forward_comm(int n, int *list, double *buf, int pbc_flag, in
buf[m++] = tab_comm[j];
// if (j < 3) printf ("[%d] %d pfc %d %d buf_send = %f \n",me,n,i,m-1,buf[m-1]);
}
return 1;
return m;
}
/* ---------------------------------------------------------------------- */
@ -3433,7 +3434,7 @@ int PairSMTBQ::pack_reverse_comm(int n, int first, double *buf)
buf[m++] = tab_comm[i];
// if (i<first+3) printf ("[%d] prc %d %d buf_send = %f \n",me,i,m-1,buf[m-1]);
}
return 1;
return m;
}
/* ---------------------------------------------------------------------- */

View File

@ -1318,45 +1318,51 @@ void Atom::data_bonus(int n, char *buf, AtomVec *avec_bonus, tagint id_offset)
void Atom::data_bodies(int n, char *buf, AtomVecBody *avec_body,
tagint id_offset)
{
int j,m,tagdata,ninteger,ndouble;
int j,m,nvalues,tagdata,ninteger,ndouble;
int maxint = 0;
int maxdouble = 0;
char **ivalues = NULL;
char **dvalues = NULL;
int *ivalues = NULL;
double *dvalues = NULL;
// loop over lines of body data
// tokenize the lines into ivalues and dvalues
// if I own atom tag, unpack its values
// if I own atom tag, tokenize lines into ivalues/dvalues, call data_body()
// else skip values
for (int i = 0; i < n; i++) {
if (i == 0) tagdata = ATOTAGINT(strtok(buf," \t\n\r\f")) + id_offset;
else tagdata = ATOTAGINT(strtok(NULL," \t\n\r\f")) + id_offset;
ninteger = atoi(strtok(NULL," \t\n\r\f"));
ndouble = atoi(strtok(NULL," \t\n\r\f"));
if (tagdata <= 0 || tagdata > map_tag_max)
error->one(FLERR,"Invalid atom ID in Bodies section of data file");
ninteger = force->inumeric(FLERR,strtok(NULL," \t\n\r\f"));
ndouble = force->inumeric(FLERR,strtok(NULL," \t\n\r\f"));
if ((m = map(tagdata)) >= 0) {
if (ninteger > maxint) {
delete [] ivalues;
maxint = ninteger;
ivalues = new char*[maxint];
ivalues = new int[maxint];
}
if (ndouble > maxdouble) {
delete [] dvalues;
maxdouble = ndouble;
dvalues = new char*[maxdouble];
dvalues = new double[maxdouble];
}
for (j = 0; j < ninteger; j++)
ivalues[j] = strtok(NULL," \t\n\r\f");
ivalues[j] = force->inumeric(FLERR,strtok(NULL," \t\n\r\f"));
for (j = 0; j < ndouble; j++)
dvalues[j] = strtok(NULL," \t\n\r\f");
dvalues[j] = force->numeric(FLERR,strtok(NULL," \t\n\r\f"));
if (tagdata <= 0 || tagdata > map_tag_max)
error->one(FLERR,"Invalid atom ID in Bodies section of data file");
if ((m = map(tagdata)) >= 0)
avec_body->data_body(m,ninteger,ndouble,ivalues,dvalues);
} else {
nvalues = ninteger + ndouble; // number of values to skip
for (j = 0; j < nvalues; j++)
strtok(NULL," \t\n\r\f");
}
}
delete [] ivalues;
@ -1588,6 +1594,12 @@ void Atom::add_molecule_atom(Molecule *onemol, int iatom,
else if (rmass_flag)
rmass[ilocal] = 4.0*MY_PI/3.0 *
radius[ilocal]*radius[ilocal]*radius[ilocal];
if (onemol->bodyflag) {
body[ilocal] = 0; // as if a body read from data file
onemol->avec_body->data_body(ilocal,onemol->nibody,onemol->ndbody,
onemol->ibodyparams,onemol->dbodyparams);
onemol->avec_body->set_quat(ilocal,onemol->quat_external);
}
if (molecular != 1) return;

View File

@ -26,6 +26,10 @@
#include "memory.h"
#include "error.h"
// debug
#include "update.h"
using namespace LAMMPS_NS;
/* ---------------------------------------------------------------------- */
@ -48,6 +52,7 @@ AtomVecBody::AtomVecBody(LAMMPS *lmp) : AtomVec(lmp)
atom->body_flag = 1;
atom->rmass_flag = 1;
atom->angmom_flag = atom->torque_flag = 1;
atom->radius_flag = 1;
nlocal_bonus = nghost_bonus = nmax_bonus = 0;
bonus = NULL;
@ -103,7 +108,7 @@ void AtomVecBody::process_args(int narg, char **arg)
// bptr values = max number of additional ivalues/dvalues from Body class
size_forward = 7 + bptr->size_forward;
size_border = 16 + bptr->size_border;
size_border = 18 + bptr->size_border;
}
/* ----------------------------------------------------------------------
@ -128,6 +133,7 @@ void AtomVecBody::grow(int n)
v = memory->grow(atom->v,nmax,3,"atom:v");
f = memory->grow(atom->f,nmax*comm->nthreads,3,"atom:f");
radius = memory->grow(atom->radius,nmax,"atom:radius");
rmass = memory->grow(atom->rmass,nmax,"atom:rmass");
angmom = memory->grow(atom->angmom,nmax,3,"atom:angmom");
torque = memory->grow(atom->torque,nmax*comm->nthreads,3,"atom:torque");
@ -147,7 +153,8 @@ void AtomVecBody::grow_reset()
tag = atom->tag; type = atom->type;
mask = atom->mask; image = atom->image;
x = atom->x; v = atom->v; f = atom->f;
rmass = atom->rmass; angmom = atom->angmom; torque = atom->torque;
radius = atom->radius; rmass = atom->rmass;
angmom = atom->angmom; torque = atom->torque;
body = atom->body;
}
@ -183,6 +190,7 @@ void AtomVecBody::copy(int i, int j, int delflag)
v[j][1] = v[i][1];
v[j][2] = v[i][2];
radius[j] = radius[i];
rmass[j] = rmass[i];
angmom[j][0] = angmom[i][0];
angmom[j][1] = angmom[i][1];
@ -571,6 +579,8 @@ int AtomVecBody::pack_border(int n, int *list, double *buf,
buf[m++] = ubuf(tag[j]).d;
buf[m++] = ubuf(type[j]).d;
buf[m++] = ubuf(mask[j]).d;
buf[m++] = radius[j];
buf[m++] = rmass[j];
if (body[j] < 0) buf[m++] = ubuf(0).d;
else {
buf[m++] = ubuf(1).d;
@ -606,6 +616,8 @@ int AtomVecBody::pack_border(int n, int *list, double *buf,
buf[m++] = ubuf(tag[j]).d;
buf[m++] = ubuf(type[j]).d;
buf[m++] = ubuf(mask[j]).d;
buf[m++] = radius[j];
buf[m++] = rmass[j];
if (body[j] < 0) buf[m++] = ubuf(0).d;
else {
buf[m++] = ubuf(1).d;
@ -651,6 +663,8 @@ int AtomVecBody::pack_border_vel(int n, int *list, double *buf,
buf[m++] = ubuf(tag[j]).d;
buf[m++] = ubuf(type[j]).d;
buf[m++] = ubuf(mask[j]).d;
buf[m++] = radius[j];
buf[m++] = rmass[j];
if (body[j] < 0) buf[m++] = ubuf(0).d;
else {
buf[m++] = ubuf(1).d;
@ -693,6 +707,8 @@ int AtomVecBody::pack_border_vel(int n, int *list, double *buf,
buf[m++] = ubuf(tag[j]).d;
buf[m++] = ubuf(type[j]).d;
buf[m++] = ubuf(mask[j]).d;
buf[m++] = radius[j];
buf[m++] = rmass[j];
if (body[j] < 0) buf[m++] = ubuf(0).d;
else {
buf[m++] = ubuf(1).d;
@ -777,6 +793,8 @@ int AtomVecBody::pack_border_hybrid(int n, int *list, double *buf)
m = 0;
for (i = 0; i < n; i++) {
j = list[i];
buf[m++] = radius[j];
buf[m++] = rmass[j];
if (body[j] < 0) buf[m++] = ubuf(0).d;
else {
buf[m++] = ubuf(1).d;
@ -814,6 +832,8 @@ void AtomVecBody::unpack_border(int n, int first, double *buf)
tag[i] = (tagint) ubuf(buf[m++]).i;
type[i] = (int) ubuf(buf[m++]).i;
mask[i] = (int) ubuf(buf[m++]).i;
radius[i] = buf[m++];
rmass[i] = buf[m++];
body[i] = (int) ubuf(buf[m++]).i;
if (body[i] == 0) body[i] = -1;
else {
@ -862,6 +882,8 @@ void AtomVecBody::unpack_border_vel(int n, int first, double *buf)
tag[i] = (tagint) ubuf(buf[m++]).i;
type[i] = (int) ubuf(buf[m++]).i;
mask[i] = (int) ubuf(buf[m++]).i;
radius[i] = buf[m++];
rmass[i] = buf[m++];
body[i] = (int) ubuf(buf[m++]).i;
if (body[i] == 0) body[i] = -1;
else {
@ -909,6 +931,8 @@ int AtomVecBody::unpack_border_hybrid(int n, int first, double *buf)
m = 0;
last = first + n;
for (i = first; i < last; i++) {
radius[i] = buf[m++];
rmass[i] = buf[m++];
body[i] = (int) ubuf(buf[m++]).i;
if (body[i] == 0) body[i] = -1;
else {
@ -954,7 +978,7 @@ int AtomVecBody::pack_exchange(int i, double *buf)
buf[m++] = ubuf(type[i]).d;
buf[m++] = ubuf(mask[i]).d;
buf[m++] = ubuf(image[i]).d;
buf[m++] = radius[i];
buf[m++] = rmass[i];
buf[m++] = angmom[i][0];
buf[m++] = angmom[i][1];
@ -1008,7 +1032,7 @@ int AtomVecBody::unpack_exchange(double *buf)
type[nlocal] = (int) ubuf(buf[m++]).i;
mask[nlocal] = (int) ubuf(buf[m++]).i;
image[nlocal] = (imageint) ubuf(buf[m++]).i;
radius[nlocal] = buf[m++];
rmass[nlocal] = buf[m++];
angmom[nlocal][0] = buf[m++];
angmom[nlocal][1] = buf[m++];
@ -1067,11 +1091,11 @@ int AtomVecBody::size_restart()
int nlocal = atom->nlocal;
for (i = 0; i < nlocal; i++)
if (body[i] >= 0) {
n += 25;
n += 26;
if (intdoubleratio == 1) n += bonus[body[i]].ninteger;
else n += (bonus[body[i]].ninteger+1)/2;
n += bonus[body[i]].ndouble;
} else n += 16;
} else n += 17;
if (atom->nextra_restart)
for (int iextra = 0; iextra < atom->nextra_restart; iextra++)
@ -1101,6 +1125,7 @@ int AtomVecBody::pack_restart(int i, double *buf)
buf[m++] = v[i][1];
buf[m++] = v[i][2];
buf[m++] = radius[i];
buf[m++] = rmass[i];
buf[m++] = angmom[i][0];
buf[m++] = angmom[i][1];
@ -1161,6 +1186,7 @@ int AtomVecBody::unpack_restart(double *buf)
v[nlocal][1] = buf[m++];
v[nlocal][2] = buf[m++];
radius[nlocal] = buf[m++];
rmass[nlocal] = buf[m++];
angmom[nlocal][0] = buf[m++];
angmom[nlocal][1] = buf[m++];
@ -1228,6 +1254,7 @@ void AtomVecBody::create_atom(int itype, double *coord)
v[nlocal][1] = 0.0;
v[nlocal][2] = 0.0;
radius[nlocal] = 0.5;
rmass[nlocal] = 1.0;
angmom[nlocal][0] = 0.0;
angmom[nlocal][1] = 0.0;
@ -1274,6 +1301,7 @@ void AtomVecBody::data_atom(double *coord, imageint imagetmp, char **values)
angmom[nlocal][0] = 0.0;
angmom[nlocal][1] = 0.0;
angmom[nlocal][2] = 0.0;
radius[nlocal] = 0.5;
atom->nlocal++;
}
@ -1302,12 +1330,12 @@ int AtomVecBody::data_atom_hybrid(int nlocal, char **values)
------------------------------------------------------------------------- */
void AtomVecBody::data_body(int m, int ninteger, int ndouble,
char **ivalues, char **dvalues)
int *ivalues, double *dvalues)
{
if (body[m]) error->one(FLERR,"Assigning body parameters to non-body atom");
if (nlocal_bonus == nmax_bonus) grow_bonus();
bptr->data_body(nlocal_bonus,ninteger,ndouble,ivalues,dvalues);
bonus[nlocal_bonus].ilocal = m;
bptr->data_body(nlocal_bonus,ninteger,ndouble,ivalues,dvalues);
body[m] = nlocal_bonus++;
}
@ -1448,6 +1476,29 @@ int AtomVecBody::write_vel_hybrid(FILE *fp, double *buf)
return 3;
}
/* ----------------------------------------------------------------------
body computes its size based on ivalues/dvalues and returns it
------------------------------------------------------------------------- */
double AtomVecBody::radius_body(int ninteger, int ndouble,
int *ivalues, double *dvalues)
{
return bptr->radius_body(ninteger,ndouble,ivalues,dvalues);
}
/* ----------------------------------------------------------------------
reset quat orientation for atom M to quat_external
called by Atom:add_molecule_atom()
------------------------------------------------------------------------- */
void AtomVecBody::set_quat(int m, double *quat_external)
{
if (body[m] < 0) error->one(FLERR,"Assigning quat to non-body atom");
double *quat = bonus[body[m]].quat;
quat[0] = quat_external[0]; quat[1] = quat_external[1];
quat[2] = quat_external[2]; quat[3] = quat_external[3];
}
/* ----------------------------------------------------------------------
return # of bytes of allocated memory
------------------------------------------------------------------------- */
@ -1464,6 +1515,7 @@ bigint AtomVecBody::memory_usage()
if (atom->memcheck("v")) bytes += memory->usage(v,nmax,3);
if (atom->memcheck("f")) bytes += memory->usage(f,nmax*comm->nthreads,3);
if (atom->memcheck("radius")) bytes += memory->usage(radius,nmax);
if (atom->memcheck("rmass")) bytes += memory->usage(rmass,nmax);
if (atom->memcheck("angmom")) bytes += memory->usage(angmom,nmax,3);
if (atom->memcheck("torque")) bytes +=

View File

@ -85,13 +85,19 @@ class AtomVecBody : public AtomVec {
// manipulate Bonus data structure for extra atom info
void clear_bonus();
void data_body(int, int, int, char **, char **);
void data_body(int, int, int, int *, double *);
// methods used by other classes to query/set body info
double radius_body(int, int, int *, double *);
void set_quat(int, double *);
private:
tagint *tag;
int *type,*mask;
imageint *image;
double **x,**v,**f;
double *radius;
double *rmass;
double **angmom,**torque;
int *body;

View File

@ -41,11 +41,13 @@ class Body : protected Pointers {
virtual int unpack_border_body(struct AtomVecBody::Bonus *,
double *) {return 0;}
virtual void data_body(int, int, int, char **, char **) = 0;
virtual void data_body(int, int, int, int*, double *) = 0;
virtual int noutrow(int) = 0;
virtual int noutcol() = 0;
virtual void output(int, int, double *) = 0;
virtual int image(int, double, double, int *&, double **&) = 0;
virtual double radius_body(int, int, int *, double *) {return 0.0;}
};
}

View File

@ -191,6 +191,7 @@ ComputeChunkAtom::ComputeChunkAtom(LAMMPS *lmp, int narg, char **arg) :
maxflag[1] = UPPER;
maxflag[2] = UPPER;
scaleflag = LATTICE;
pbcflag = 0;
while (iarg < narg) {
if (strcmp(arg[iarg],"region") == 0) {
@ -269,6 +270,12 @@ ComputeChunkAtom::ComputeChunkAtom(LAMMPS *lmp, int narg, char **arg) :
else if (strcmp(arg[iarg+1],"reduced") == 0) scaleflag = REDUCED;
else error->all(FLERR,"Illegal compute chunk/atom command");
iarg += 2;
} else if (strcmp(arg[iarg],"pbc") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal compute chunk/atom command");
if (strcmp(arg[iarg+1],"no") == 0) pbcflag = 0;
else if (strcmp(arg[iarg+1],"yes") == 0) pbcflag = 1;
else error->all(FLERR,"Illegal compute chunk/atom command");
iarg += 2;
} else error->all(FLERR,"Illegal compute chunk/atom command");
}
@ -572,13 +579,14 @@ void ComputeChunkAtom::init()
strcpy(id_fix,id);
strcat(id_fix,"_COMPUTE_STORE");
char **newarg = new char*[5];
char **newarg = new char*[6];
newarg[0] = id_fix;
newarg[1] = group->names[igroup];
newarg[2] = (char *) "STORE";
newarg[3] = (char *) "1";
newarg[3] = (char *) "peratom";
newarg[4] = (char *) "1";
modify->add_fix(5,newarg);
newarg[5] = (char *) "1";
modify->add_fix(6,newarg);
fixstore = (FixStore *) modify->fix[modify->nfix-1];
delete [] newarg;
}
@ -1265,7 +1273,7 @@ int ComputeChunkAtom::setup_xyz_bins()
int ComputeChunkAtom::setup_sphere_bins()
{
// convert sorigin_user to sorigin
// sorigin is always in box units, for orthogonal or triclinic domains
// sorigin,srad are always in box units, for orthogonal or triclinic domains
// lamda2x works for either orthogonal or triclinic
if (scaleflag == REDUCED) {
@ -1280,6 +1288,23 @@ int ComputeChunkAtom::setup_sphere_bins()
sradmax = sradmax_user;
}
// if pbcflag set, sradmax must be < 1/2 box in any periodic dim
// treat orthongonal and triclinic the same
// check every time bins are created
if (pbcflag) {
double *prd_half = domain->prd_half;
int *periodicity = domain->periodicity;
int flag = 0;
if (periodicity[0] && sradmax > prd_half[0]) flag = 1;
if (periodicity[1] && sradmax > prd_half[1]) flag = 1;
if (domain->dimension == 3 &&
periodicity[2] && sradmax > prd_half[2]) flag = 1;
if (flag)
error->all(FLERR,"Compute chunk/atom bin/sphere radius "
"is too large for periodic box");
}
sinvrad = nsbin / (sradmax-sradmin);
// allocate and set bin coordinates
@ -1328,6 +1353,21 @@ int ComputeChunkAtom::setup_cylinder_bins()
cradmax = cradmax_user;
}
// if pbcflag set, sradmax must be < 1/2 box in any periodic non-axis dim
// treat orthongonal and triclinic the same
// check every time bins are created
if (pbcflag) {
double *prd_half = domain->prd_half;
int *periodicity = domain->periodicity;
int flag = 0;
if (periodicity[cdim1] && sradmax > prd_half[cdim1]) flag = 1;
if (periodicity[cdim2] && sradmax > prd_half[cdim2]) flag = 1;
if (flag)
error->all(FLERR,"Compute chunk/atom bin/cylinder radius "
"is too large for periodic box");
}
cinvrad = ncbin / (cradmax-cradmin);
// allocate and set radial bin coordinates
@ -1743,6 +1783,7 @@ void ComputeChunkAtom::atom2binsphere()
double *boxlo = domain->boxlo;
double *boxhi = domain->boxhi;
double *prd = domain->prd;
double *prd_half = domain->prd_half;
int *periodicity = domain->periodicity;
// remap each atom's relevant coords back into box via PBC if necessary
@ -1773,6 +1814,33 @@ void ComputeChunkAtom::atom2binsphere()
dx = xremap - sorigin[0];
dy = yremap - sorigin[1];
dz = zremap - sorigin[2];
// if requested, apply PBC to distance from sphere center
// treat orthogonal and triclinic the same
// with dx,dy,dz = lengths independent of each other
// so do not use domain->minimum_image() which couples for triclinic
if (pbcflag) {
if (periodicity[0]) {
if (fabs(dx) > prd_half[0]) {
if (dx < 0.0) dx += prd[0];
else dx -= prd[0];
}
}
if (periodicity[1]) {
if (fabs(dy) > prd_half[1]) {
if (dy < 0.0) dy += prd[1];
else dy -= prd[1];
}
}
if (periodicity[2]) {
if (fabs(dz) > prd_half[2]) {
if (dz < 0.0) dz += prd[2];
else dz -= prd[2];
}
}
}
r = sqrt(dx*dx + dy*dy + dz*dz);
ibin = static_cast<int> ((r - sradmin) * sinvrad);
@ -1792,7 +1860,6 @@ void ComputeChunkAtom::atom2binsphere()
/* ----------------------------------------------------------------------
assign each atom to a cylindrical bin
------------------------------------------------------------------------- */
void ComputeChunkAtom::atom2bincylinder()
@ -1812,6 +1879,7 @@ void ComputeChunkAtom::atom2bincylinder()
double *boxlo = domain->boxlo;
double *boxhi = domain->boxhi;
double *prd = domain->prd;
double *prd_half = domain->prd_half;
int *periodicity = domain->periodicity;
// remap each atom's relevant coords back into box via PBC if necessary
@ -1837,6 +1905,26 @@ void ComputeChunkAtom::atom2bincylinder()
d1 = remap1 - corigin[cdim1];
d2 = remap2 - corigin[cdim2];
// if requested, apply PBC to distance from cylinder axis
// treat orthogonal and triclinic the same
// with d1,d2 = lengths independent of each other
if (pbcflag) {
if (periodicity[cdim1]) {
if (fabs(d1) > prd_half[cdim1]) {
if (d1 < 0.0) d1 += prd[cdim1];
else d1 -= prd[cdim1];
}
}
if (periodicity[cdim2]) {
if (fabs(d2) > prd_half[cdim2]) {
if (d2 < 0.0) d2 += prd[cdim2];
else d2 -= prd[cdim2];
}
}
}
r = sqrt(d1*d1 + d2*d2);
rbin = static_cast<int> ((r - cradmin) * cinvrad);

View File

@ -51,7 +51,7 @@ class ComputeChunkAtom : public Compute {
int which,binflag;
int regionflag,nchunksetflag,nchunkflag,discard;
int limit,limitstyle,limitfirst;
int scaleflag;
int scaleflag,pbcflag;
double xscale,yscale,zscale;
int argindex;
char *cfvid;

View File

@ -45,13 +45,14 @@ ComputeDisplaceAtom::ComputeDisplaceAtom(LAMMPS *lmp, int narg, char **arg) :
strcpy(id_fix,id);
strcat(id_fix,"_COMPUTE_STORE");
char **newarg = new char*[5];
char **newarg = new char*[6];
newarg[0] = id_fix;
newarg[1] = group->names[igroup];
newarg[2] = (char *) "STORE";
newarg[3] = (char *) "1";
newarg[4] = (char *) "3";
modify->add_fix(5,newarg);
newarg[3] = (char *) "peratom";
newarg[4] = (char *) "1";
newarg[5] = (char *) "3";
modify->add_fix(6,newarg);
fix = (FixStore *) modify->fix[modify->nfix-1];
delete [] newarg;

View File

@ -65,13 +65,14 @@ ComputeMSD::ComputeMSD(LAMMPS *lmp, int narg, char **arg) :
strcpy(id_fix,id);
strcat(id_fix,"_COMPUTE_STORE");
char **newarg = new char*[5];
char **newarg = new char*[6];
newarg[0] = id_fix;
newarg[1] = group->names[igroup];
newarg[2] = (char *) "STORE";
newarg[3] = (char *) "1";
newarg[4] = (char *) "3";
modify->add_fix(5,newarg);
newarg[3] = (char *) "peratom";
newarg[4] = (char *) "1";
newarg[5] = (char *) "3";
modify->add_fix(6,newarg);
fix = (FixStore *) modify->fix[modify->nfix-1];
delete [] newarg;

View File

@ -42,13 +42,14 @@ ComputeVACF::ComputeVACF(LAMMPS *lmp, int narg, char **arg) :
strcpy(id_fix,id);
strcat(id_fix,"_COMPUTE_STORE");
char **newarg = new char*[5];
char **newarg = new char*[6];
newarg[0] = id_fix;
newarg[1] = group->names[igroup];
newarg[2] = (char *) "STORE";
newarg[3] = (char *) "1";
newarg[4] = (char *) "3";
modify->add_fix(5,newarg);
newarg[3] = (char *) "peratom";
newarg[4] = (char *) "1";
newarg[5] = (char *) "3";
modify->add_fix(6,newarg);
fix = (FixStore *) modify->fix[modify->nfix-1];
delete [] newarg;

View File

@ -347,6 +347,13 @@ void CreateAtoms::command(int narg, char **arg)
}
}
// clear ghost count and any ghost bonus data internal to AtomVec
// same logic as beginning of Comm::exchange()
// do it now b/c creating atoms will overwrite ghost atoms
atom->nghost = 0;
atom->avec->clear_bonus();
// add atoms/molecules in one of 3 ways
bigint natoms_previous = atom->natoms;
@ -400,11 +407,10 @@ void CreateAtoms::command(int narg, char **arg)
if (atom->tag_enable) atom->tag_extend();
atom->tag_check();
// create global mapping of atoms
// zero nghost in case are adding new atoms to existing atoms
// if global map exists, reset it
// invoke map_init() b/c atom count has grown
if (atom->map_style) {
atom->nghost = 0;
atom->map_init();
atom->map_set();
}
@ -595,7 +601,7 @@ void CreateAtoms::add_single()
if (mode == ATOM) atom->avec->create_atom(ntype,xone);
else if (quatone[0] == 0.0 && quatone[1] == 0.0 && quatone[2] == 0.0)
add_molecule(xone);
else add_molecule(xone,&quatone[0]);
else add_molecule(xone,quatone);
}
}
@ -807,6 +813,10 @@ void CreateAtoms::add_molecule(double *center, double *quat_user)
int n;
double r[3],rotmat[3][3],quat[4],xnew[3];
if (quat_user) {
quat[0] = quat_user[0]; quat[1] = quat_user[1];
quat[2] = quat_user[2]; quat[3] = quat_user[3];
} else {
if (domain->dimension == 3) {
r[0] = ranmol->uniform() - 0.5;
r[1] = ranmol->uniform() - 0.5;
@ -815,16 +825,13 @@ void CreateAtoms::add_molecule(double *center, double *quat_user)
r[0] = r[1] = 0.0;
r[2] = 1.0;
}
if (quat_user) {
quat[0] = quat_user[0]; quat[1] = quat_user[1];
quat[2] = quat_user[2]; quat[3] = quat_user[3];
} else {
double theta = ranmol->uniform() * MY_2PI;
MathExtra::norm3(r);
double theta = ranmol->uniform() * MY_2PI;
MathExtra::axisangle_to_quat(r,theta,quat);
}
MathExtra::quat_to_mat(quat,rotmat);
onemol->quat_external = quat;
// create atoms in molecule with atom ID = 0 and mol ID = 0
// reset in caller after all moleclues created by all procs

View File

@ -210,11 +210,12 @@ void FixAdapt::post_constructor()
id_fix_diam = NULL;
id_fix_chg = NULL;
char **newarg = new char*[5];
char **newarg = new char*[6];
newarg[1] = group->names[igroup];
newarg[2] = (char *) "STORE";
newarg[3] = (char *) "1";
newarg[3] = (char *) "peratom";
newarg[4] = (char *) "1";
newarg[5] = (char *) "1";
if (diamflag) {
int n = strlen(id) + strlen("_FIX_STORE_DIAM") + 1;
@ -222,7 +223,7 @@ void FixAdapt::post_constructor()
strcpy(id_fix_diam,id);
strcat(id_fix_diam,"_FIX_STORE_DIAM");
newarg[0] = id_fix_diam;
modify->add_fix(5,newarg);
modify->add_fix(6,newarg);
fix_diam = (FixStore *) modify->fix[modify->nfix-1];
if (fix_diam->restart_reset) fix_diam->restart_reset = 0;
@ -245,7 +246,7 @@ void FixAdapt::post_constructor()
strcpy(id_fix_chg,id);
strcat(id_fix_chg,"_FIX_STORE_CHG");
newarg[0] = id_fix_chg;
modify->add_fix(5,newarg);
modify->add_fix(6,newarg);
fix_chg = (FixStore *) modify->fix[modify->nfix-1];
if (fix_chg->restart_reset) fix_chg->restart_reset = 0;

View File

@ -15,49 +15,91 @@
#include <string.h>
#include "fix_store.h"
#include "atom.h"
#include "comm.h"
#include "force.h"
#include "memory.h"
#include "error.h"
#include "force.h"
using namespace LAMMPS_NS;
using namespace FixConst;
enum{GLOBAL,PERATOM};
/* ---------------------------------------------------------------------- */
FixStore::FixStore(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg)
{
if (narg != 5) error->all(FLERR,"Illegal fix store command");
if (narg != 6) error->all(FLERR,"Illegal fix store command");
// syntax: id group style 0/1 nvalue
// 0/1 flag = not-store or store values in restart file
// 4th arg determines GLOBAL vs PERATOM values
// syntax: id group style peratom 0/1 nvalue
// 0/1 flag = not-store or store peratom values in restart file
// nvalue = # of peratom values, N=1 is vector, N>1 is array
// syntax: id group style global nrow ncol
// Nrow by Ncol array of global values
// Ncol=1 is vector, Nrow>1 is array
restart_peratom = force->inumeric(FLERR,arg[3]);
nvalues = force->inumeric(FLERR,arg[4]);
if (strcmp(arg[3],"peratom") == 0) flavor = PERATOM;
else if (strcmp(arg[3],"global") == 0) flavor = GLOBAL;
else error->all(FLERR,"Invalid fix store command");
// GLOBAL values are always written to restart file
// PERATOM restart_peratom is set by caller
if (flavor == GLOBAL) {
restart_global = 1;
nrow = force->inumeric(FLERR,arg[4]);
ncol = force->inumeric(FLERR,arg[5]);
if (nrow <= 0 || ncol <= 0)
error->all(FLERR,"Invalid fix store command");
vecflag = 0;
if (ncol == 1) vecflag = 1;
} else {
restart_peratom = force->inumeric(FLERR,arg[4]);
nvalues = force->inumeric(FLERR,arg[5]);
if (restart_peratom < 0 or restart_peratom > 1 || nvalues <= 0)
error->all(FLERR,"Invalid fix store command");
vecflag = 0;
if (nvalues == 1) vecflag = 1;
// perform initial allocation of atom-based array
// register with Atom class
}
vstore = NULL;
astore = NULL;
// allocate vector or array
// for PERATOM, register with Atom class
if (flavor == GLOBAL) {
if (vecflag) memory->create(vstore,nrow,"fix/store:vstore");
else memory->create(astore,nrow,ncol,"fix/store:astore");
}
if (flavor == PERATOM) {
grow_arrays(atom->nmax);
atom->add_callback(0);
if (restart_peratom) atom->add_callback(1);
}
// zero the storage
// since may be exchanged before filled by caller
// PERATOM may be exchanged before filled by caller
if (flavor == GLOBAL) {
if (vecflag)
for (int i = 0; i < nrow; i++) vstore[i] = 0.0;
else
for (int i = 0; i < nrow; i++)
for (int j = 0; j < ncol; j++)
astore[i][j] = 0.0;
}
if (flavor == PERATOM) {
int nlocal = atom->nlocal;
if (vecflag)
for (int i = 0; i < nlocal; i++)
vstore[i] = 0.0;
for (int i = 0; i < nlocal; i++) vstore[i] = 0.0;
else
for (int i = 0; i < nlocal; i++)
for (int j = 0; j < nvalues; j++)
astore[i][j] = 0.0;
}
}
/* ---------------------------------------------------------------------- */
@ -65,8 +107,10 @@ FixStore::~FixStore()
{
// unregister callbacks to this fix from Atom class
if (flavor == PERATOM) {
atom->delete_callback(id,0);
if (restart_peratom) atom->delete_callback(id,1);
}
memory->destroy(vstore);
memory->destroy(astore);
@ -81,13 +125,32 @@ int FixStore::setmask()
}
/* ----------------------------------------------------------------------
memory usage of local atom-based array
write global array to restart file
------------------------------------------------------------------------- */
double FixStore::memory_usage()
void FixStore::write_restart(FILE *fp)
{
double bytes = atom->nmax*nvalues * sizeof(double);
return bytes;
int n = nrow*ncol;
if (comm->me == 0) {
int size = n * sizeof(double);
fwrite(&size,sizeof(int),1,fp);
if (vecflag) fwrite(vstore,sizeof(double),n,fp);
else fwrite(&astore[0][0],sizeof(double),n,fp);
}
}
/* ----------------------------------------------------------------------
use global array from restart file to restart the Fix
------------------------------------------------------------------------- */
void FixStore::restart(char *buf)
{
// HOWTO insure size of buf is the same
int n = nrow*ncol;
double *dbuf = (double *) buf;
if (vecflag) memcpy(vstore,dbuf,n*sizeof(double));
else memcpy(&astore[0][0],dbuf,n*sizeof(double));
}
/* ----------------------------------------------------------------------
@ -189,3 +252,15 @@ int FixStore::size_restart(int nlocal)
{
return nvalues+1;
}
/* ----------------------------------------------------------------------
memory usage of local atom-based array
------------------------------------------------------------------------- */
double FixStore::memory_usage()
{
double bytes;
if (flavor == GLOBAL) bytes += nrow*ncol * sizeof(double);
if (flavor == PERATOM) bytes += atom->nmax*nvalues * sizeof(double);
return bytes;
}

View File

@ -20,20 +20,25 @@ FixStyle(STORE,FixStore)
#ifndef LMP_FIX_STORE_H
#define LMP_FIX_STORE_H
#include <stdio.h>
#include "fix.h"
namespace LAMMPS_NS {
class FixStore : public Fix {
public:
double *vstore; // vector storage if nvalues = 1
double **astore; // array storage if nvalues > 1
int nrow,ncol; // size of global data array
int nvalues; // number of per-atom values
double *vstore; // vector storage for GLOBAL or PERATOM
double **astore; // array storage for GLOBAL or PERATOM
FixStore(class LAMMPS *, int, char **);
~FixStore();
int setmask();
double memory_usage();
void write_restart(FILE *);
void restart(char *);
void grow_arrays(int);
void copy_arrays(int, int, int);
int pack_exchange(int, double *);
@ -43,9 +48,11 @@ class FixStore : public Fix {
int size_restart(int);
int maxsize_restart();
double memory_usage();
private:
int nvalues; // total # of values per atom
int vecflag; // 1 if nvalues = 1
int flavor; // GLOBAL or PERATOM
int vecflag; // 1 if ncol=1 or nvalues=1
};
}

View File

@ -105,6 +105,19 @@ FixStoreState::FixStoreState(LAMMPS *lmp, int narg, char **arg) :
if (domain->triclinic)
pack_choice[nvalues++] = &FixStoreState::pack_zu_triclinic;
else pack_choice[nvalues++] = &FixStoreState::pack_zu;
} else if (strcmp(arg[iarg],"xsu") == 0) {
if (domain->triclinic)
pack_choice[nvalues++] = &FixStoreState::pack_xsu_triclinic;
else pack_choice[nvalues++] = &FixStoreState::pack_xsu;
} else if (strcmp(arg[iarg],"ysu") == 0) {
if (domain->triclinic)
pack_choice[nvalues++] = &FixStoreState::pack_ysu_triclinic;
else pack_choice[nvalues++] = &FixStoreState::pack_ysu;
} else if (strcmp(arg[iarg],"zsu") == 0) {
if (domain->triclinic)
pack_choice[nvalues++] = &FixStoreState::pack_zsu_triclinic;
else pack_choice[nvalues++] = &FixStoreState::pack_zsu;
} else if (strcmp(arg[iarg],"ix") == 0) {
pack_choice[nvalues++] = &FixStoreState::pack_ix;
} else if (strcmp(arg[iarg],"iy") == 0) {
@ -145,12 +158,22 @@ FixStoreState::FixStoreState(LAMMPS *lmp, int narg, char **arg) :
error->all(FLERR,
"Fix store/state for atom property that isn't allocated");
pack_choice[nvalues++] = &FixStoreState::pack_muz;
} else if (strcmp(arg[iarg],"mu") == 0) {
if (!atom->mu_flag)
error->all(FLERR,
"Fix store/state for atom property that isn't allocated");
pack_choice[nvalues++] = &FixStoreState::pack_mu;
} else if (strcmp(arg[iarg],"radius") == 0) {
if (!atom->radius_flag)
error->all(FLERR,
"Fix store/state for atom property that isn't allocated");
pack_choice[nvalues++] = &FixStoreState::pack_radius;
} else if (strcmp(arg[iarg],"diameter") == 0) {
if (!atom->radius_flag)
error->all(FLERR,
"Fix store/state for atom property that isn't allocated");
pack_choice[nvalues++] = &FixStoreState::pack_diameter;
} else if (strcmp(arg[iarg],"omegax") == 0) {
if (!atom->omega_flag)
error->all(FLERR,
@ -989,6 +1012,128 @@ void FixStoreState::pack_zu_triclinic(int n)
/* ---------------------------------------------------------------------- */
void FixStoreState::pack_xsu(int n)
{
double **x = atom->x;
imageint *image = atom->image;
int *mask = atom->mask;
int nlocal = atom->nlocal;
double boxxlo = domain->boxlo[0];
double invxprd = 1.0/domain->xprd;
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit)
vbuf[n] = (x[i][0]-boxxlo)*invxprd + ((image[i] & IMGMASK) - IMGMAX);
else vbuf[n] = 0.0;
n += nvalues;
}
}
/* ---------------------------------------------------------------------- */
void FixStoreState::pack_ysu(int n)
{
double **x = atom->x;
imageint *image = atom->image;
int *mask = atom->mask;
int nlocal = atom->nlocal;
double boxylo = domain->boxlo[1];
double invyprd = 1.0/domain->yprd;
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit)
vbuf[n] = (x[i][1]-boxylo)*invyprd + (image[i] >> IMGBITS & IMGMASK) - IMGMAX;
else vbuf[n] = 0.0;
n += nvalues;
}
}
/* ---------------------------------------------------------------------- */
void FixStoreState::pack_zsu(int n)
{
double **x = atom->x;
imageint *image = atom->image;
int *mask = atom->mask;
int nlocal = atom->nlocal;
double boxzlo = domain->boxlo[2];
double invzprd = 1.0/domain->zprd;
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit)
vbuf[n] = (x[i][2]-boxzlo)*invzprd + (image[i] >> IMG2BITS) - IMGMAX;
else vbuf[n] = 0.0;
n += nvalues;
}
}
/* ---------------------------------------------------------------------- */
void DumpCustom::pack_xsu_triclinic(int n)
{
double **x = atom->x;
imageint *image = atom->image;
int *mask = atom->mask;
int nlocal = atom->nlocal;
double *boxlo = domain->boxlo;
double *h_inv = domain->h_inv;
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit)
vbuf[n] = h_inv[0]*(x[i][0]-boxlo[0]) + h_inv[5]*(x[i][1]-boxlo[1]) +
h_inv[4]*(x[i][2]-boxlo[2]) + (image[i] & IMGMASK) - IMGMAX;
else vbuf[n] = 0.0;
n += nvalues;
}
}
/* ---------------------------------------------------------------------- */
void DumpCustom::pack_ysu_triclinic(int n)
{
double **x = atom->x;
imageint *image = atom->image;
int *mask = atom->mask;
int nlocal = atom->nlocal;
double *boxlo = domain->boxlo;
double *h_inv = domain->h_inv;
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit)
vbuf[n] = h_inv[1]*(x[i][1]-boxlo[1]) + h_inv[3]*(x[i][2]-boxlo[2]) +
(image[i] >> IMGBITS & IMGMASK) - IMGMAX;
else vbuf[n] = 0.0;
n += nvalues;
}
}
/* ---------------------------------------------------------------------- */
void DumpCustom::pack_zsu_triclinic(int n)
{
double **x = atom->x;
imageint *image = atom->image;
int *mask = atom->mask;
int nlocal = atom->nlocal;
double *boxlo = domain->boxlo;
double *h_inv = domain->h_inv;
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit)
buf[n] = h_inv[2]*(x[i][2]-boxlo[2]) + (image[i] >> IMG2BITS) - IMGMAX;
else vbuf[n] = 0.0;
n += nvalues;
}
}
/* ---------------------------------------------------------------------- */
void FixStoreState::pack_ix(int n)
{
imageint *image = atom->image;
@ -1184,6 +1329,21 @@ void FixStoreState::pack_muz(int n)
/* ---------------------------------------------------------------------- */
void FixStoreState::pack_mu(int n)
{
double **mu = atom->mu;
int *mask = atom->mask;
int nlocal = atom->nlocal;
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) vbuf[n] = mu[i][3];
else vbuf[n] = 0.0;
n += nvalues;
}
}
/* ---------------------------------------------------------------------- */
void FixStoreState::pack_radius(int n)
{
double *radius = atom->radius;
@ -1199,6 +1359,21 @@ void FixStoreState::pack_radius(int n)
/* ---------------------------------------------------------------------- */
void FixStoreState::pack_diameter(int n)
{
double *radius = atom->radius;
int *mask = atom->mask;
int nlocal = atom->nlocal;
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) vbuf[n] = 2.0*radius[i];
else vbuf[n] = 0.0;
n += nvalues;
}
}
/* ---------------------------------------------------------------------- */
void FixStoreState::pack_omegax(int n)
{
double **omega = atom->omega;

View File

@ -79,6 +79,13 @@ class FixStoreState : public Fix {
void pack_xu_triclinic(int);
void pack_yu_triclinic(int);
void pack_zu_triclinic(int);
void pack_xsu(int);
void pack_ysu(int);
void pack_zsu(int);
void pack_xsu_triclinic(int);
void pack_ysu_triclinic(int);
void pack_zsu_triclinic(int);
void pack_ix(int);
void pack_iy(int);
void pack_iz(int);
@ -93,7 +100,9 @@ class FixStoreState : public Fix {
void pack_mux(int);
void pack_muy(int);
void pack_muz(int);
void pack_mu(int);
void pack_radius(int);
void pack_diameter(int);
void pack_omegax(int);
void pack_omegay(int);
void pack_omegaz(int);

View File

@ -16,6 +16,7 @@
#include "molecule.h"
#include "atom.h"
#include "atom_vec.h"
#include "atom_vec_body.h"
#include "force.h"
#include "comm.h"
#include "domain.h"
@ -163,7 +164,8 @@ Molecule::~Molecule()
compute center = geometric center of molecule
also compute:
dx = displacement of each atom from center
molradius = radius of molecule from center including finite-size particles
molradius = radius of molecule from center
including finite-size particles or body particles
------------------------------------------------------------------------- */
void Molecule::compute_center()
@ -446,6 +448,13 @@ void Molecule::read(int flag)
itensor[4] *= sizescale*sizescale*sizescale*sizescale*sizescale;
itensor[5] *= sizescale*sizescale*sizescale*sizescale*sizescale;
}
else if (strstr(line,"body")) {
bodyflag = 1;
avec_body = (AtomVecBody *) atom->style_match("body");
if (!avec_body)
error->all(FLERR,"Molecule file requires atom style body");
sscanf(line,"%d %d",&nibody,&ndbody);
}
else break;
}
@ -542,6 +551,19 @@ void Molecule::read(int flag)
if (flag) shaketype_read(line);
else skip_lines(natoms,line);
} else if (strcmp(keyword,"Body Integers") == 0) {
if (bodyflag == 0 || nibody == 0)
error->all(FLERR,"Molecule file has body params "
"but no setting for them");
ibodyflag = 1;
body(flag,0,line);
} else if (strcmp(keyword,"Body Doubles") == 0) {
if (bodyflag == 0 || ndbody == 0)
error->all(FLERR,"Molecule file has body params "
"but no setting for them");
dbodyflag = 1;
body(flag,1,line);
} else error->one(FLERR,"Unknown section in molecule file");
parse_keyword(1,line,keyword);
@ -560,6 +582,10 @@ void Molecule::read(int flag)
error->all(FLERR,"Molecule file has special flags but no bonds");
if ((shakeflagflag || shakeatomflag || shaketypeflag) && !shakeflag)
error->all(FLERR,"Molecule file shake info is incomplete");
if (bodyflag && nibody && ibodyflag == 0)
error->all(FLERR,"Molecule file has no Body Integers section");
if (bodyflag && ndbody && dbodyflag == 0)
error->all(FLERR,"Molecule file has no Body Doubles section");
}
// auto-generate special bonds
@ -570,6 +596,21 @@ void Molecule::read(int flag)
maxspecial = atom->maxspecial;
if (flag) special_generate();
}
// body particle must have natom = 1
// set radius by having body class compute its own radius
if (bodyflag) {
radiusflag = 1;
if (natoms != 1)
error->all(FLERR,"Molecule natoms must be 1 for body particle");
if (sizescale != 1.0)
error->all(FLERR,"Molecule sizescale must be 1.0 for body particle");
if (flag) {
radius[0] = avec_body->radius_body(nibody,ndbody,ibodyparams,dbodyparams);
maxradius = radius[0];
}
}
}
/* ----------------------------------------------------------------------
@ -1244,6 +1285,44 @@ void Molecule::shaketype_read(char *line)
}
}
/* ----------------------------------------------------------------------
read body params from file
pflag = 0/1 for integer/double params
------------------------------------------------------------------------- */
void Molecule::body(int flag, int pflag, char *line)
{
int i,ncount;
int nparam = nibody;
if (pflag) nparam = ndbody;
int nword = 0;
while (nword < nparam) {
readline(line);
ncount = atom->count_words(line);
if (ncount == 0)
error->one(FLERR,"Too few values in body section of molecule file");
if (nword+ncount > nparam)
error->all(FLERR,"Too many values in body section of molecule file");
if (flag) {
if (pflag == 0) {
ibodyparams[nword++] = force->inumeric(FLERR,strtok(line," \t\n\r\f"));
for (i = 1; i < ncount; i++)
ibodyparams[nword++] =
force->inumeric(FLERR,strtok(NULL," \t\n\r\f"));
} else {
dbodyparams[nword++] = force->numeric(FLERR,strtok(line," \t\n\r\f"));
for (i = 1; i < ncount; i++)
dbodyparams[nword++] =
force->numeric(FLERR,strtok(NULL," \t\n\r\f"));
}
} else nword += ncount;
}
}
/* ----------------------------------------------------------------------
error check molecule attributes and topology against system settings
flag = 0, just check this molecule
@ -1318,6 +1397,7 @@ void Molecule::initialize()
nbonds = nangles = ndihedrals = nimpropers = 0;
ntypes = 0;
nbondtypes = nangletypes = ndihedraltypes = nimpropertypes = 0;
nibody = ndbody = 0;
bond_per_atom = angle_per_atom = dihedral_per_atom = improper_per_atom = 0;
maxspecial = 0;
@ -1326,6 +1406,7 @@ void Molecule::initialize()
bondflag = angleflag = dihedralflag = improperflag = 0;
nspecialflag = specialflag = 0;
shakeflag = shakeflagflag = shakeatomflag = shaketypeflag = 0;
bodyflag = ibodyflag = dbodyflag = 0;
centerflag = massflag = comflag = inertiaflag = 0;
tag_require = 0;
@ -1359,6 +1440,9 @@ void Molecule::initialize()
shake_atom = NULL;
shake_type = NULL;
ibodyparams = NULL;
dbodyparams = NULL;
dx = NULL;
dxcom = NULL;
dxbody = NULL;
@ -1445,6 +1529,11 @@ void Molecule::allocate()
memory->create(shake_atom,natoms,4,"molecule:shake_flag");
memory->create(shake_type,natoms,3,"molecule:shake_flag");
}
if (bodyflag) {
if (nibody) memory->create(ibodyparams,nibody,"molecule:ibodyparams");
if (ndbody) memory->create(dbodyparams,ndbody,"molecule:dbodyparams");
}
}
/* ----------------------------------------------------------------------
@ -1493,6 +1582,9 @@ void Molecule::deallocate()
memory->destroy(dx);
memory->destroy(dxcom);
memory->destroy(dxbody);
memory->destroy(ibodyparams);
memory->destroy(dbodyparams);
}
/* ----------------------------------------------------------------------

View File

@ -26,11 +26,13 @@ class Molecule : protected Pointers {
int last; // 1 if last molecule in set, else 0
// number of atoms,bonds,etc in molecule
// nibody,ndbody = # of integer/double fields in body
int natoms;
int nbonds,nangles,ndihedrals,nimpropers;
int ntypes;
int nbondtypes,nangletypes,ndihedraltypes,nimpropertypes;
int nibody,ndbody;
// max bond,angle,etc per atom
@ -43,6 +45,7 @@ class Molecule : protected Pointers {
int bondflag,angleflag,dihedralflag,improperflag;
int nspecialflag,specialflag;
int shakeflag,shakeflagflag,shakeatomflag,shaketypeflag;
int bodyflag,ibodyflag,dbodyflag;
// 1 if attribute defined or computed, 0 if not
@ -83,6 +86,10 @@ class Molecule : protected Pointers {
tagint **shake_atom;
int **shake_type;
class AtomVecBody *avec_body;
int *ibodyparams; // integer and double body params
double *dbodyparams;
double center[3]; // geometric center of molecule
double masstotal; // total mass of molecule
double com[3]; // center of mass of molecule
@ -102,6 +109,9 @@ class Molecule : protected Pointers {
double **dxbody; // displacement of each atom relative to COM
// in body frame (diagonalized interia tensor)
double *quat_external; // orientation imposed by external class
// e.g. FixPour or CreateAtoms
Molecule(class LAMMPS *, int, char **, int &);
~Molecule();
void compute_center();
@ -134,6 +144,7 @@ class Molecule : protected Pointers {
void shakeflag_read(char *);
void shakeatom_read(char *);
void shaketype_read(char *);
void body(int, int, char *);
void initialize();
void allocate();

View File

@ -1469,7 +1469,7 @@ void ReadData::bonus(bigint nbonus, AtomVec *ptr, const char *type)
void ReadData::bodies(int firstpass)
{
int i,m,nchunk,nline,nmax,ninteger,ndouble,nword,onebody,tmp;
int i,m,nchunk,nline,nmax,ninteger,ndouble,nword,ncount,onebody,tmp;
char *eof;
int mapflag = 0;
@ -1500,8 +1500,7 @@ void ReadData::bodies(int firstpass)
sscanf(&buffer[m],"%d %d %d",&tmp,&ninteger,&ndouble);
m += strlen(&buffer[m]);
// read lines one at a time into buffer
// make copy of line and count words
// read lines one at a time into buffer and count words
// count to ninteger and ndouble until have enough lines
onebody = 0;
@ -1510,23 +1509,29 @@ void ReadData::bodies(int firstpass)
while (nword < ninteger) {
eof = fgets(&buffer[m],MAXLINE,fp);
if (eof == NULL) error->one(FLERR,"Unexpected end of data file");
nword += atom->count_words(&buffer[m],copy);
ncount = atom->count_words(&buffer[m],copy);
if (ncount == 0)
error->one(FLERR,"Too few values in body lines in data file");
nword += ncount;
m += strlen(&buffer[m]);
onebody++;
}
if (nword > ninteger)
error->one(FLERR,"Too many value in body lines in data file");
error->one(FLERR,"Too many values in body lines in data file");
nword = 0;
while (nword < ndouble) {
eof = fgets(&buffer[m],MAXLINE,fp);
if (eof == NULL) error->one(FLERR,"Unexpected end of data file");
nword += atom->count_words(&buffer[m],copy);
ncount = atom->count_words(&buffer[m],copy);
if (ncount == 0)
error->one(FLERR,"Too few values in body lines in data file");
nword += ncount;
m += strlen(&buffer[m]);
onebody++;
}
if (nword > ndouble)
error->one(FLERR,"Too many value in body lines in data file");
error->one(FLERR,"Too many values in body lines in data file");
if (onebody+1 > MAXBODY)
error->one(FLERR,

View File

@ -21,6 +21,7 @@
#include "atom_vec_ellipsoid.h"
#include "atom_vec_line.h"
#include "atom_vec_tri.h"
#include "atom_vec_body.h"
#include "domain.h"
#include "region.h"
#include "group.h"
@ -41,7 +42,7 @@ using namespace MathConst;
enum{ATOM_SELECT,MOL_SELECT,TYPE_SELECT,GROUP_SELECT,REGION_SELECT};
enum{TYPE,TYPE_FRACTION,MOLECULE,X,Y,Z,CHARGE,MASS,SHAPE,LENGTH,TRI,
DIPOLE,DIPOLE_RANDOM,QUAT,QUAT_RANDOM,THETA,ANGMOM,OMEGA,
DIPOLE,DIPOLE_RANDOM,QUAT,QUAT_RANDOM,THETA,THETA_RANDOM,ANGMOM,OMEGA,
DIAMETER,DENSITY,VOLUME,IMAGE,BOND,ANGLE,DIHEDRAL,IMPROPER,
MESO_E,MESO_CV,MESO_RHO,SMD_MASS_DENSITY,SMD_CONTACT_RADIUS,INAME,DNAME};
@ -221,7 +222,7 @@ void Set::command(int narg, char **arg)
else zvalue = force->numeric(FLERR,arg[iarg+3]);
if (strstr(arg[iarg+4],"v_") == arg[iarg+4]) varparse(arg[iarg+4],4);
else wvalue = force->numeric(FLERR,arg[iarg+4]);
if (!atom->ellipsoid_flag && !atom->tri_flag)
if (!atom->ellipsoid_flag && !atom->tri_flag && !atom->body_flag)
error->all(FLERR,"Cannot set this attribute for this atom style");
set(QUAT);
iarg += 5;
@ -229,7 +230,7 @@ void Set::command(int narg, char **arg)
} else if (strcmp(arg[iarg],"quat/random") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal set command");
ivalue = force->inumeric(FLERR,arg[iarg+1]);
if (!atom->ellipsoid_flag && !atom->tri_flag)
if (!atom->ellipsoid_flag && !atom->tri_flag && !atom->body_flag)
error->all(FLERR,"Cannot set this attribute for this atom style");
if (ivalue <= 0)
error->all(FLERR,"Invalid random number seed in set command");
@ -248,6 +249,16 @@ void Set::command(int narg, char **arg)
set(THETA);
iarg += 2;
} else if (strcmp(arg[iarg],"theta/random") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal set command");
ivalue = force->inumeric(FLERR,arg[iarg+1]);
if (!atom->line_flag)
error->all(FLERR,"Cannot set this attribute for this atom style");
if (ivalue <= 0)
error->all(FLERR,"Invalid random number seed in set command");
set(THETA_RANDOM);
iarg += 2;
} else if (strcmp(arg[iarg],"angmom") == 0) {
if (iarg+4 > narg) error->all(FLERR,"Illegal set command");
if (strstr(arg[iarg+1],"v_") == arg[iarg+1]) varparse(arg[iarg+1],1);
@ -256,7 +267,7 @@ void Set::command(int narg, char **arg)
else yvalue = force->numeric(FLERR,arg[iarg+2]);
if (strstr(arg[iarg+3],"v_") == arg[iarg+3]) varparse(arg[iarg+3],3);
else zvalue = force->numeric(FLERR,arg[iarg+3]);
if (!atom->ellipsoid_flag && !atom->tri_flag)
if (!atom->angmom_flag)
error->all(FLERR,"Cannot set this attribute for this atom style");
set(ANGMOM);
iarg += 4;
@ -269,7 +280,7 @@ void Set::command(int narg, char **arg)
else yvalue = force->numeric(FLERR,arg[iarg+2]);
if (strstr(arg[iarg+3],"v_") == arg[iarg+3]) varparse(arg[iarg+3],3);
else zvalue = force->numeric(FLERR,arg[iarg+3]);
if (!atom->sphere_flag)
if (!atom->omega_flag)
error->all(FLERR,"Cannot set this attribute for this atom style");
set(OMEGA);
iarg += 4;
@ -558,6 +569,7 @@ void Set::set(int keyword)
(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");
int nlocal = atom->nlocal;
for (int i = 0; i < nlocal; i++) {
@ -604,7 +616,8 @@ void Set::set(int keyword)
else if (keyword == MESO_E) atom->e[i] = dvalue;
else if (keyword == MESO_CV) atom->cv[i] = dvalue;
else if (keyword == MESO_RHO) atom->rho[i] = dvalue;
else if (keyword == SMD_MASS_DENSITY) { // set mass from volume and supplied mass density
else if (keyword == SMD_MASS_DENSITY) {
// set mass from volume and supplied mass density
atom->rmass[i] = atom->vfrac[i] * dvalue;
}
else if (keyword == SMD_CONTACT_RADIUS) atom->contact_radius[i] = dvalue;
@ -678,7 +691,8 @@ void Set::set(int keyword)
mu[i][2]*mu[i][2]);
}
// set quaternion orientation of ellipsoid or tri particle
// set quaternion orientation of ellipsoid or tri or body particle
// enforce quat rotation vector in z dir for 2d systems
else if (keyword == QUAT) {
double *quat;
@ -686,8 +700,13 @@ void Set::set(int keyword)
quat = avec_ellipsoid->bonus[atom->ellipsoid[i]].quat;
else if (avec_tri && atom->tri[i] >= 0)
quat = avec_tri->bonus[atom->tri[i]].quat;
else if (avec_body && atom->body[i] >= 0)
quat = avec_body->bonus[atom->body[i]].quat;
else
error->one(FLERR,"Cannot set quaternion for atom that has none");
if (domain->dimension == 2 && (xvalue != 0.0 || yvalue != 0.0))
error->one(FLERR,"Cannot set quaternion with xy components "
"for 2d system");
double theta2 = MY_PI2 * wvalue/180.0;
double sintheta2 = sin(theta2);
@ -706,7 +725,7 @@ void Set::set(int keyword)
avec_line->bonus[atom->line[i]].theta = dvalue;
}
// set angmom of ellipsoidal or tri particle
// set angmom or omega of particle
else if (keyword == ANGMOM) {
atom->angmom[i][0] = xvalue;
@ -720,7 +739,6 @@ void Set::set(int keyword)
atom->omega[i][2] = zvalue;
}
// reset any or all of 3 image flags
else if (keyword == IMAGE) {
@ -768,7 +786,9 @@ void Set::setrandom(int keyword)
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");
RanPark *random = new RanPark(lmp,1);
double **x = atom->x;
@ -828,7 +848,7 @@ void Set::setrandom(int keyword)
}
}
// set quaternions to random orientations in 3d or 2d
// set quaternions to random orientations in 3d and 2d
} else if (keyword == QUAT_RANDOM) {
int nlocal = atom->nlocal;
@ -842,6 +862,8 @@ void Set::setrandom(int keyword)
quat = avec_ellipsoid->bonus[atom->ellipsoid[i]].quat;
else if (avec_tri && atom->tri[i] >= 0)
quat = avec_tri->bonus[atom->tri[i]].quat;
else if (avec_body && atom->body[i] >= 0)
quat = avec_body->bonus[atom->body[i]].quat;
else
error->one(FLERR,"Cannot set quaternion for atom that has none");
@ -864,6 +886,8 @@ void Set::setrandom(int keyword)
if (select[i]) {
if (avec_ellipsoid && atom->ellipsoid[i] >= 0)
quat = avec_ellipsoid->bonus[atom->ellipsoid[i]].quat;
else if (avec_body && atom->body[i] >= 0)
quat = avec_body->bonus[atom->body[i]].quat;
else
error->one(FLERR,"Cannot set quaternion for atom that has none");
@ -876,6 +900,21 @@ void Set::setrandom(int keyword)
count++;
}
}
// set theta to random orientation in 2d
} else if (keyword == QUAT_RANDOM) {
int nlocal = atom->nlocal;
double theta;
for (i = 0; i < nlocal; i++) {
if (select[i]) {
if (atom->line[i] < 0)
error->one(FLERR,"Cannot set theta for atom that is not a line");
random->reset(seed,x[i]);
avec_line->bonus[atom->line[i]].theta = MY_2PI*random->uniform();
count++;
}
}
}
delete random;

View File

@ -4480,13 +4480,14 @@ VarReader::VarReader(LAMMPS *lmp, char *name, char *file, int flag) :
strcpy(id_fix,name);
strcat(id_fix,"_VARIABLE_STORE");
char **newarg = new char*[5];
char **newarg = new char*[6];
newarg[0] = id_fix;
newarg[1] = (char *) "all";
newarg[2] = (char *) "STORE";
newarg[3] = (char *) "0";
newarg[4] = (char *) "1";
modify->add_fix(5,newarg);
newarg[3] = (char *) "peratom";
newarg[4] = (char *) "0";
newarg[5] = (char *) "1";
modify->add_fix(6,newarg);
fixstore = (FixStore *) modify->fix[modify->nfix-1];
delete [] newarg;

View File

@ -1 +1 @@
#define LAMMPS_VERSION "24 Dec 2015"
#define LAMMPS_VERSION "11 Jan 2016"