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

Merge branch 'lammps-icms' of https://github.com/lammps/lammps into weighted_load_balance

Browse Source
This commit is contained in:
Iain Bethune
2016-07-29 15:05:05 +01:00
parent fd8794f52a ba0147f747
commit 48cf74f49e
52 changed files with 1101 additions and 311 deletions
Download Patch File Download Diff File Expand all files Collapse all files

2
doc/html/Manual.html
View File

@ -135,7 +135,7 @@
<H1></H1><div class="section" id="lammps-icms-documentation">
<h1>LAMMPS-ICMS Documentation</h1>
<div class="section" id="jul-2016-version">
<h2>22 Jul 2016 version</h2>
<h2>27 Jul 2016 version</h2>
</div>
<div class="section" id="version-info">
<h2>Version info:</h2>

2
doc/html/_sources/Manual.txt
View File

@ -5,7 +5,7 @@
LAMMPS Documentation
====================
22 Jul 2016 version
27 Jul 2016 version
-------------------
Version info:

73
doc/html/_sources/compute_centro_atom.txt
View File

@ -8,11 +8,20 @@ Syntax
.. parsed-literal::
compute ID group-ID centro/atom lattice
compute ID group-ID centro/atom lattice keyword value ...
* ID, group-ID are documented in :doc:`compute <compute>` command
* centro/atom = style name of this compute command
* lattice = *fcc* or *bcc* or N = # of neighbors per atom to include
centro/atom = style name of this compute command
lattice = *fcc* or *bcc* or N = # of neighbors per atom to include
* zero or more keyword/value pairs may be appended
* keyword = *axes*
.. parsed-literal::
*axes* value = *no* or *yes*
*no* = do not calulate 3 symmetry axes
*yes* = calulate 3 symmetry axes
Examples
""""""""
@ -29,11 +38,12 @@ Description
"""""""""""
Define a computation that calculates the centro-symmetry parameter for
each atom in the group. In solid-state systems the centro-symmetry
parameter is a useful measure of the local lattice disorder around an
atom and can be used to characterize whether the atom is part of a
perfect lattice, a local defect (e.g. a dislocation or stacking
fault), or at a surface.
each atom in the group, for either FCC or BCC lattices, depending on
the choice of the *lattice* argument. In solid-state systems the
centro-symmetry parameter is a useful measure of the local lattice
disorder around an atom and can be used to characterize whether the
atom is part of a perfect lattice, a local defect (e.g. a dislocation
or stacking fault), or at a surface.
The value of the centro-symmetry parameter will be 0.0 for atoms not
in the specified compute group.
@ -44,7 +54,7 @@ This parameter is computed using the following formula from
.. image:: Eqs/centro_symmetry.jpg
:align: center
where the *N* nearest neighbors or each atom are identified and Ri and
where the *N* nearest neighbors of each atom are identified and Ri and
Ri+N/2 are vectors from the central atom to a particular pair of
nearest neighbors. There are N*(N-1)/2 possible neighbor pairs that
can contribute to this formula. The quantity in the sum is computed
@ -67,6 +77,21 @@ positive parameter. If the atom does not have *N* neighbors (within
the potential cutoff), then its centro-symmetry parameter is set to
0.0.
If the keyword *axes* has the setting *yes*\ , then this compute also
estimates three symmetry axes for each atom's local neighborhood. The
first two of these are the vectors joining the two pairs of neighbor
atoms with smallest contributions to the centrosymmetry parameter,
i.e. the two most symmetric pairs of atoms. The third vector is
normal to the first two by the right-hand rule. All three vectors are
normalized to unit length. For FCC crystals, the first two vectors
will lie along a <110> direction, while the third vector will lie
along either a <100> or <111> direction. For HCP crystals, the first
two vectors will lie along <1000> directions, while the third vector
will lie along <0001>. This provides a simple way to measure local
orientation in HCP structures. In general, the *axes* keyword can be
used to estimate the orientation of symmetry axes in the neighborhood
of any atom.
Only atoms within the cutoff of the pairwise neighbor list are
considered as possible neighbors. Atoms not in the compute group are
included in the *N* neighbors used in this calculation.
@ -79,16 +104,23 @@ too frequently or to have multiple compute/dump commands, each with a
**Output info:**
This compute calculates a per-atom vector, which can be accessed by
any command that uses per-atom values from a compute as input. See
:ref:`Section_howto 15 <howto_15>` for an overview of
LAMMPS output options.
By default, this compute calculates the centrosymmetry value for each
atom as a per-atom vector, which can be accessed by any command that
uses per-atom values from a compute as input. See :ref:`Section_howto 15 <howto_15>` for an overview of LAMMPS output
options.
The per-atom vector values are unitless values >= 0.0. Their
magnitude depends on the lattice style due to the number of
contibuting neighbor pairs in the summation in the formula above. And
it depends on the local defects surrounding the central atom, as
described above.
If the *axes* keyword setting is *yes*\ , then a per-atom array is
calculated. The first column is the centrosymmetry parameter. The
next three columns are the x, y, and z components of the first
symmetry axis, followed by the second, and third symmetry axes in
columns 5-7 and 8-10.
The centrosymmetry values are unitless values >= 0.0. Their magnitude
depends on the lattice style due to the number of contibuting neighbor
pairs in the summation in the formula above. And it depends on the
local defects surrounding the central atom, as described above. For
the *axes yes* case, the vector components are also unitless, since
they represent spatial directions.
Here are typical centro-symmetry values, from a a nanoindentation
simulation into gold (FCC). These were provided by Jon Zimmerman
@ -124,7 +156,10 @@ Related commands
:doc:`compute cna/atom <compute_cna_atom>`
**Default:** none
Default
"""""""
The default value for the optional keyword is axes = no.
----------

26
doc/html/_sources/fix_nh.txt
View File

@ -76,7 +76,9 @@ Syntax
*flip* value = *yes* or *no* = allow or disallow box flips when it becomes highly skewed
*fixedpoint* values = x y z
x,y,z = perform barostat dilation/contraction around this point (distance units)
*update* value = *dipole* update dipole orientation (only for sphere variants)
*update* value = *dipole* or *dipole/dlm*
dipole = update dipole orientation (only for sphere variants)
dipole/dlm = use DLM integrator to update dipole orientation (only for sphere variants)
@ -376,6 +378,13 @@ during the time integration. This option should be used for models
where a dipole moment is assigned to finite-size particles,
e.g. spheroids via use of the :doc:`atom_style hybrid sphere dipole <atom_style>` command.
The default dipole orientation integrator can be changed to the
Dullweber-Leimkuhler-McLachlan integration scheme
:ref:`(Dullweber) <nh-Dullweber>` when using *update* with the value
*dipole/dlm*\ . This integrator is symplectic and time-reversible,
giving better energy conservation and allows slightly longer timesteps
at only a small additional computational cost.
----------
@ -479,7 +488,7 @@ thermal degrees of freedom, and the bias is added back in.
----------
These fixes can be used with either the *verlet* or *respa*
These fixes can be used with either the *verlet* or *respa*
:doc:`integrators <run_style>`. When using one of the barostat fixes
with *respa*\ , LAMMPS uses an integrator constructed
according to the following factorization of the Liouville propagator
@ -505,7 +514,7 @@ of the underlying non-Hamiltonian equations of motion.
up to machine precision under NVT dynamics. Under NPT dynamics,
for a system with zero initial total linear momentum, the total
momentum fluctuates close to zero. It may occasionally undergo brief
excursions to non-negligible values, before returning close to zero.
excursions to non-negligible values, before returning close to zero.
Over long simulations, this has the effect of causing the center-of-mass
to undergo a slow random walk. This can be mitigated by resetting
the momentum at infrequent intervals using the
@ -517,7 +526,7 @@ of the underlying non-Hamiltonian equations of motion.
up to machine precision under NVT dynamics. Under NPT dynamics,
for a system with zero initial total linear momentum, the total
momentum fluctuates close to zero. It may occasionally undergo brief
excursions to non-negligible values, before returning close to zero.
excursions to non-negligible values, before returning close to zero.
Over long simulations, this has the effect of causing the center-of-mass
to undergo a slow random walk. This can be mitigated by resetting
the momentum at infrequent intervals using the
@ -685,7 +694,7 @@ Default
The keyword defaults are tchain = 3, pchain = 3, mtk = yes, tloop =
ploop = 1, nreset = 0, drag = 0.0, dilate = all, couple = none,
scaleyz = scalexz = scalexy = yes if periodic in 2nd dimension and
scaleyz = scalexz = scalexy = yes if periodic in 2nd dimension and
not coupled to barostat, otherwise no.
@ -717,6 +726,13 @@ Martyna, J Phys A: Math Gen, 39, 5629 (2006).
**(Shinoda)** Shinoda, Shiga, and Mikami, Phys Rev B, 69, 134103 (2004).
.. _nh-Dullweber:
**(Dullweber)** Dullweber, Leimkuhler and McLachlan, J Chem Phys, 107,
5840 (1997).
.. _lws: http://lammps.sandia.gov
.. _ld: Manual.html

24
doc/html/_sources/fix_nve_sphere.txt
View File

@ -16,11 +16,12 @@ Syntax
* ID, group-ID are documented in :doc:`fix <fix>` command
* nve/sphere = style name of this fix command
* zero or more keyword/value pairs may be appended
* keyword = *update*
.. parsed-literal::
*update* value = *dipole*
keyword = *update*
*update* value = *dipole* or *dipole/dlm*
dipole = update orientation of dipole moment during integration
dipole/dlm = use DLM integrator to update dipole orientation
@ -31,6 +32,7 @@ Examples
fix 1 all nve/sphere
fix 1 all nve/sphere update dipole
fix 1 all nve/sphere update dipole/dlm
Description
"""""""""""
@ -49,6 +51,13 @@ during the time integration. This option should be used for models
where a dipole moment is assigned to finite-size particles,
e.g. spheroids via use of the :doc:`atom_style hybrid sphere dipole <atom_style>` command.
The default dipole orientation integrator can be changed to the
Dullweber-Leimkuhler-McLachlan integration scheme
:ref:`(Dullweber) <nh-Dullweber>` when using *update* with the value
*dipole/dlm*\ . This integrator is symplectic and time-reversible,
giving better energy conservation and allows slightly longer timesteps
at only a small additional computational cost.
----------
@ -106,6 +115,17 @@ Related commands
**Default:** none
----------
.. _nve-Dullweber:
**(Dullweber)** Dullweber, Leimkuhler and McLachlan, J Chem Phys, 107,
5840 (1997).
.. _lws: http://lammps.sandia.gov
.. _ld: Manual.html
.. _lc: Section_commands.html#comm

71
doc/html/compute_centro_atom.html
View File

@ -128,14 +128,21 @@
<span id="index-0"></span><h1>compute centro/atom command</h1>
<div class="section" id="syntax">
<h2>Syntax</h2>
<div class="highlight-default"><div class="highlight"><pre><span></span><span class="n">compute</span> <span class="n">ID</span> <span class="n">group</span><span class="o">-</span><span class="n">ID</span> <span class="n">centro</span><span class="o">/</span><span class="n">atom</span> <span class="n">lattice</span>
<div class="highlight-default"><div class="highlight"><pre><span></span><span class="n">compute</span> <span class="n">ID</span> <span class="n">group</span><span class="o">-</span><span class="n">ID</span> <span class="n">centro</span><span class="o">/</span><span class="n">atom</span> <span class="n">lattice</span> <span class="n">keyword</span> <span class="n">value</span> <span class="o">...</span>
</pre></div>
</div>
<ul class="simple">
<li>ID, group-ID are documented in <a class="reference internal" href="compute.html"><span class="doc">compute</span></a> command</li>
<li>centro/atom = style name of this compute command</li>
<li>lattice = <em>fcc</em> or <em>bcc</em> or N = # of neighbors per atom to include</li>
<li>ID, group-ID are documented in <a class="reference internal" href="compute.html"><span class="doc">compute</span></a> command
centro/atom = style name of this compute command
lattice = <em>fcc</em> or <em>bcc</em> or N = # of neighbors per atom to include</li>
<li>zero or more keyword/value pairs may be appended</li>
<li>keyword = <em>axes</em></li>
</ul>
<pre class="literal-block">
<em>axes</em> value = <em>no</em> or <em>yes</em>
<em>no</em> = do not calulate 3 symmetry axes
<em>yes</em> = calulate 3 symmetry axes
</pre>
</div>
<div class="section" id="examples">
<h2>Examples</h2>
@ -149,17 +156,18 @@
<div class="section" id="description">
<h2>Description</h2>
<p>Define a computation that calculates the centro-symmetry parameter for
each atom in the group. In solid-state systems the centro-symmetry
parameter is a useful measure of the local lattice disorder around an
atom and can be used to characterize whether the atom is part of a
perfect lattice, a local defect (e.g. a dislocation or stacking
fault), or at a surface.</p>
each atom in the group, for either FCC or BCC lattices, depending on
the choice of the <em>lattice</em> argument. In solid-state systems the
centro-symmetry parameter is a useful measure of the local lattice
disorder around an atom and can be used to characterize whether the
atom is part of a perfect lattice, a local defect (e.g. a dislocation
or stacking fault), or at a surface.</p>
<p>The value of the centro-symmetry parameter will be 0.0 for atoms not
in the specified compute group.</p>
<p>This parameter is computed using the following formula from
<a class="reference internal" href="#kelchner"><span class="std std-ref">(Kelchner)</span></a></p>
<img alt="_images/centro_symmetry.jpg" class="align-center" src="_images/centro_symmetry.jpg" />
<p>where the <em>N</em> nearest neighbors or each atom are identified and Ri and
<p>where the <em>N</em> nearest neighbors of each atom are identified and Ri and
Ri+N/2 are vectors from the central atom to a particular pair of
nearest neighbors. There are N*(N-1)/2 possible neighbor pairs that
can contribute to this formula. The quantity in the sum is computed
@ -179,6 +187,20 @@ larger positive value. An atom at a surface will have a large
positive parameter. If the atom does not have <em>N</em> neighbors (within
the potential cutoff), then its centro-symmetry parameter is set to
0.0.</p>
<p>If the keyword <em>axes</em> has the setting <em>yes</em>, then this compute also
estimates three symmetry axes for each atom&#8217;s local neighborhood. The
first two of these are the vectors joining the two pairs of neighbor
atoms with smallest contributions to the centrosymmetry parameter,
i.e. the two most symmetric pairs of atoms. The third vector is
normal to the first two by the right-hand rule. All three vectors are
normalized to unit length. For FCC crystals, the first two vectors
will lie along a &lt;110&gt; direction, while the third vector will lie
along either a &lt;100&gt; or &lt;111&gt; direction. For HCP crystals, the first
two vectors will lie along &lt;1000&gt; directions, while the third vector
will lie along &lt;0001&gt;. This provides a simple way to measure local
orientation in HCP structures. In general, the <em>axes</em> keyword can be
used to estimate the orientation of symmetry axes in the neighborhood
of any atom.</p>
<p>Only atoms within the cutoff of the pairwise neighbor list are
considered as possible neighbors. Atoms not in the compute group are
included in the <em>N</em> neighbors used in this calculation.</p>
@ -188,15 +210,21 @@ is dumped). Thus it can be inefficient to compute/dump this quantity
too frequently or to have multiple compute/dump commands, each with a
<em>centro/atom</em> style.</p>
<p><strong>Output info:</strong></p>
<p>This compute calculates a per-atom vector, which can be accessed by
any command that uses per-atom values from a compute as input. See
<a class="reference internal" href="Section_howto.html#howto-15"><span class="std std-ref">Section_howto 15</span></a> for an overview of
LAMMPS output options.</p>
<p>The per-atom vector values are unitless values &gt;= 0.0. Their
magnitude depends on the lattice style due to the number of
contibuting neighbor pairs in the summation in the formula above. And
it depends on the local defects surrounding the central atom, as
described above.</p>
<p>By default, this compute calculates the centrosymmetry value for each
atom as a per-atom vector, which can be accessed by any command that
uses per-atom values from a compute as input. See <a class="reference internal" href="Section_howto.html#howto-15"><span class="std std-ref">Section_howto 15</span></a> for an overview of LAMMPS output
options.</p>
<p>If the <em>axes</em> keyword setting is <em>yes</em>, then a per-atom array is
calculated. The first column is the centrosymmetry parameter. The
next three columns are the x, y, and z components of the first
symmetry axis, followed by the second, and third symmetry axes in
columns 5-7 and 8-10.</p>
<p>The centrosymmetry values are unitless values &gt;= 0.0. Their magnitude
depends on the lattice style due to the number of contibuting neighbor
pairs in the summation in the formula above. And it depends on the
local defects surrounding the central atom, as described above. For
the <em>axes yes</em> case, the vector components are also unitless, since
they represent spatial directions.</p>
<p>Here are typical centro-symmetry values, from a a nanoindentation
simulation into gold (FCC). These were provided by Jon Zimmerman
(Sandia):</p>
@ -226,7 +254,10 @@ of 12.</p>
<div class="section" id="related-commands">
<h2>Related commands</h2>
<p><a class="reference internal" href="compute_cna_atom.html"><span class="doc">compute cna/atom</span></a></p>
<p><strong>Default:</strong> none</p>
</div>
<div class="section" id="default">
<h2>Default</h2>
<p>The default value for the optional keyword is axes = no.</p>
<hr class="docutils" />
<p id="kelchner"><strong>(Kelchner)</strong> Kelchner, Plimpton, Hamilton, Phys Rev B, 58, 11085 (1998).</p>
</div>

