jianmu
/
lammps
forked from lijiext/lammps
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Merge branch 'integration' into doc-corrections-and-updates

This commit is contained in:
Axel Kohlmeyer 2016-09-06 20:19:02 -04:00
parent ace5dc3c7c 0252347d43
commit 1c17b98500
16 changed files with 376 additions and 241 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
doc/html/Manual.html
View File

@ -134,8 +134,8 @@
<H1></H1><div class="section" id="lammps-documentation">
<h1>LAMMPS Documentation</h1>
<div class="section" id="aug-2016-version">
<h2>27 Aug 2016 version</h2>
<div class="section" id="sep-2016-version">
<h2>7 Sep 2016 version</h2>
</div>
<div class="section" id="version-info">
<h2>Version info:</h2>

4
doc/html/Section_howto.html
View File

@ -1773,9 +1773,9 @@ translation velocities but also rotational velocities for spherical
and aspherical particles.</p>
<p>All of the barostatting fixes use the <a class="reference internal" href="compute_pressure.html"><span class="doc">compute pressure</span></a> compute to calculate a current
pressure. By default, this compute is created with a simple <a class="reference internal" href="compute_temp.html"><span class="doc">compute temp</span></a> (see the last argument of the <a class="reference internal" href="compute_pressure.html"><span class="doc">compute pressure</span></a> command), which is used to calculated
the kinetic componenet of the pressure. The barostatting fixes can
the kinetic component of the pressure. The barostatting fixes can
also use temperature computes that remove bias for the purpose of
computing the kinetic componenet which contributes to the current
computing the kinetic component which contributes to the current
pressure. See the doc pages for the individual fixes and for the
<a class="reference internal" href="fix_modify.html"><span class="doc">fix_modify</span></a> command for instructions on how to assign
a temperature or pressure compute to a barostatting fix.</p>

122
doc/html/Section_modify.html
View File

@ -1057,11 +1057,35 @@ then be accessed by variables) was discussed
</div>
<div class="section" id="submitting-new-features-for-inclusion-in-lammps">
<span id="mod-15"></span><h2>10.15. Submitting new features for inclusion in LAMMPS</h2>
<p>We encourage users to submit new features to <a class="reference external" href="http://lammps.sandia.gov/authors.html">the developers</a> that they add to
LAMMPS, especially if you think they will be of interest to other
users. The preferred way to do this is via GitHub. Once you have
prepared the content described below, see <a class="reference internal" href="tutorial_github.html"><span class="doc">this tutorial</span></a> for instructions on how to submit
your changes or new files.</p>
<p>We encourage users to submit new features or modifications for
LAMMPS to <a class="reference external" href="http://lammps.sandia.gov/authors.html">the core developers</a>
so they can be added to the LAMMPS distribution. The preferred way to
manage and coordinate this is as of Fall 2016 via the LAMMPS project on
<a class="reference external" href="https://github.com/lammps/lammps">GitHub</a>. An alternative is to contact
the LAMMPS developers or the indicated developer of a package or feature
directly and send in your contribution via e-mail.</p>
<p>For any larger modifications or programming project, you are encouraged
to contact the LAMMPS developers ahead of time, in order to discuss
implementation strategies and coding guidelines, that will make it
easier to integrate your contribution and result in less work for
everybody involved. You are also encouraged to search through the list
of <a class="reference external" href="https://github.com/lammps/lammps/issues">open issues on GitHub</a> and
submit a new issue for a planned feature, so you would not duplicate
the work of others (and possibly get scooped by them) or have your work
duplicated by others.</p>
<p>How quickly your contribution will be integrated
depends largely on how much effort it will cause to integrate and test
it, how much it requires changes to the core codebase, and of how much
interest it is to the larger LAMMPS community. Please see below for a
checklist of typical requirements. Once you have prepared everything,
see <a class="reference internal" href="tutorial_github.html"><span class="doc">this tutorial</span></a> for instructions on how to
submit your changes or new files through a GitHub pull request. If you
prefer to submit patches or full files, you should first make certain,
that your code works correctly with the latest patch-level version of
LAMMPS and contains all bugfixes from it. Then create a gzipped tar
file of all changed or added files or a corresponding patch file using
&#8216;diff -u&#8217; or &#8216;diff -c&#8217; and compress it with gzip. Please only use
gzip compression, as this works well on all platforms.</p>
<p>If the new features/files are broadly useful we may add them as core
files to LAMMPS or as part of a <a class="reference internal" href="Section_start.html#start-3"><span class="std std-ref">standard package</span></a>. Else we will add them as a
user-contributed file or package. Examples of user packages are in
@ -1104,20 +1128,34 @@ changes in core LAMMPS files, you&#8217;ll need to <a class="reference external"
not want to make those changes. An example of a trivial change is
making a parent-class method &#8220;virtual&#8221; when you derive a new child
class from it.</p>
<p>Here are the steps you need to follow to submit a single file or user
package for our consideration. Following these steps will save both
you and us time. See existing files in packages in the src dir for
examples.</p>
<p>Here is a checklist of steps you need to follow to submit a single file
or user package for our consideration. Following these steps will save
both you and us time. See existing files in packages in the src dir for
examples. If you are uncertain, please ask.</p>
<ul class="simple">
<li>All source files you provide must compile with the most current
version of LAMMPS.</li>
<li>If you want your file(s) to be added to main LAMMPS or one of its
standard packages, then it needs to be written in a style compatible
with other LAMMPS source files. This is so the developers can
understand it and hopefully maintain it. This basically means that
the code accesses data structures, performs its operations, and is
formatted similar to other LAMMPS source files, including the use of
the error class for error and warning messages.</li>
version of LAMMPS with multiple configurations. In particular you
need to test compiling LAMMPS from scratch with -DLAMMPS_BIGBIG
set in addition to the default -DLAMMPS_SMALLBIG setting. Your code
will need to work correctly in serial and in parallel using MPI.</li>
<li>For consistency with the rest of LAMMPS and especially, if you want
your contribution(s) to be added to main LAMMPS code or one of its
standard packages, it needs to be written in a style compatible with
other LAMMPS source files. This means: 2-character indentation per
level, <strong>no tabs</strong>, no lines over 80 characters. I/O is done via
the C-style stdio library, class header files should not import any
system headers outside &lt;stdio.h&gt;, STL containers should be avoided
in headers, and forward declarations used where possible or needed.
All added code should be placed into the LAMMPS_NS namespace or a
sub-namespace; global or static variables should be avoided, as they
conflict with the modular nature of LAMMPS and the C++ class structure.
Header files must <strong>not</strong> import namespaces with <em>using</em>.
This all is so the developers can more easily understand, integrate,
and maintain your contribution and reduce conflicts with other parts
of LAMMPS. This basically means that the code accesses data
structures, performs its operations, and is formatted similar to other
LAMMPS source files, including the use of the error class for error
and warning messages.</li>
<li>If you want your contribution to be added as a user-contributed
feature, and it&#8217;s a single file (actually a <a href="#id8"><span class="problematic" id="id9">*</span></a>.cpp and <a href="#id10"><span class="problematic" id="id11">*</span></a>.h file) it can
rapidly be added to the USER-MISC directory. Send us the one-line
@ -1144,13 +1182,21 @@ coding style (see above). I.e. the files do not need to be in the
same stylistic format and syntax as other LAMMPS files, though that
would be nice for developers as well as users who try to read your
code.</li>
<li>You must also create a documentation file for each new command or
style you are adding to LAMMPS. This will be one file for a
single-file feature. For a package, it might be several files. These
are simple text files which we auto-convert to HTML. Thus they must
be in the same format as other <em>.txt files in the lammps/doc directory
for similar commands and styles; use one or more of them as a starting
point. As appropriate, the text files can include links to equations
<li>You <strong>must</strong> also create a <strong>documentation</strong> file for each new command or
style you are adding to LAMMPS. For simplicity and convenience, the
documentation of groups of closely related commands or styles may be
combined into a single file. This will be one file for a single-file
feature. For a package, it might be several files. These are simple
text files with a specific markup language, that are then auto-converted
to HTML and PDF. The tools for this conversion are included in the
source distribution, and the translation can be as simple as doing
&#8220;make html pdf&#8221; in the doc folder.
Thus the documentation source files must be in the same format and
style as other <em>.txt files in the lammps/doc/src directory for similar
commands and styles; use one or more of them as a starting point.
A description of the markup can also be found in
lammps/doc/utils/txt2html/README.html
As appropriate, the text files can include links to equations
(see doc/Eqs/</em>.tex for examples, we auto-create the associated JPG
files), or figures (see doc/JPG for examples), or even additional PDF
files with further details (see doc/PDF for examples). The doc page
@ -1159,14 +1205,19 @@ bottom of doc/fix_nh.txt for examples and the earlier part of the same
file for how to format the cite itself. The &#8220;Restrictions&#8221; section of
the doc page should indicate that your command is only available if
LAMMPS is built with the appropriate USER-MISC or USER-FOO package.
See other user package doc files for examples of how to do this. The
txt2html tool we use to convert to HTML can be downloaded from <a class="reference external" href="http://www.sandia.gov/~sjplimp/download.html">this site</a>, so you can perform
the HTML conversion yourself to proofread your doc page.</li>
<li>For a new package (or even a single command) you can include one or
more example scripts. These should run in no more than 1 minute, even
on a single processor, and not require large data files as input. See
directories under examples/USER for examples of input scripts other
users provided for their packages.</li>
See other user package doc files for examples of how to do this. The
prerequiste for building the HTML format files are Python 3.x and
virtualenv, the requirement for generating the PDF format manual
is the <a class="reference external" href="http://www.htmldoc.org/">htmldoc</a> software. Please run at least
&#8220;make html&#8221; and carefully inspect and proofread the resuling HTML format
doc page before submitting your code.</li>
<li>For a new package (or even a single command) you should include one or
more example scripts demonstrating its use. These should run in no
more than a couple minutes, even on a single processor, and not require
large data files as input. See directories under examples/USER for
examples of input scripts other users provided for their packages.
These example inputs are also required for validating memory accesses
and testing for memory leaks with valgrind</li>
<li>If there is a paper of yours describing your feature (either the
algorithm/science behind the feature itself, or its initial usage, or
its implementation in LAMMPS), you can add the citation to the <a href="#id12"><span class="problematic" id="id13">*</span></a>.cpp
@ -1183,10 +1234,9 @@ usage. That kind of citation should just be in the doc page you
provide.</li>
</ul>
<p>Finally, as a general rule-of-thumb, the more clear and
self-explanatory you make your doc and README files, and the easier
you make it for people to get started, e.g. by providing example
scripts, the more likely it is that users will try out your new
feature.</p>
self-explanatory you make your documentation and README files, and the
easier you make it for people to get started, e.g. by providing example
scripts, the more likely it is that users will try out your new feature.</p>
<p id="foo"><strong>(Foo)</strong> Foo, Morefoo, and Maxfoo, J of Classic Potentials, 75, 345 (1997).</p>
</div>
</div>