12
doc/html/fix_nh.html
View File

@ -198,7 +198,9 @@ keyword = <em>temp</em> or <em>iso</em> or <em>aniso</em> or <em>tri</em> or <em
<em>flip</em> value = <em>yes</em> or <em>no</em> = allow or disallow box flips when it becomes highly skewed
<em>fixedpoint</em> values = x y z
x,y,z = perform barostat dilation/contraction around this point (distance units)
<em>update</em> value = <em>dipole</em> update dipole orientation (only for sphere variants)
<em>update</em> value = <em>dipole</em> or <em>dipole/dlm</em>
dipole = update dipole orientation (only for sphere variants)
dipole/dlm = use DLM integrator to update dipole orientation (only for sphere variants)
</pre>
</div>
<div class="section" id="examples">
@ -451,6 +453,12 @@ orientation of the dipole moment of each particle is also updated
during the time integration. This option should be used for models
where a dipole moment is assigned to finite-size particles,
e.g. spheroids via use of the <a class="reference internal" href="atom_style.html"><span class="doc">atom_style hybrid sphere dipole</span></a> command.</p>
<p>The default dipole orientation integrator can be changed to the
Dullweber-Leimkuhler-McLachlan integration scheme
<a class="reference internal" href="#nh-dullweber"><span class="std std-ref">(Dullweber)</span></a> when using <em>update</em> with the value
<em>dipole/dlm</em>. This integrator is symplectic and time-reversible,
giving better energy conservation and allows slightly longer timesteps
at only a small additional computational cost.</p>
<hr class="docutils" />
<div class="admonition note">
<p class="first admonition-title">Note</p>
@ -715,6 +723,8 @@ not coupled to barostat, otherwise no.</p>
<p id="nh-tuckerman"><strong>(Tuckerman)</strong> Tuckerman, Alejandre, Lopez-Rendon, Jochim, and
Martyna, J Phys A: Math Gen, 39, 5629 (2006).</p>
<p id="nh-shinoda"><strong>(Shinoda)</strong> Shinoda, Shiga, and Mikami, Phys Rev B, 69, 134103 (2004).</p>
<p id="nh-dullweber"><strong>(Dullweber)</strong> Dullweber, Leimkuhler and McLachlan, J Chem Phys, 107,
5840 (1997).</p>
</div>
</div>

17
doc/html/fix_nve_sphere.html
View File

@ -138,17 +138,19 @@
<li>ID, group-ID are documented in <a class="reference internal" href="fix.html"><span class="doc">fix</span></a> command</li>
<li>nve/sphere = style name of this fix command</li>
<li>zero or more keyword/value pairs may be appended</li>
<li>keyword = <em>update</em></li>
</ul>
<pre class="literal-block">
<em>update</em> value = <em>dipole</em>
dipole = update orientation of dipole moment during integration
keyword = <em>update</em>
<em>update</em> value = <em>dipole</em> or <em>dipole/dlm</em>
dipole = update orientation of dipole moment during integration
dipole/dlm = use DLM integrator to update dipole orientation
</pre>
</div>
<div class="section" id="examples">
<h2>Examples</h2>
<div class="highlight-default"><div class="highlight"><pre><span></span><span class="n">fix</span> <span class="mi">1</span> <span class="nb">all</span> <span class="n">nve</span><span class="o">/</span><span class="n">sphere</span>
<span class="n">fix</span> <span class="mi">1</span> <span class="nb">all</span> <span class="n">nve</span><span class="o">/</span><span class="n">sphere</span> <span class="n">update</span> <span class="n">dipole</span>
<span class="n">fix</span> <span class="mi">1</span> <span class="nb">all</span> <span class="n">nve</span><span class="o">/</span><span class="n">sphere</span> <span class="n">update</span> <span class="n">dipole</span><span class="o">/</span><span class="n">dlm</span>
</pre></div>
</div>
</div>
@ -165,6 +167,12 @@ orientation of the dipole moment of each particle is also updated
during the time integration. This option should be used for models
where a dipole moment is assigned to finite-size particles,
e.g. spheroids via use of the <a class="reference internal" href="atom_style.html"><span class="doc">atom_style hybrid sphere dipole</span></a> command.</p>
<p>The default dipole orientation integrator can be changed to the
Dullweber-Leimkuhler-McLachlan integration scheme
<a class="reference internal" href="fix_nh.html#nh-dullweber"><span class="std std-ref">(Dullweber)</span></a> when using <em>update</em> with the value
<em>dipole/dlm</em>. This integrator is symplectic and time-reversible,
giving better energy conservation and allows slightly longer timesteps
at only a small additional computational cost.</p>
<hr class="docutils" />
<p>Styles with a <em>gpu</em>, <em>intel</em>, <em>kk</em>, <em>omp</em>, or <em>opt</em> suffix are
functionally the same as the corresponding style without the suffix.
@ -205,6 +213,9 @@ be point particles.</p>
<h2>Related commands</h2>
<p><a class="reference internal" href="fix_nve.html"><span class="doc">fix nve</span></a>, <a class="reference internal" href="fix_nve_asphere.html"><span class="doc">fix nve/asphere</span></a></p>
<p><strong>Default:</strong> none</p>
<hr class="docutils" />
<p id="nve-dullweber"><strong>(Dullweber)</strong> Dullweber, Leimkuhler and McLachlan, J Chem Phys, 107,
5840 (1997).</p>
</div>
</div>

2
doc/html/searchindex.js
View File

File diff suppressed because one or more lines are too long

4
doc/src/Manual.txt
View File

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

68
doc/src/compute_centro_atom.txt
View File

@ -10,11 +10,17 @@ compute centro/atom command :h3
[Syntax:]
compute ID group-ID centro/atom lattice :pre
compute ID group-ID centro/atom lattice keyword value ... :pre
ID, group-ID are documented in "compute"_compute.html command
centro/atom = style name of this compute command
lattice = {fcc} or {bcc} or N = # of neighbors per atom to include :ul
lattice = {fcc} or {bcc} or N = # of neighbors per atom to include :l
zero or more keyword/value pairs may be appended :l
keyword = {axes} :l
{axes} value = {no} or {yes}
{no} = do not calulate 3 symmetry axes
{yes} = calulate 3 symmetry axes :pre
:ule
[Examples:]
@ -24,11 +30,12 @@ compute 1 all centro/atom 8 :pre
[Description:]
Define a computation that calculates the centro-symmetry parameter for
each atom in the group. In solid-state systems the centro-symmetry
parameter is a useful measure of the local lattice disorder around an
atom and can be used to characterize whether the atom is part of a
perfect lattice, a local defect (e.g. a dislocation or stacking
fault), or at a surface.
each atom in the group, for either FCC or BCC lattices, depending on
the choice of the {lattice} argument. In solid-state systems the
centro-symmetry parameter is a useful measure of the local lattice
disorder around an atom and can be used to characterize whether the
atom is part of a perfect lattice, a local defect (e.g. a dislocation
or stacking fault), or at a surface.
The value of the centro-symmetry parameter will be 0.0 for atoms not
in the specified compute group.
@ -38,7 +45,7 @@ This parameter is computed using the following formula from
:c,image(Eqs/centro_symmetry.jpg)
where the {N} nearest neighbors or each atom are identified and Ri and
where the {N} nearest neighbors of each atom are identified and Ri and
Ri+N/2 are vectors from the central atom to a particular pair of
nearest neighbors. There are N*(N-1)/2 possible neighbor pairs that
can contribute to this formula. The quantity in the sum is computed
@ -61,6 +68,21 @@ positive parameter. If the atom does not have {N} neighbors (within
the potential cutoff), then its centro-symmetry parameter is set to
0.0.
If the keyword {axes} has the setting {yes}, then this compute also
estimates three symmetry axes for each atom's local neighborhood. The
first two of these are the vectors joining the two pairs of neighbor
atoms with smallest contributions to the centrosymmetry parameter,
i.e. the two most symmetric pairs of atoms. The third vector is
normal to the first two by the right-hand rule. All three vectors are
normalized to unit length. For FCC crystals, the first two vectors
will lie along a <110> direction, while the third vector will lie
along either a <100> or <111> direction. For HCP crystals, the first
two vectors will lie along <1000> directions, while the third vector
will lie along <0001>. This provides a simple way to measure local
orientation in HCP structures. In general, the {axes} keyword can be
used to estimate the orientation of symmetry axes in the neighborhood
of any atom.
Only atoms within the cutoff of the pairwise neighbor list are
considered as possible neighbors. Atoms not in the compute group are
included in the {N} neighbors used in this calculation.
@ -73,16 +95,24 @@ too frequently or to have multiple compute/dump commands, each with a
[Output info:]
This compute calculates a per-atom vector, which can be accessed by
any command that uses per-atom values from a compute as input. See
"Section_howto 15"_Section_howto.html#howto_15 for an overview of
LAMMPS output options.
By default, this compute calculates the centrosymmetry value for each
atom as a per-atom vector, which can be accessed by any command that
uses per-atom values from a compute as input. See "Section_howto
15"_Section_howto.html#howto_15 for an overview of LAMMPS output
options.
The per-atom vector values are unitless values >= 0.0. Their
magnitude depends on the lattice style due to the number of
contibuting neighbor pairs in the summation in the formula above. And
it depends on the local defects surrounding the central atom, as
described above.
If the {axes} keyword setting is {yes}, then a per-atom array is
calculated. The first column is the centrosymmetry parameter. The
next three columns are the x, y, and z components of the first
symmetry axis, followed by the second, and third symmetry axes in
columns 5-7 and 8-10.
The centrosymmetry values are unitless values >= 0.0. Their magnitude
depends on the lattice style due to the number of contibuting neighbor
pairs in the summation in the formula above. And it depends on the
local defects surrounding the central atom, as described above. For
the {axes yes} case, the vector components are also unitless, since
they represent spatial directions.
Here are typical centro-symmetry values, from a a nanoindentation
simulation into gold (FCC). These were provided by Jon Zimmerman
@ -111,7 +141,9 @@ of 12.
"compute cna/atom"_compute_cna_atom.html
[Default:] none
[Default:]
The default value for the optional keyword is axes = no.
:line

16
doc/src/fix_nh.txt
View File