2
doc/html/Section_packages.html
View File

@ -493,7 +493,7 @@ name links to a sub-section below with more details.</p>
<td>python</td>
<td>lib/python</td>
</tr>
<tr class="row-even"><td><a class="reference internal" href="Section_tools.html#reax"><span class="std std-ref">REAX</span></a></td>
<tr class="row-even"><td><a class="reference internal" href="#reax"><span class="std std-ref">REAX</span></a></td>
<td>ReaxFF potential</td>
<td>Aidan Thompson (Sandia)</td>
<td><a class="reference internal" href="pair_reax.html"><span class="doc">pair_style reax</span></a></td>

2
doc/html/compute_rdf.html
View File

@ -227,7 +227,7 @@ atoms of type jtypeN.</p>
to a file is to use the <a class="reference internal" href="fix_ave_time.html"><span class="doc">fix ave/time</span></a> command, for
example:</p>
<div class="highlight-default"><div class="highlight"><pre><span></span><span class="n">compute</span> <span class="n">myRDF</span> <span class="nb">all</span> <span class="n">rdf</span> <span class="mi">50</span>
<span class="n">fix</span> <span class="mi">1</span> <span class="nb">all</span> <span class="n">ave</span><span class="o">/</span><span class="n">time</span> <span class="mi">100</span> <span class="mi">1</span> <span class="mi">100</span> <span class="n">c_myRDF</span> <span class="n">file</span> <span class="n">tmp</span><span class="o">.</span><span class="n">rdf</span> <span class="n">mode</span> <span class="n">vector</span>
<span class="n">fix</span> <span class="mi">1</span> <span class="nb">all</span> <span class="n">ave</span><span class="o">/</span><span class="n">time</span> <span class="mi">100</span> <span class="mi">1</span> <span class="mi">100</span> <span class="n">c_myRDF</span><span class="p">[</span><span class="o">*</span><span class="p">]</span> <span class="n">file</span> <span class="n">tmp</span><span class="o">.</span><span class="n">rdf</span> <span class="n">mode</span> <span class="n">vector</span>
</pre></div>
</div>
<p><strong>Output info:</strong></p>