@ -56,8 +56,9 @@ keyword = {temp} or {iso} or {aniso} or {tri} or {x} or {y} or {z} or {xy} or {y
{flip} value = {yes} or {no} = allow or disallow box flips when it becomes highly skewed
{fixedpoint} values = x y z
x,y,z = perform barostat dilation/contraction around this point (distance units)
{update} value = {dipole} update dipole orientation (only for sphere variants) :pre
{update} value = {dipole} or {dipole/dlm}
dipole = update dipole orientation (only for sphere variants)
dipole/dlm = use DLM integrator to update dipole orientation (only for sphere variants) :pre
:ule
[Examples:]
@ -335,6 +336,13 @@ where a dipole moment is assigned to finite-size particles,
e.g. spheroids via use of the "atom_style hybrid sphere
dipole"_atom_style.html command.
The default dipole orientation integrator can be changed to the
Dullweber-Leimkuhler-McLachlan integration scheme
"(Dullweber)"_#nh-Dullweber when using {update} with the value
{dipole/dlm}. This integrator is symplectic and time-reversible,
giving better energy conservation and allows slightly longer timesteps
at only a small additional computational cost.
:line
NOTE: Using a barostat coupled to tilt dimensions {xy}, {xz}, {yz} can
@ -638,3 +646,7 @@ Martyna, J Phys A: Math Gen, 39, 5629 (2006).
:link(nh-Shinoda)
[(Shinoda)] Shinoda, Shiga, and Mikami, Phys Rev B, 69, 134103 (2004).
:link(nh-Dullweber)
[(Dullweber)] Dullweber, Leimkuhler and McLachlan, J Chem Phys, 107,
5840 (1997).

23
doc/src/fix_nve_sphere.txt
View File

@ -16,15 +16,17 @@ fix ID group-ID nve/sphere :pre
ID, group-ID are documented in "fix"_fix.html command :ulb,l
nve/sphere = style name of this fix command :l
zero or more keyword/value pairs may be appended :l
keyword = {update} :l
{update} value = {dipole}
dipole = update orientation of dipole moment during integration :pre
keyword = {update}
{update} value = {dipole} or {dipole/dlm}
dipole = update orientation of dipole moment during integration
dipole/dlm = use DLM integrator to update dipole orientation :pre
:ule
[Examples:]
fix 1 all nve/sphere
fix 1 all nve/sphere update dipole :pre
fix 1 all nve/sphere update dipole
fix 1 all nve/sphere update dipole/dlm :pre
[Description:]
@ -43,6 +45,13 @@ where a dipole moment is assigned to finite-size particles,
e.g. spheroids via use of the "atom_style hybrid sphere
dipole"_atom_style.html command.
The default dipole orientation integrator can be changed to the
Dullweber-Leimkuhler-McLachlan integration scheme
"(Dullweber)"_#nh-Dullweber when using {update} with the value
{dipole/dlm}. This integrator is symplectic and time-reversible,
giving better energy conservation and allows slightly longer timesteps
at only a small additional computational cost.
:line
Styles with a {gpu}, {intel}, {kk}, {omp}, or {opt} suffix are
@ -94,3 +103,9 @@ be point particles.
"fix nve"_fix_nve.html, "fix nve/asphere"_fix_nve_asphere.html
[Default:] none
:line
:link(nve-Dullweber)
[(Dullweber)] Dullweber, Leimkuhler and McLachlan, J Chem Phys, 107,
5840 (1997).

4
doc/src/min_style.txt
View File

@ -69,8 +69,8 @@ defined via the "timestep"_timestep.html command. Often they will
converge more quickly if you use a timestep about 10x larger than you
would normally use for dynamics simulations.
NOTE: The {quickmin} and {fire} styles do not yet support the use of
the "fix box/relax"_fix_box_relax.html command or minimizations
NOTE: The {quickmin}, {fire} and {hftn} styles do not yet support the
use of the "fix box/relax"_fix_box_relax.html command or minimizations
involving the electron radius in "eFF"_pair_eff.html models.
[Restrictions:] none

29
src/KOKKOS/comm_kokkos.cpp
View File

@ -43,27 +43,18 @@ enum{SINGLE,MULTI};
CommKokkos::CommKokkos(LAMMPS *lmp) : CommBrick(lmp)
{
sendlist = NULL; // need to free this since parent allocated?
if (sendlist) for (int i = 0; i < maxswap; i++) memory->destroy(sendlist[i]);
memory->sfree(sendlist);
sendlist = NULL;
k_sendlist = ArrayTypes<LMPDeviceType>::tdual_int_2d();
// error check for disallow of OpenMP threads?
// initialize comm buffers & exchange memory
// maxsend = BUFMIN;
// k_buf_send = ArrayTypes<LMPDeviceType>::
// tdual_xfloat_2d("comm:k_buf_send",(maxsend+BUFEXTRA+5)/6,6);
// buf_send = k_buf_send.view<LMPHostType>().ptr_on_device();
maxsend = 0;
memory->destroy(buf_send);
buf_send = NULL;
// maxrecv = BUFMIN;
// k_buf_recv = ArrayTypes<LMPDeviceType>::
// tdual_xfloat_2d("comm:k_buf_recv",(maxrecv+5)/6,6);
// buf_recv = k_buf_recv.view<LMPHostType>().ptr_on_device();
maxrecv = 0;
memory->destroy(buf_recv);
buf_recv = NULL;
k_exchange_sendlist = ArrayTypes<LMPDeviceType>::
@ -73,8 +64,8 @@ CommKokkos::CommKokkos(LAMMPS *lmp) : CommBrick(lmp)
k_count = ArrayTypes<LMPDeviceType>::tdual_int_1d("comm:k_count",1);
k_sendflag = ArrayTypes<LMPDeviceType>::tdual_int_1d("comm:k_sendflag",100);
// next line is bogus?
memory->destroy(maxsendlist);
maxsendlist = NULL;
memory->create(maxsendlist,maxswap,"comm:maxsendlist");
for (int i = 0; i < maxswap; i++) {
maxsendlist[i] = BUFMIN;
@ -87,14 +78,20 @@ CommKokkos::CommKokkos(LAMMPS *lmp) : CommBrick(lmp)
CommKokkos::~CommKokkos()
{
memory->destroy_kokkos(k_sendlist,sendlist);
sendlist = NULL;
memory->destroy_kokkos(k_buf_send,buf_send);
buf_send = NULL;
memory->destroy_kokkos(k_buf_recv,buf_recv);
buf_recv = NULL;
}
/* ---------------------------------------------------------------------- */
void CommKokkos::init()
{
grow_send_kokkos(maxsend+bufextra,0,Host);
grow_recv_kokkos(maxrecv,Host);
atomKK = (AtomKokkos *) atom;
exchange_comm_classic = lmp->kokkos->exchange_comm_classic;
forward_comm_classic = lmp->kokkos->forward_comm_classic;

80
src/KSPACE/pair_lj_charmm_coul_long.cpp
View File

@ -172,9 +172,9 @@ void PairLJCharmmCoulLong::compute(int eflag, int vflag)
forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]);
if (rsq > cut_lj_innersq) {
switch1 = (cut_ljsq-rsq) * (cut_ljsq-rsq) *
(cut_ljsq + 2.0*rsq - 3.0*cut_lj_innersq) / denom_lj;
(cut_ljsq + 2.0*rsq - 3.0*cut_lj_innersq) * denom_lj_inv;
switch2 = 12.0*rsq * (cut_ljsq-rsq) *
(rsq-cut_lj_innersq) / denom_lj;
(rsq-cut_lj_innersq) * denom_lj_inv;
philj = r6inv * (lj3[itype][jtype]*r6inv - lj4[itype][jtype]);
forcelj = forcelj*switch1 + philj*switch2;
}
@ -206,7 +206,7 @@ void PairLJCharmmCoulLong::compute(int eflag, int vflag)
evdwl = r6inv*(lj3[itype][jtype]*r6inv-lj4[itype][jtype]);
if (rsq > cut_lj_innersq) {
switch1 = (cut_ljsq-rsq) * (cut_ljsq-rsq) *
(cut_ljsq + 2.0*rsq - 3.0*cut_lj_innersq) / denom_lj;
(cut_ljsq + 2.0*rsq - 3.0*cut_lj_innersq) * denom_lj_inv;
evdwl *= switch1;
}
evdwl *= factor_lj;
@ -247,13 +247,6 @@ void PairLJCharmmCoulLong::compute_inner()
numneigh = listinner->numneigh;
firstneigh = listinner->firstneigh;
double cut_out_on = cut_respa[0];
double cut_out_off = cut_respa[1];
double cut_out_diff = cut_out_off - cut_out_on;
double cut_out_on_sq = cut_out_on*cut_out_on;
double cut_out_off_sq = cut_out_off*cut_out_off;
// loop over neighbors of my atoms
for (ii = 0; ii < inum; ii++) {
@ -289,7 +282,7 @@ void PairLJCharmmCoulLong::compute_inner()
fpair = (forcecoul + factor_lj*forcelj) * r2inv;
if (rsq > cut_out_on_sq) {
rsw = (sqrt(rsq) - cut_out_on)/cut_out_diff;
rsw = (sqrt(rsq) - cut_out_on)*cut_out_diff_inv;
fpair *= 1.0 + rsw*rsw*(2.0*rsw-3.0);
}
@ -332,18 +325,6 @@ void PairLJCharmmCoulLong::compute_middle()
numneigh = listmiddle->numneigh;
firstneigh = listmiddle->firstneigh;
double cut_in_off = cut_respa[0];
double cut_in_on = cut_respa[1];
double cut_out_on = cut_respa[2];
double cut_out_off = cut_respa[3];
double cut_in_diff = cut_in_on - cut_in_off;
double cut_out_diff = cut_out_off - cut_out_on;
double cut_in_off_sq = cut_in_off*cut_in_off;
double cut_in_on_sq = cut_in_on*cut_in_on;
double cut_out_on_sq = cut_out_on*cut_out_on;
double cut_out_off_sq = cut_out_off*cut_out_off;
// loop over neighbors of my atoms
for (ii = 0; ii < inum; ii++) {
@ -378,20 +359,20 @@ void PairLJCharmmCoulLong::compute_middle()
forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]);
if (rsq > cut_lj_innersq) {
switch1 = (cut_ljsq-rsq) * (cut_ljsq-rsq) *
(cut_ljsq + 2.0*rsq - 3.0*cut_lj_innersq) / denom_lj;
(cut_ljsq + 2.0*rsq - 3.0*cut_lj_innersq) * denom_lj_inv;
switch2 = 12.0*rsq * (cut_ljsq-rsq) *
(rsq-cut_lj_innersq) / denom_lj;
(rsq-cut_lj_innersq) * denom_lj_inv;
philj = r6inv * (lj3[itype][jtype]*r6inv - lj4[itype][jtype]);
forcelj = forcelj*switch1 + philj*switch2;
}
fpair = (forcecoul + factor_lj*forcelj) * r2inv;
if (rsq < cut_in_on_sq) {
rsw = (sqrt(rsq) - cut_in_off)/cut_in_diff;
rsw = (sqrt(rsq) - cut_in_off)*cut_in_diff_inv;
fpair *= rsw*rsw*(3.0 - 2.0*rsw);
}
if (rsq > cut_out_on_sq) {
rsw = (sqrt(rsq) - cut_out_on)/cut_out_diff;
rsw = (sqrt(rsq) - cut_out_on)*cut_out_diff_inv;
fpair *= 1.0 + rsw*rsw*(2.0*rsw - 3.0);
}
@ -441,13 +422,6 @@ void PairLJCharmmCoulLong::compute_outer(int eflag, int vflag)
numneigh = listouter->numneigh;
firstneigh = listouter->firstneigh;
double cut_in_off = cut_respa[2];
double cut_in_on = cut_respa[3];
double cut_in_diff = cut_in_on - cut_in_off;
double cut_in_off_sq = cut_in_off*cut_in_off;
double cut_in_on_sq = cut_in_on*cut_in_on;
// loop over neighbors of my atoms
for (ii = 0; ii < inum; ii++) {
@ -518,9 +492,9 @@ void PairLJCharmmCoulLong::compute_outer(int eflag, int vflag)
forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]);
if (rsq > cut_lj_innersq) {
switch1 = (cut_ljsq-rsq) * (cut_ljsq-rsq) *
(cut_ljsq + 2.0*rsq - 3.0*cut_lj_innersq) / denom_lj;
(cut_ljsq + 2.0*rsq - 3.0*cut_lj_innersq) * denom_lj_inv;
switch2 = 12.0*rsq * (cut_ljsq-rsq) *
(rsq-cut_lj_innersq) / denom_lj;
(rsq-cut_lj_innersq) * denom_lj_inv;
philj = r6inv * (lj3[itype][jtype]*r6inv - lj4[itype][jtype]);
forcelj = forcelj*switch1 + philj*switch2;
}
@ -562,7 +536,7 @@ void PairLJCharmmCoulLong::compute_outer(int eflag, int vflag)
evdwl = r6inv*(lj3[itype][jtype]*r6inv-lj4[itype][jtype]);
if (rsq > cut_lj_innersq) {
switch1 = (cut_ljsq-rsq) * (cut_ljsq-rsq) *
(cut_ljsq + 2.0*rsq - 3.0*cut_lj_innersq) / denom_lj;
(cut_ljsq + 2.0*rsq - 3.0*cut_lj_innersq) * denom_lj_inv;
evdwl *= switch1;
}
evdwl *= factor_lj;
@ -590,9 +564,9 @@ void PairLJCharmmCoulLong::compute_outer(int eflag, int vflag)
forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]);
if (rsq > cut_lj_innersq) {
switch1 = (cut_ljsq-rsq) * (cut_ljsq-rsq) *
(cut_ljsq + 2.0*rsq - 3.0*cut_lj_innersq) / denom_lj;
(cut_ljsq + 2.0*rsq - 3.0*cut_lj_innersq) * denom_lj_inv;
switch2 = 12.0*rsq * (cut_ljsq-rsq) *
(rsq-cut_lj_innersq) / denom_lj;
(rsq-cut_lj_innersq) * denom_lj_inv;
philj = r6inv * (lj3[itype][jtype]*r6inv - lj4[itype][jtype]);
forcelj = forcelj*switch1 + philj*switch2;
}
@ -601,9 +575,9 @@ void PairLJCharmmCoulLong::compute_outer(int eflag, int vflag)
forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]);
if (rsq > cut_lj_innersq) {
switch1 = (cut_ljsq-rsq) * (cut_ljsq-rsq) *
(cut_ljsq + 2.0*rsq - 3.0*cut_lj_innersq) / denom_lj;
(cut_ljsq + 2.0*rsq - 3.0*cut_lj_innersq) * denom_lj_inv;
switch2 = 12.0*rsq * (cut_ljsq-rsq) *
(rsq-cut_lj_innersq) / denom_lj;
(rsq-cut_lj_innersq) * denom_lj_inv;
philj = r6inv * (lj3[itype][jtype]*r6inv - lj4[itype][jtype]);
forcelj = forcelj*switch1 + philj*switch2;
}
@ -759,14 +733,28 @@ void PairLJCharmmCoulLong::init_style()
cut_coulsq = cut_coul * cut_coul;
cut_bothsq = MAX(cut_ljsq,cut_coulsq);
denom_lj = (cut_ljsq-cut_lj_innersq) * (cut_ljsq-cut_lj_innersq) *
(cut_ljsq-cut_lj_innersq);
denom_lj = ( (cut_ljsq-cut_lj_innersq) * (cut_ljsq-cut_lj_innersq) *
(cut_ljsq-cut_lj_innersq) );
denom_lj_inv = 1.0 / denom_lj;
// set & error check interior rRESPA cutoffs
if (strstr(update->integrate_style,"respa") &&
((Respa *) update->integrate)->level_inner >= 0) {
cut_respa = ((Respa *) update->integrate)->cutoff;
cut_in_off = cut_respa[0];
cut_in_on = cut_respa[1];
cut_out_on = cut_respa[2];
cut_out_off = cut_respa[3];
cut_in_diff = cut_in_on - cut_in_off;
cut_out_diff = cut_out_off - cut_out_on;
cut_in_diff_inv = 1.0 / (cut_in_diff);
cut_out_diff_inv = 1.0 / (cut_out_diff);
cut_in_off_sq = cut_in_off*cut_in_off;
cut_in_on_sq = cut_in_on*cut_in_on;
cut_out_on_sq = cut_out_on*cut_out_on;
cut_out_off_sq = cut_out_off*cut_out_off;
if (MIN(cut_lj,cut_coul) < cut_respa[3])
error->all(FLERR,"Pair cutoff < Respa interior cutoff");
if (cut_lj_inner < cut_respa[1])
@ -992,9 +980,9 @@ double PairLJCharmmCoulLong::single(int i, int j, int itype, int jtype,
forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]);
if (rsq > cut_lj_innersq) {
switch1 = (cut_ljsq-rsq) * (cut_ljsq-rsq) *
(cut_ljsq + 2.0*rsq - 3.0*cut_lj_innersq) / denom_lj;
(cut_ljsq + 2.0*rsq - 3.0*cut_lj_innersq) * denom_lj_inv;
switch2 = 12.0*rsq * (cut_ljsq-rsq) *
(rsq-cut_lj_innersq) / denom_lj;
(rsq-cut_lj_innersq) * denom_lj_inv;
philj = r6inv * (lj3[itype][jtype]*r6inv - lj4[itype][jtype]);
forcelj = forcelj*switch1 + philj*switch2;
}
@ -1017,7 +1005,7 @@ double PairLJCharmmCoulLong::single(int i, int j, int itype, int jtype,
philj = r6inv*(lj3[itype][jtype]*r6inv-lj4[itype][jtype]);
if (rsq > cut_lj_innersq) {
switch1 = (cut_ljsq-rsq) * (cut_ljsq-rsq) *
(cut_ljsq + 2.0*rsq - 3.0*cut_lj_innersq) / denom_lj;
(cut_ljsq + 2.0*rsq - 3.0*cut_lj_innersq) * denom_lj_inv;
philj *= switch1;
}
eng += factor_lj*philj;

6
src/KSPACE/pair_lj_charmm_coul_long.h
View File

@ -54,7 +54,11 @@ class PairLJCharmmCoulLong : public Pair {
double cut_lj_innersq,cut_ljsq;
double cut_coul,cut_coulsq;
double cut_bothsq;
double denom_lj;
double cut_in_off, cut_in_on, cut_out_off, cut_out_on;
double cut_in_diff, cut_out_diff;
double cut_in_diff_inv, cut_out_diff_inv;
double cut_in_off_sq, cut_in_on_sq, cut_out_off_sq, cut_out_on_sq;
double denom_lj, denom_lj_inv;
double **epsilon,**sigma,**eps14,**sigma14;
double **lj1,**lj2,**lj3,**lj4,**offset;
double **lj14_1,**lj14_2,**lj14_3,**lj14_4;

26
src/KSPACE/pair_lj_charmm_coul_msm.cpp
View File

@ -161,9 +161,9 @@ void PairLJCharmmCoulMSM::compute(int eflag, int vflag)
forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]);
if (rsq > cut_lj_innersq) {
switch1 = (cut_ljsq-rsq) * (cut_ljsq-rsq) *
(cut_ljsq + 2.0*rsq - 3.0*cut_lj_innersq) / denom_lj;
(cut_ljsq + 2.0*rsq - 3.0*cut_lj_innersq) * denom_lj_inv;
switch2 = 12.0*rsq * (cut_ljsq-rsq) *
(rsq-cut_lj_innersq) / denom_lj;
(rsq-cut_lj_innersq) * denom_lj_inv;
philj = r6inv * (lj3[itype][jtype]*r6inv - lj4[itype][jtype]);
forcelj = forcelj*switch1 + philj*switch2;
}
@ -221,7 +221,7 @@ void PairLJCharmmCoulMSM::compute(int eflag, int vflag)
evdwl = r6inv*(lj3[itype][jtype]*r6inv-lj4[itype][jtype]);
if (rsq > cut_lj_innersq) {
switch1 = (cut_ljsq-rsq) * (cut_ljsq-rsq) *
(cut_ljsq + 2.0*rsq - 3.0*cut_lj_innersq) / denom_lj;
(cut_ljsq + 2.0*rsq - 3.0*cut_lj_innersq) * denom_lj_inv;
evdwl *= switch1;
}
evdwl *= factor_lj;
@ -358,9 +358,9 @@ void PairLJCharmmCoulMSM::compute_outer(int eflag, int vflag)
forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]);
if (rsq > cut_lj_innersq) {
switch1 = (cut_ljsq-rsq) * (cut_ljsq-rsq) *
(cut_ljsq + 2.0*rsq - 3.0*cut_lj_innersq) / denom_lj;
(cut_ljsq + 2.0*rsq - 3.0*cut_lj_innersq) * denom_lj_inv;
switch2 = 12.0*rsq * (cut_ljsq-rsq) *
(rsq-cut_lj_innersq) / denom_lj;
(rsq-cut_lj_innersq) * denom_lj_inv;
philj = r6inv * (lj3[itype][jtype]*r6inv - lj4[itype][jtype]);
forcelj = forcelj*switch1 + philj*switch2;
}
@ -402,7 +402,7 @@ void PairLJCharmmCoulMSM::compute_outer(int eflag, int vflag)
evdwl = r6inv*(lj3[itype][jtype]*r6inv-lj4[itype][jtype]);
if (rsq > cut_lj_innersq) {
switch1 = (cut_ljsq-rsq) * (cut_ljsq-rsq) *
(cut_ljsq + 2.0*rsq - 3.0*cut_lj_innersq) / denom_lj;
(cut_ljsq + 2.0*rsq - 3.0*cut_lj_innersq) * denom_lj_inv;
evdwl *= switch1;
}
evdwl *= factor_lj;
@ -430,9 +430,9 @@ void PairLJCharmmCoulMSM::compute_outer(int eflag, int vflag)
forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]);
if (rsq > cut_lj_innersq) {
switch1 = (cut_ljsq-rsq) * (cut_ljsq-rsq) *
(cut_ljsq + 2.0*rsq - 3.0*cut_lj_innersq) / denom_lj;
(cut_ljsq + 2.0*rsq - 3.0*cut_lj_innersq) * denom_lj_inv;
switch2 = 12.0*rsq * (cut_ljsq-rsq) *
(rsq-cut_lj_innersq) / denom_lj;
(rsq-cut_lj_innersq) * denom_lj_inv;
philj = r6inv * (lj3[itype][jtype]*r6inv - lj4[itype][jtype]);
forcelj = forcelj*switch1 + philj*switch2;
}
@ -441,9 +441,9 @@ void PairLJCharmmCoulMSM::compute_outer(int eflag, int vflag)
forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]);
if (rsq > cut_lj_innersq) {
switch1 = (cut_ljsq-rsq) * (cut_ljsq-rsq) *
(cut_ljsq + 2.0*rsq - 3.0*cut_lj_innersq) / denom_lj;
(cut_ljsq + 2.0*rsq - 3.0*cut_lj_innersq) * denom_lj_inv;
switch2 = 12.0*rsq * (cut_ljsq-rsq) *
(rsq-cut_lj_innersq) / denom_lj;
(rsq-cut_lj_innersq) * denom_lj_inv;
philj = r6inv * (lj3[itype][jtype]*r6inv - lj4[itype][jtype]);
forcelj = forcelj*switch1 + philj*switch2;
}
@ -499,9 +499,9 @@ double PairLJCharmmCoulMSM::single(int i, int j, int itype, int jtype,
forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]);
if (rsq > cut_lj_innersq) {
switch1 = (cut_ljsq-rsq) * (cut_ljsq-rsq) *
(cut_ljsq + 2.0*rsq - 3.0*cut_lj_innersq) / denom_lj;
(cut_ljsq + 2.0*rsq - 3.0*cut_lj_innersq) * denom_lj_inv;
switch2 = 12.0*rsq * (cut_ljsq-rsq) *
(rsq-cut_lj_innersq) / denom_lj;
(rsq-cut_lj_innersq) * denom_lj_inv;
philj = r6inv * (lj3[itype][jtype]*r6inv - lj4[itype][jtype]);
forcelj = forcelj*switch1 + philj*switch2;
}
@ -524,7 +524,7 @@ double PairLJCharmmCoulMSM::single(int i, int j, int itype, int jtype,
philj = r6inv*(lj3[itype][jtype]*r6inv-lj4[itype][jtype]);
if (rsq > cut_lj_innersq) {
switch1 = (cut_ljsq-rsq) * (cut_ljsq-rsq) *
(cut_ljsq + 2.0*rsq - 3.0*cut_lj_innersq) / denom_lj;
(cut_ljsq + 2.0*rsq - 3.0*cut_lj_innersq) * denom_lj_inv;
philj *= switch1;
}
eng += factor_lj*philj;

3
src/QEQ/fix_qeq.cpp
View File

@ -282,7 +282,8 @@ void FixQEq::setup_pre_force(int vflag)
if (force->newton_pair == 0)
error->all(FLERR,"QEQ with 'newton pair off' not supported");
neighbor->build_one(list);
// should not be needed
// neighbor->build_one(list);
deallocate_storage();
allocate_storage();

9
src/RIGID/fix_shake.cpp
View File

@ -446,6 +446,15 @@ void FixShake::setup(int vflag)
else
respa = 1;
if (!respa) {
dtv = update->dt;
dtfsq = 0.5 * update->dt * update->dt * force->ftm2v;
if (!rattle) dtfsq = update->dt * update->dt * force->ftm2v;
} else {
dtv = step_respa[0];
dtf_innerhalf = 0.5 * step_respa[0] * force->ftm2v;
dtf_inner = dtf_innerhalf;
}
// correct geometry of cluster if necessary

8
src/STUBS/mpi.c
View File

@ -190,7 +190,7 @@ int MPI_Type_size(MPI_Datatype datatype, int *size)
/* ---------------------------------------------------------------------- */
int MPI_Send(void *buf, int count, MPI_Datatype datatype,
int MPI_Send(const void *buf, int count, MPI_Datatype datatype,
int dest, int tag, MPI_Comm comm)
{
printf("MPI Stub WARNING: Should not send message to self\n");
@ -199,7 +199,7 @@ int MPI_Send(void *buf, int count, MPI_Datatype datatype,
/* ---------------------------------------------------------------------- */
int MPI_Isend(void *buf, int count, MPI_Datatype datatype,
int MPI_Isend(const void *buf, int count, MPI_Datatype datatype,
int source, int tag, MPI_Comm comm, MPI_Request *request)
{
printf("MPI Stub WARNING: Should not send message to self\n");
@ -208,7 +208,7 @@ int MPI_Isend(void *buf, int count, MPI_Datatype datatype,
/* ---------------------------------------------------------------------- */
int MPI_Rsend(void *buf, int count, MPI_Datatype datatype,
int MPI_Rsend(const void *buf, int count, MPI_Datatype datatype,
int dest, int tag, MPI_Comm comm)
{
printf("MPI Stub WARNING: Should not rsend message to self\n");
@ -260,7 +260,7 @@ int MPI_Waitany(int count, MPI_Request *request, int *index,
/* ---------------------------------------------------------------------- */
int MPI_Sendrecv(void *sbuf, int scount, MPI_Datatype sdatatype,
int MPI_Sendrecv(const void *sbuf, int scount, MPI_Datatype sdatatype,
int dest, int stag, void *rbuf, int rcount,
MPI_Datatype rdatatype, int source, int rtag,
MPI_Comm comm, MPI_Status *status)

8
src/STUBS/mpi.h
View File

@ -90,11 +90,11 @@ double MPI_Wtime();
int MPI_Type_size(int, int *);
int MPI_Send(void *buf, int count, MPI_Datatype datatype,
int MPI_Send(const void *buf, int count, MPI_Datatype datatype,
int dest, int tag, MPI_Comm comm);
int MPI_Isend(void *buf, int count, MPI_Datatype datatype,
int MPI_Isend(const void *buf, int count, MPI_Datatype datatype,
int source, int tag, MPI_Comm comm, MPI_Request *request);
int MPI_Rsend(void *buf, int count, MPI_Datatype datatype,
int MPI_Rsend(const void *buf, int count, MPI_Datatype datatype,
int dest, int tag, MPI_Comm comm);
int MPI_Recv(void *buf, int count, MPI_Datatype datatype,
int source, int tag, MPI_Comm comm, MPI_Status *status);
@ -104,7 +104,7 @@ int MPI_Wait(MPI_Request *request, MPI_Status *status);
int MPI_Waitall(int n, MPI_Request *request, MPI_Status *status);
int MPI_Waitany(int count, MPI_Request *request, int *index,
MPI_Status *status);
int MPI_Sendrecv(void *sbuf, int scount, MPI_Datatype sdatatype,
int MPI_Sendrecv(const void *sbuf, int scount, MPI_Datatype sdatatype,
int dest, int stag, void *rbuf, int rcount,
MPI_Datatype rdatatype, int source, int rtag,
MPI_Comm comm, MPI_Status *status);

14
src/USER-MISC/pair_lj_sf_dipole_sf.cpp
View File

@ -133,8 +133,9 @@ void PairLJSFDipoleSF::compute(int eflag, int vflag)
if (rsq < cut_coulsq[itype][jtype]) {
rcutcoul2inv=1.0/cut_coulsq[itype][jtype];
if (qtmp != 0.0 && q[j] != 0.0) {
pre1 = qtmp*q[j]*rinv*(r2inv-1.0/cut_coulsq[itype][jtype]);
pre1 = qtmp*q[j]*rinv*(r2inv-rcutcoul2inv);
forcecoulx += pre1*delx;
forcecouly += pre1*dely;
@ -144,7 +145,6 @@ void PairLJSFDipoleSF::compute(int eflag, int vflag)
if (mu[i][3] > 0.0 && mu[j][3] > 0.0) {
r3inv = r2inv*rinv;
r5inv = r3inv*r2inv;
rcutcoul2inv=1.0/cut_coulsq[itype][jtype];
pdotp = mu[i][0]*mu[j][0] + mu[i][1]*mu[j][1] + mu[i][2]*mu[j][2];
pidotr = mu[i][0]*delx + mu[i][1]*dely + mu[i][2]*delz;
@ -156,7 +156,7 @@ void PairLJSFDipoleSF::compute(int eflag, int vflag)
aforcecouly = pre1*dely;
aforcecoulz = pre1*delz;
bfac = 1.0 - 4.0*rsq*sqrt(rsq)*rcutcoul2inv*sqrt(rcutcoul2inv) +
bfac = 1.0 - 4.0*rsq*sqrt(rsq*rcutcoul2inv)*rcutcoul2inv +
3.0*rsq*rsq*rcutcoul2inv*rcutcoul2inv;
presf = 2.0 * r2inv * pidotr * pjdotr;
bforcecoulx = bfac * (pjdotr*mu[i][0]+pidotr*mu[j][0]-presf*delx);
@ -187,10 +187,9 @@ void PairLJSFDipoleSF::compute(int eflag, int vflag)
r3inv = r2inv*rinv;
r5inv = r3inv*r2inv;
pidotr = mu[i][0]*delx + mu[i][1]*dely + mu[i][2]*delz;
rcutcoul2inv=1.0/cut_coulsq[itype][jtype];
pre1 = 3.0 * q[j] * r5inv * pidotr * (1-rsq*rcutcoul2inv);
pqfac = 1.0 - 3.0*rsq*rcutcoul2inv +
2.0*rsq*sqrt(rsq)*rcutcoul2inv*sqrt(rcutcoul2inv);
2.0*rsq*sqrt(rsq*rcutcoul2inv)*rcutcoul2inv;
pre2 = q[j] * r3inv * pqfac;
forcecoulx += pre2*mu[i][0] - pre1*delx;
@ -205,10 +204,9 @@ void PairLJSFDipoleSF::compute(int eflag, int vflag)
r3inv = r2inv*rinv;
r5inv = r3inv*r2inv;
pjdotr = mu[j][0]*delx + mu[j][1]*dely + mu[j][2]*delz;
rcutcoul2inv=1.0/cut_coulsq[itype][jtype];
pre1 = 3.0 * qtmp * r5inv * pjdotr * (1-rsq*rcutcoul2inv);
qpfac = 1.0 - 3.0*rsq*rcutcoul2inv +
2.0*rsq*sqrt(rsq)*rcutcoul2inv*sqrt(rcutcoul2inv);
2.0*rsq*sqrt(rsq*rcutcoul2inv)*rcutcoul2inv;
pre2 = qtmp * r3inv * qpfac;
forcecoulx += pre1*delx - pre2*mu[j][0];
@ -261,7 +259,7 @@ void PairLJSFDipoleSF::compute(int eflag, int vflag)
if (eflag) {
if (rsq < cut_coulsq[itype][jtype]) {
ecoul = (1.0-sqrt(rsq)/sqrt(cut_coulsq[itype][jtype]));
ecoul = (1.0-sqrt(rsq/cut_coulsq[itype][jtype]));
ecoul *= ecoul;
ecoul *= qtmp * q[j] * rinv;
if (mu[i][3] > 0.0 && mu[j][3] > 0.0)

3
src/USER-OMP/fix_nh_sphere_omp.cpp
View File

@ -55,7 +55,7 @@ void FixNHSphereOMP::init()
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit)
if (radius[i] == 0.0)
error->one(FLERR,"Fix nvt/sphere requires extended particles");
error->one(FLERR,"Fix nvt/npt/nph/sphere/omp require extended particles");
FixNHOMP::init();
}
@ -97,6 +97,7 @@ void FixNHSphereOMP::nve_v()
v[i].x += dtfm*f[i].x;
v[i].y += dtfm*f[i].y;
v[i].z += dtfm*f[i].z;
const double dtirotate = dtfrotate / (radius[i]*radius[i]*rmass[i]);
omega[i].x += dtirotate*torque[i].x;
omega[i].y += dtirotate*torque[i].y;

135
src/USER-OMP/fix_nve_sphere_omp.cpp
View File

@ -21,13 +21,17 @@
#include "respa.h"
#include "force.h"
#include "error.h"
#include "math_vector.h"
#include "math_extra.h"
using namespace LAMMPS_NS;
using namespace FixConst;
using namespace MathExtra;
#define INERTIA 0.4 // moment of inertia prefactor for sphere
enum{NONE,DIPOLE};
enum{NODLM,DLM};
/* ---------------------------------------------------------------------- */
@ -77,21 +81,128 @@ void FixNVESphereOMP::initial_integrate(int vflag)
if (extra == DIPOLE) {
double * const * const mu = atom->mu;
if (dlm == NODLM) {
#if defined(_OPENMP)
#pragma omp parallel for private(i) default(none)
#endif
for (i = 0; i < nlocal; i++) {
double g0,g1,g2,msq,scale;
if (mask[i] & groupbit) {
if (mu[i][3] > 0.0) {
g0 = mu[i][0] + dtv * (omega[i][1]*mu[i][2]-omega[i][2]*mu[i][1]);
g1 = mu[i][1] + dtv * (omega[i][2]*mu[i][0]-omega[i][0]*mu[i][2]);
g2 = mu[i][2] + dtv * (omega[i][0]*mu[i][1]-omega[i][1]*mu[i][0]);
msq = g0*g0 + g1*g1 + g2*g2;
scale = mu[i][3]/sqrt(msq);
mu[i][0] = g0*scale;
mu[i][1] = g1*scale;
mu[i][2] = g2*scale;
for (i = 0; i < nlocal; i++) {
double g0,g1,g2,msq,scale;
if (mask[i] & groupbit) {
if (mu[i][3] > 0.0) {
g0 = mu[i][0] + dtv * (omega[i][1]*mu[i][2]-omega[i][2]*mu[i][1]);
g1 = mu[i][1] + dtv * (omega[i][2]*mu[i][0]-omega[i][0]*mu[i][2]);
g2 = mu[i][2] + dtv * (omega[i][0]*mu[i][1]-omega[i][1]*mu[i][0]);
msq = g0*g0 + g1*g1 + g2*g2;
scale = mu[i][3]/sqrt(msq);
mu[i][0] = g0*scale;
mu[i][1] = g1*scale;
mu[i][2] = g2*scale;
}
}
}
} else {
#if defined(_OPENMP)
#pragma omp parallel for private(i) default(none)
#endif
// Integrate orientation following Dullweber-Leimkuhler-Maclachlan scheme
for (i = 0; i < nlocal; i++) {
vector w, w_temp, a;
matrix Q, Q_temp, R;
if (mask[i] & groupbit && mu[i][3] > 0.0) {
// Construct Q from dipole:
// Q is the rotation matrix from space frame to body frame
// i.e. v_b = Q.v_s
// Define mu to lie along the z axis in the body frame
// We take the unit dipole to avoid getting a scaling matrix
const double inv_len_mu = 1.0/mu[i][3];
a[0] = mu[i][0]*inv_len_mu;
a[1] = mu[i][1]*inv_len_mu;
a[2] = mu[i][2]*inv_len_mu;
// v = a x [0 0 1] - cross product of mu in space and body frames
// s = |v|
// c = a.[0 0 1] = a[2]
// vx = [ 0 -v[2] v[1]
// v[2] 0 -v[0]
// -v[1] v[0] 0 ]
// then
// Q = I + vx + vx^2 * (1-c)/s^2
const double s2 = a[0]*a[0] + a[1]*a[1];
if (s2 != 0.0){ // i.e. the vectors are not parallel
const double scale = (1.0 - a[2])/s2;
Q[0][0] = 1.0 - scale*a[0]*a[0]; Q[0][1] = -scale*a[0]*a[1]; Q[0][2] = -a[0];
Q[1][0] = -scale*a[0]*a[1]; Q[1][1] = 1.0 - scale*a[1]*a[1]; Q[1][2] = -a[1];
Q[2][0] = a[0]; Q[2][1] = a[1]; Q[2][2] = 1.0 - scale*(a[0]*a[0] + a[1]*a[1]);
} else { // if parallel then we just have I or -I
Q[0][0] = 1.0/a[2]; Q[0][1] = 0.0; Q[0][2] = 0.0;
Q[1][0] = 0.0; Q[1][1] = 1.0/a[2]; Q[1][2] = 0.0;
Q[2][0] = 0.0; Q[2][1] = 0.0; Q[2][2] = 1.0/a[2];
}
// Local copy of this particle's angular velocity (in space frame)
w[0] = omega[i][0]; w[1] = omega[i][1]; w[2] = omega[i][2];
// Transform omega into body frame: w_temp= Q.w
matvec(Q,w,w_temp);
// Construct rotation R1
BuildRxMatrix(R, dtf/force->ftm2v*w_temp[0]);
// Apply R1 to w: w = R.w_temp
matvec(R,w_temp,w);
// Apply R1 to Q: Q_temp = R^T.Q
transpose_times3(R,Q,Q_temp);
// Construct rotation R2
BuildRyMatrix(R, dtf/force->ftm2v*w[1]);
// Apply R2 to w: w_temp = R.w
matvec(R,w,w_temp);
// Apply R2 to Q: Q = R^T.Q_temp
transpose_times3(R,Q_temp,Q);
// Construct rotation R3
BuildRzMatrix(R, 2.0*dtf/force->ftm2v*w_temp[2]);
// Apply R3 to w: w = R.w_temp
matvec(R,w_temp,w);
// Apply R3 to Q: Q_temp = R^T.Q
transpose_times3(R,Q,Q_temp);
// Construct rotation R4
BuildRyMatrix(R, dtf/force->ftm2v*w[1]);
// Apply R4 to w: w_temp = R.w
matvec(R,w,w_temp);
// Apply R4 to Q: Q = R^T.Q_temp
transpose_times3(R,Q_temp,Q);
// Construct rotation R5
BuildRxMatrix(R, dtf/force->ftm2v*w_temp[0]);
// Apply R5 to w: w = R.w_temp
matvec(R,w_temp,w);
// Apply R5 to Q: Q_temp = R^T.Q
transpose_times3(R,Q,Q_temp);
// Transform w back into space frame w_temp = Q^T.w
transpose_matvec(Q_temp,w,w_temp);
omega[i][0] = w_temp[0]; omega[i][1] = w_temp[1]; omega[i][2] = w_temp[2];
// Set dipole according to updated Q: mu = Q^T.[0 0 1] * |mu|
mu[i][0] = Q_temp[2][0] * mu[i][3];
mu[i][1] = Q_temp[2][1] * mu[i][3];
mu[i][2] = Q_temp[2][2] * mu[i][3];
}
}
}

5
src/USER-OMP/fix_omp.cpp
View File

@ -105,6 +105,7 @@ FixOMP::FixOMP(LAMMPS *lmp, int narg, char **arg)
// print summary of settings
if (comm->me == 0) {
#if defined(_OPENMP)
const char * const nmode = _neighbor ? "multi-threaded" : "serial";
if (screen) {
@ -118,6 +119,10 @@ FixOMP::FixOMP(LAMMPS *lmp, int narg, char **arg)
fprintf(logfile,"set %d OpenMP thread(s) per MPI task\n", nthreads);
fprintf(logfile,"using %s neighbor list subroutines\n", nmode);
}
#else
error->warning(FLERR,"OpenMP support not enabled during compilation; "
"using 1 thread only.");
#endif
}
// allocate list for per thread accumulator manager class instances

16
src/USER-OMP/pair_lj_sf_dipole_sf_omp.cpp
View File

@ -123,8 +123,8 @@ void PairLJSFDipoleSFOMP::eval(int iifrom, int iito, ThrData * const thr)
for (jj = 0; jj < jnum; jj++) {
j = jlist[jj];
factor_coul = special_coul[sbmask(j)];
factor_lj = special_lj[sbmask(j)];
factor_coul = special_coul[sbmask(j)];
j &= NEIGHMASK;
delx = xtmp - x[j].x;
@ -137,8 +137,6 @@ void PairLJSFDipoleSFOMP::eval(int iifrom, int iito, ThrData * const thr)
r2inv = 1.0/rsq;
rinv = sqrt(r2inv);
// atom can have both a charge and dipole
// i,j = charge-charge, dipole-dipole, dipole-charge, or charge-dipole
// atom can have both a charge and dipole
// i,j = charge-charge, dipole-dipole, dipole-charge, or charge-dipole
@ -148,8 +146,9 @@ void PairLJSFDipoleSFOMP::eval(int iifrom, int iito, ThrData * const thr)
if (rsq < cut_coulsq[itype][jtype]) {
rcutcoul2inv=1.0/cut_coulsq[itype][jtype];
if (qtmp != 0.0 && q[j] != 0.0) {
pre1 = qtmp*q[j]*rinv*(r2inv-1.0/cut_coulsq[itype][jtype]);
pre1 = qtmp*q[j]*rinv*(r2inv-rcutcoul2inv);
forcecoulx += pre1*delx;
forcecouly += pre1*dely;
@ -159,7 +158,6 @@ void PairLJSFDipoleSFOMP::eval(int iifrom, int iito, ThrData * const thr)
if (mu[i].w > 0.0 && mu[j].w > 0.0) {
r3inv = r2inv*rinv;
r5inv = r3inv*r2inv;
rcutcoul2inv=1.0/cut_coulsq[itype][jtype];
pdotp = mu[i].x*mu[j].x + mu[i].y*mu[j].y + mu[i].z*mu[j].z;
pidotr = mu[i].x*delx + mu[i].y*dely + mu[i].z*delz;
@ -171,7 +169,7 @@ void PairLJSFDipoleSFOMP::eval(int iifrom, int iito, ThrData * const thr)
aforcecouly = pre1*dely;
aforcecoulz = pre1*delz;
bfac = 1.0 - 4.0*rsq*sqrt(rsq)*rcutcoul2inv*sqrt(rcutcoul2inv) +
bfac = 1.0 - 4.0*rsq*sqrt(rsq*rcutcoul2inv)*rcutcoul2inv +
3.0*rsq*rsq*rcutcoul2inv*rcutcoul2inv;
presf = 2.0 * r2inv * pidotr * pjdotr;
bforcecoulx = bfac * (pjdotr*mu[i].x+pidotr*mu[j].x-presf*delx);
@ -202,10 +200,9 @@ void PairLJSFDipoleSFOMP::eval(int iifrom, int iito, ThrData * const thr)
r3inv = r2inv*rinv;
r5inv = r3inv*r2inv;
pidotr = mu[i].x*delx + mu[i].y*dely + mu[i].z*delz;
rcutcoul2inv=1.0/cut_coulsq[itype][jtype];
pre1 = 3.0 * q[j] * r5inv * pidotr * (1-rsq*rcutcoul2inv);
pqfac = 1.0 - 3.0*rsq*rcutcoul2inv +
2.0*rsq*sqrt(rsq)*rcutcoul2inv*sqrt(rcutcoul2inv);
2.0*rsq*sqrt(rsq*rcutcoul2inv)*rcutcoul2inv;
pre2 = q[j] * r3inv * pqfac;
forcecoulx += pre2*mu[i].x - pre1*delx;
@ -220,10 +217,9 @@ void PairLJSFDipoleSFOMP::eval(int iifrom, int iito, ThrData * const thr)
r3inv = r2inv*rinv;
r5inv = r3inv*r2inv;
pjdotr = mu[j].x*delx + mu[j].y*dely + mu[j].z*delz;
rcutcoul2inv=1.0/cut_coulsq[itype][jtype];
pre1 = 3.0 * qtmp * r5inv * pjdotr * (1-rsq*rcutcoul2inv);
qpfac = 1.0 - 3.0*rsq*rcutcoul2inv +
2.0*rsq*sqrt(rsq)*rcutcoul2inv*sqrt(rcutcoul2inv);
2.0*rsq*sqrt(rsq*rcutcoul2inv)*rcutcoul2inv;
pre2 = qtmp * r3inv * qpfac;
forcecoulx += pre1*delx - pre2*mu[j].x;

4
src/USER-QMMM/fix_qmmm.cpp
View File

@ -410,8 +410,8 @@ void FixQMMM::exchange_positions()
const double * const mass = atom->mass;
const tagint * const tag = atom->tag;
const int nlocal = atom->nlocal;
const int ntypes = atom->ntypes;
const int natoms = (int) atom->natoms;
const int ntypes = (int) atom->ntypes;
if (atom->natoms > MAXSMALLINT)
error->all(FLERR,"Too many MM atoms for fix qmmmm");
@ -444,7 +444,7 @@ void FixQMMM::exchange_positions()
isend_buf[1] = num_qm;
isend_buf[2] = num_mm;
isend_buf[3] = ntypes;
MPI_Send(isend_buf, 4, MPI_INTEGER,1, QMMM_TAG_SIZE, qm_comm);
MPI_Send(isend_buf, 4, MPI_INT, 1, QMMM_TAG_SIZE, qm_comm);
MPI_Send(celldata, 9, MPI_DOUBLE, 1, QMMM_TAG_CELL, qm_comm);
}
if (verbose > 0) {

3
src/USER-REAXC/fix_qeq_reax.cpp
View File

@ -416,7 +416,8 @@ void FixQEqReax::init_taper()
void FixQEqReax::setup_pre_force(int vflag)
{
neighbor->build_one(list);
// should not be needed
// neighbor->build_one(list);
deallocate_storage();
allocate_storage();

3
src/comm.cpp
View File

@ -87,7 +87,8 @@ Comm::Comm(LAMMPS *lmp) : Pointers(lmp)
} else if (getenv("OMP_NUM_THREADS") == NULL) {
nthreads = 1;
if (me == 0)
error->warning(FLERR,"OMP_NUM_THREADS environment is not set.");
error->warning(FLERR,"OMP_NUM_THREADS environment is not set. "
"Defaulting to 1 thread.");
} else {
nthreads = omp_get_max_threads();
}

176
src/compute_centro_atom.cpp
View File

@ -27,6 +27,7 @@
#include "force.h"
#include "pair.h"
#include "comm.h"
#include "math_extra.h"
#include "memory.h"
#include "error.h"
@ -37,23 +38,42 @@ using namespace LAMMPS_NS;
ComputeCentroAtom::ComputeCentroAtom(LAMMPS *lmp, int narg, char **arg) :
Compute(lmp, narg, arg)
{
if (narg != 4) error->all(FLERR,"Illegal compute centro/atom command");
if (narg < 4 || narg > 6) error->all(FLERR,"Illegal compute centro/atom command");
if (strcmp(arg[3],"fcc") == 0) nnn = 12;
else if (strcmp(arg[3],"bcc") == 0) nnn = 8;
else nnn = force->inumeric(FLERR,arg[3]);
// default values
axes_flag = 0;
// optional keywords
int iarg = 4;
while (iarg < narg) {
if (strcmp(arg[iarg],"axes") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal compute centro/atom command3");
if (strcmp(arg[iarg+1],"yes") == 0) axes_flag = 1;
else if (strcmp(arg[iarg+1],"no") == 0) axes_flag = 0;
else error->all(FLERR,"Illegal compute centro/atom command2");
iarg += 2;
} else error->all(FLERR,"Illegal compute centro/atom command1");
}
if (nnn <= 0 || nnn % 2)
error->all(FLERR,"Illegal neighbor value for compute centro/atom command");
peratom_flag = 1;
size_peratom_cols = 0;
if (!axes_flag) size_peratom_cols = 0;
else size_peratom_cols = 10;
nmax = 0;
centro = NULL;
maxneigh = 0;
distsq = NULL;
nearest = NULL;
array_atom = NULL;
}
/* ---------------------------------------------------------------------- */
@ -63,6 +83,7 @@ ComputeCentroAtom::~ComputeCentroAtom()
memory->destroy(centro);
memory->destroy(distsq);
memory->destroy(nearest);
if (axes_flag) memory->destroy(array_atom);
}
/* ---------------------------------------------------------------------- */
@ -106,12 +127,21 @@ void ComputeCentroAtom::compute_peratom()
invoked_peratom = update->ntimestep;
// grow centro array if necessary
// grow array_atom array if axes_flag set
if (atom->nmax > nmax) {
memory->destroy(centro);
nmax = atom->nmax;
memory->create(centro,nmax,"centro/atom:centro");
vector_atom = centro;
if (!axes_flag) {
memory->destroy(centro);
nmax = atom->nmax;
memory->create(centro,nmax,"centro/atom:centro");
vector_atom = centro;
} else {
memory->destroy(centro);
memory->destroy(array_atom);
nmax = atom->nmax;
memory->create(centro,nmax,"centro/atom:centro");
memory->create(array_atom,nmax,size_peratom_cols,"centro/atom:array_atom");
}
}
// invoke full neighbor list (will copy or build if necessary)
@ -174,45 +204,123 @@ void ComputeCentroAtom::compute_peratom()
}
}
// if not nnn neighbors, centro = 0.0
// check whether to include local crystal symmetry axes
if (!axes_flag) {
// if not nnn neighbors, centro = 0.0
if (n < nnn) {
centro[i] = 0.0;
continue;
if (n < nnn) {
centro[i] = 0.0;
continue;
}
// store nnn nearest neighs in 1st nnn locations of distsq and nearest
select2(nnn,n,distsq,nearest);
// R = Ri + Rj for each of npairs i,j pairs among nnn neighbors
// pairs = squared length of each R
n = 0;
for (j = 0; j < nnn; j++) {
jj = nearest[j];
for (k = j+1; k < nnn; k++) {
kk = nearest[k];
delx = x[jj][0] + x[kk][0] - 2.0*xtmp;
dely = x[jj][1] + x[kk][1] - 2.0*ytmp;
delz = x[jj][2] + x[kk][2] - 2.0*ztmp;
pairs[n++] = delx*delx + dely*dely + delz*delz;
}
}
} else {
// calculate local crystal symmetry axes
// rsq1, rsq2 are two smallest values of R^2
// R1, R2 are corresponding vectors Ri - Rj
// R3 is normal to R1, R2
double rsq1,rsq2;
double* r1 = &array_atom[i][1];
double* r2 = &array_atom[i][4];
double* r3 = &array_atom[i][7];
if (n < nnn) {
centro[i] = 0.0;
MathExtra::zero3(r1);
MathExtra::zero3(r2);
MathExtra::zero3(r3);
continue;
}
// store nnn nearest neighs in 1st nnn locations of distsq and nearest
select2(nnn,n,distsq,nearest);
n = 0;
rsq1 = rsq2 = cutsq;
for (j = 0; j < nnn; j++) {
jj = nearest[j];
for (k = j+1; k < nnn; k++) {
kk = nearest[k];
delx = x[jj][0] + x[kk][0] - 2.0*xtmp;
dely = x[jj][1] + x[kk][1] - 2.0*ytmp;
delz = x[jj][2] + x[kk][2] - 2.0*ztmp;
double rsq = delx*delx + dely*dely + delz*delz;
pairs[n++] = rsq;
if (rsq < rsq2) {
if (rsq < rsq1) {
rsq2 = rsq1;
MathExtra::copy3(r1, r2);
rsq1 = rsq;
MathExtra::sub3(x[jj],x[kk],r1);
} else {
rsq2 = rsq;
MathExtra::sub3(x[jj],x[kk],r2);
}
}
}
}
MathExtra::cross3(r1,r2,r3);
MathExtra::norm3(r1);
MathExtra::norm3(r2);
MathExtra::norm3(r3);
}
// store nnn nearest neighs in 1st nnn locations of distsq and nearest
select2(nnn,n,distsq,nearest);
// R = Ri + Rj for each of npairs i,j pairs among nnn neighbors
// pairs = squared length of each R
n = 0;
for (j = 0; j < nnn; j++) {
jj = nearest[j];
for (k = j+1; k < nnn; k++) {
kk = nearest[k];
delx = x[jj][0] + x[kk][0] - 2.0*xtmp;
dely = x[jj][1] + x[kk][1] - 2.0*ytmp;
delz = x[jj][2] + x[kk][2] - 2.0*ztmp;
pairs[n++] = delx*delx + dely*dely + delz*delz;
}
}
// store nhalf smallest pair distances in 1st nhalf locations of pairs
select(nhalf,npairs,pairs);
// centrosymmetry = sum of nhalf smallest squared values
value = 0.0;
for (j = 0; j < nhalf; j++) value += pairs[j];
centro[i] = value;
} else centro[i] = 0.0;
} else {
centro[i] = 0.0;
if (axes_flag) {
MathExtra::zero3(&array_atom[i][1]);
MathExtra::zero3(&array_atom[i][4]);
MathExtra::zero3(&array_atom[i][7]);
}
}
}
delete [] pairs;
if (axes_flag)
for (ii = 0; ii < inum; ii++) {
i = ilist[ii];
if (mask[i] & groupbit)
array_atom[i][0] = centro[i];
}
}
/* ----------------------------------------------------------------------

3
src/compute_centro_atom.h
View File

@ -39,7 +39,8 @@ class ComputeCentroAtom : public Compute {
int *nearest;
class NeighList *list;
double *centro;
int axes_flag;
void select(int, int, double *);
void select2(int, int, double *, int *);
};

1
src/compute_orientorder_atom.cpp
View File

@ -124,6 +124,7 @@ ComputeOrientOrderAtom::~ComputeOrientOrderAtom()
{
memory->destroy(qnarray);
memory->destroy(distsq);
memory->destroy(rlist);
memory->destroy(nearest);
memory->destroy(qlist);
memory->destroy(qnm_r);

8
src/fix_nh.cpp
View File

@ -79,6 +79,7 @@ FixNH::FixNH(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg)
etap_mass_flag = 0;
flipflag = 1;
dipole_flag = 0;
dlm_flag = 0;
tcomputeflag = 0;
pcomputeflag = 0;
@ -330,7 +331,10 @@ FixNH::FixNH(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg)
} else if (strcmp(arg[iarg],"update") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal fix nvt/npt/nph command");
if (strcmp(arg[iarg+1],"dipole") == 0) dipole_flag = 1;
else error->all(FLERR,"Illegal fix nvt/npt/nph command");
else if (strcmp(arg[iarg+1],"dipole/dlm") == 0) {
dipole_flag = 1;
dlm_flag = 1;
} else error->all(FLERR,"Illegal fix nvt/npt/nph command");
iarg += 2;
} else if (strcmp(arg[iarg],"fixedpoint") == 0) {
if (iarg+4 > narg) error->all(FLERR,"Illegal fix nvt/npt/nph command");
@ -434,7 +438,7 @@ FixNH::FixNH(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg)
if (!atom->mu_flag)
error->all(FLERR,"Using update dipole flag requires atom attribute mu");
}
if ((tstat_flag && t_period <= 0.0) ||
(p_flag[0] && p_period[0] <= 0.0) ||
(p_flag[1] && p_period[1] <= 0.0) ||

5
src/fix_nh.h
View File

@ -111,6 +111,7 @@ class FixNH : public Fix {
int omega_mass_flag; // 1 if omega_mass updated, 0 if not.
int etap_mass_flag; // 1 if etap_mass updated, 0 if not.
int dipole_flag; // 1 if dipole is updated, 0 if not.
int dlm_flag; // 1 if using the DLM rotational integrator, 0 if not
int scaleyz; // 1 if yz scaled with lz
int scalexz; // 1 if xz scaled with lz
@ -219,6 +220,10 @@ Self-explanatory.
E: Using update dipole flag requires atom attribute mu
Self-explanatory.
E: The dlm flag must be used with update dipole
Self-explanatory.
E: Fix nvt/npt/nph damping parameters must be > 0.0

143
src/fix_nh_sphere.cpp
View File

@ -21,9 +21,13 @@
#include "atom_vec.h"
#include "group.h"
#include "error.h"
#include "force.h"
#include "math_vector.h"
#include "math_extra.h"
using namespace LAMMPS_NS;
using namespace FixConst;
using namespace MathExtra;
#define INERTIA 0.4 // moment of inertia prefactor for sphere
@ -50,7 +54,7 @@ void FixNHSphere::init()
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit)
if (radius[i] == 0.0)
error->one(FLERR,"Fix nvt/sphere requires extended particles");
error->one(FLERR,"Fix nvt/npt/nph/sphere require extended particles");
FixNH::init();
}
@ -102,28 +106,133 @@ void FixNHSphere::nve_x()
FixNH::nve_x();
// update mu for dipoles
// d_mu/dt = omega cross mu
// renormalize mu to dipole length
if (dipole_flag) {
double msq,scale,g[3];
double **mu = atom->mu;
double **omega = atom->omega;
int *mask = atom->mask;
int nlocal = atom->nlocal;
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit)
if (mu[i][3] > 0.0) {
g[0] = mu[i][0] + dtv * (omega[i][1]*mu[i][2]-omega[i][2]*mu[i][1]);
g[1] = mu[i][1] + dtv * (omega[i][2]*mu[i][0]-omega[i][0]*mu[i][2]);
g[2] = mu[i][2] + dtv * (omega[i][0]*mu[i][1]-omega[i][1]*mu[i][0]);
msq = g[0]*g[0] + g[1]*g[1] + g[2]*g[2];
scale = mu[i][3]/sqrt(msq);
mu[i][0] = g[0]*scale;
mu[i][1] = g[1]*scale;
mu[i][2] = g[2]*scale;
if (dlm_flag == 0){
// d_mu/dt = omega cross mu
// renormalize mu to dipole length
double msq,scale,g[3];
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit)
if (mu[i][3] > 0.0) {
g[0] = mu[i][0] + dtv * (omega[i][1]*mu[i][2]-omega[i][2]*mu[i][1]);
g[1] = mu[i][1] + dtv * (omega[i][2]*mu[i][0]-omega[i][0]*mu[i][2]);
g[2] = mu[i][2] + dtv * (omega[i][0]*mu[i][1]-omega[i][1]*mu[i][0]);
msq = g[0]*g[0] + g[1]*g[1] + g[2]*g[2];
scale = mu[i][3]/sqrt(msq);
mu[i][0] = g[0]*scale;
mu[i][1] = g[1]*scale;
mu[i][2] = g[2]*scale;
}
} else {
// Integrate orientation following Dullweber-Leimkuhler-Maclachlan scheme
vector w, w_temp, a;
matrix Q, Q_temp, R;
double scale,s2,inv_len_mu;
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit && mu[i][3] > 0.0) {
// Construct Q from dipole:
// Q is the rotation matrix from space frame to body frame
// i.e. v_b = Q.v_s
// Define mu to lie along the z axis in the body frame
// We take the unit dipole to avoid getting a scaling matrix
inv_len_mu = 1.0/mu[i][3];
a[0] = mu[i][0]*inv_len_mu;
a[1] = mu[i][1]*inv_len_mu;
a[2] = mu[i][2]*inv_len_mu;
// v = a x [0 0 1] - cross product of mu in space and body frames
// s = |v|
// c = a.[0 0 1] = a[2]
// vx = [ 0 -v[2] v[1]
// v[2] 0 -v[0]
// -v[1] v[0] 0 ]
// then
// Q = I + vx + vx^2 * (1-c)/s^2
s2 = a[0]*a[0] + a[1]*a[1];
if (s2 != 0.0){ // i.e. the vectors are not parallel
scale = (1.0 - a[2])/s2;
Q[0][0] = 1.0 - scale*a[0]*a[0]; Q[0][1] = -scale*a[0]*a[1]; Q[0][2] = -a[0];
Q[1][0] = -scale*a[0]*a[1]; Q[1][1] = 1.0 - scale*a[1]*a[1]; Q[1][2] = -a[1];
Q[2][0] = a[0]; Q[2][1] = a[1]; Q[2][2] = 1.0 - scale*(a[0]*a[0] + a[1]*a[1]);
} else { // if parallel then we just have I or -I
Q[0][0] = 1.0/a[2]; Q[0][1] = 0.0; Q[0][2] = 0.0;
Q[1][0] = 0.0; Q[1][1] = 1.0/a[2]; Q[1][2] = 0.0;
Q[2][0] = 0.0; Q[2][1] = 0.0; Q[2][2] = 1.0/a[2];
}
// Local copy of this particle's angular velocity (in space frame)
w[0] = omega[i][0]; w[1] = omega[i][1]; w[2] = omega[i][2];
// Transform omega into body frame: w_temp= Q.w
matvec(Q,w,w_temp);
// Construct rotation R1
BuildRxMatrix(R, dtf/force->ftm2v*w_temp[0]);
// Apply R1 to w: w = R.w_temp
matvec(R,w_temp,w);
// Apply R1 to Q: Q_temp = R^T.Q
transpose_times3(R,Q,Q_temp);
// Construct rotation R2
BuildRyMatrix(R, dtf/force->ftm2v*w[1]);
// Apply R2 to w: w_temp = R.w
matvec(R,w,w_temp);
// Apply R2 to Q: Q = R^T.Q_temp
transpose_times3(R,Q_temp,Q);
// Construct rotation R3
BuildRzMatrix(R, 2.0*dtf/force->ftm2v*w_temp[2]);
// Apply R3 to w: w = R.w_temp
matvec(R,w_temp,w);
// Apply R3 to Q: Q_temp = R^T.Q
transpose_times3(R,Q,Q_temp);
// Construct rotation R4
BuildRyMatrix(R, dtf/force->ftm2v*w[1]);
// Apply R4 to w: w_temp = R.w
matvec(R,w,w_temp);
// Apply R4 to Q: Q = R^T.Q_temp
transpose_times3(R,Q_temp,Q);
// Construct rotation R5
BuildRxMatrix(R, dtf/force->ftm2v*w_temp[0]);
// Apply R5 to w: w = R.w_temp
matvec(R,w_temp,w);
// Apply R5 to Q: Q_temp = R^T.Q
transpose_times3(R,Q,Q_temp);
// Transform w back into space frame w_temp = Q^T.w
transpose_matvec(Q_temp,w,w_temp);
omega[i][0] = w_temp[0]; omega[i][1] = w_temp[1]; omega[i][2] = w_temp[2];
// Set dipole according to updated Q: mu = Q^T.[0 0 1] * |mu|
mu[i][0] = Q_temp[2][0] * mu[i][3];
mu[i][1] = Q_temp[2][1] * mu[i][3];
mu[i][2] = Q_temp[2][2] * mu[i][3];
}
}
}
}
}

149
src/fix_nve_sphere.cpp
View File

@ -21,13 +21,17 @@
#include "respa.h"
#include "force.h"
#include "error.h"
#include "math_vector.h"
#include "math_extra.h"
using namespace LAMMPS_NS;
using namespace FixConst;
using namespace MathExtra;
#define INERTIA 0.4 // moment of inertia prefactor for sphere
enum{NONE,DIPOLE};
enum{NODLM,DLM};
/* ---------------------------------------------------------------------- */
@ -41,13 +45,17 @@ FixNVESphere::FixNVESphere(LAMMPS *lmp, int narg, char **arg) :
// process extra keywords
extra = NONE;
dlm = NODLM;
int iarg = 3;
while (iarg < narg) {
if (strcmp(arg[iarg],"update") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal fix nve/sphere command");
if (strcmp(arg[iarg+1],"dipole") == 0) extra = DIPOLE;
else error->all(FLERR,"Illegal fix nve/sphere command");
else if (strcmp(arg[iarg+1],"dipole/dlm") == 0) {
extra = DIPOLE;
dlm = DLM;
} else error->all(FLERR,"Illegal fix nve/sphere command");
iarg += 2;
} else error->all(FLERR,"Illegal fix nve/sphere command");
}
@ -57,7 +65,7 @@ FixNVESphere::FixNVESphere(LAMMPS *lmp, int narg, char **arg) :
if (!atom->sphere_flag)
error->all(FLERR,"Fix nve/sphere requires atom style sphere");
if (extra == DIPOLE && !atom->mu_flag)
error->all(FLERR,"Fix nve/sphere dipole requires atom attribute mu");
error->all(FLERR,"Fix nve/sphere update dipole requires atom attribute mu");
}
/* ---------------------------------------------------------------------- */
@ -83,8 +91,10 @@ void FixNVESphere::init()
void FixNVESphere::initial_integrate(int vflag)
{
double dtfm,dtirotate,msq,scale;
double dtfm,dtirotate,msq,scale,s2,inv_len_mu;
double g[3];
vector w, w_temp, a;
matrix Q, Q_temp, R;
double **x = atom->x;
double **v = atom->v;
@ -120,25 +130,127 @@ void FixNVESphere::initial_integrate(int vflag)
omega[i][2] += dtirotate * torque[i][2];
}
}
// update mu for dipoles
// d_mu/dt = omega cross mu
// renormalize mu to dipole length
if (extra == DIPOLE) {
double **mu = atom->mu;
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit)
if (mu[i][3] > 0.0) {
g[0] = mu[i][0] + dtv * (omega[i][1]*mu[i][2]-omega[i][2]*mu[i][1]);
g[1] = mu[i][1] + dtv * (omega[i][2]*mu[i][0]-omega[i][0]*mu[i][2]);
g[2] = mu[i][2] + dtv * (omega[i][0]*mu[i][1]-omega[i][1]*mu[i][0]);
msq = g[0]*g[0] + g[1]*g[1] + g[2]*g[2];
scale = mu[i][3]/sqrt(msq);
mu[i][0] = g[0]*scale;
mu[i][1] = g[1]*scale;
mu[i][2] = g[2]*scale;
if (dlm == NODLM) {
// d_mu/dt = omega cross mu
// renormalize mu to dipole length
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit)
if (mu[i][3] > 0.0) {
g[0] = mu[i][0] + dtv * (omega[i][1]*mu[i][2]-omega[i][2]*mu[i][1]);
g[1] = mu[i][1] + dtv * (omega[i][2]*mu[i][0]-omega[i][0]*mu[i][2]);
g[2] = mu[i][2] + dtv * (omega[i][0]*mu[i][1]-omega[i][1]*mu[i][0]);
msq = g[0]*g[0] + g[1]*g[1] + g[2]*g[2];
scale = mu[i][3]/sqrt(msq);
mu[i][0] = g[0]*scale;
mu[i][1] = g[1]*scale;
mu[i][2] = g[2]*scale;
}
} else {
// Integrate orientation following Dullweber-Leimkuhler-Maclachlan scheme
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit && mu[i][3] > 0.0) {
// Construct Q from dipole:
// Q is the rotation matrix from space frame to body frame
// i.e. v_b = Q.v_s
// Define mu to lie along the z axis in the body frame
// We take the unit dipole to avoid getting a scaling matrix
inv_len_mu = 1.0/mu[i][3];
a[0] = mu[i][0]*inv_len_mu;
a[1] = mu[i][1]*inv_len_mu;
a[2] = mu[i][2]*inv_len_mu;
// v = a x [0 0 1] - cross product of mu in space and body frames
// s = |v|
// c = a.[0 0 1] = a[2]
// vx = [ 0 -v[2] v[1]
// v[2] 0 -v[0]
// -v[1] v[0] 0 ]
// then
// Q = I + vx + vx^2 * (1-c)/s^2
s2 = a[0]*a[0] + a[1]*a[1];
if (s2 != 0.0){ // i.e. the vectors are not parallel
scale = (1.0 - a[2])/s2;
Q[0][0] = 1.0 - scale*a[0]*a[0]; Q[0][1] = -scale*a[0]*a[1]; Q[0][2] = -a[0];
Q[1][0] = -scale*a[0]*a[1]; Q[1][1] = 1.0 - scale*a[1]*a[1]; Q[1][2] = -a[1];
Q[2][0] = a[0]; Q[2][1] = a[1]; Q[2][2] = 1.0 - scale*(a[0]*a[0] + a[1]*a[1]);
} else { // if parallel then we just have I or -I
Q[0][0] = 1.0/a[2]; Q[0][1] = 0.0; Q[0][2] = 0.0;
Q[1][0] = 0.0; Q[1][1] = 1.0/a[2]; Q[1][2] = 0.0;
Q[2][0] = 0.0; Q[2][1] = 0.0; Q[2][2] = 1.0/a[2];
}
// Local copy of this particle's angular velocity (in space frame)
w[0] = omega[i][0]; w[1] = omega[i][1]; w[2] = omega[i][2];
// Transform omega into body frame: w_temp= Q.w
matvec(Q,w,w_temp);
// Construct rotation R1
BuildRxMatrix(R, dtf/force->ftm2v*w_temp[0]);
// Apply R1 to w: w = R.w_temp
matvec(R,w_temp,w);
// Apply R1 to Q: Q_temp = R^T.Q
transpose_times3(R,Q,Q_temp);
// Construct rotation R2
BuildRyMatrix(R, dtf/force->ftm2v*w[1]);
// Apply R2 to w: w_temp = R.w
matvec(R,w,w_temp);
// Apply R2 to Q: Q = R^T.Q_temp
transpose_times3(R,Q_temp,Q);
// Construct rotation R3
BuildRzMatrix(R, 2.0*dtf/force->ftm2v*w_temp[2]);
// Apply R3 to w: w = R.w_temp
matvec(R,w_temp,w);
// Apply R3 to Q: Q_temp = R^T.Q
transpose_times3(R,Q,Q_temp);
// Construct rotation R4
BuildRyMatrix(R, dtf/force->ftm2v*w[1]);
// Apply R4 to w: w_temp = R.w
matvec(R,w,w_temp);
// Apply R4 to Q: Q = R^T.Q_temp
transpose_times3(R,Q_temp,Q);
// Construct rotation R5
BuildRxMatrix(R, dtf/force->ftm2v*w_temp[0]);
// Apply R5 to w: w = R.w_temp
matvec(R,w_temp,w);
// Apply R5 to Q: Q_temp = R^T.Q
transpose_times3(R,Q,Q_temp);
// Transform w back into space frame w_temp = Q^T.w
transpose_matvec(Q_temp,w,w_temp);
omega[i][0] = w_temp[0]; omega[i][1] = w_temp[1]; omega[i][2] = w_temp[2];
// Set dipole according to updated Q: mu = Q^T.[0 0 1] * |mu|
mu[i][0] = Q_temp[2][0] * mu[i][3];
mu[i][1] = Q_temp[2][1] * mu[i][3];
mu[i][2] = Q_temp[2][2] * mu[i][3];
}
}
}
}
}
@ -165,6 +277,7 @@ void FixNVESphere::final_integrate()
// update v,omega for all particles
// d_omega/dt = torque / inertia
double rke = 0.0;
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
dtfm = dtf / rmass[i];
@ -176,5 +289,7 @@ void FixNVESphere::final_integrate()
omega[i][0] += dtirotate * torque[i][0];
omega[i][1] += dtirotate * torque[i][1];
omega[i][2] += dtirotate * torque[i][2];
rke += (omega[i][0]*omega[i][0] + omega[i][1]*omega[i][1] + omega[i][2]*omega[i][2])*radius[i]*radius[i]*rmass[i];
}
}

8
src/fix_nve_sphere.h
View File

@ -34,6 +34,7 @@ class FixNVESphere : public FixNVE {
protected:
int extra;
int dlm;
};
}
@ -53,12 +54,17 @@ E: Fix nve/sphere requires atom style sphere
Self-explanatory.
E: Fix nve/sphere dipole requires atom attribute mu
E: Fix nve/sphere update dipole requires atom attribute mu
An atom style with this attribute is needed.
E: Fix nve/sphere requires extended particles
This fix can only be used for particles of a finite size.
E: Fix nve/sphere dlm must be used with update dipole
The DLM algorithm can only be used in conjunction with update dipole.
*/

3
src/info.cpp
View File

@ -264,7 +264,8 @@ void Info::command(int narg, char **arg)
}
fprintf(out,"Nprocs = %d, Nthreads = %d\n",
comm->nprocs, comm->nthreads);
fprintf(out,"Processor grid = %d x %d x %d\n",comm->procgrid[0],
if (domain->box_exist)
fprintf(out,"Processor grid = %d x %d x %d\n",comm->procgrid[0],
comm->procgrid[1], comm->procgrid[2]);
}

45
src/math_extra.cpp
View File

@ -607,6 +607,51 @@ void inertia_triangle(double *idiag, double *quat, double mass,
inertia[5] = tensor[0][1];
}
/* ----------------------------------------------------------------------
Build rotation matrix for a small angle rotation around the X axis
------------------------------------------------------------------------- */
void BuildRxMatrix(double R[3][3], const double angle)
{
const double angleSq = angle * angle;
const double cosAngle = (1.0 - angleSq * 0.25) / (1.0 + angleSq * 0.25);
const double sinAngle = angle / (1.0 + angleSq * 0.25);
R[0][0] = 1.0; R[0][1] = 0.0; R[0][2] = 0.0;
R[1][0] = 0.0; R[1][1] = cosAngle; R[1][2] = -sinAngle;
R[2][0] = 0.0; R[2][1] = sinAngle; R[2][2] = cosAngle;
}
/* ----------------------------------------------------------------------
Build rotation matrix for a small angle rotation around the Y axis
------------------------------------------------------------------------- */
void BuildRyMatrix(double R[3][3], const double angle)
{
const double angleSq = angle * angle;
const double cosAngle = (1.0 - angleSq * 0.25) / (1.0 + angleSq * 0.25);
const double sinAngle = angle / (1.0 + angleSq * 0.25);
R[0][0] = cosAngle; R[0][1] = 0.0; R[0][2] = sinAngle;
R[1][0] = 0.0; R[1][1] = 1.0; R[1][2] = 0.0;
R[2][0] = -sinAngle; R[2][1] = 0.0; R[2][2] = cosAngle;
}
/* ----------------------------------------------------------------------
Build rotation matrix for a small angle rotation around the Y axis
------------------------------------------------------------------------- */
void BuildRzMatrix(double R[3][3], const double angle)
{
const double angleSq = angle * angle;
const double cosAngle = (1.0 - angleSq * 0.25) / (1.0 + angleSq * 0.25);
const double sinAngle = angle / (1.0 + angleSq * 0.25);
R[0][0] = cosAngle; R[0][1] = -sinAngle; R[0][2] = 0.0;
R[1][0] = sinAngle; R[1][1] = cosAngle; R[1][2] = 0.0;
R[2][0] = 0.0; R[2][1] = 0.0; R[2][2] = 1.0;
}
/* ---------------------------------------------------------------------- */
}

40
src/math_extra.h
View File

@ -27,12 +27,15 @@ namespace MathExtra {
// 3 vector operations
inline void copy3(const double *v, double *ans);
inline void zero3(double *v);
inline void norm3(double *v);
inline void normalize3(const double *v, double *ans);
inline void snormalize3(const double, const double *v, double *ans);
inline void negate3(double *v);
inline void scale3(double s, double *v);
inline void add3(const double *v1, const double *v2, double *ans);
inline void scaleadd3(double s, const double *v1, const double *v2, double *ans);
inline void sub3(const double *v1, const double *v2, double *ans);
inline double len3(const double *v);
inline double lensq3(const double *v);
@ -111,6 +114,10 @@ namespace MathExtra {
inline void rotation_generator_x(const double m[3][3], double ans[3][3]);
inline void rotation_generator_y(const double m[3][3], double ans[3][3]);
inline void rotation_generator_z(const double m[3][3], double ans[3][3]);
void BuildRxMatrix(double R[3][3], const double angle);
void BuildRyMatrix(double R[3][3], const double angle);
void BuildRzMatrix(double R[3][3], const double angle);
// moment of inertia operations
@ -124,6 +131,28 @@ namespace MathExtra {
double *inertia);
}
/* ----------------------------------------------------------------------
copy a vector, return in ans
------------------------------------------------------------------------- */
inline void MathExtra::copy3(const double *v, double *ans)
{
ans[0] = v[0];
ans[1] = v[1];
ans[2] = v[2];
}
/* ----------------------------------------------------------------------
set vector equal to zero
------------------------------------------------------------------------- */
inline void MathExtra::zero3(double *v)
{
v[0] = 0.0;
v[1] = 0.0;
v[2] = 0.0;
}
/* ----------------------------------------------------------------------
normalize a vector in place
------------------------------------------------------------------------- */
@ -194,6 +223,17 @@ inline void MathExtra::add3(const double *v1, const double *v2, double *ans)
ans[2] = v1[2] + v2[2];
}
/* ----------------------------------------------------------------------
ans = s*v1 + v2
------------------------------------------------------------------------- */
inline void MathExtra::scaleadd3(double s, const double *v1, const double *v2, double *ans)
{
ans[0] = s*v1[0] + v2[0];
ans[1] = s*v1[1] + v2[1];
ans[2] = s*v1[2] + v2[2];
}
/* ----------------------------------------------------------------------
ans = v1 - v2
------------------------------------------------------------------------- */

21
src/min.cpp
View File

@ -241,12 +241,21 @@ void Min::setup()
// remove these restriction eventually
if (nextra_global && searchflag == 0)
error->all(FLERR,
"Cannot use a damped dynamics min style with fix box/relax");
if (nextra_atom && searchflag == 0)
error->all(FLERR,
"Cannot use a damped dynamics min style with per-atom DOF");
if (searchflag == 0) {
if (nextra_global)
error->all(FLERR,
"Cannot use a damped dynamics min style with fix box/relax");
if (nextra_atom)
error->all(FLERR,
"Cannot use a damped dynamics min style with per-atom DOF");
}
if (strcmp(update->minimize_style,"hftn") == 0) {
if (nextra_global)
error->all(FLERR, "Cannot use hftn min style with fix box/relax");
if (nextra_atom)
error->all(FLERR, "Cannot use hftn min style with per-atom DOF");
}
// atoms may have migrated in comm->exchange()

5
src/min_hftn.cpp
View File

@ -138,8 +138,11 @@ void MinHFTN::setup_style()
//---- ALLOCATE MEMORY FOR EXTRA GLOBAL DEGREES OF FREEDOM.
//---- THE FIX MODULE TAKES CARE OF THE FIRST VECTOR, X0 (XK).
if (nextra_global) {
for (int i = 1; i < NUM_HFTN_ATOM_BASED_VECTORS; i++)
for (int i = 1; i < NUM_HFTN_ATOM_BASED_VECTORS; i++) {
_daExtraGlobal[i] = new double[nextra_global];
for (int j = 0; j < nextra_global; j++)
_daExtraGlobal[i][j] = 0.0;
}
}
//---- ALLOCATE MEMORY FOR EXTRA PER-ATOM DEGREES OF FREEDOM.

24
src/min_hftn.h
View File

@ -24,17 +24,6 @@ MinimizeStyle(hftn,MinHFTN)
namespace LAMMPS_NS
{
//---- THE ALGORITHM NEEDS TO STORE THIS MANY ATOM-BASED VECTORS,
//---- IN ADDITION TO ATOM POSITIONS AND THE FORCE VECTOR.
static const int NUM_HFTN_ATOM_BASED_VECTORS = 7;
static const int VEC_XK = 0; //-- ATOM POSITIONS AT SUBITER START
static const int VEC_CG_P = 1; //-- STEP p IN CG SUBITER
static const int VEC_CG_D = 2; //-- DIRECTION d IN CG SUBITER
static const int VEC_CG_HD = 3; //-- HESSIAN-VECTOR PRODUCT Hd
static const int VEC_CG_R = 4; //-- RESIDUAL r IN CG SUBITER
static const int VEC_DIF1 = 5; //-- FOR FINITE DIFFERENCING
static const int VEC_DIF2 = 6; //-- FOR FINITE DIFFERENCING
class MinHFTN : public Min
{
@ -49,6 +38,19 @@ class MinHFTN : public Min
private:
//---- THE ALGORITHM NEEDS TO STORE THIS MANY ATOM-BASED VECTORS,
//---- IN ADDITION TO ATOM POSITIONS AND THE FORCE VECTOR.
enum {
VEC_XK=0, //-- ATOM POSITIONS AT SUBITER START
VEC_CG_P, //-- STEP p IN CG SUBITER
VEC_CG_D, //-- DIRECTION d IN CG SUBITER
VEC_CG_HD, //-- HESSIAN-VECTOR PRODUCT Hd
VEC_CG_R, //-- RESIDUAL r IN CG SUBITER
VEC_DIF1, //-- FOR FINITE DIFFERENCING
VEC_DIF2, //-- FOR FINITE DIFFERENCING
NUM_HFTN_ATOM_BASED_VECTORS
};
//---- ATOM-BASED STORAGE VECTORS.
double * _daAVectors[NUM_HFTN_ATOM_BASED_VECTORS];
double ** _daExtraAtom[NUM_HFTN_ATOM_BASED_VECTORS];

1
src/neigh_request.h
View File

@ -105,6 +105,7 @@ class NeighRequest : protected Pointers {
int skip; // 1 if this list skips atom types from another list
int *iskip; // iskip[i] if atoms of type I are not in list
int **ijskip; // ijskip[i][j] if pairs of type I,J are not in list
int off2on; // 1 if this is newton on list, but skips from off list
int otherlist; // index of other list to copy or skip from

8
src/pair_dpd.cpp
View File

@ -164,12 +164,13 @@ void PairDPD::compute(int eflag, int vflag)
void PairDPD::allocate()
{
int i,j;
allocated = 1;
int n = atom->ntypes;
memory->create(setflag,n+1,n+1,"pair:setflag");
for (int i = 1; i <= n; i++)
for (int j = i; j <= n; j++)
for (i = 1; i <= n; i++)
for (j = i; j <= n; j++)
setflag[i][j] = 0;
memory->create(cutsq,n+1,n+1,"pair:cutsq");
@ -178,6 +179,9 @@ void PairDPD::allocate()
memory->create(a0,n+1,n+1,"pair:a0");
memory->create(gamma,n+1,n+1,"pair:gamma");
memory->create(sigma,n+1,n+1,"pair:sigma");
for (i = 0; i <= atom->ntypes; i++)
for (j = 0; j <= atom->ntypes; j++)
sigma[i][j] = gamma[i][j] = 0.0;
}
/* ----------------------------------------------------------------------

14
src/pair_hybrid.cpp
View File

@ -684,9 +684,10 @@ void PairHybrid::modify_requests()
}
// adjustments to newly added granular parent requests (gran = 1)
// parent newton = 2 if has children with granonesided = 0 and 1
// set parent newton = 2 if has children with granonesided = 0 and 1
// else newton = 0 = setting of children
// parent gran onesided = 0 if has children with granonesided = 0 and 1
// if 2, also set child off2on for both granonesided kinds of children
// set parent gran onesided = 0 if has children with granonesided = 0 and 1
// else onesided = setting of children
for (i = nrequest_original; i < neighbor->nrequest; i++) {
@ -698,6 +699,7 @@ void PairHybrid::modify_requests()
for (j = 0; j < nrequest_original; j++) {
if (!neighbor->requests[j]->pair) continue;
if (!neighbor->requests[j]->gran) continue;
if (neighbor->requests[j]->otherlist != i) continue;
jrq = neighbor->requests[j];
if (onesided < 0) onesided = jrq->granonesided;
@ -708,6 +710,14 @@ void PairHybrid::modify_requests()
if (onesided == 2) {
irq->newton = 2;
irq->granonesided = 0;
for (j = 0; j < nrequest_original; j++) {
if (!neighbor->requests[j]->pair) continue;
if (!neighbor->requests[j]->gran) continue;
if (neighbor->requests[j]->otherlist != i) continue;
jrq = neighbor->requests[j];
jrq->off2on = 1;
}
}
}
}

14
src/thermo.cpp
View File

@ -573,7 +573,7 @@ void Thermo::allocate()
int n = nfield_initial + 1;
keyword = new char*[n];
for (int i = 0; i < n; i++) keyword[i] = new char[32];
for (int i = 0; i < n; i++) keyword[i] = NULL;
vfunc = new FnPtr[n];
vtype = new int[n];
@ -821,7 +821,6 @@ void Thermo::parse_fields(char *str)
// compute value = c_ID, fix value = f_ID, variable value = v_ID
// count trailing [] and store int arguments
// copy = at most 8 chars of ID to pass to addfield
} else if ((strncmp(word,"c_",2) == 0) || (strncmp(word,"f_",2) == 0) ||
(strncmp(word,"v_",2) == 0)) {
@ -829,9 +828,6 @@ void Thermo::parse_fields(char *str)
int n = strlen(word);
char *id = new char[n];
strcpy(id,&word[2]);
char copy[9];
strncpy(copy,id,8);
copy[8] = '\0';
// parse zero or one or two trailing brackets from ID
// argindex1,argindex2 = int inside each bracket pair, 0 if no bracket
@ -878,7 +874,7 @@ void Thermo::parse_fields(char *str)
field2index[nfield] = add_compute(id,VECTOR);
else
field2index[nfield] = add_compute(id,ARRAY);
addfield(copy,&Thermo::compute_compute,FLOAT);
addfield(word,&Thermo::compute_compute,FLOAT);
} else if (word[0] == 'f') {
n = modify->find_fix(id);
@ -903,7 +899,7 @@ void Thermo::parse_fields(char *str)
}
field2index[nfield] = add_fix(id);
addfield(copy,&Thermo::compute_fix,FLOAT);
addfield(word,&Thermo::compute_fix,FLOAT);
} else if (word[0] == 'v') {
n = input->variable->find(id);
@ -919,7 +915,7 @@ void Thermo::parse_fields(char *str)
error->all(FLERR,"Thermo custom variable cannot have two indices");
field2index[nfield] = add_variable(id);
addfield(copy,&Thermo::compute_variable,FLOAT);
addfield(word,&Thermo::compute_variable,FLOAT);
}
delete [] id;
@ -936,6 +932,8 @@ void Thermo::parse_fields(char *str)
void Thermo::addfield(const char *key, FnPtr func, int typeflag)
{
int n = strlen(key) + 1;
keyword[nfield] = new char[n];
strcpy(keyword[nfield],key);
vfunc[nfield] = func;
vtype[nfield] = typeflag;

2
src/version.h
View File

@ -1 +1 @@
#define LAMMPS_VERSION "22 Jul 2016"
#define LAMMPS_VERSION "27 Jul 2016"

12
src/write_coeff.cpp
View File

@ -33,7 +33,7 @@ using namespace LAMMPS_NS;
void WriteCoeff::command(int narg, char **arg)
{
if (domain->box_exist == 0)
error->all(FLERR,"write_coeff command before simulation box is defined");
error->all(FLERR,"Write_coeff command before simulation box is defined");
if (narg != 1) error->all(FLERR,"Illegal write_coeff command");
@ -43,6 +43,9 @@ void WriteCoeff::command(int narg, char **arg)
strcpy(file,"tmp.");
strcat(file,arg[0]);
// initialize relevant styles
force->init();
if (comm->me == 0) {
char str[256], coeff[256];
FILE *one = fopen(file,"wb+");
@ -52,7 +55,6 @@ void WriteCoeff::command(int narg, char **arg)
}
if (force->pair && force->pair->writedata) {
force->pair->init();
fprintf(one,"# pair_style %s\npair_coeff\n",force->pair_style);
force->pair->write_data_all(one);
fprintf(one,"end\n");
@ -68,12 +70,14 @@ void WriteCoeff::command(int narg, char **arg)
fprintf(one,"end\n");
}
if (force->dihedral && force->dihedral->writedata) {
fprintf(one,"# dihedral_style %s\ndihedral_coeff\n",force->dihedral_style);
fprintf(one,"# dihedral_style %s\ndihedral_coeff\n",
force->dihedral_style);
force->dihedral->write_data(one);
fprintf(one,"end\n");
}
if (force->improper && force->improper->writedata) {
fprintf(one,"# improper_style %s\nimproper_coeff\n",force->improper_style);
fprintf(one,"# improper_style %s\nimproper_coeff\n",
force->improper_style);
force->improper->write_data(one);
fprintf(one,"end\n");
}

6
test/02_commands/in.min-007
View File

@ -17,8 +17,10 @@ pair_coeff 1 1 1.0 1.0 2.5
neighbor 0.3 bin
displace_atoms all random 0.1 0.1 0.1 12314
min_style hftn
fix 1 all box/relax iso 0.0
min_style sd
fix 1 all box/relax aniso 0.0
thermo 100
minimize 0.0 0.0 100 1000
min_style cg
minimize 0.0 0.0 100 1000

28
test/02_commands/in.min-009 Normal file
View File

@ -0,0 +1,28 @@
timer off
units lj
atom_style atomic
lattice fcc 0.8442
region box block 0 10 0 10 0 10
create_box 1 box
create_atoms 1 box
mass 1 1.0
change_box all triclinic
velocity all create 3.0 87287
pair_style lj/cut 2.5
pair_coeff 1 1 1.0 1.0 2.5
neighbor 0.3 bin
displace_atoms all random 0.1 0.1 0.1 12314
min_style cg
fix 1 all box/relax tri 0.0 fixedpoint 5.0 5.0 5.0
thermo 100
minimize 0.0 0.0 100 1000
min_style sd
minimize 0.0 0.0 100 1000