2
doc/html/fix_lb_pc.html
View File

@ -146,7 +146,7 @@
<h2>Description</h2>
<p>Update the positions and velocities of the individual particles
described by <em>group-ID</em>, experiencing velocity-dependent hydrodynamic
forces, using the integration algorithm described in <a class="reference internal" href="fix_lb_viscous.html#mackay"><span class="std std-ref">Mackay et al.</span></a>. This integration algorithm should only be used if a
forces, using the integration algorithm described in <a class="reference internal" href="#mackay"><span class="std std-ref">Mackay et al.</span></a>. This integration algorithm should only be used if a
user-specified value for the force-coupling constant used in <a class="reference internal" href="fix_lb_fluid.html"><span class="doc">fix lb/fluid</span></a> has been set; do not use this integration
algorithm if the force coupling constant has been set by default.</p>
</div>

4
doc/html/pair_thole.html
View File

@ -189,8 +189,8 @@ short distances by a function</p>
\frac{s_{ij} r_{ij} }{2} \right)
\exp \left( - s_{ij} r_{ij} \right) \end{equation}\]</div>
<p>This function results from an adaptation to point charges
<a class="reference internal" href="tutorial_drude.html#noskov"><span class="std std-ref">(Noskov)</span></a> of the dipole screening scheme originally proposed
by <a class="reference internal" href="tutorial_drude.html#thole"><span class="std std-ref">Thole</span></a>. The scaling coefficient <span class="math">\(s_{ij}\)</span> is determined
<a class="reference internal" href="#noskov"><span class="std std-ref">(Noskov)</span></a> of the dipole screening scheme originally proposed
by <a class="reference internal" href="#thole"><span class="std std-ref">Thole</span></a>. The scaling coefficient <span class="math">\(s_{ij}\)</span> is determined
by the polarizability of the atoms, <span class="math">\(\alpha_i\)</span>, and by a Thole
damping parameter <span class="math">\(a\)</span>. This Thole damping parameter usually takes
a value of 2.6, but in certain force fields the value can depend upon

4
doc/src/Manual.txt
View File

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

4
doc/src/Section_howto.txt
View File

@ -1698,9 +1698,9 @@ pressure"_compute_pressure.html compute to calculate a current
pressure. By default, this compute is created with a simple "compute
temp"_compute_temp.html (see the last argument of the "compute
pressure"_compute_pressure.html command), which is used to calculated
the kinetic componenet of the pressure. The barostatting fixes can
the kinetic component of the pressure. The barostatting fixes can
also use temperature computes that remove bias for the purpose of
computing the kinetic componenet which contributes to the current
computing the kinetic component which contributes to the current
pressure. See the doc pages for the individual fixes and for the
"fix_modify"_fix_modify.html command for instructions on how to assign
a temperature or pressure compute to a barostatting fix.

130
doc/src/Section_modify.txt
View File

@ -621,13 +621,37 @@ then be accessed by variables) was discussed
10.15 Submitting new features for inclusion in LAMMPS :link(mod_15),h4
We encourage users to submit new features to "the
developers"_http://lammps.sandia.gov/authors.html that they add to
LAMMPS, especially if you think they will be of interest to other
users. The preferred way to do this is via GitHub. Once you have
prepared the content described below, see "this
tutorial"_tutorial_github.html for instructions on how to submit
your changes or new files.
We encourage users to submit new features or modifications for
LAMMPS to "the core developers"_http://lammps.sandia.gov/authors.html
so they can be added to the LAMMPS distribution. The preferred way to
manage and coordinate this is as of Fall 2016 via the LAMMPS project on
"GitHub"_https://github.com/lammps/lammps. An alternative is to contact
the LAMMPS developers or the indicated developer of a package or feature
directly and send in your contribution via e-mail.
For any larger modifications or programming project, you are encouraged
to contact the LAMMPS developers ahead of time, in order to discuss
implementation strategies and coding guidelines, that will make it
easier to integrate your contribution and result in less work for
everybody involved. You are also encouraged to search through the list
of "open issues on GitHub"_https://github.com/lammps/lammps/issues and
submit a new issue for a planned feature, so you would not duplicate
the work of others (and possibly get scooped by them) or have your work
duplicated by others.
How quickly your contribution will be integrated
depends largely on how much effort it will cause to integrate and test
it, how much it requires changes to the core codebase, and of how much
interest it is to the larger LAMMPS community. Please see below for a
checklist of typical requirements. Once you have prepared everything,
see "this tutorial"_tutorial_github.html for instructions on how to
submit your changes or new files through a GitHub pull request. If you
prefer to submit patches or full files, you should first make certain,
that your code works correctly with the latest patch-level version of
LAMMPS and contains all bugfixes from it. Then create a gzipped tar
file of all changed or added files or a corresponding patch file using
'diff -u' or 'diff -c' and compress it with gzip. Please only use
gzip compression, as this works well on all platforms.
If the new features/files are broadly useful we may add them as core
files to LAMMPS or as part of a "standard
@ -677,21 +701,35 @@ not want to make those changes. An example of a trivial change is
making a parent-class method "virtual" when you derive a new child
class from it.
Here are the steps you need to follow to submit a single file or user
package for our consideration. Following these steps will save both
you and us time. See existing files in packages in the src dir for
examples.
Here is a checklist of steps you need to follow to submit a single file
or user package for our consideration. Following these steps will save
both you and us time. See existing files in packages in the src dir for
examples. If you are uncertain, please ask.
All source files you provide must compile with the most current
version of LAMMPS. :ulb,l
version of LAMMPS with multiple configurations. In particular you
need to test compiling LAMMPS from scratch with -DLAMMPS_BIGBIG
set in addition to the default -DLAMMPS_SMALLBIG setting. Your code
will need to work correctly in serial and in parallel using MPI. :ulb,l
If you want your file(s) to be added to main LAMMPS or one of its
standard packages, then it needs to be written in a style compatible
with other LAMMPS source files. This is so the developers can
understand it and hopefully maintain it. This basically means that
the code accesses data structures, performs its operations, and is
formatted similar to other LAMMPS source files, including the use of
the error class for error and warning messages. :l
For consistency with the rest of LAMMPS and especially, if you want
your contribution(s) to be added to main LAMMPS code or one of its
standard packages, it needs to be written in a style compatible with
other LAMMPS source files. This means: 2-character indentation per
level, [no tabs], no lines over 80 characters. I/O is done via
the C-style stdio library, class header files should not import any
system headers outside <stdio.h>, STL containers should be avoided
in headers, and forward declarations used where possible or needed.
All added code should be placed into the LAMMPS_NS namespace or a
sub-namespace; global or static variables should be avoided, as they
conflict with the modular nature of LAMMPS and the C++ class structure.
Header files must [not] import namespaces with {using}.
This all is so the developers can more easily understand, integrate,
and maintain your contribution and reduce conflicts with other parts
of LAMMPS. This basically means that the code accesses data
structures, performs its operations, and is formatted similar to other
LAMMPS source files, including the use of the error class for error
and warning messages. :l
If you want your contribution to be added as a user-contributed
feature, and it's a single file (actually a *.cpp and *.h file) it can
@ -722,13 +760,21 @@ same stylistic format and syntax as other LAMMPS files, though that
would be nice for developers as well as users who try to read your
code. :l
You must also create a documentation file for each new command or
style you are adding to LAMMPS. This will be one file for a
single-file feature. For a package, it might be several files. These
are simple text files which we auto-convert to HTML. Thus they must
be in the same format as other *.txt files in the lammps/doc directory
for similar commands and styles; use one or more of them as a starting
point. As appropriate, the text files can include links to equations
You [must] also create a [documentation] file for each new command or
style you are adding to LAMMPS. For simplicity and convenience, the
documentation of groups of closely related commands or styles may be
combined into a single file. This will be one file for a single-file
feature. For a package, it might be several files. These are simple
text files with a specific markup language, that are then auto-converted
to HTML and PDF. The tools for this conversion are included in the
source distribution, and the translation can be as simple as doing
"make html pdf" in the doc folder.
Thus the documentation source files must be in the same format and
style as other *.txt files in the lammps/doc/src directory for similar
commands and styles; use one or more of them as a starting point.
A description of the markup can also be found in
lammps/doc/utils/txt2html/README.html
As appropriate, the text files can include links to equations
(see doc/Eqs/*.tex for examples, we auto-create the associated JPG
files), or figures (see doc/JPG for examples), or even additional PDF
files with further details (see doc/PDF for examples). The doc page
@ -737,16 +783,20 @@ bottom of doc/fix_nh.txt for examples and the earlier part of the same
file for how to format the cite itself. The "Restrictions" section of
the doc page should indicate that your command is only available if
LAMMPS is built with the appropriate USER-MISC or USER-FOO package.
See other user package doc files for examples of how to do this. The
txt2html tool we use to convert to HTML can be downloaded from "this
site"_http://www.sandia.gov/~sjplimp/download.html, so you can perform
the HTML conversion yourself to proofread your doc page. :l
See other user package doc files for examples of how to do this. The
prerequiste for building the HTML format files are Python 3.x and
virtualenv, the requirement for generating the PDF format manual
is the "htmldoc"_http://www.htmldoc.org/ software. Please run at least
"make html" and carefully inspect and proofread the resuling HTML format
doc page before submitting your code. :l
For a new package (or even a single command) you can include one or
more example scripts. These should run in no more than 1 minute, even
on a single processor, and not require large data files as input. See
directories under examples/USER for examples of input scripts other
users provided for their packages. :l
For a new package (or even a single command) you should include one or
more example scripts demonstrating its use. These should run in no
more than a couple minutes, even on a single processor, and not require
large data files as input. See directories under examples/USER for
examples of input scripts other users provided for their packages.
These example inputs are also required for validating memory accesses
and testing for memory leaks with valgrind :l
If there is a paper of yours describing your feature (either the
algorithm/science behind the feature itself, or its initial usage, or
@ -761,13 +811,13 @@ should only use this for a paper you or your group authored.
E.g. adding a cite in the code for a paper by Nose and Hoover if you
write a fix that implements their integrator is not the intended
usage. That kind of citation should just be in the doc page you
provide. :l,ule
provide. :l
:ule
Finally, as a general rule-of-thumb, the more clear and
self-explanatory you make your doc and README files, and the easier
you make it for people to get started, e.g. by providing example
scripts, the more likely it is that users will try out your new
feature.
self-explanatory you make your documentation and README files, and the
easier you make it for people to get started, e.g. by providing example
scripts, the more likely it is that users will try out your new feature.
:line
:line

2
doc/src/compute_rdf.txt
View File

@ -110,7 +110,7 @@ to a file is to use the "fix ave/time"_fix_ave_time.html command, for
example:
compute myRDF all rdf 50
fix 1 all ave/time 100 1 100 c_myRDF file tmp.rdf mode vector :pre
fix 1 all ave/time 100 1 100 c_myRDF\[*\] file tmp.rdf mode vector :pre
[Output info:]

135
src/KOKKOS/pair_eam_kokkos.cpp
View File

@ -340,81 +340,64 @@ void PairEAMKokkos<DeviceType>::array2spline()
rdr = 1.0/dr;
rdrho = 1.0/drho;
tdual_ffloat4 k_frho_spline_a = tdual_ffloat4("pair:frho_a",nfrho,nrho+1);
tdual_ffloat4 k_rhor_spline_a = tdual_ffloat4("pair:rhor_a",nrhor,nr+1);
tdual_ffloat4 k_z2r_spline_a = tdual_ffloat4("pair:z2r_a",nz2r,nr+1);
tdual_ffloat4 k_frho_spline_b = tdual_ffloat4("pair:frho_b",nfrho,nrho+1);
tdual_ffloat4 k_rhor_spline_b = tdual_ffloat4("pair:rhor_b",nrhor,nr+1);
tdual_ffloat4 k_z2r_spline_b = tdual_ffloat4("pair:z2r_b",nz2r,nr+1);
tdual_ffloat_2d_n7 k_frho_spline = tdual_ffloat_2d_n7("pair:frho",nfrho,nrho+1);
tdual_ffloat_2d_n7 k_rhor_spline = tdual_ffloat_2d_n7("pair:rhor",nrhor,nr+1);
tdual_ffloat_2d_n7 k_z2r_spline = tdual_ffloat_2d_n7("pair:z2r",nz2r,nr+1);
t_host_ffloat4 h_frho_spline_a = k_frho_spline_a.h_view;
t_host_ffloat4 h_rhor_spline_a = k_rhor_spline_a.h_view;
t_host_ffloat4 h_z2r_spline_a = k_z2r_spline_a.h_view;
t_host_ffloat4 h_frho_spline_b = k_frho_spline_b.h_view;
t_host_ffloat4 h_rhor_spline_b = k_rhor_spline_b.h_view;
t_host_ffloat4 h_z2r_spline_b = k_z2r_spline_b.h_view;
t_host_ffloat_2d_n7 h_frho_spline = k_frho_spline.h_view;
t_host_ffloat_2d_n7 h_rhor_spline = k_rhor_spline.h_view;
t_host_ffloat_2d_n7 h_z2r_spline = k_z2r_spline.h_view;
for (int i = 0; i < nfrho; i++)
interpolate(nrho,drho,frho[i],h_frho_spline_a,h_frho_spline_b,i);
k_frho_spline_a.template modify<LMPHostType>();
k_frho_spline_a.template sync<DeviceType>();
k_frho_spline_b.template modify<LMPHostType>();
k_frho_spline_b.template sync<DeviceType>();
interpolate(nrho,drho,frho[i],h_frho_spline,i);
k_frho_spline.template modify<LMPHostType>();
k_frho_spline.template sync<DeviceType>();
for (int i = 0; i < nrhor; i++)
interpolate(nr,dr,rhor[i],h_rhor_spline_a,h_rhor_spline_b,i);
k_rhor_spline_a.template modify<LMPHostType>();
k_rhor_spline_a.template sync<DeviceType>();
k_rhor_spline_b.template modify<LMPHostType>();
k_rhor_spline_b.template sync<DeviceType>();
interpolate(nr,dr,rhor[i],h_rhor_spline,i);
k_rhor_spline.template modify<LMPHostType>();
k_rhor_spline.template sync<DeviceType>();
for (int i = 0; i < nz2r; i++)
interpolate(nr,dr,z2r[i],h_z2r_spline_a,h_z2r_spline_b,i);
k_z2r_spline_a.template modify<LMPHostType>();
k_z2r_spline_a.template sync<DeviceType>();
k_z2r_spline_b.template modify<LMPHostType>();
k_z2r_spline_b.template sync<DeviceType>();
d_frho_spline_a = k_frho_spline_a.d_view;
d_rhor_spline_a = k_rhor_spline_a.d_view;
d_z2r_spline_a = k_z2r_spline_a.d_view;
d_frho_spline_b = k_frho_spline_b.d_view;
d_rhor_spline_b = k_rhor_spline_b.d_view;
d_z2r_spline_b = k_z2r_spline_b.d_view;
interpolate(nr,dr,z2r[i],h_z2r_spline,i);
k_z2r_spline.template modify<LMPHostType>();
k_z2r_spline.template sync<DeviceType>();
d_frho_spline = k_frho_spline.d_view;
d_rhor_spline = k_rhor_spline.d_view;
d_z2r_spline = k_z2r_spline.d_view;
}
/* ---------------------------------------------------------------------- */
template<class DeviceType>
void PairEAMKokkos<DeviceType>::interpolate(int n, double delta, double *f, t_host_ffloat4 h_spline_a, t_host_ffloat4 h_spline_b, int i)
void PairEAMKokkos<DeviceType>::interpolate(int n, double delta, double *f, t_host_ffloat_2d_n7 h_spline, int i)
{
for (int m = 1; m <= n; m++) h_spline_b(i,m).w = f[m];
for (int m = 1; m <= n; m++) h_spline(i,m,6) = f[m];
h_spline_b(i,1).z = h_spline_b(i,2).w - h_spline_b(i,1).w;
h_spline_b(i,2).z = 0.5 * (h_spline_b(i,3).w-h_spline_b(i,1).w);
h_spline_b(i,n-1).z = 0.5 * (h_spline_b(i,n).w-h_spline_b(i,n-2).w);
h_spline_b(i,n).z = h_spline_b(i,n).w - h_spline_b(i,n-1).w;
h_spline(i,1,5) = h_spline(i,2,6) - h_spline(i,1,6);
h_spline(i,2,5) = 0.5 * (h_spline(i,3,6)-h_spline(i,1,6));
h_spline(i,n-1,5) = 0.5 * (h_spline(i,n,6)-h_spline(i,n-2,6));
h_spline(i,n,5) = h_spline(i,n,6) - h_spline(i,n-1,6);
for (int m = 3; m <= n-2; m++)
h_spline_b(i,m).z = ((h_spline_b(i,m-2).w-h_spline_b(i,m+2).w) +
8.0*(h_spline_b(i,m+1).w-h_spline_b(i,m-1).w)) / 12.0;
h_spline(i,m,5) = ((h_spline(i,m-2,6)-h_spline(i,m+2,6)) +
8.0*(h_spline(i,m+1,6)-h_spline(i,m-1,6))) / 12.0;
for (int m = 1; m <= n-1; m++) {
h_spline_b(i,m).y = 3.0*(h_spline_b(i,m+1).w-h_spline_b(i,m).w) -
2.0*h_spline_b(i,m).z - h_spline_b(i,m+1).z;
h_spline_b(i,m).x = h_spline_b(i,m).z + h_spline_b(i,m+1).z -
2.0*(h_spline_b(i,m+1).w-h_spline_b(i,m).w);
h_spline(i,m,4) = 3.0*(h_spline(i,m+1,6)-h_spline(i,m,6)) -
2.0*h_spline(i,m,5) - h_spline(i,m+1,5);
h_spline(i,m,3) = h_spline(i,m,5) + h_spline(i,m+1,5) -
2.0*(h_spline(i,m+1,6)-h_spline(i,m,6));
}
h_spline_b(i,n).y = 0.0;
h_spline_b(i,n).x = 0.0;
h_spline(i,n,4) = 0.0;
h_spline(i,n,3) = 0.0;
for (int m = 1; m <= n; m++) {
h_spline_a(i,m).z = h_spline_b(i,m).z/delta;
h_spline_a(i,m).y = 2.0*h_spline_b(i,m).y/delta;
h_spline_a(i,m).x = 3.0*h_spline_b(i,m).x/delta;
h_spline(i,m,2) = h_spline(i,m,5)/delta;
h_spline(i,m,1) = 2.0*h_spline(i,m,4)/delta;
h_spline(i,m,0) = 3.0*h_spline(i,m,3)/delta;
}
}
@ -558,13 +541,12 @@ void PairEAMKokkos<DeviceType>::operator()(TagPairEAMKernelA<NEIGHFLAG,NEWTON_PA
p -= m;
p = MIN(p,1.0);
const int d_type2rhor_ji = d_type2rhor(jtype,itype);
const F_FLOAT4 rhor = d_rhor_spline_b(d_type2rhor_ji,m);
rhotmp += ((rhor.x*p + rhor.y)*p + rhor.z)*p + rhor.w;
rhotmp += ((d_rhor_spline(d_type2rhor_ji,m,3)*p + d_rhor_spline(d_type2rhor_ji,m,4))*p +
d_rhor_spline(d_type2rhor_ji,m,5))*p + d_rhor_spline(d_type2rhor_ji,m,6);
if (NEWTON_PAIR || j < nlocal) {
const int d_type2rhor_ij = d_type2rhor(itype,jtype);
const F_FLOAT4 rhor = d_rhor_spline_b(d_type2rhor_ij,m);
rho[j] += ((rhor.x*p + rhor.y)*p + rhor.z)*p + rhor.w;
rho[j] += ((d_rhor_spline(d_type2rhor_ij,m,3)*p + d_rhor_spline(d_type2rhor_ij,m,4))*p +
d_rhor_spline(d_type2rhor_ij,m,5))*p + d_rhor_spline(d_type2rhor_ij,m,6);
}
}
@ -593,11 +575,10 @@ void PairEAMKokkos<DeviceType>::operator()(TagPairEAMKernelB<EFLAG>, const int &
p -= m;
p = MIN(p,1.0);
const int d_type2frho_i = d_type2frho[itype];
const F_FLOAT4 frho = d_frho_spline_a(d_type2frho_i,m);
d_fp[i] = (frho.x*p + frho.y)*p + frho.z;
d_fp[i] = (d_frho_spline(d_type2frho_i,m,0)*p + d_frho_spline(d_type2frho_i,m,1))*p + d_frho_spline(d_type2frho_i,m,2);
if (EFLAG) {
const F_FLOAT4 frho_b = d_frho_spline_b(d_type2frho_i,m);
F_FLOAT phi = ((frho_b.x*p + frho_b.y)*p + frho_b.z)*p + frho_b.w;
F_FLOAT phi = ((d_frho_spline(d_type2frho_i,m,3)*p + d_frho_spline(d_type2frho_i,m,4))*p +
d_frho_spline(d_type2frho_i,m,5))*p + d_frho_spline(d_type2frho_i,m,6);
if (d_rho[i] > rhomax) phi += d_fp[i] * (d_rho[i]-rhomax);
if (eflag_global) ev.evdwl += phi;
if (eflag_atom) d_eatom[i] += phi;
@ -651,8 +632,8 @@ void PairEAMKokkos<DeviceType>::operator()(TagPairEAMKernelAB<EFLAG>, const int
p -= m;
p = MIN(p,1.0);
const int d_type2rhor_ji = d_type2rhor(jtype,itype);
const F_FLOAT4 rhor = d_rhor_spline_b(d_type2rhor_ji,m);
rhotmp += ((rhor.x*p + rhor.y)*p + rhor.z)*p + rhor.w;
rhotmp += ((d_rhor_spline(d_type2rhor_ji,m,3)*p + d_rhor_spline(d_type2rhor_ji,m,4))*p +
d_rhor_spline(d_type2rhor_ji,m,5))*p + d_rhor_spline(d_type2rhor_ji,m,6);
}
}
@ -669,11 +650,10 @@ void PairEAMKokkos<DeviceType>::operator()(TagPairEAMKernelAB<EFLAG>, const int
p -= m;
p = MIN(p,1.0);
const int d_type2frho_i = d_type2frho[itype];
const F_FLOAT4 frho = d_frho_spline_a(d_type2frho_i,m);
d_fp[i] = (frho.x*p + frho.y)*p + frho.z;
d_fp[i] = (d_frho_spline(d_type2frho_i,m,0)*p + d_frho_spline(d_type2frho_i,m,1))*p + d_frho_spline(d_type2frho_i,m,2);
if (EFLAG) {
const F_FLOAT4 frho_b = d_frho_spline_b(d_type2frho_i,m);
F_FLOAT phi = ((frho_b.x*p + frho_b.y)*p + frho_b.z)*p + frho_b.w;
F_FLOAT phi = ((d_frho_spline(d_type2frho_i,m,3)*p + d_frho_spline(d_type2frho_i,m,4))*p +
d_frho_spline(d_type2frho_i,m,5))*p + d_frho_spline(d_type2frho_i,m,6);
if (d_rho[i] > rhomax) phi += d_fp[i] * (d_rho[i]-rhomax);
if (eflag_global) ev.evdwl += phi;
if (eflag_atom) d_eatom[i] += phi;
@ -742,18 +722,16 @@ void PairEAMKokkos<DeviceType>::operator()(TagPairEAMKernelC<NEIGHFLAG,NEWTON_PA
// hence embed' = Fi(sum rho_ij) rhojp + Fj(sum rho_ji) rhoip
const int d_type2rhor_ij = d_type2rhor(itype,jtype);
const F_FLOAT4 rhor_ij = d_rhor_spline_a(d_type2rhor_ij,m);
const F_FLOAT rhoip = (rhor_ij.x*p + rhor_ij.y)*p + rhor_ij.z;
const F_FLOAT rhoip = (d_rhor_spline(d_type2rhor_ij,m,0)*p + d_rhor_spline(d_type2rhor_ij,m,1))*p +
d_rhor_spline(d_type2rhor_ij,m,2);
const int d_type2rhor_ji = d_type2rhor(jtype,itype);
const F_FLOAT4 rhor_ji = d_rhor_spline_a(d_type2rhor_ji,m);
const F_FLOAT rhojp = (rhor_ji.x*p + rhor_ji.y)*p + rhor_ji.z;
const F_FLOAT rhojp = (d_rhor_spline(d_type2rhor_ji,m,0)*p + d_rhor_spline(d_type2rhor_ji,m,1))*p +
d_rhor_spline(d_type2rhor_ji,m,2);
const int d_type2z2r_ij = d_type2z2r(itype,jtype);
const F_FLOAT4 z2r_a = d_z2r_spline_a(d_type2z2r_ij,m);
const F_FLOAT z2p = (z2r_a.x*p + z2r_a.y)*p + z2r_a.z;
const F_FLOAT4 z2r_b = d_z2r_spline_b(d_type2z2r_ij,m);
const F_FLOAT z2 = ((z2r_b.x*p + z2r_b.y)*p + z2r_b.z)*p + z2r_b.w;
const F_FLOAT z2p = (d_z2r_spline(d_type2z2r_ij,m,0)*p + d_z2r_spline(d_type2z2r_ij,m,1))*p +
d_z2r_spline(d_type2z2r_ij,m,2);
const F_FLOAT z2 = ((d_z2r_spline(d_type2z2r_ij,m,3)*p + d_z2r_spline(d_type2z2r_ij,m,4))*p +
d_z2r_spline(d_type2z2r_ij,m,5))*p + d_z2r_spline(d_type2z2r_ij,m,6);
const F_FLOAT recip = 1.0/r;
const F_FLOAT phi = z2*recip;
@ -894,5 +872,4 @@ template class PairEAMKokkos<LMPDeviceType>;
#ifdef KOKKOS_HAVE_CUDA
template class PairEAMKokkos<LMPHostType>;
#endif
}
}

14
src/KOKKOS/pair_eam_kokkos.h
View File

@ -136,14 +136,14 @@ class PairEAMKokkos : public PairEAM {
DAT::t_int_2d_randomread d_type2rhor;
DAT::t_int_2d_randomread d_type2z2r;
typedef Kokkos::DualView<F_FLOAT4**,Kokkos::LayoutLeft,DeviceType> tdual_ffloat4;
typedef typename tdual_ffloat4::t_dev_const_randomread t_ffloat4_randomread;
typedef typename tdual_ffloat4::t_host t_host_ffloat4;
typedef Kokkos::DualView<F_FLOAT**[7],Kokkos::LayoutRight,DeviceType> tdual_ffloat_2d_n7;
typedef typename tdual_ffloat_2d_n7::t_dev_const_randomread t_ffloat_2d_n7_randomread;
typedef typename tdual_ffloat_2d_n7::t_host t_host_ffloat_2d_n7;
t_ffloat4_randomread d_frho_spline_a, d_frho_spline_b;
t_ffloat4_randomread d_rhor_spline_a, d_rhor_spline_b;
t_ffloat4_randomread d_z2r_spline_a, d_z2r_spline_b;
void interpolate(int, double, double *, t_host_ffloat4, t_host_ffloat4, int);
t_ffloat_2d_n7_randomread d_frho_spline;
t_ffloat_2d_n7_randomread d_rhor_spline;
t_ffloat_2d_n7_randomread d_z2r_spline;
void interpolate(int, double, double *, t_host_ffloat_2d_n7, int);
virtual void file2array();
void array2spline();

115
src/compute_omega_chunk.cpp
View File

@ -18,11 +18,14 @@
#include "modify.h"
#include "compute_chunk_atom.h"
#include "domain.h"
#include "math_extra.h"
#include "memory.h"
#include "error.h"
using namespace LAMMPS_NS;
#define EPSILON 1.0e-6
/* ---------------------------------------------------------------------- */
ComputeOmegaChunk::ComputeOmegaChunk(LAMMPS *lmp, int narg, char **arg) :
@ -192,49 +195,97 @@ void ComputeOmegaChunk::compute_array()
MPI_Allreduce(&angmom[0][0],&angmomall[0][0],3*nchunk,
MPI_DOUBLE,MPI_SUM,world);
// compute omega for each chunk from L = Iw, inverting I to solve for w
// compute omega for each chunk
double ione[3][3],inverse[3][3];
double determinant,invdeterminant;
double idiag[3],ex[3],ey[3],ez[3],cross[3];
double ione[3][3],inverse[3][3],evectors[3][3];
double *iall,*mall;
for (m = 0; m < nchunk; m++) {
ione[0][0] = inertiaall[m][0];
ione[1][1] = inertiaall[m][1];
ione[2][2] = inertiaall[m][2];
ione[0][1] = inertiaall[m][3];
ione[1][2] = inertiaall[m][4];
ione[0][2] = inertiaall[m][5];
ione[1][0] = ione[0][1];
ione[2][1] = ione[1][2];
ione[2][0] = ione[0][2];
inverse[0][0] = ione[1][1]*ione[2][2] - ione[1][2]*ione[2][1];
inverse[0][1] = -(ione[0][1]*ione[2][2] - ione[0][2]*ione[2][1]);
inverse[0][2] = ione[0][1]*ione[1][2] - ione[0][2]*ione[1][1];
// determinant = triple product of rows of inertia matrix
inverse[1][0] = -(ione[1][0]*ione[2][2] - ione[1][2]*ione[2][0]);
inverse[1][1] = ione[0][0]*ione[2][2] - ione[0][2]*ione[2][0];
inverse[1][2] = -(ione[0][0]*ione[1][2] - ione[0][2]*ione[1][0]);
iall = &inertiaall[m][0];
determinant = iall[0] * (iall[1]*iall[2] - iall[4]*iall[4]) +
iall[3] * (iall[4]*iall[5] - iall[3]*iall[2]) +
iall[5] * (iall[3]*iall[4] - iall[1]*iall[5]);
inverse[2][0] = ione[1][0]*ione[2][1] - ione[1][1]*ione[2][0];
inverse[2][1] = -(ione[0][0]*ione[2][1] - ione[0][1]*ione[2][0]);
inverse[2][2] = ione[0][0]*ione[1][1] - ione[0][1]*ione[1][0];
ione[0][0] = iall[0];
ione[1][1] = iall[1];
ione[2][2] = iall[2];
ione[0][1] = ione[1][0] = iall[3];
ione[1][2] = ione[2][1] = iall[4];
ione[0][2] = ione[2][0] = iall[5];
double determinant = ione[0][0]*ione[1][1]*ione[2][2] +
ione[0][1]*ione[1][2]*ione[2][0] + ione[0][2]*ione[1][0]*ione[2][1] -
ione[0][0]*ione[1][2]*ione[2][1] - ione[0][1]*ione[1][0]*ione[2][2] -
ione[2][0]*ione[1][1]*ione[0][2];
// non-singular I matrix
// use L = Iw, inverting I to solve for w
if (determinant > 0.0)
if (determinant > EPSILON) {
inverse[0][0] = ione[1][1]*ione[2][2] - ione[1][2]*ione[2][1];
inverse[0][1] = -(ione[0][1]*ione[2][2] - ione[0][2]*ione[2][1]);
inverse[0][2] = ione[0][1]*ione[1][2] - ione[0][2]*ione[1][1];
inverse[1][0] = -(ione[1][0]*ione[2][2] - ione[1][2]*ione[2][0]);
inverse[1][1] = ione[0][0]*ione[2][2] - ione[0][2]*ione[2][0];
inverse[1][2] = -(ione[0][0]*ione[1][2] - ione[0][2]*ione[1][0]);
inverse[2][0] = ione[1][0]*ione[2][1] - ione[1][1]*ione[2][0];
inverse[2][1] = -(ione[0][0]*ione[2][1] - ione[0][1]*ione[2][0]);
inverse[2][2] = ione[0][0]*ione[1][1] - ione[0][1]*ione[1][0];
invdeterminant = 1.0/determinant;
for (i = 0; i < 3; i++)
for (j = 0; j < 3; j++)
inverse[i][j] /= determinant;
inverse[i][j] *= invdeterminant;
mall = &angmomall[m][0];
omega[m][0] = inverse[0][0]*mall[0] + inverse[0][1]*mall[1] +
inverse[0][2]*mall[2];
omega[m][1] = inverse[1][0]*mall[0] + inverse[1][1]*mall[1] +
inverse[1][2]*mall[2];
omega[m][2] = inverse[2][0]*mall[0] + inverse[2][1]*mall[1] +
inverse[2][2]*mall[2];
omega[m][0] = inverse[0][0]*angmom[m][0] + inverse[0][1]*angmom[m][1] +
inverse[0][2]*angmom[m][2];
omega[m][1] = inverse[1][0]*angmom[m][0] + inverse[1][1]*angmom[m][1] +
inverse[1][2]*angmom[m][2];
omega[m][2] = inverse[2][0]*angmom[m][0] + inverse[2][1]*angmom[m][1] +
inverse[2][2]*angmom[i][2];
// handle each (nearly) singular I matrix
// due to 2-atom chunk or linear molecule
// use jacobi() and angmom_to_omega() to calculate valid omega
} else {
int ierror = MathExtra::jacobi(ione,idiag,evectors);
if (ierror) error->all(FLERR,
"Insufficient Jacobi rotations for omega/chunk");
ex[0] = evectors[0][0];
ex[1] = evectors[1][0];
ex[2] = evectors[2][0];
ey[0] = evectors[0][1];
ey[1] = evectors[1][1];
ey[2] = evectors[2][1];
ez[0] = evectors[0][2];
ez[1] = evectors[1][2];
ez[2] = evectors[2][2];
// enforce 3 evectors as a right-handed coordinate system
// flip 3rd vector if needed
MathExtra::cross3(ex,ey,cross);
if (MathExtra::dot3(cross,ez) < 0.0) MathExtra::negate3(ez);
// if any principal moment < scaled EPSILON, set to 0.0
double max;
max = MAX(idiag[0],idiag[1]);
max = MAX(max,idiag[2]);
if (idiag[0] < EPSILON*max) idiag[0] = 0.0;
if (idiag[1] < EPSILON*max) idiag[1] = 0.0;
if (idiag[2] < EPSILON*max) idiag[2] = 0.0;
// calculate omega using diagonalized inertia matrix
MathExtra::angmom_to_omega(&angmomall[m][0],ex,ey,ez,idiag,&omega[m][0]);
}
}
}

71
src/group.cpp
View File

@ -29,6 +29,7 @@
#include "input.h"
#include "variable.h"
#include "dump.h"
#include "math_extra.h"
#include "memory.h"
#include "error.h"
@ -37,6 +38,7 @@
using namespace LAMMPS_NS;
#define MAX_GROUP 32
#define EPSILON 1.0e-6
enum{TYPE,MOLECULE,ID};
enum{LT,LE,GT,GE,EQ,NEQ,BETWEEN};
@ -1662,42 +1664,47 @@ void Group::inertia(int igroup, double *cm, double itensor[3][3], int iregion)
}
/* ----------------------------------------------------------------------
compute angular velocity omega from L = Iw, inverting I to solve for w
compute angular velocity omega from L and I
really not a group/region operation, but L,I were computed for a group/region
diagonalize I instead of inverting it, to allow for a singular matrix
------------------------------------------------------------------------- */
void Group::omega(double *angmom, double inertia[3][3], double *w)
{
double inverse[3][3];
double idiag[3],ex[3],ey[3],ez[3],cross[3];
double evectors[3][3];
int ierror = MathExtra::jacobi(inertia,idiag,evectors);
if (ierror) error->all(FLERR,
"Insufficient Jacobi rotations for group::omega");
inverse[0][0] = inertia[1][1]*inertia[2][2] - inertia[1][2]*inertia[2][1];
inverse[0][1] = -(inertia[0][1]*inertia[2][2] - inertia[0][2]*inertia[2][1]);
inverse[0][2] = inertia[0][1]*inertia[1][2] - inertia[0][2]*inertia[1][1];
inverse[1][0] = -(inertia[1][0]*inertia[2][2] - inertia[1][2]*inertia[2][0]);
inverse[1][1] = inertia[0][0]*inertia[2][2] - inertia[0][2]*inertia[2][0];
inverse[1][2] = -(inertia[0][0]*inertia[1][2] - inertia[0][2]*inertia[1][0]);
inverse[2][0] = inertia[1][0]*inertia[2][1] - inertia[1][1]*inertia[2][0];
inverse[2][1] = -(inertia[0][0]*inertia[2][1] - inertia[0][1]*inertia[2][0]);
inverse[2][2] = inertia[0][0]*inertia[1][1] - inertia[0][1]*inertia[1][0];
double determinant = inertia[0][0]*inertia[1][1]*inertia[2][2] +
inertia[0][1]*inertia[1][2]*inertia[2][0] +
inertia[0][2]*inertia[1][0]*inertia[2][1] -
inertia[0][0]*inertia[1][2]*inertia[2][1] -
inertia[0][1]*inertia[1][0]*inertia[2][2] -
inertia[2][0]*inertia[1][1]*inertia[0][2];
if (determinant > 0.0)
for (int i = 0; i < 3; i++)
for (int j = 0; j < 3; j++)
inverse[i][j] /= determinant;
w[0] = inverse[0][0]*angmom[0] + inverse[0][1]*angmom[1] +
inverse[0][2]*angmom[2];
w[1] = inverse[1][0]*angmom[0] + inverse[1][1]*angmom[1] +
inverse[1][2]*angmom[2];
w[2] = inverse[2][0]*angmom[0] + inverse[2][1]*angmom[1] +
inverse[2][2]*angmom[2];
ex[0] = evectors[0][0];
ex[1] = evectors[1][0];
ex[2] = evectors[2][0];
ey[0] = evectors[0][1];
ey[1] = evectors[1][1];
ey[2] = evectors[2][1];
ez[0] = evectors[0][2];
ez[1] = evectors[1][2];
ez[2] = evectors[2][2];
// enforce 3 evectors as a right-handed coordinate system
// flip 3rd vector if needed
MathExtra::cross3(ex,ey,cross);
if (MathExtra::dot3(cross,ez) < 0.0) MathExtra::negate3(ez);
// if any principal moment < scaled EPSILON, set to 0.0
double max;
max = MAX(idiag[0],idiag[1]);
max = MAX(max,idiag[2]);
if (idiag[0] < EPSILON*max) idiag[0] = 0.0;
if (idiag[1] < EPSILON*max) idiag[1] = 0.0;
if (idiag[2] < EPSILON*max) idiag[2] = 0.0;
// calculate omega using diagonalized inertia matrix
MathExtra::angmom_to_omega(angmom,ex,ey,ez,idiag,w);
}

2
src/version.h
View File

@ -1 +1 @@
#define LAMMPS_VERSION "27 Aug 2016"
#define LAMMPS_VERSION "7 Sep 2016"