Merge branch 'master' into replicate-template-fixes

This commit is contained in:
Axel Kohlmeyer
2020-11-18 10:37:57 -05:00
88 changed files with 4805 additions and 1799 deletions

View File

@ -108,7 +108,7 @@ For bug reports, the next step is that one of the core LAMMPS developers will se
For submitting pull requests, there is a [detailed tutorial](https://lammps.sandia.gov/doc/Howto_github.html) in the LAMMPS manual. Thus only a brief breakdown of the steps is presented here. Please note, that the LAMMPS developers are still reviewing and trying to improve the process. If you are unsure about something, do not hesitate to post a question on the lammps-users mailing list or contact one fo the core LAMMPS developers. For submitting pull requests, there is a [detailed tutorial](https://lammps.sandia.gov/doc/Howto_github.html) in the LAMMPS manual. Thus only a brief breakdown of the steps is presented here. Please note, that the LAMMPS developers are still reviewing and trying to improve the process. If you are unsure about something, do not hesitate to post a question on the lammps-users mailing list or contact one fo the core LAMMPS developers.
Immediately after the submission, the LAMMPS continuing integration server at ci.lammps.org will download your submitted branch and perform a simple compilation test, i.e. will test whether your submitted code can be compiled under various conditions. It will also do a check on whether your included documentation translates cleanly. Whether these tests are successful or fail will be recorded. If a test fails, please inspect the corresponding output on the CI server and take the necessary steps, if needed, so that the code can compile cleanly again. The test will be re-run each the pull request is updated with a push to the remote branch on GitHub. Immediately after the submission, the LAMMPS continuing integration server at ci.lammps.org will download your submitted branch and perform a simple compilation test, i.e. will test whether your submitted code can be compiled under various conditions. It will also do a check on whether your included documentation translates cleanly. Whether these tests are successful or fail will be recorded. If a test fails, please inspect the corresponding output on the CI server and take the necessary steps, if needed, so that the code can compile cleanly again. The test will be re-run each the pull request is updated with a push to the remote branch on GitHub.
Next a LAMMPS core developer will self-assign and do an overall technical assessment of the submission. If you are not yet registered as a LAMMPS collaborator, you will receive an invitation for that. As part of the assesment, the pull request will be categorized with labels. There are two special labels: `needs_work` (indicates that work from the submitter of the pull request is needed) and `work_in_progress` (indicates, that the assigned LAMMPS developer will make changes, if not done by the contributor who made the submit). Next a LAMMPS core developer will self-assign and do an overall technical assessment of the submission. If you are not yet registered as a LAMMPS collaborator, you will receive an invitation for that. As part of the assessment, the pull request will be categorized with labels. There are two special labels: `needs_work` (indicates that work from the submitter of the pull request is needed) and `work_in_progress` (indicates, that the assigned LAMMPS developer will make changes, if not done by the contributor who made the submit).
You may also receive comments and suggestions on the overall submission or specific details and on occasion specific requests for changes as part of the review. If permitted, also additional changes may be pushed into your pull request branch or a pull request may be filed in your LAMMPS fork on GitHub to include those changes. You may also receive comments and suggestions on the overall submission or specific details and on occasion specific requests for changes as part of the review. If permitted, also additional changes may be pushed into your pull request branch or a pull request may be filed in your LAMMPS fork on GitHub to include those changes.
The LAMMPS developer may then decide to assign the pull request to another developer (e.g. when that developer is more knowledgeable about the submitted feature or enhancement or has written the modified code). It may also happen, that additional developers are requested to provide a review and approve the changes. For submissions, that may change the general behavior of LAMMPS, or where a possibility of unwanted side effects exists, additional tests may be requested by the assigned developer. The LAMMPS developer may then decide to assign the pull request to another developer (e.g. when that developer is more knowledgeable about the submitted feature or enhancement or has written the modified code). It may also happen, that additional developers are requested to provide a review and approve the changes. For submissions, that may change the general behavior of LAMMPS, or where a possibility of unwanted side effects exists, additional tests may be requested by the assigned developer.
If the assigned developer is satisfied and considers the submission ready for inclusion into LAMMPS, the pull request will receive approvals and be merged into the master branch by one of the core LAMMPS developers. After the pull request is merged, you may delete the feature branch used for the pull request in your personal LAMMPS fork. If the assigned developer is satisfied and considers the submission ready for inclusion into LAMMPS, the pull request will receive approvals and be merged into the master branch by one of the core LAMMPS developers. After the pull request is merged, you may delete the feature branch used for the pull request in your personal LAMMPS fork.

32
cmake/Modules/YAML.cmake Normal file
View File

@ -0,0 +1,32 @@
message(STATUS "Downloading and building YAML library")
include(ExternalProject)
set(YAML_URL "https://pyyaml.org/download/libyaml/yaml-0.2.5.tar.gz" CACHE STRING "URL for libyaml tarball")
mark_as_advanced(YAML_URL)
ExternalProject_Add(libyaml
URL ${YAML_URL}
URL_MD5 bb15429d8fb787e7d3f1c83ae129a999
SOURCE_DIR "${CMAKE_BINARY_DIR}/yaml-src"
BINARY_DIR "${CMAKE_BINARY_DIR}/yaml-build"
CONFIGURE_COMMAND <SOURCE_DIR>/configure ${CONFIGURE_REQUEST_PIC}
CXX=${CMAKE_CXX_COMPILER}
CC=${CMAKE_C_COMPILER}
--prefix=<INSTALL_DIR> --disable-shared
BUILD_BYPRODUCTS <INSTALL_DIR>/lib/${CMAKE_FIND_LIBRARY_PREFIXES}yaml.a
TEST_COMMAND "")
ExternalProject_Get_Property(libyaml INSTALL_DIR)
set(YAML_INCLUDE_DIR ${INSTALL_DIR}/include)
set(YAML_LIBRARY_DIR ${INSTALL_DIR}/lib)
# workaround for CMake 3.10 on ubuntu 18.04
file(MAKE_DIRECTORY ${YAML_INCLUDE_DIR})
file(MAKE_DIRECTORY ${YAML_LIBRARY_DIR})
set(YAML_LIBRARY_PATH ${INSTALL_DIR}/lib/${CMAKE_FIND_LIBRARY_PREFIXES}yaml.a)
add_library(Yaml::Yaml UNKNOWN IMPORTED)
set_target_properties(Yaml::Yaml PROPERTIES
IMPORTED_LOCATION ${YAML_LIBRARY_PATH}
INTERFACE_INCLUDE_DIRECTORIES ${YAML_INCLUDE_DIR})
add_dependencies(Yaml::Yaml libyaml)

View File

@ -95,7 +95,7 @@ on the pull request discussion page on GitHub, so that other developers
can later review the entire discussion after the fact and understand the can later review the entire discussion after the fact and understand the
rationale behind choices made. Exceptions to this policy are technical rationale behind choices made. Exceptions to this policy are technical
discussions, that are centered on tools or policies themselves discussions, that are centered on tools or policies themselves
(git, github, c++) rather than on the content of the pull request. (git, GitHub, c++) rather than on the content of the pull request.
### Checklist for Pull Requests ### Checklist for Pull Requests

View File

@ -111,8 +111,10 @@ error margin). The status of this automated testing can be viewed on
The unit testing facility is integrated into the CMake build process The unit testing facility is integrated into the CMake build process
of the LAMMPS source code distribution itself. It can be enabled by of the LAMMPS source code distribution itself. It can be enabled by
setting ``-D ENABLE_TESTING=on`` during the CMake configuration step. setting ``-D ENABLE_TESTING=on`` during the CMake configuration step.
It requires the `PyYAML <http://pyyaml.org/>`_ library and development It requires the `YAML <http://pyyaml.org/>`_ library and development
headers to compile and will download and compile a recent version of the headers (if not found locally a recent version will be downloaded
and compiled transparently) to compile and will download and compile
a specific recent version of the
`Googletest <https://github.com/google/googletest/>`_ C++ test framework `Googletest <https://github.com/google/googletest/>`_ C++ test framework
for implementing the tests. for implementing the tests.

View File

@ -241,6 +241,7 @@ OPT.
* :doc:`spin/dipole/long <pair_spin_dipole>` * :doc:`spin/dipole/long <pair_spin_dipole>`
* :doc:`spin/dmi <pair_spin_dmi>` * :doc:`spin/dmi <pair_spin_dmi>`
* :doc:`spin/exchange <pair_spin_exchange>` * :doc:`spin/exchange <pair_spin_exchange>`
* :doc:`spin/exchange/biquadratic <pair_spin_exchange>`
* :doc:`spin/magelec <pair_spin_magelec>` * :doc:`spin/magelec <pair_spin_magelec>`
* :doc:`spin/neel <pair_spin_neel>` * :doc:`spin/neel <pair_spin_neel>`
* :doc:`srp <pair_srp>` * :doc:`srp <pair_srp>`

View File

@ -119,7 +119,6 @@ Doc page with :doc:`ERROR messages <Errors_messages>`
:doc:`pair style zero <pair_zero>` with a suitable cutoff or use :doc:`comm_modify cutoff <comm_modify>`. :doc:`pair style zero <pair_zero>` with a suitable cutoff or use :doc:`comm_modify cutoff <comm_modify>`.
*Communication cutoff is shorter than a bond length based estimate. This may lead to errors.* *Communication cutoff is shorter than a bond length based estimate. This may lead to errors.*
Since LAMMPS stores topology data with individual atoms, all atoms Since LAMMPS stores topology data with individual atoms, all atoms
comprising a bond, angle, dihedral or improper must be present on any comprising a bond, angle, dihedral or improper must be present on any
sub-domain that "owns" the atom with the information, either as a sub-domain that "owns" the atom with the information, either as a

View File

@ -72,7 +72,7 @@ explained in more detail here: `feature branch workflow <https://www.atlassian.c
**Feature branches** **Feature branches**
First of all, create a clone of your version on github on your local First of all, create a clone of your version on GitHub on your local
machine via HTTPS: machine via HTTPS:
.. code-block:: bash .. code-block:: bash
@ -155,7 +155,7 @@ useful message that explains the change.
.. code-block:: bash .. code-block:: bash
$ git commit -m 'Finally updated the github tutorial' $ git commit -m 'Finally updated the GitHub tutorial'
After the commit, the changes can be pushed to the same branch on GitHub: After the commit, the changes can be pushed to the same branch on GitHub:

View File

@ -67,5 +67,5 @@ rotate.
The only frictional idealized walls currently in LAMMPS are flat or The only frictional idealized walls currently in LAMMPS are flat or
curved surfaces specified by the :doc:`fix wall/gran <fix_wall_gran>` curved surfaces specified by the :doc:`fix wall/gran <fix_wall_gran>`
command. At some point we plan to allow regoin surfaces to be used as command. At some point we plan to allow region surfaces to be used as
frictional walls, as well as triangulated surfaces. frictional walls, as well as triangulated surfaces.

View File

@ -24,13 +24,15 @@ DOI for the LAMMPS code
LAMMPS developers use the `Zenodo service at CERN LAMMPS developers use the `Zenodo service at CERN
<https://zenodo.org/>`_ to create digital object identifies (DOI) for <https://zenodo.org/>`_ to create digital object identifies (DOI) for
stable releases of the LAMMPS code. There are two types of DOIs for the stable releases of the LAMMPS code. There are two types of DOIs for the
LAMMPS source code: 1) the canonical DOI for **all** versions of LAMMPS, LAMMPS source code: the canonical DOI for **all** versions of LAMMPS,
which will always point to the latest stable release version is: which will always point to the **latest** stable release version is:
`DOI: 10.5281/zenodo.3726416 <https://dx.doi/org/10.5281/zenodo.3726416>`_ - DOI: `10.5281/zenodo.3726416 <https://dx.doi.org/10.5281/zenodo.3726416>`_
In addition there are DOIs for individual stable releases starting with In addition there are DOIs for individual stable releases. Currently there are:
the `3 March 2020 version, DOI:10.5281/zenodo.3726417 <https://dx.doi/org/10.5281/zenodo.3726416>`_
- 3 March 2020 version: `DOI:10.5281/zenodo.3726417 <https://dx.doi.org/10.5281/zenodo.3726417>`_
- 29 October 2020 version: `DOI:10.5281/zenodo.4157471 <https://dx.doi.org/10.5281/zenodo.4157471>`_
Home page Home page

View File

@ -1036,9 +1036,11 @@ the usual manner via MD. Various pair, fix, and compute styles.
* :doc:`pair_style spin/dipole/long <pair_spin_dipole>` * :doc:`pair_style spin/dipole/long <pair_spin_dipole>`
* :doc:`pair_style spin/dmi <pair_spin_dmi>` * :doc:`pair_style spin/dmi <pair_spin_dmi>`
* :doc:`pair_style spin/exchange <pair_spin_exchange>` * :doc:`pair_style spin/exchange <pair_spin_exchange>`
* :doc:`pair_style spin/exchange/biquadratic <pair_spin_exchange>`
* :doc:`pair_style spin/magelec <pair_spin_magelec>` * :doc:`pair_style spin/magelec <pair_spin_magelec>`
* :doc:`pair_style spin/neel <pair_spin_neel>` * :doc:`pair_style spin/neel <pair_spin_neel>`
* :doc:`fix nve/spin <fix_nve_spin>` * :doc:`fix nve/spin <fix_nve_spin>`
* :doc:`fix langevin/spin <fix_langevin_spin>`
* :doc:`fix precession/spin <fix_precession_spin>` * :doc:`fix precession/spin <fix_precession_spin>`
* :doc:`compute spin <compute_spin>` * :doc:`compute spin <compute_spin>`
* :doc:`neb/spin <neb_spin>` * :doc:`neb/spin <neb_spin>`

View File

@ -14,7 +14,7 @@ Syntax
* AtC fixID = ID of :doc:`fix atc <fix_atc>` instance * AtC fixID = ID of :doc:`fix atc <fix_atc>` instance
* *output* or *output index* = name of the AtC sub-command * *output* or *output index* = name of the AtC sub-command
* filename_prefix = prefix for data files (for *output*) * filename_prefix = prefix for data files (for *output*)
* frequency = frequency of output in time-steps (for *output*) * frequency = frequency of output in timesteps (for *output*)
* optional keywords for *output*: * optional keywords for *output*:
- text = creates text output of index, step and nodal variable values for unique nodes - text = creates text output of index, step and nodal variable values for unique nodes

View File

@ -56,7 +56,7 @@ is slightly modified only for the computation of long-range forces. A
good cluster decomposition constitutes in building clusters which good cluster decomposition constitutes in building clusters which
contain the fastest covalent bonds inside clusters. contain the fastest covalent bonds inside clusters.
If the clusters are chosen suitably, the :doc:`run_style respa <run_style>` is stable for outer time-steps of at least 8fs. If the clusters are chosen suitably, the :doc:`run_style respa <run_style>` is stable for outer timesteps of at least 8fs.
---------- ----------

View File

@ -62,7 +62,7 @@ with:
The field value in Tesla is multiplied by the gyromagnetic The field value in Tesla is multiplied by the gyromagnetic
ratio, :math:`g \cdot \mu_B/\hbar`, converting it into a precession frequency in ratio, :math:`g \cdot \mu_B/\hbar`, converting it into a precession frequency in
rad.THz (in metal units and with :math:`\mu_B = 5.788 eV/T`). rad.THz (in metal units and with :math:`\mu_B = 5.788\cdot 10^{-5}` eV/T).
As a comparison, the figure below displays the simulation of a As a comparison, the figure below displays the simulation of a
single spin (of norm :math:`\mu_i = 1.0`) submitted to an external single spin (of norm :math:`\mu_i = 1.0`) submitted to an external

View File

@ -90,10 +90,10 @@ accepted, *h* is increased by a proportional amount, and the next ODE step is be
Otherwise, *h* is shrunk and the ODE step is repeated. Otherwise, *h* is shrunk and the ODE step is repeated.
Run-time diagnostics are available for the rkf45 ODE solver. The frequency Run-time diagnostics are available for the rkf45 ODE solver. The frequency
(in time-steps) that diagnostics are reported is controlled by the last (optional) (in timesteps) that diagnostics are reported is controlled by the last (optional)
12th argument. A negative frequency means that diagnostics are reported once at the 12th argument. A negative frequency means that diagnostics are reported once at the
end of each run. A positive value N means that the diagnostics are reported once end of each run. A positive value N means that the diagnostics are reported once
per N time-steps. per N timesteps.
The diagnostics report the average # of integrator steps and RHS function evaluations The diagnostics report the average # of integrator steps and RHS function evaluations
and run-time per ODE as well as the average/RMS/min/max per process. If the and run-time per ODE as well as the average/RMS/min/max per process. If the

View File

@ -143,7 +143,7 @@ combinations, else an error will result.
Mixing, shift, table, tail correction, restart, rRESPA info Mixing, shift, table, tail correction, restart, rRESPA info
""""""""""""""""""""""""""""""""""""""""""""""""""""""""""" """""""""""""""""""""""""""""""""""""""""""""""""""""""""""
This pair styles do not support the :doc:`pair_modify <pair_modify>` This pair style do not support the :doc:`pair_modify <pair_modify>`
mix, shift, table, and tail options. mix, shift, table, and tail options.
This pair style writes its information to :doc:`binary restart files This pair style writes its information to :doc:`binary restart files

View File

@ -117,7 +117,7 @@ global Coulombic cutoff is allowed.
Mixing, shift, table, tail correction, restart, rRESPA info Mixing, shift, table, tail correction, restart, rRESPA info
""""""""""""""""""""""""""""""""""""""""""""""""""""""""""" """""""""""""""""""""""""""""""""""""""""""""""""""""""""""
This pair styles does not support mixing. Thus, coefficients for all This pair style does not support mixing. Thus, coefficients for all
I,J pairs must be specified explicitly. I,J pairs must be specified explicitly.
This pair style supports the :doc:`pair_modify <pair_modify>` shift This pair style supports the :doc:`pair_modify <pair_modify>` shift

View File

@ -160,7 +160,7 @@ For atom type pairs I,J and I != J, the epsilon and sigma coefficients
and cutoff distance for this pair style can be mixed. The default mix and cutoff distance for this pair style can be mixed. The default mix
value is *geometric*\ . See the "pair_modify" command for details. value is *geometric*\ . See the "pair_modify" command for details.
This pair styles supports the :doc:`pair_modify <pair_modify>` shift This pair style supports the :doc:`pair_modify <pair_modify>` shift
option for the energy of the Lennard-Jones portion of the pair option for the energy of the Lennard-Jones portion of the pair
interaction, but only for sphere-sphere interactions. There is no interaction, but only for sphere-sphere interactions. There is no
shifting performed for ellipsoidal interactions due to the anisotropic shifting performed for ellipsoidal interactions due to the anisotropic

View File

@ -75,14 +75,15 @@ This pair style can only be used via the *pair* keyword of the
Restrictions Restrictions
"""""""""""" """"""""""""
This pair styles is part of the MANYBODY package. It is only enabled This pair style is part of the MANYBODY package. It is only enabled
if LAMMPS was built with that package. See the :doc:`Build package <Build_package>` doc page for more info. if LAMMPS was built with that package.
See the :doc:`Build package <Build_package>` doc page for more info.
This pair potential requires the :doc:`newton <newton>` setting to be This pair potential requires the :doc:`newton <newton>` setting to be
"on" for pair interactions. "on" for pair interactions.
The C.lcbop potential file provided with LAMMPS (see the potentials The ``C.lcbop`` potential file provided with LAMMPS (see the potentials
directory) is parameterized for metal :doc:`units <units>`. You can use directory) is parameterized for :doc:`metal units <units>`. You can use
the LCBOP potential with any LAMMPS units, but you would need to the LCBOP potential with any LAMMPS units, but you would need to
create your own LCBOP potential file with coefficients listed in the create your own LCBOP potential file with coefficients listed in the
appropriate units if your simulation does not use "metal" units. appropriate units if your simulation does not use "metal" units.

View File

@ -298,7 +298,7 @@ described above. For each of the F functions, nx values are listed.
Mixing, shift, table, tail correction, restart, rRESPA info Mixing, shift, table, tail correction, restart, rRESPA info
""""""""""""""""""""""""""""""""""""""""""""""""""""""""""" """""""""""""""""""""""""""""""""""""""""""""""""""""""""""
This pair styles does not support the :doc:`pair_modify <pair_modify>` This pair style does not support the :doc:`pair_modify <pair_modify>`
shift, table, and tail options. shift, table, and tail options.
This pair style does not write their information to :doc:`binary restart This pair style does not write their information to :doc:`binary restart

View File

@ -173,7 +173,7 @@ equation for the Hamaker constant presented here. Mixing of sigma and
epsilon followed by calculation of the energy prefactors using the epsilon followed by calculation of the energy prefactors using the
equations above is recommended. equations above is recommended.
This pair styles supports the :doc:`pair_modify <pair_modify>` shift This pair style supports the :doc:`pair_modify <pair_modify>` shift
option for the energy of the Lennard-Jones portion of the pair option for the energy of the Lennard-Jones portion of the pair
interaction, but only for sphere-sphere interactions. There is no interaction, but only for sphere-sphere interactions. There is no
shifting performed for ellipsoidal interactions due to the anisotropic shifting performed for ellipsoidal interactions due to the anisotropic

View File

@ -1,14 +1,19 @@
.. index:: pair_style spin/exchange .. index:: pair_style spin/exchange
.. index:: pair_style spin/exchange/biquadratic
pair_style spin/exchange command pair_style spin/exchange command
================================ ================================
pair_style spin/exchange/biquadratic command
============================================
Syntax Syntax
"""""" """"""
.. code-block:: LAMMPS .. code-block:: LAMMPS
pair_style spin/exchange cutoff pair_style spin/exchange cutoff
pair_style spin/exchange/biquadratic cutoff
* cutoff = global cutoff pair (distance in metal units) * cutoff = global cutoff pair (distance in metal units)
@ -19,7 +24,11 @@ Examples
pair_style spin/exchange 4.0 pair_style spin/exchange 4.0
pair_coeff * * exchange 4.0 0.0446928 0.003496 1.4885 pair_coeff * * exchange 4.0 0.0446928 0.003496 1.4885
pair_coeff 1 2 exchange 6.0 -0.01575 0.0 1.965 pair_coeff 1 2 exchange 6.0 -0.01575 0.0 1.965 offset yes
pair_style spin/exchange/biquadratic 4.0
pair_coeff * * biquadratic 4.0 0.05 0.03 1.48 0.05 0.03 1.48 offset no
pair_coeff 1 2 biquadratic 6.0 -0.01 0.0 1.9 0.0 0.1 19
Description Description
""""""""""" """""""""""
@ -31,69 +40,163 @@ pairs of magnetic spins:
H_{ex} = -\sum_{i,j}^N J_{ij} (r_{ij}) \,\vec{s}_i \cdot \vec{s}_j H_{ex} = -\sum_{i,j}^N J_{ij} (r_{ij}) \,\vec{s}_i \cdot \vec{s}_j
where :math:`\vec{s}_i` and :math:`\vec{s}_j` are two neighboring magnetic spins of two particles, where :math:`\vec{s}_i` and :math:`\vec{s}_j` are two unit vectors representing
:math:`r_{ij} = \vert \vec{r}_i - \vec{r}_j \vert` is the inter-atomic distance between the two the magnetic spins of two particles (usually atoms), and
particles. The summation is over pairs of nearest neighbors. :math:`r_{ij} = \vert \vec{r}_i - \vec{r}_j \vert` is the inter-atomic distance
:math:`J(r_{ij})` is a function defining the intensity and the sign of the exchange between those two particles. The summation is over pairs of nearest neighbors.
interaction for different neighboring shells. This function is defined as: :math:`J(r_{ij})` is a function defining the intensity and the sign of the
exchange interaction for different neighboring shells.
Style *spin/exchange/biquadratic* computes a biquadratic exchange interaction
between pairs of magnetic spins:
.. math:: .. math::
{J}\left( r_{ij} \right) = 4 a \left( \frac{r_{ij}}{d} \right)^2 \left( 1 - b \left( \frac{r_{ij}}{d} \right)^2 \right) e^{-\left( \frac{r_{ij}}{d} \right)^2 }\Theta (R_c - r_{ij}) H_{bi} = -\sum_{i, j}^{N} {J}_{ij} \left(r_{ij} \right)\,
\vec{s}_{i}\cdot \vec{s}_{j}
-\sum_{i, j}^{N} {K}_{ij} \left(r_{ij} \right)\,
\left(\vec{s}_{i}\cdot
\vec{s}_{j}\right)^2
where :math:`a`, :math:`b` and :math:`d` are the three constant coefficients defined in the associated where :math:`\vec{s}_i`, :math:`\vec{s}_j`, :math:`r_{ij}` and
"pair_coeff" command, and :math:`R_c` is the radius cutoff associated to :math:`J(r_{ij})` have the same definitions as above, and :math:`K(r_{ij})` is
the pair interaction (see below for more explanations). a second function, defining the intensity and the sign of the biquadratic term.
The coefficients :math:`a`, :math:`b`, and :math:`d` need to be fitted so that the function above matches with The interatomic dependence of :math:`J(r_{ij})` and :math:`K(r_{ij})` in both
the value of the exchange interaction for the :math:`N` neighbor shells taken into account. interactions above is defined by the following function:
Examples and more explanations about this function and its parameterization are reported
in :ref:`(Tranchida) <Tranchida3>`. .. math::
{f}\left( r_{ij} \right) = 4 a \left( \frac{r_{ij}}{d} \right)^2
\left( 1 - b \left( \frac{r_{ij}}{d} \right)^2 \right)
e^{-\left( \frac{r_{ij}}{d} \right)^2 }\Theta (R_c - r_{ij})
where :math:`a`, :math:`b` and :math:`d` are the three constant coefficients
defined in the associated "pair_coeff" command, and :math:`R_c` is the radius
cutoff associated to the pair interaction (see below for more explanations).
The coefficients :math:`a`, :math:`b`, and :math:`d` need to be fitted so that
the function above matches with the value of the exchange interaction for the
:math:`N` neighbor shells taken into account.
Examples and more explanations about this function and its parameterization
are reported in :ref:`(Tranchida) <Tranchida3>`.
When a *spin/exchange/biquadratic* pair style is defined, six coefficients
(three for :math:`J(r_{ij})`, and three for :math:`K(r_{ij})`) have to be
fitted.
From this exchange interaction, each spin :math:`i` will be submitted From this exchange interaction, each spin :math:`i` will be submitted
to a magnetic torque :math:`\vec{\omega}`, and its associated atom can be submitted to a to a magnetic torque :math:`\vec{\omega}_{i}`, and its associated atom can be
force :math:`\vec{F}` for spin-lattice calculations (see :doc:`fix nve/spin <fix_nve_spin>`), submitted to a force :math:`\vec{F}_{i}` for spin-lattice calculations (see
such as: :doc:`fix nve/spin <fix_nve_spin>`), such as:
.. math:: .. math::
\vec{\omega}_{i} = \frac{1}{\hbar} \sum_{j}^{Neighb} {J} \vec{\omega}_{i} = \frac{1}{\hbar} \sum_{j}^{Neighb} {J}
\left(r_{ij} \right)\,\vec{s}_{j} \left(r_{ij} \right)\,\vec{s}_{j}
~~{\rm and}~~ ~~{\rm and}~~
\vec{F}_{i} = \sum_{j}^{Neighb} \frac{\partial {J} \left(r_{ij} \right)}{ \partial r_{ij}} \left( \vec{s}_{i}\cdot \vec{s}_{j} \right) \vec{e}_{ij} \vec{F}_{i} = \sum_{j}^{Neighb} \frac{\partial {J} \left(r_{ij} \right)}{
\partial r_{ij}} \left( \vec{s}_{i}\cdot \vec{s}_{j} \right) \vec{e}_{ij}
with :math:`\hbar` the Planck constant (in metal units), and :math:`\vec{e}_{ij} = \frac{\vec{r}_i - \vec{r}_j}{\vert \vec{r}_i-\vec{r}_j \vert}` the unit with :math:`\hbar` the Planck constant (in metal units), and :math:`\vec{e}_{ij}
= \frac{\vec{r}_i - \vec{r}_j}{\vert \vec{r}_i-\vec{r}_j \vert}` the unit
vector between sites :math:`i` and :math:`j`. vector between sites :math:`i` and :math:`j`.
Equivalent forces and magnetic torques are generated for the biquadratic term
when a *spin/exchange/biquadratic* pair style is defined.
More details about the derivation of these torques/forces are reported in More details about the derivation of these torques/forces are reported in
:ref:`(Tranchida) <Tranchida3>`. :ref:`(Tranchida) <Tranchida3>`.
For the *spin/exchange* pair style, the following coefficients must be defined For the *spin/exchange* and *spin/exchange/biquadratic* pair styles, the
for each pair of atoms types via the :doc:`pair_coeff <pair_coeff>` command as in following coefficients must be defined for each pair of atoms types via the
the examples above, or in the data file or restart files read by the :doc:`pair_coeff <pair_coeff>` command as in the examples above, or in the data
:doc:`read_data <read_data>` or :doc:`read_restart <read_restart>` commands, and file or restart files read by the :doc:`read_data <read_data>` or
set in the following order: :doc:`read_restart <read_restart>` commands, and set in the following order:
* :math:`R_c` (distance units) * :math:`R_c` (distance units)
* :math:`a` (energy units) * :math:`a` (energy units)
* :math:`b` (adim parameter) * :math:`b` (adim parameter)
* :math:`d` (distance units) * :math:`d` (distance units)
Note that :math:`R_c` is the radius cutoff of the considered exchange interaction, for the *spin/exchange* pair style, and:
and :math:`a`, :math:`b` and :math:`d` are the three coefficients performing the parameterization
of the function :math:`J(r_{ij})` defined above. * :math:`R_c` (distance units)
* :math:`a_j` (energy units)
* :math:`b_j` (adim parameter)
* :math:`d_j` (distance units)
* :math:`a_k` (energy units)
* :math:`b_k` (adim parameter)
* :math:`d_k` (distance units)
for the *spin/exchange/biquadratic* pair style.
Note that :math:`R_c` is the radius cutoff of the considered exchange
interaction, and :math:`a`, :math:`b` and :math:`d` are the three coefficients
performing the parameterization of the function :math:`J(r_{ij})` defined
above (in the *biquadratic* style, :math:`a_j`, :math:`b_j`, :math:`d_j` and
:math:`a_k`, :math:`b_k`, :math:`d_k` are the coefficients of :math:`J(r_{ij})`
and :math:`K(r_{ij})` respectively).
None of those coefficients is optional. If not specified, the None of those coefficients is optional. If not specified, the
*spin/exchange* pair style cannot be used. *spin/exchange* pair style cannot be used.
---------- ----------
**Offsetting magnetic forces and energies**\ :
For spin-lattice simulation, it can be useful to offset the
mechanical forces and energies generated by the exchange
interaction.
The *offset* keyword allows to apply this offset.
By setting *offset* to *yes*, the energy definitions above are
replaced by:
.. math::
H_{ex} = -\sum_{i,j}^N J_{ij} (r_{ij}) \,[ \vec{s}_i \cdot \vec{s}_j-1 ]
for the *spin/exchange* pair style, and:
.. math::
H_{bi} = -\sum_{i, j}^{N} {J}_{ij} \left(r_{ij} \right)\,
[ \vec{s}_{i}\cdot \vec{s}_{j} -1 ]
-\sum_{i, j}^{N} {K}_{ij} \left(r_{ij} \right)\,
[ \left(\vec{s}_{i}\cdot
\vec{s}_{j}\right)^2 -1]
for the *spin/exchange/biquadratic* pair style.
Note that this offset only affects the calculation of the energy
and mechanical forces. It does not modify the calculation of the
precession vectors (and thus does no impact the purely magnetic
properties).
This ensures that when all spins are aligned, the magnetic energy
and the associated mechanical forces (and thus the pressure
generated by the magnetic potential) are null.
.. note::
This offset term can be very important when calculations such as
equations of state (energy vs volume, or energy vs pressure) are
being performed. Indeed, setting the *offset* term ensures that
at the ground state of the crystal and at the equilibrium magnetic
configuration (typically ferromagnetic), the pressure is null,
as expected.
Otherwise, magnetic forces could generate a residual pressure.
When the *offset* option is set to *no*, no offset is applied
(also corresponding to the default option).
----------
Restrictions Restrictions
"""""""""""" """"""""""""
All the *pair/spin* styles are part of the SPIN package. These styles All the *pair/spin* styles are part of the SPIN package. These styles
are only enabled if LAMMPS was built with this package, and if the are only enabled if LAMMPS was built with this package, and if the
atom_style "spin" was declared. See the :doc:`Build package <Build_package>` doc page for more info. atom_style "spin" was declared.
See the :doc:`Build package <Build_package>` doc page for more info.
Related commands Related commands
"""""""""""""""" """"""""""""""""
@ -105,7 +208,7 @@ Default
""""""" """""""
none The default *offset* keyword value is *no*.
---------- ----------

View File

@ -124,7 +124,7 @@ at the cutoff distance :math:`r_c`.
Mixing, shift, table, tail correction, restart, rRESPA info Mixing, shift, table, tail correction, restart, rRESPA info
""""""""""""""""""""""""""""""""""""""""""""""""""""""""""" """""""""""""""""""""""""""""""""""""""""""""""""""""""""""
This pair styles does not support mixing. This pair style does not support mixing.
This pair style does not support the :doc:`pair_modify <pair_modify>` This pair style does not support the :doc:`pair_modify <pair_modify>`
shift option for the energy of the pair interaction. Note that as shift option for the energy of the pair interaction. Note that as

View File

@ -305,6 +305,7 @@ accelerated styles exist.
* :doc:`spin/dipole/long <pair_spin_dipole>` - * :doc:`spin/dipole/long <pair_spin_dipole>` -
* :doc:`spin/dmi <pair_spin_dmi>` - * :doc:`spin/dmi <pair_spin_dmi>` -
* :doc:`spin/exchange <pair_spin_exchange>` - * :doc:`spin/exchange <pair_spin_exchange>` -
* :doc:`spin/exchange/biquadratic <pair_spin_exchange>` -
* :doc:`spin/magelec <pair_spin_magelec>` - * :doc:`spin/magelec <pair_spin_magelec>` -
* :doc:`spin/neel <pair_spin_neel>` - * :doc:`spin/neel <pair_spin_neel>` -
* :doc:`srp <pair_srp>` - * :doc:`srp <pair_srp>` -

View File

@ -370,6 +370,8 @@ needed to generate absolute, unscaled coordinates.
Restrictions Restrictions
"""""""""""" """"""""""""
The *native* dump file reader does not support binary .bin dump files.
To read gzipped dump files, you must compile LAMMPS with the To read gzipped dump files, you must compile LAMMPS with the
-DLAMMPS_GZIP option. See the :doc:`Build settings <Build_settings>` -DLAMMPS_GZIP option. See the :doc:`Build settings <Build_settings>`
doc page for details. doc page for details.

View File

@ -99,14 +99,15 @@ files do not match the specified output frequency.
---------- ----------
If more than one dump file is specified, the dump files are read one If more than one dump file is specified, the dump files are read one
after the other. It is assumed that snapshot timesteps will be in after the other in the order specified. It is assumed that snapshot
ascending order. If a snapshot is encountered that is not in timesteps will be in ascending order. If a snapshot is encountered that
ascending order, it will skip the snapshot until it reads one that is. is not in ascending order, it will skip the snapshot until it reads one
that is.
This allows skipping of a duplicate snapshot (same timestep), This allows skipping of a duplicate snapshot (same timestep),
e.g. that appeared at the end of one file and beginning of the next. e.g. that appeared at the end of one file and beginning of the next.
However if you specify a series of dump files in an incorrect order However if you specify a series of dump files in an incorrect order
(with respect to the timesteps they contain), you may skip large (with respect to the timesteps they contain), you may skip large
numbers of snapshots numbers of snapshots.
Note that the dump files specified as part of the *dump* keyword can be Note that the dump files specified as part of the *dump* keyword can be
parallel files, i.e. written as multiple files either per processor parallel files, i.e. written as multiple files either per processor
@ -118,17 +119,24 @@ and write parallel dump files.
The *first*\ , *last*\ , *every*\ , *skip* keywords determine which The *first*\ , *last*\ , *every*\ , *skip* keywords determine which
snapshots are read from the dump file(s). Snapshots are skipped until snapshots are read from the dump file(s). Snapshots are skipped until
they have a timestamp >= *Nfirst*\ . When a snapshot with a timestamp > they have a timestep >= *Nfirst*\ . When a snapshot with a timestep >
*Nlast* is encountered, the rerun command finishes. Note below that *Nlast* is encountered, the rerun command finishes. Note that
the defaults for *first* and *last* are to read all snapshots. If the the defaults for *first* and *last* are to read all snapshots. If the
*every* keyword is set to a value > 0, then only snapshots with *every* keyword is set to a value > 0, then only snapshots with
timestamps that are a multiple of *Nevery* are read (the first timesteps that are a multiple of *Nevery* are read (the first
snapshot is always read). If *Nevery* = 0, then this criterion is snapshot is always read). If *Nevery* = 0, then this criterion is
ignored, i.e. every snapshot is read that meets the other criteria. ignored, i.e. every snapshot is read that meets the other criteria.
If the *skip* keyword is used, then after the first snapshot is read, If the *skip* keyword is used, then after the first snapshot is read,
every Nth snapshot is read, where N = *Nskip*\ . E.g. if *Nskip* = 3, every Nth snapshot is read, where N = *Nskip*\ . E.g. if *Nskip* = 3,
then only 1 out of every 3 snapshots is read, assuming the snapshot then only 1 out of every 3 snapshots is read, assuming the snapshot
timestamp is also consistent with the other criteria. timestep is also consistent with the other criteria.
.. note::
Not all dump formats contain the timestep and not all dump readers
support reading it. In that case individual snapshots are assigned
consecutive timestep numbers starting at 1.
The *start* and *stop* keywords do not affect which snapshots are read The *start* and *stop* keywords do not affect which snapshots are read
from the dump file(s). Rather, they have the same meaning that they from the dump file(s). Rather, they have the same meaning that they
@ -205,9 +213,8 @@ thermodynamic output or new dump file output.
Restrictions Restrictions
"""""""""""" """"""""""""
To read gzipped dump files, you must compile LAMMPS with the The *rerun* command is subject to all restrictions of
-DLAMMPS_GZIP option. See the :doc:`Build settings <Build_settings>` the :doc:`read_dump <read_dump>` command.
doc page for details.
Related commands Related commands
"""""""""""""""" """"""""""""""""

View File

@ -240,6 +240,7 @@ bigint
Bij Bij
bilayer bilayer
bilayers bilayers
biquadratic
binsize binsize
binstyle binstyle
binutils binutils
@ -2614,7 +2615,6 @@ Ree
refactored refactored
refactoring refactoring
reflectionstyle reflectionstyle
regoin
Reinders Reinders
reinit reinit
relaxbox relaxbox

View File

@ -26,7 +26,7 @@ velocity all create 100 4928459 rot yes dist gaussian
#pair_style hybrid/overlay eam/alloy spin/exchange 4.0 spin/neel 4.0 #pair_style hybrid/overlay eam/alloy spin/exchange 4.0 spin/neel 4.0
pair_style hybrid/overlay eam/alloy spin/exchange 4.0 pair_style hybrid/overlay eam/alloy spin/exchange 4.0
pair_coeff * * eam/alloy Co_PurjaPun_2012.eam.alloy Co pair_coeff * * eam/alloy Co_PurjaPun_2012.eam.alloy Co
pair_coeff * * spin/exchange exchange 4.0 -0.3593 1.135028015e-05 1.064568567 pair_coeff * * spin/exchange exchange 4.0 -0.3593 1.135028015e-05 1.0645 offset yes
#pair_coeff * * spin/neel neel 4.0 0.0048 0.234 1.168 2.6905 0.705 0.652 #pair_coeff * * spin/neel neel 4.0 0.0048 0.234 1.168 2.6905 0.705 0.652
neighbor 0.1 bin neighbor 0.1 bin

View File

@ -25,7 +25,7 @@ velocity all create 100 4928459 rot yes dist gaussian
pair_style hybrid/overlay eam/alloy spin/exchange 3.5 pair_style hybrid/overlay eam/alloy spin/exchange 3.5
pair_coeff * * eam/alloy Fe_Mishin2006.eam.alloy Fe pair_coeff * * eam/alloy Fe_Mishin2006.eam.alloy Fe
pair_coeff * * spin/exchange exchange 3.4 0.02726 0.2171 1.841 pair_coeff * * spin/exchange exchange 3.4 0.02726 0.2171 1.841 offset yes
neighbor 0.1 bin neighbor 0.1 bin
neigh_modify every 10 check yes delay 20 neigh_modify every 10 check yes delay 20

View File

@ -21,9 +21,9 @@ mass 1 55.845
set group all spin 2.2 -1.0 0.0 0.0 set group all spin 2.2 -1.0 0.0 0.0
velocity all create 100 4928459 rot yes dist gaussian velocity all create 100 4928459 rot yes dist gaussian
pair_style hybrid/overlay eam/alloy spin/exchange 3.5 pair_style hybrid/overlay eam/alloy spin/exchange/biquadratic 3.5
pair_coeff * * eam/alloy Fe_Mishin2006.eam.alloy Fe pair_coeff * * eam/alloy Fe_Mishin2006.eam.alloy Fe
pair_coeff * * spin/exchange exchange 3.4 0.02726 0.2171 1.841 pair_coeff * * spin/exchange/biquadratic biquadratic 3.4 0.02726 0.2171 1.841 0.0 0.0 2.0 offset yes
neighbor 0.1 bin neighbor 0.1 bin
neigh_modify every 10 check yes delay 20 neigh_modify every 10 check yes delay 20

View File

@ -6,9 +6,17 @@ import matplotlib.pyplot as plt
import mpmath as mp import mpmath as mp
hbar=0.658212 # Planck's constant (eV.fs/rad) hbar=0.658212 # Planck's constant (eV.fs/rad)
J0=0.05 # per-neighbor exchange interaction (eV) # J0=0.05 # per-neighbor exchange interaction (eV)
# exchange interaction parameters
J1 = 11.254 # in eV
J2 = 0.0 # adim
J3 = 1.0 # in Ang.
# initial spins
S1 = np.array([1.0, 0.0, 0.0]) S1 = np.array([1.0, 0.0, 0.0])
S2 = np.array([0.0, 1.0, 0.0]) S2 = np.array([0.0, 1.0, 0.0])
alpha=0.01 # damping coefficient alpha=0.01 # damping coefficient
pi=math.pi pi=math.pi
@ -30,6 +38,14 @@ def rotation_matrix(axis, theta):
[2 * (bc - ad), aa + cc - bb - dd, 2 * (cd + ab)], [2 * (bc - ad), aa + cc - bb - dd, 2 * (cd + ab)],
[2 * (bd + ac), 2 * (cd - ab), aa + dd - bb - cc]]) [2 * (bd + ac), 2 * (cd - ab), aa + dd - bb - cc]])
#Definition of the Bethe-Slater function
def func_BS(x,a,b,c):
return 4*a*((x/c)**2)*(1-b*(x/c)**2)*np.exp(-(x/c)**2)
#Definition of the derivative of the Bethe-Slater function
def func_dBS(x,a,b,c):
return 4*a*((x/c)**2)*(1-b*(x/c)**2)*np.exp(-(x/c)**2)
# calculating precession field of spin Sr # calculating precession field of spin Sr
def calc_rot_vector(Sr,Sf): def calc_rot_vector(Sr,Sf):
rot = (J0/hbar)*(Sf-alpha*np.cross(Sf,Sr))/(1.0+alpha**2) rot = (J0/hbar)*(Sf-alpha*np.cross(Sf,Sr))/(1.0+alpha**2)
@ -65,6 +81,6 @@ for t in range (0,N):
# calc. average magnetization # calc. average magnetization
Sm = (S1+S2)*0.5 Sm = (S1+S2)*0.5
# calc. energy # calc. energy
en = -2.0*J0*(np.dot(S1,S2)) en = -J0*(np.dot(S1,S2))
# print res. in ps for comparison with LAMMPS # print res. in ps for comparison with LAMMPS
print(t*dt/1000.0,Sm[0],Sm[1],Sm[2],en) print(t*dt/1000.0,Sm[0],Sm[1],Sm[2],en)

View File

@ -13,7 +13,7 @@ en="$(echo "$en-$in" | bc -l)"
tail -n +$in log.lammps | head -n $en > res_lammps.dat tail -n +$in log.lammps | head -n $en > res_lammps.dat
# compute Langevin # compute Langevin
python3 -m llg_exchange.py > res_llg.dat python3 llg_exchange.py > res_llg.dat
# plot results # plot results
python3 -m plot_precession.py res_lammps.dat res_llg.dat python3 plot_precession.py res_lammps.dat res_llg.dat

View File

@ -5,22 +5,24 @@ atom_style spin
atom_modify map array atom_modify map array
boundary f f f boundary f f f
read_data two_spins.data atom_modify map array
lattice sc 3.0
region box block 0 2 0 1 0 1
create_box 1 box
create_atoms 1 box
mass 1 55.845
set atom 1 spin 2.0 1.0 0.0 0.0
set atom 2 spin 2.0 0.0 1.0 0.0
pair_style spin/exchange 3.1 pair_style spin/exchange 3.1
pair_coeff * * exchange 3.1 11.254 0.0 1.0 pair_coeff * * exchange 3.1 11.254 0.0 1.0
group bead type 1
variable H equal 0.0
variable Kan equal 0.0
variable Temperature equal 0.0 variable Temperature equal 0.0
variable RUN equal 30000 variable RUN equal 30000
fix 1 all nve/spin lattice no fix 1 all nve/spin lattice frozen
fix 2 all precession/spin zeeman ${H} 0.0 0.0 1.0 anisotropy ${Kan} 0.0 0.0 1.0 fix 2 all langevin/spin ${Temperature} 0.01 12345
fix_modify 2 energy yes
fix 3 all langevin/spin ${Temperature} 0.01 12345
compute out_mag all spin compute out_mag all spin
compute out_pe all pe compute out_pe all pe
@ -34,6 +36,9 @@ variable emag equal c_out_mag[5]
thermo_style custom step time v_magx v_magy v_magz v_emag pe etotal thermo_style custom step time v_magx v_magy v_magz v_emag pe etotal
thermo 10 thermo 10
compute outsp all property/atom spx spy spz sp fmx fmy fmz
dump 1 all custom 10 dump.data type x y z c_outsp[1] c_outsp[2] c_outsp[3] fx fy fz
timestep 0.0001 timestep 0.0001
run ${RUN} run ${RUN}

View File

@ -1,22 +0,0 @@
LAMMPS data file via write_data, version 19 Sep 2019, timestep = 0
2 atoms
1 atom types
0.0 6.0 xlo xhi
0.0 3.0 ylo yhi
0.0 3.0 zlo zhi
Masses
1 1
Atoms # spin
1 1 2.0 0.0 0.0 0.0 1.0 0.0 0.0 0 0 0
2 1 2.0 3.0 0.0 0.0 0.0 1.0 0.0 0 0 0
Velocities
1 0.0 0.0 0.0
2 0.0 0.0 0.0

View File

@ -13,4 +13,4 @@ en="$(echo "$en-$in" | bc -l)"
tail -n +$in log.lammps | head -n $en > res_lammps.dat tail -n +$in log.lammps | head -n $en > res_lammps.dat
# plot results # plot results
python3 -m plot_nve.py res_lammps.dat res_llg.dat python3 plot_nve.py res_lammps.dat res_llg.dat

View File

@ -30,7 +30,7 @@ neighbor 0.1 bin
neigh_modify every 10 check yes delay 20 neigh_modify every 10 check yes delay 20
fix 1 all precession/spin zeeman 0.0 0.0 0.0 1.0 fix 1 all precession/spin zeeman 0.0 0.0 0.0 1.0
fix 2 all langevin 200.0 200.0 10.0 48279 fix 2 all langevin 200.0 200.0 1.0 48279
fix 3 all langevin/spin 0.0 0.00001 321 fix 3 all langevin/spin 0.0 0.00001 321
fix 4 all nve/spin lattice moving fix 4 all nve/spin lattice moving
timestep 0.001 timestep 0.001

View File

@ -29,7 +29,7 @@ neighbor 0.1 bin
neigh_modify every 10 check yes delay 20 neigh_modify every 10 check yes delay 20
fix 1 all precession/spin zeeman 0.0 0.0 0.0 1.0 fix 1 all precession/spin zeeman 0.0 0.0 0.0 1.0
fix 2 all langevin/spin 200.0 0.1 321 fix 2 all langevin/spin 200.0 0.01 321
fix 3 all nve/spin lattice moving fix 3 all nve/spin lattice moving
timestep 0.001 timestep 0.001

View File

@ -39,5 +39,5 @@ plt.xlabel('Time (in ps)')
plt.legend() plt.legend()
plt.show() plt.show()
fig.savefig(os.path.join(os.getcwd(), "nve_spin_lattice.pdf"), bbox_inches="tight") fig.savefig(os.path.join(os.getcwd(), "nvt_spin_lattice.pdf"), bbox_inches="tight")
plt.close(fig) plt.close(fig)

View File

@ -41,7 +41,10 @@ processors * * 1
region boxreg block 0 ${boxx} 0 ${boxy} 0 ${boxz} region boxreg block 0 ${boxx} 0 ${boxy} 0 ${boxz}
create_box 2 boxreg create_box 2 boxreg
change_box all boundary p p f change_box all boundary p p f
comm_modify vel yes
pair_style granular
pair_coeff 1 * hertz/material 1e5 0.2 0.3 tangential mindlin NULL 1.0 0.5 damping tsuji
pair_coeff 2 2 jkr 1e5 0.1 0.3 50 tangential mindlin NULL 1.0 0.5 rolling sds 1e3 1e3 0.1 twisting marshall damping tsuji
variable theta equal 0 variable theta equal 0
@ -61,10 +64,6 @@ comm_modify vel yes
neighbor ${skin} bin neighbor ${skin} bin
neigh_modify delay 0 every 1 check yes neigh_modify delay 0 every 1 check yes
pair_style granular
pair_coeff 1 * hertz/material 1e5 0.2 0.3 tangential mindlin NULL 1.0 0.5 damping tsuji
pair_coeff 2 2 jkr 1e5 0.1 0.3 50 tangential mindlin NULL 1.0 0.5 rolling sds 1e3 1e3 0.1 twisting marshall damping tsuji
fix 3 all wall/gran/region granular hertz/material 1e5 0.1 0.3 tangential mindlin NULL 1.0 0.5 damping tsuji region curved_wall fix 3 all wall/gran/region granular hertz/material 1e5 0.1 0.3 tangential mindlin NULL 1.0 0.5 damping tsuji region curved_wall
fix 4 all wall/gran/region granular hertz/material 1e5 0.1 0.3 tangential mindlin NULL 1.0 0.5 damping tsuji region bottom_wall fix 4 all wall/gran/region granular hertz/material 1e5 0.1 0.3 tangential mindlin NULL 1.0 0.5 damping tsuji region bottom_wall

View File

@ -35,6 +35,10 @@ region boxreg block 0 ${boxx} 0 ${boxy} 0 ${boxz}
create_box 2 boxreg create_box 2 boxreg
change_box all boundary p p f change_box all boundary p p f
pair_style granular
pair_coeff 1 * jkr 1000.0 50.0 0.3 10 tangential mindlin 800.0 1.0 0.5 rolling sds 500.0 200.0 0.5 twisting marshall
pair_coeff 2 2 hertz 200.0 20.0 tangential linear_history 300.0 1.0 0.1 rolling sds 200.0 100.0 0.1 twisting marshall
comm_modify vel yes comm_modify vel yes
region insreg1 cylinder z ${xc1} ${yc} 5 15 ${boxz} region insreg1 cylinder z ${xc1} ${yc} 5 15 ${boxz}
@ -45,18 +49,12 @@ fix grav all gravity 10.0 vector 0 0 -1
fix ins1 all pour 1500 1 3123 region insreg1 diam range ${dlo} ${dhi} dens ${dens} ${dens} fix ins1 all pour 1500 1 3123 region insreg1 diam range ${dlo} ${dhi} dens ${dens} ${dens}
fix ins2 all pour 1500 2 3123 region insreg2 diam range ${dlo} ${dhi} dens ${dens} ${dens} fix ins2 all pour 1500 2 3123 region insreg2 diam range ${dlo} ${dhi} dens ${dens} ${dens}
comm_modify vel yes
neighbor ${skin} bin neighbor ${skin} bin
neigh_modify delay 0 every 1 check yes neigh_modify delay 0 every 1 check yes
pair_style granular
pair_coeff 1 * jkr 1000.0 50.0 0.3 10 tangential mindlin 800.0 1.0 0.5 rolling sds 500.0 200.0 0.5 twisting marshall
pair_coeff 2 2 hertz 200.0 20.0 tangential linear_history 300.0 1.0 0.1 rolling sds 200.0 100.0 0.1 twisting marshall
fix 3 all wall/gran granular hertz/material 1e5 1e3 0.3 tangential mindlin NULL 1.0 0.5 zplane 0 NULL fix 3 all wall/gran granular hertz/material 1e5 1e3 0.3 tangential mindlin NULL 1.0 0.5 zplane 0 NULL
thermo_style custom step cpu atoms ke thermo_style custom step atoms ke
thermo_modify lost warn thermo_modify lost warn
thermo 100 thermo 100

View File

@ -1,5 +1,4 @@
LAMMPS (29 Mar 2019) LAMMPS (29 Oct 2020)
OMP_NUM_THREADS environment is not set. Defaulting to 1 thread. (../comm.cpp:88)
using 1 OpenMP thread(s) per MPI task using 1 OpenMP thread(s) per MPI task
# pour two types of particles (cohesive and non-cohesive) into cylinder # pour two types of particles (cohesive and non-cohesive) into cylinder
# 'turn' cylinder by changing direction of gravity, then rotate it. # 'turn' cylinder by changing direction of gravity, then rotate it.
@ -55,10 +54,14 @@ region boxreg block 0 30 0 ${boxy} 0 ${boxz}
region boxreg block 0 30 0 30 0 ${boxz} region boxreg block 0 30 0 30 0 ${boxz}
region boxreg block 0 30 0 30 0 50 region boxreg block 0 30 0 30 0 50
create_box 2 boxreg create_box 2 boxreg
Created orthogonal box = (0 0 0) to (30 30 50) Created orthogonal box = (0.0000000 0.0000000 0.0000000) to (30.000000 30.000000 50.000000)
1 by 1 by 1 MPI processor grid 1 by 1 by 1 MPI processor grid
change_box all boundary p p f change_box all boundary p p f
comm_modify vel yes Changing box ...
pair_style granular
pair_coeff 1 * hertz/material 1e5 0.2 0.3 tangential mindlin NULL 1.0 0.5 damping tsuji
pair_coeff 2 2 jkr 1e5 0.1 0.3 50 tangential mindlin NULL 1.0 0.5 rolling sds 1e3 1e3 0.1 twisting marshall damping tsuji
variable theta equal 0 variable theta equal 0
@ -104,10 +107,6 @@ neighbor ${skin} bin
neighbor 0.2 bin neighbor 0.2 bin
neigh_modify delay 0 every 1 check yes neigh_modify delay 0 every 1 check yes
pair_style granular
pair_coeff 1 * hertz/material 1e5 0.2 0.3 tangential mindlin NULL 1.0 0.5 damping tsuji
pair_coeff 2 2 jkr 1e5 0.1 0.3 50 tangential mindlin NULL 1.0 0.5 rolling sds 1e3 1e3 0.1 twisting marshall damping tsuji
fix 3 all wall/gran/region granular hertz/material 1e5 0.1 0.3 tangential mindlin NULL 1.0 0.5 damping tsuji region curved_wall fix 3 all wall/gran/region granular hertz/material 1e5 0.1 0.3 tangential mindlin NULL 1.0 0.5 damping tsuji region curved_wall
fix 4 all wall/gran/region granular hertz/material 1e5 0.1 0.3 tangential mindlin NULL 1.0 0.5 damping tsuji region bottom_wall fix 4 all wall/gran/region granular hertz/material 1e5 0.1 0.3 tangential mindlin NULL 1.0 0.5 damping tsuji region bottom_wall
@ -117,8 +116,7 @@ thermo 100
timestep 0.001 timestep 0.001
dump 1 all custom 100 ${name}.dump id type radius mass x y z #dump 1 all custom 100 ${name}.dump id type radius mass x y z
dump 1 all custom 100 rotating_drum_two_types.dump id type radius mass x y z
#For removal later #For removal later
compute 1 all property/atom radius compute 1 all property/atom radius
@ -140,7 +138,7 @@ Neighbor list info ...
pair build: half/size/bin/newton pair build: half/size/bin/newton
stencil: half/bin/3d/newton stencil: half/bin/3d/newton
bin: standard bin: standard
Per MPI rank memory allocation (min/avg/max) = 13.02 | 13.02 | 13.02 Mbytes Per MPI rank memory allocation (min/avg/max) = 13.03 | 13.03 | 13.03 Mbytes
Step Atoms KinEng v_theta Step Atoms KinEng v_theta
0 0 -0 0 0 0 -0 0
100 4000 -0 0 100 4000 -0 0
@ -163,32 +161,32 @@ Step Atoms KinEng v_theta
1800 4000 -0 0 1800 4000 -0 0
1900 4000 -0 0 1900 4000 -0 0
2000 4000 -0 0 2000 4000 -0 0
Loop time of 3.54461 on 1 procs for 2000 steps with 4000 atoms Loop time of 10.5178 on 1 procs for 2000 steps with 4000 atoms
Performance: 48750.057 tau/day, 564.237 timesteps/s Performance: 16429.309 tau/day, 190.154 timesteps/s
99.5% CPU use with 1 MPI tasks x 1 OpenMP threads 99.8% CPU use with 1 MPI tasks x 1 OpenMP threads
MPI task timing breakdown: MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total Section | min time | avg time | max time |%varavg| %total
--------------------------------------------------------------- ---------------------------------------------------------------
Pair | 0.61949 | 0.61949 | 0.61949 | 0.0 | 17.48 Pair | 1.0701 | 1.0701 | 1.0701 | 0.0 | 10.17
Neigh | 1.2492 | 1.2492 | 1.2492 | 0.0 | 35.24 Neigh | 4.2135 | 4.2135 | 4.2135 | 0.0 | 40.06
Comm | 0.046404 | 0.046404 | 0.046404 | 0.0 | 1.31 Comm | 0.38276 | 0.38276 | 0.38276 | 0.0 | 3.64
Output | 0.15901 | 0.15901 | 0.15901 | 0.0 | 4.49 Output | 0.0013647 | 0.0013647 | 0.0013647 | 0.0 | 0.01
Modify | 1.4165 | 1.4165 | 1.4165 | 0.0 | 39.96 Modify | 4.7076 | 4.7076 | 4.7076 | 0.0 | 44.76
Other | | 0.05391 | | | 1.52 Other | | 0.1424 | | | 1.35
Nlocal: 4000 ave 4000 max 4000 min Nlocal: 4000.00 ave 4000 max 4000 min
Histogram: 1 0 0 0 0 0 0 0 0 0 Histogram: 1 0 0 0 0 0 0 0 0 0
Nghost: 166 ave 166 max 166 min Nghost: 171.000 ave 171 max 171 min
Histogram: 1 0 0 0 0 0 0 0 0 0 Histogram: 1 0 0 0 0 0 0 0 0 0
Neighs: 8195 ave 8195 max 8195 min Neighs: 8093.00 ave 8093 max 8093 min
Histogram: 1 0 0 0 0 0 0 0 0 0 Histogram: 1 0 0 0 0 0 0 0 0 0
Total # of neighbors = 8195 Total # of neighbors = 8093
Ave neighs/atom = 2.04875 Ave neighs/atom = 2.0232500
Neighbor list builds = 1004 Neighbor list builds = 1004
Dangerous builds = 3 Dangerous builds = 4
#Remove any particles that are above z > 0.5*drum_height #Remove any particles that are above z > 0.5*drum_height
delete_atoms group delgroup delete_atoms group delgroup
@ -205,67 +203,69 @@ fix 5 all wall/gran/region granular hertz/material 1e5 0.1 0.3 tangential mindl
# 'Turn' drum by switching the direction of gravity # 'Turn' drum by switching the direction of gravity
unfix grav unfix grav
unfix ins1
unfix ins2
fix grav all gravity 10 vector 0 -1 0 fix grav all gravity 10 vector 0 -1 0
variable theta equal 2*PI*elapsed/20000.0 variable theta equal 2*PI*elapsed/20000.0
run 3000 run 3000
Per MPI rank memory allocation (min/avg/max) = 24.81 | 24.81 | 24.81 Mbytes Per MPI rank memory allocation (min/avg/max) = 19.37 | 19.37 | 19.37 Mbytes
Step Atoms KinEng v_theta Step Atoms KinEng v_theta
2000 4000 64.333531 0 2000 4000 65.647582 0
2100 4000 106.69182 0.031415927 2100 4000 105.60001 0.031415927
2200 4000 121.8461 0.062831853 2200 4000 112.27573 0.062831853
2300 4000 88.767952 0.09424778 2300 4000 92.758671 0.09424778
2400 4000 82.850721 0.12566371 2400 4000 88.925835 0.12566371
2500 4000 91.683284 0.15707963 2500 4000 81.369163 0.15707963
2600 4000 31.56344 0.18849556 2600 4000 32.046943 0.18849556
2700 4000 4.5697672 0.21991149 2700 4000 4.1926368 0.21991149
2800 4000 3.9879051 0.25132741 2800 4000 3.9933453 0.25132741
2900 4000 4.4394235 0.28274334 2900 4000 4.5062193 0.28274334
3000 4000 5.1212931 0.31415927 3000 4000 5.3409521 0.31415927
3100 4000 5.8608892 0.34557519 3100 4000 6.0165991 0.34557519
3200 4000 6.600714 0.37699112 3200 4000 6.606767 0.37699112
3300 4000 7.3497851 0.40840704 3300 4000 7.3997751 0.40840704
3400 4000 8.0490988 0.43982297 3400 4000 8.1098807 0.43982297
3500 4000 8.6712396 0.4712389 3500 4000 8.6552424 0.4712389
3600 4000 9.1328667 0.50265482 3600 4000 9.8445204 0.50265482
3700 4000 9.4683561 0.53407075 3700 4000 10.098753 0.53407075
3800 4000 9.5878145 0.56548668 3800 4000 10.039489 0.56548668
3900 4000 9.387745 0.5969026 3900 4000 9.6376278 0.5969026
4000 4000 8.9117631 0.62831853 4000 4000 9.2598836 0.62831853
4100 4000 8.2344368 0.65973446 4100 4000 8.7116037 0.65973446
4200 4000 7.5335088 0.69115038 4200 4000 8.1274117 0.69115038
4300 4000 6.8426179 0.72256631 4300 4000 7.1487627 0.72256631
4400 4000 6.0567247 0.75398224 4400 4000 6.2253778 0.75398224
4500 4000 5.4166132 0.78539816 4500 4000 5.3061398 0.78539816
4600 4000 4.6012409 0.81681409 4600 4000 4.4319316 0.81681409
4700 4000 3.8314982 0.84823002 4700 4000 4.205607 0.84823002
4800 4000 3.1916415 0.87964594 4800 4000 3.2112987 0.87964594
4900 4000 2.7833964 0.91106187 4900 4000 2.6449777 0.91106187
5000 4000 2.5051362 0.9424778 5000 4000 2.3475497 0.9424778
Loop time of 11.9545 on 1 procs for 3000 steps with 4000 atoms Loop time of 32.4926 on 1 procs for 3000 steps with 4000 atoms
Performance: 21682.142 tau/day, 250.951 timesteps/s Performance: 7977.205 tau/day, 92.329 timesteps/s
99.7% CPU use with 1 MPI tasks x 1 OpenMP threads 99.9% CPU use with 1 MPI tasks x 1 OpenMP threads
MPI task timing breakdown: MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total Section | min time | avg time | max time |%varavg| %total
--------------------------------------------------------------- ---------------------------------------------------------------
Pair | 4.8291 | 4.8291 | 4.8291 | 0.0 | 40.40 Pair | 8.0124 | 8.0124 | 8.0124 | 0.0 | 24.66
Neigh | 2.7489 | 2.7489 | 2.7489 | 0.0 | 22.99 Neigh | 10.993 | 10.993 | 10.993 | 0.0 | 33.83
Comm | 0.071249 | 0.071249 | 0.071249 | 0.0 | 0.60 Comm | 0.86697 | 0.86697 | 0.86697 | 0.0 | 2.67
Output | 0.20547 | 0.20547 | 0.20547 | 0.0 | 1.72 Output | 0.0021827 | 0.0021827 | 0.0021827 | 0.0 | 0.01
Modify | 4.0179 | 4.0179 | 4.0179 | 0.0 | 33.61 Modify | 12.367 | 12.367 | 12.367 | 0.0 | 38.06
Other | | 0.0819 | | | 0.69 Other | | 0.2515 | | | 0.77
Nlocal: 4000 ave 4000 max 4000 min Nlocal: 4000.00 ave 4000 max 4000 min
Histogram: 1 0 0 0 0 0 0 0 0 0 Histogram: 1 0 0 0 0 0 0 0 0 0
Nghost: 322 ave 322 max 322 min Nghost: 318.000 ave 318 max 318 min
Histogram: 1 0 0 0 0 0 0 0 0 0 Histogram: 1 0 0 0 0 0 0 0 0 0
Neighs: 14849 ave 14849 max 14849 min Neighs: 14807.0 ave 14807 max 14807 min
Histogram: 1 0 0 0 0 0 0 0 0 0 Histogram: 1 0 0 0 0 0 0 0 0 0
Total # of neighbors = 14849 Total # of neighbors = 14807
Ave neighs/atom = 3.71225 Ave neighs/atom = 3.7017500
Neighbor list builds = 1290 Neighbor list builds = 2189
Dangerous builds = 672 Dangerous builds = 1536
Total wall time: 0:00:15 Total wall time: 0:00:43

View File

@ -1,5 +1,4 @@
LAMMPS (29 Mar 2019) LAMMPS (29 Oct 2020)
OMP_NUM_THREADS environment is not set. Defaulting to 1 thread. (../comm.cpp:88)
using 1 OpenMP thread(s) per MPI task using 1 OpenMP thread(s) per MPI task
# pour two types of particles (cohesive and non-cohesive) into cylinder # pour two types of particles (cohesive and non-cohesive) into cylinder
# 'turn' cylinder by changing direction of gravity, then rotate it. # 'turn' cylinder by changing direction of gravity, then rotate it.
@ -55,10 +54,14 @@ region boxreg block 0 30 0 ${boxy} 0 ${boxz}
region boxreg block 0 30 0 30 0 ${boxz} region boxreg block 0 30 0 30 0 ${boxz}
region boxreg block 0 30 0 30 0 50 region boxreg block 0 30 0 30 0 50
create_box 2 boxreg create_box 2 boxreg
Created orthogonal box = (0 0 0) to (30 30 50) Created orthogonal box = (0.0000000 0.0000000 0.0000000) to (30.000000 30.000000 50.000000)
2 by 2 by 1 MPI processor grid 2 by 2 by 1 MPI processor grid
change_box all boundary p p f change_box all boundary p p f
comm_modify vel yes Changing box ...
pair_style granular
pair_coeff 1 * hertz/material 1e5 0.2 0.3 tangential mindlin NULL 1.0 0.5 damping tsuji
pair_coeff 2 2 jkr 1e5 0.1 0.3 50 tangential mindlin NULL 1.0 0.5 rolling sds 1e3 1e3 0.1 twisting marshall damping tsuji
variable theta equal 0 variable theta equal 0
@ -104,10 +107,6 @@ neighbor ${skin} bin
neighbor 0.2 bin neighbor 0.2 bin
neigh_modify delay 0 every 1 check yes neigh_modify delay 0 every 1 check yes
pair_style granular
pair_coeff 1 * hertz/material 1e5 0.2 0.3 tangential mindlin NULL 1.0 0.5 damping tsuji
pair_coeff 2 2 jkr 1e5 0.1 0.3 50 tangential mindlin NULL 1.0 0.5 rolling sds 1e3 1e3 0.1 twisting marshall damping tsuji
fix 3 all wall/gran/region granular hertz/material 1e5 0.1 0.3 tangential mindlin NULL 1.0 0.5 damping tsuji region curved_wall fix 3 all wall/gran/region granular hertz/material 1e5 0.1 0.3 tangential mindlin NULL 1.0 0.5 damping tsuji region curved_wall
fix 4 all wall/gran/region granular hertz/material 1e5 0.1 0.3 tangential mindlin NULL 1.0 0.5 damping tsuji region bottom_wall fix 4 all wall/gran/region granular hertz/material 1e5 0.1 0.3 tangential mindlin NULL 1.0 0.5 damping tsuji region bottom_wall
@ -117,8 +116,7 @@ thermo 100
timestep 0.001 timestep 0.001
dump 1 all custom 100 ${name}.dump id type radius mass x y z #dump 1 all custom 100 ${name}.dump id type radius mass x y z
dump 1 all custom 100 rotating_drum_two_types.dump id type radius mass x y z
#For removal later #For removal later
compute 1 all property/atom radius compute 1 all property/atom radius
@ -163,32 +161,32 @@ Step Atoms KinEng v_theta
1800 4000 -0 0 1800 4000 -0 0
1900 4000 -0 0 1900 4000 -0 0
2000 4000 -0 0 2000 4000 -0 0
Loop time of 2.0709 on 4 procs for 2000 steps with 4000 atoms Loop time of 3.86825 on 4 procs for 2000 steps with 4000 atoms
Performance: 83442.024 tau/day, 965.764 timesteps/s Performance: 44671.398 tau/day, 517.030 timesteps/s
97.7% CPU use with 4 MPI tasks x 1 OpenMP threads 96.7% CPU use with 4 MPI tasks x 1 OpenMP threads
MPI task timing breakdown: MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total Section | min time | avg time | max time |%varavg| %total
--------------------------------------------------------------- ---------------------------------------------------------------
Pair | 0.24679 | 0.26336 | 0.28853 | 3.0 | 12.72 Pair | 0.26114 | 0.27918 | 0.28728 | 2.0 | 7.22
Neigh | 0.52279 | 0.5332 | 0.53858 | 0.9 | 25.75 Neigh | 1.2044 | 1.2414 | 1.3105 | 3.7 | 32.09
Comm | 0.17418 | 0.20253 | 0.23266 | 4.7 | 9.78 Comm | 0.38592 | 0.47065 | 0.51052 | 7.4 | 12.17
Output | 0.092897 | 0.093531 | 0.09515 | 0.3 | 4.52 Output | 0.0007236 | 0.0013456 | 0.0024846 | 1.8 | 0.03
Modify | 0.88151 | 0.89571 | 0.90582 | 0.9 | 43.25 Modify | 1.6217 | 1.6723 | 1.7801 | 5.0 | 43.23
Other | | 0.08257 | | | 3.99 Other | | 0.2034 | | | 5.26
Nlocal: 1000 ave 1001 max 999 min Nlocal: 1000.00 ave 1012 max 988 min
Histogram: 2 0 0 0 0 0 0 0 0 2 Histogram: 2 0 0 0 0 0 0 0 0 2
Nghost: 267.75 ave 276 max 262 min Nghost: 269.250 ave 278 max 256 min
Histogram: 1 0 1 0 1 0 0 0 0 1 Histogram: 1 0 0 0 0 0 1 1 0 1
Neighs: 2031.5 ave 2091 max 1958 min Neighs: 2060.50 ave 2156 max 1921 min
Histogram: 1 0 0 0 1 0 0 1 0 1 Histogram: 1 0 0 1 0 0 0 0 0 2
Total # of neighbors = 8126 Total # of neighbors = 8242
Ave neighs/atom = 2.0315 Ave neighs/atom = 2.0605000
Neighbor list builds = 1004 Neighbor list builds = 1004
Dangerous builds = 3 Dangerous builds = 4
#Remove any particles that are above z > 0.5*drum_height #Remove any particles that are above z > 0.5*drum_height
delete_atoms group delgroup delete_atoms group delgroup
@ -205,67 +203,69 @@ fix 5 all wall/gran/region granular hertz/material 1e5 0.1 0.3 tangential mindl
# 'Turn' drum by switching the direction of gravity # 'Turn' drum by switching the direction of gravity
unfix grav unfix grav
unfix ins1
unfix ins2
fix grav all gravity 10 vector 0 -1 0 fix grav all gravity 10 vector 0 -1 0
variable theta equal 2*PI*elapsed/20000.0 variable theta equal 2*PI*elapsed/20000.0
run 3000 run 3000
Per MPI rank memory allocation (min/avg/max) = 21.6 | 22.6 | 23.82 Mbytes Per MPI rank memory allocation (min/avg/max) = 18.55 | 18.55 | 18.55 Mbytes
Step Atoms KinEng v_theta Step Atoms KinEng v_theta
2000 4000 64.255821 0 2000 4000 65.819213 0
2100 4000 106.47082 0.031415927 2100 4000 105.02389 0.031415927
2200 4000 121.52634 0.062831853 2200 4000 112.02469 0.062831853
2300 4000 87.748818 0.09424778 2300 4000 92.271262 0.09424778
2400 4000 82.712784 0.12566371 2400 4000 89.369506 0.12566371
2500 4000 90.618713 0.15707963 2500 4000 80.910925 0.15707963
2600 4000 30.096031 0.18849556 2600 4000 31.620722 0.18849556
2700 4000 4.0838611 0.21991149 2700 4000 4.3019937 0.21991149
2800 4000 3.7485959 0.25132741 2800 4000 3.9913967 0.25132741
2900 4000 4.2159774 0.28274334 2900 4000 4.5203726 0.28274334
3000 4000 4.8730048 0.31415927 3000 4000 5.484886 0.31415927
3100 4000 5.6109465 0.34557519 3100 4000 6.1085958 0.34557519
3200 4000 6.4290528 0.37699112 3200 4000 6.7085635 0.37699112
3300 4000 7.2699677 0.40840704 3300 4000 7.4787777 0.40840704
3400 4000 8.0895944 0.43982297 3400 4000 8.2116413 0.43982297
3500 4000 8.7222781 0.4712389 3500 4000 8.7979302 0.4712389
3600 4000 9.133205 0.50265482 3600 4000 9.871649 0.50265482
3700 4000 9.3404584 0.53407075 3700 4000 10.012426 0.53407075
3800 4000 9.3359844 0.56548668 3800 4000 9.9067754 0.56548668
3900 4000 9.0916854 0.5969026 3900 4000 9.725458 0.5969026
4000 4000 8.5596424 0.62831853 4000 4000 9.3350056 0.62831853
4100 4000 7.9734883 0.65973446 4100 4000 8.8337295 0.65973446
4200 4000 7.2154383 0.69115038 4200 4000 8.2712493 0.69115038
4300 4000 6.7039232 0.72256631 4300 4000 6.9609934 0.72256631
4400 4000 6.1542738 0.75398224 4400 4000 6.0120294 0.75398224
4500 4000 5.4049454 0.78539816 4500 4000 5.0490036 0.78539816
4600 4000 4.4603192 0.81681409 4600 4000 4.2796544 0.81681409
4700 4000 3.6197985 0.84823002 4700 4000 4.1736483 0.84823002
4800 4000 2.9895571 0.87964594 4800 4000 3.0860106 0.87964594
4900 4000 2.5314553 0.91106187 4900 4000 2.6670909 0.91106187
5000 4000 2.2645533 0.9424778 5000 4000 2.2901814 0.9424778
Loop time of 6.64209 on 4 procs for 3000 steps with 4000 atoms Loop time of 10.7627 on 4 procs for 3000 steps with 4000 atoms
Performance: 39023.861 tau/day, 451.665 timesteps/s Performance: 24083.252 tau/day, 278.741 timesteps/s
96.6% CPU use with 4 MPI tasks x 1 OpenMP threads 97.9% CPU use with 4 MPI tasks x 1 OpenMP threads
MPI task timing breakdown: MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total Section | min time | avg time | max time |%varavg| %total
--------------------------------------------------------------- ---------------------------------------------------------------
Pair | 1.8376 | 2.126 | 2.3131 | 12.6 | 32.01 Pair | 1.6731 | 2.0701 | 2.3327 | 18.9 | 19.23
Neigh | 0.97762 | 1.0518 | 1.1337 | 5.4 | 15.84 Neigh | 2.7389 | 3.1706 | 3.5146 | 15.7 | 29.46
Comm | 0.53699 | 0.84265 | 1.2325 | 27.6 | 12.69 Comm | 0.93507 | 1.5441 | 2.1182 | 39.1 | 14.35
Output | 0.13922 | 0.14159 | 0.14388 | 0.4 | 2.13 Output | 0.0021682 | 0.0044412 | 0.006026 | 2.2 | 0.04
Modify | 1.8815 | 2.1026 | 2.3368 | 11.2 | 31.66 Modify | 3.0031 | 3.4223 | 3.9262 | 18.3 | 31.80
Other | | 0.3774 | | | 5.68 Other | | 0.5511 | | | 5.12
Nlocal: 1000 ave 1256 max 744 min Nlocal: 1000.00 ave 1277 max 723 min
Histogram: 2 0 0 0 0 0 0 0 0 2 Histogram: 2 0 0 0 0 0 0 0 0 2
Nghost: 579.5 ave 789 max 498 min Nghost: 569.750 ave 809 max 454 min
Histogram: 2 1 0 0 0 0 0 0 0 1 Histogram: 1 2 0 0 0 0 0 0 0 1
Neighs: 3696.25 ave 4853 max 2590 min Neighs: 3690.50 ave 4937 max 2426 min
Histogram: 2 0 0 0 0 0 0 0 1 1 Histogram: 1 1 0 0 0 0 0 0 0 2
Total # of neighbors = 14785 Total # of neighbors = 14762
Ave neighs/atom = 3.69625 Ave neighs/atom = 3.6905000
Neighbor list builds = 1230 Neighbor list builds = 2187
Dangerous builds = 676 Dangerous builds = 1610
Total wall time: 0:00:08 Total wall time: 0:00:14

View File

@ -1,5 +1,4 @@
LAMMPS (29 Mar 2019) LAMMPS (29 Oct 2020)
OMP_NUM_THREADS environment is not set. Defaulting to 1 thread. (../comm.cpp:88)
using 1 OpenMP thread(s) per MPI task using 1 OpenMP thread(s) per MPI task
# pour two types of particles (cohesive and non-cohesive) on flat wall # pour two types of particles (cohesive and non-cohesive) on flat wall
@ -45,9 +44,14 @@ region boxreg block 0 20 0 ${boxy} 0 ${boxz}
region boxreg block 0 20 0 20 0 ${boxz} region boxreg block 0 20 0 20 0 ${boxz}
region boxreg block 0 20 0 20 0 30 region boxreg block 0 20 0 20 0 30
create_box 2 boxreg create_box 2 boxreg
Created orthogonal box = (0 0 0) to (20 20 30) Created orthogonal box = (0.0000000 0.0000000 0.0000000) to (20.000000 20.000000 30.000000)
1 by 1 by 1 MPI processor grid 1 by 1 by 1 MPI processor grid
change_box all boundary p p f change_box all boundary p p f
Changing box ...
pair_style granular
pair_coeff 1 * jkr 1000.0 50.0 0.3 10 tangential mindlin 800.0 1.0 0.5 rolling sds 500.0 200.0 0.5 twisting marshall
pair_coeff 2 2 hertz 200.0 20.0 tangential linear_history 300.0 1.0 0.1 rolling sds 200.0 100.0 0.1 twisting marshall
comm_modify vel yes comm_modify vel yes
@ -75,19 +79,13 @@ fix ins2 all pour 1500 2 3123 region insreg2 diam range 0.5 1 dens 1 ${dens}
fix ins2 all pour 1500 2 3123 region insreg2 diam range 0.5 1 dens 1 1 fix ins2 all pour 1500 2 3123 region insreg2 diam range 0.5 1 dens 1 1
Particle insertion: 562 every 346 steps, 1500 by step 693 Particle insertion: 562 every 346 steps, 1500 by step 693
comm_modify vel yes
neighbor ${skin} bin neighbor ${skin} bin
neighbor 0.15 bin neighbor 0.15 bin
neigh_modify delay 0 every 1 check yes neigh_modify delay 0 every 1 check yes
pair_style granular
pair_coeff 1 * jkr 1000.0 50.0 0.3 10 tangential mindlin 800.0 1.0 0.5 rolling sds 500.0 200.0 0.5 twisting marshall
pair_coeff 2 2 hertz 200.0 20.0 tangential linear_history 300.0 1.0 0.1 rolling sds 200.0 100.0 0.1 twisting marshall
fix 3 all wall/gran granular hertz/material 1e5 1e3 0.3 tangential mindlin NULL 1.0 0.5 zplane 0 NULL fix 3 all wall/gran granular hertz/material 1e5 1e3 0.3 tangential mindlin NULL 1.0 0.5 zplane 0 NULL
thermo_style custom step cpu atoms ke thermo_style custom step atoms ke
thermo_modify lost warn thermo_modify lost warn
thermo 100 thermo 100
@ -109,26 +107,82 @@ Neighbor list info ...
stencil: half/bin/3d/newton stencil: half/bin/3d/newton
bin: standard bin: standard
Per MPI rank memory allocation (min/avg/max) = 12.22 | 12.22 | 12.22 Mbytes Per MPI rank memory allocation (min/avg/max) = 12.22 | 12.22 | 12.22 Mbytes
Step CPU Atoms KinEng Step Atoms KinEng
0 0 0 -0 0 0 -0
100 3.8153191 855 -0 100 926 -0
200 4.195287 855 -0 200 926 -0
300 4.5890362 855 -0 300 926 -0
400 10.636087 1500 -0 400 1498 -0
500 11.306909 1500 -0 500 1498 -0
600 11.968198 1500 -0 600 1498 -0
700 22.631892 2288 -0 700 2275 -0
800 23.711387 2288 -0 800 2275 -0
900 24.754344 2288 -0 900 2275 -0
1000 25.811778 2288 -0 1000 2275 -0
1100 35.368869 2845 -0 1100 2954 -0
1200 37.149843 2845 -0 1200 2954 -0
1300 39.026458 2845 -0 1300 2954 -0
1400 41.757583 3000 -0 1400 3000 -0
1500 45.155503 3000 -0 1500 3000 -0
1600 48.570241 3000 -0 1600 3000 -0
1700 52.839322 3000 -0 1700 3000 -0
1800 59.772697 3000 -0 1800 3000 -0
1900 69.493305 3000 -0 1900 3000 -0
2000 114.61886 3000 -0 2000 3000 -0
2100 152.89232 3000 -0 2100 3000 -0
2200 3000 -0
2300 3000 -0
2400 3000 -0
2500 3000 -0
2600 3000 -0
2700 3000 -0
2800 3000 -0
2900 3000 -0
3000 3000 -0
3100 3000 -0
3200 3000 -0
3300 3000 -0
3400 3000 -0
3500 3000 -0
3600 3000 -0
3700 3000 -0
3800 3000 -0
3900 3000 -0
4000 3000 -0
4100 3000 -0
4200 3000 -0
4300 3000 -0
4400 3000 -0
4500 3000 -0
4600 3000 -0
4700 3000 -0
4800 3000 -0
4900 3000 -0
5000 3000 -0
Loop time of 24.3889 on 1 procs for 5000 steps with 3000 atoms
Performance: 17713.003 tau/day, 205.012 timesteps/s
99.9% CPU use with 1 MPI tasks x 1 OpenMP threads
MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total
---------------------------------------------------------------
Pair | 14.362 | 14.362 | 14.362 | 0.0 | 58.89
Neigh | 3.3483 | 3.3483 | 3.3483 | 0.0 | 13.73
Comm | 0.42893 | 0.42893 | 0.42893 | 0.0 | 1.76
Output | 0.0025065 | 0.0025065 | 0.0025065 | 0.0 | 0.01
Modify | 6.059 | 6.059 | 6.059 | 0.0 | 24.84
Other | | 0.1876 | | | 0.77
Nlocal: 3000.00 ave 3000 max 3000 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Nghost: 462.000 ave 462 max 462 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Neighs: 17694.0 ave 17694 max 17694 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Total # of neighbors = 17694
Ave neighs/atom = 5.8980000
Neighbor list builds = 1133
Dangerous builds = 0
Total wall time: 0:00:24

View File

@ -1,5 +1,4 @@
LAMMPS (29 Mar 2019) LAMMPS (29 Oct 2020)
OMP_NUM_THREADS environment is not set. Defaulting to 1 thread. (../comm.cpp:88)
using 1 OpenMP thread(s) per MPI task using 1 OpenMP thread(s) per MPI task
# pour two types of particles (cohesive and non-cohesive) on flat wall # pour two types of particles (cohesive and non-cohesive) on flat wall
@ -45,9 +44,14 @@ region boxreg block 0 20 0 ${boxy} 0 ${boxz}
region boxreg block 0 20 0 20 0 ${boxz} region boxreg block 0 20 0 20 0 ${boxz}
region boxreg block 0 20 0 20 0 30 region boxreg block 0 20 0 20 0 30
create_box 2 boxreg create_box 2 boxreg
Created orthogonal box = (0 0 0) to (20 20 30) Created orthogonal box = (0.0000000 0.0000000 0.0000000) to (20.000000 20.000000 30.000000)
2 by 2 by 1 MPI processor grid 2 by 2 by 1 MPI processor grid
change_box all boundary p p f change_box all boundary p p f
Changing box ...
pair_style granular
pair_coeff 1 * jkr 1000.0 50.0 0.3 10 tangential mindlin 800.0 1.0 0.5 rolling sds 500.0 200.0 0.5 twisting marshall
pair_coeff 2 2 hertz 200.0 20.0 tangential linear_history 300.0 1.0 0.1 rolling sds 200.0 100.0 0.1 twisting marshall
comm_modify vel yes comm_modify vel yes
@ -75,26 +79,19 @@ fix ins2 all pour 1500 2 3123 region insreg2 diam range 0.5 1 dens 1 ${dens}
fix ins2 all pour 1500 2 3123 region insreg2 diam range 0.5 1 dens 1 1 fix ins2 all pour 1500 2 3123 region insreg2 diam range 0.5 1 dens 1 1
Particle insertion: 562 every 346 steps, 1500 by step 693 Particle insertion: 562 every 346 steps, 1500 by step 693
comm_modify vel yes
neighbor ${skin} bin neighbor ${skin} bin
neighbor 0.15 bin neighbor 0.15 bin
neigh_modify delay 0 every 1 check yes neigh_modify delay 0 every 1 check yes
pair_style granular
pair_coeff 1 * jkr 1000.0 50.0 0.3 10 tangential mindlin 800.0 1.0 0.5 rolling sds 500.0 200.0 0.5 twisting marshall
pair_coeff 2 2 hertz 200.0 20.0 tangential linear_history 300.0 1.0 0.1 rolling sds 200.0 100.0 0.1 twisting marshall
fix 3 all wall/gran granular hertz/material 1e5 1e3 0.3 tangential mindlin NULL 1.0 0.5 zplane 0 NULL fix 3 all wall/gran granular hertz/material 1e5 1e3 0.3 tangential mindlin NULL 1.0 0.5 zplane 0 NULL
thermo_style custom step cpu atoms ke thermo_style custom step atoms ke
thermo_modify lost warn thermo_modify lost warn
thermo 100 thermo 100
timestep 0.001 timestep 0.001
dump 1 all custom 100 ${name}.dump id type radius mass x y z #dump 1 all custom 100 ${name}.dump id type radius mass x y z
dump 1 all custom 100 pour_two_types.dump id type radius mass x y z
run 5000 run 5000
Neighbor list info ... Neighbor list info ...
@ -110,82 +107,82 @@ Neighbor list info ...
stencil: half/bin/3d/newton stencil: half/bin/3d/newton
bin: standard bin: standard
Per MPI rank memory allocation (min/avg/max) = 11.98 | 11.98 | 11.98 Mbytes Per MPI rank memory allocation (min/avg/max) = 11.98 | 11.98 | 11.98 Mbytes
Step CPU Atoms KinEng Step Atoms KinEng
0 0 0 -0 0 0 -0
100 0.11584234 855 -0 100 926 -0
200 0.12743592 855 -0 200 926 -0
300 0.13925815 855 -0 300 926 -0
400 0.35203671 1500 -0 400 1498 -0
500 0.37055922 1500 -0 500 1498 -0
600 0.38671875 1500 -0 600 1498 -0
700 0.71736908 2288 -0 700 2275 -0
800 0.74506783 2288 -0 800 2275 -0
900 0.77112222 2288 -0 900 2275 -0
1000 0.79632139 2288 -0 1000 2275 -0
1100 1.0384252 2845 -0 1100 2954 -0
1200 1.08093 2845 -0 1200 2954 -0
1300 1.1224561 2845 -0 1300 2954 -0
1400 1.1811485 3000 -0 1400 3000 -0
1500 1.2414908 3000 -0 1500 3000 -0
1600 1.3105879 3000 -0 1600 3000 -0
1700 1.390928 3000 -0 1700 3000 -0
1800 1.4869275 3000 -0 1800 3000 -0
1900 1.5958266 3000 -0 1900 3000 -0
2000 1.7172487 3000 -0 2000 3000 -0
2100 1.851155 3000 -0 2100 3000 -0
2200 1.9957182 3000 -0 2200 3000 -0
2300 2.1593764 3000 -0 2300 3000 -0
2400 2.3433132 3000 -0 2400 3000 -0
2500 2.532742 3000 -0 2500 3000 -0
2600 2.7376895 3000 -0 2600 3000 -0
2700 2.9463468 3000 -0 2700 3000 -0
2800 3.1645725 3000 -0 2800 3000 -0
2900 3.3879526 3000 -0 2900 3000 -0
3000 3.6152103 3000 -0 3000 3000 -0
3100 3.8467371 3000 -0 3100 3000 -0
3200 4.0787683 3000 -0 3200 3000 -0
3300 4.3097105 3000 -0 3300 3000 -0
3400 4.5423617 3000 -0 3400 3000 -0
3500 4.7773693 3000 -0 3500 3000 -0
3600 5.0127218 3000 -0 3600 3000 -0
3700 5.2519271 3000 -0 3700 3000 -0
3800 5.4951298 3000 -0 3800 3000 -0
3900 5.7210469 3000 -0 3900 3000 -0
4000 5.9432652 3000 -0 4000 3000 -0
4100 6.1687591 3000 -0 4100 3000 -0
4200 6.3942792 3000 -0 4200 3000 -0
4300 6.6331475 3000 -0 4300 3000 -0
4400 6.8632154 3000 -0 4400 3000 -0
4500 7.0979366 3000 -0 4500 3000 -0
4600 7.3305347 3000 -0 4600 3000 -0
4700 7.5670528 3000 -0 4700 3000 -0
4800 7.8086057 3000 -0 4800 3000 -0
4900 8.0407174 3000 -0 4900 3000 -0
5000 8.2765219 3000 -0 5000 3000 -0
Loop time of 8.27669 on 4 procs for 5000 steps with 3000 atoms Loop time of 12.1982 on 4 procs for 5000 steps with 3000 atoms
Performance: 52194.788 tau/day, 604.106 timesteps/s Performance: 35414.923 tau/day, 409.895 timesteps/s
97.7% CPU use with 4 MPI tasks x 1 OpenMP threads 97.0% CPU use with 4 MPI tasks x 1 OpenMP threads
MPI task timing breakdown: MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total Section | min time | avg time | max time |%varavg| %total
--------------------------------------------------------------- ---------------------------------------------------------------
Pair | 1.6106 | 3.4073 | 5.4191 | 95.7 | 41.17 Pair | 1.7141 | 3.8131 | 6.2143 | 107.3 | 31.26
Neigh | 0.51456 | 0.64572 | 0.81542 | 16.6 | 7.80 Neigh | 0.77648 | 0.96585 | 1.1892 | 18.3 | 7.92
Comm | 0.2808 | 2.5222 | 4.4998 | 121.9 | 30.47 Comm | 0.7427 | 3.5566 | 5.9731 | 128.4 | 29.16
Output | 0.15695 | 0.15919 | 0.16502 | 0.8 | 1.92 Output | 0.0067544 | 0.0086352 | 0.011408 | 1.8 | 0.07
Modify | 1.3517 | 1.4192 | 1.4904 | 4.9 | 17.15 Modify | 3.3476 | 3.5826 | 3.8235 | 11.5 | 29.37
Other | | 0.123 | | | 1.49 Other | | 0.2715 | | | 2.23
Nlocal: 750 ave 1036 max 482 min Nlocal: 750.000 ave 1033 max 463 min
Histogram: 2 0 0 0 0 0 0 0 1 1 Histogram: 2 0 0 0 0 0 0 0 0 2
Nghost: 429.75 ave 475 max 386 min Nghost: 435.000 ave 492 max 378 min
Histogram: 2 0 0 0 0 0 0 0 0 2
Neighs: 4434.50 ave 7028 max 1967 min
Histogram: 2 0 0 0 0 0 0 0 0 2 Histogram: 2 0 0 0 0 0 0 0 0 2
Neighs: 4051.75 ave 6274 max 2057 min
Histogram: 2 0 0 0 0 0 0 0 1 1
Total # of neighbors = 16207 Total # of neighbors = 17738
Ave neighs/atom = 5.40233 Ave neighs/atom = 5.9126667
Neighbor list builds = 1165 Neighbor list builds = 1139
Dangerous builds = 0 Dangerous builds = 0
Total wall time: 0:00:08 Total wall time: 0:00:12

View File

@ -18,7 +18,7 @@ profiling and debugging tools (https://github.com/kokkos/kokkos-tools).
A programming guide can be found on the Wiki, the API reference is under development. A programming guide can be found on the Wiki, the API reference is under development.
For questions find us on Slack: https://kokkosteam.slack.com or open a github issue. For questions find us on Slack: https://kokkosteam.slack.com or open a GitHub issue.
For non-public questions send an email to For non-public questions send an email to
crtrott(at)sandia.gov crtrott(at)sandia.gov
@ -44,7 +44,7 @@ To learn more about Kokkos consider watching one of our presentations:
We are open and try to encourage contributions from external developers. We are open and try to encourage contributions from external developers.
To do so please first open an issue describing the contribution and then issue To do so please first open an issue describing the contribution and then issue
a pull request against the develop branch. For larger features it may be good a pull request against the develop branch. For larger features it may be good
to get guidance from the core development team first through the github issue. to get guidance from the core development team first through the GitHub issue.
Note that Kokkos Core is licensed under standard 3-clause BSD terms of use. Note that Kokkos Core is licensed under standard 3-clause BSD terms of use.
Which means contributing to Kokkos allows anyone else to use your contributions Which means contributing to Kokkos allows anyone else to use your contributions

View File

@ -17,7 +17,7 @@ Building LAMMPS with QUIP support:
1) Building QUIP 1) Building QUIP
1.1) Obtaining QUIP 1.1) Obtaining QUIP
The most current release of QUIP can be obtained from github: The most current release of QUIP can be obtained from GitHub:
$ git clone https://github.com/libAtoms/QUIP.git QUIP $ git clone https://github.com/libAtoms/QUIP.git QUIP

View File

@ -3,7 +3,7 @@ is required to use the KSPACE scafacos and its kspace_style
scafacos command in a LAMMPS input script. scafacos command in a LAMMPS input script.
The ScaFaCoS library is available at http://scafacos.de or The ScaFaCoS library is available at http://scafacos.de or
on github at https://github.com/scafacos, the library was on GitHub at https://github.com/scafacos, the library was
developed by a consortium of different universities in developed by a consortium of different universities in
Germany (Bonn, Chemnitz, Stuttgart, Wuppertal) and Germany (Bonn, Chemnitz, Stuttgart, Wuppertal) and
the Research Centre Juelich (Juelich Supercomputing Centre). the Research Centre Juelich (Juelich Supercomputing Centre).

View File

@ -1102,15 +1102,8 @@ void PairGranular::init_style()
// this is so its order in the fix list is preserved // this is so its order in the fix list is preserved
if (use_history && fix_history == nullptr) { if (use_history && fix_history == nullptr) {
char dnumstr[16]; modify->replace_fix("NEIGH_HISTORY_GRANULAR_DUMMY","NEIGH_HISTORY_GRANULAR"
sprintf(dnumstr,"%d",size_history); " all NEIGH_HISTORY " + std::to_string(size_history),1);
char **fixarg = new char*[4];
fixarg[0] = (char *) "NEIGH_HISTORY_GRANULAR";
fixarg[1] = (char *) "all";
fixarg[2] = (char *) "NEIGH_HISTORY";
fixarg[3] = dnumstr;
modify->replace_fix("NEIGH_HISTORY_GRANULAR_DUMMY",4,fixarg,1);
delete [] fixarg;
int ifix = modify->find_fix("NEIGH_HISTORY_GRANULAR"); int ifix = modify->find_fix("NEIGH_HISTORY_GRANULAR");
fix_history = (FixNeighHistory *) modify->fix[ifix]; fix_history = (FixNeighHistory *) modify->fix[ifix];
fix_history->pair = this; fix_history->pair = this;
@ -1405,11 +1398,15 @@ double PairGranular::single(int i, int j, int itype, int jtype,
int *jlist; int *jlist;
double *history,*allhistory; double *history,*allhistory;
int nall = atom->nlocal + atom->nghost;
if ((i >= nall) || (j >= nall))
error->all(FLERR,"Not enough atoms for pair granular single function");
double *radius = atom->radius; double *radius = atom->radius;
radi = radius[i]; radi = radius[i];
radj = radius[j]; radj = radius[j];
radsum = radi + radj; radsum = radi + radj;
Reff = radi*radj/radsum; Reff = (radsum > 0.0) ? radi*radj/radsum : 0.0;
bool touchflag; bool touchflag;
E = normal_coeffs[itype][jtype][0]; E = normal_coeffs[itype][jtype][0];
@ -1536,6 +1533,8 @@ double PairGranular::single(int i, int j, int itype, int jtype,
jlist = list->firstneigh[i]; jlist = list->firstneigh[i];
if (use_history) { if (use_history) {
if ((fix_history == nullptr) || (fix_history->firstvalue == nullptr))
error->one(FLERR,"Pair granular single computation needs history");
allhistory = fix_history->firstvalue[i]; allhistory = fix_history->firstvalue[i];
for (int jj = 0; jj < jnum; jj++) { for (int jj = 0; jj < jnum; jj++) {
neighprev++; neighprev++;

View File

@ -228,7 +228,7 @@ void FFT3dKokkos<DeviceType>::fft_3d_kokkos(typename FFT_AT::t_FFT_DATA_1d d_in,
cufftExec(plan->plan_fast,d_data.data(),d_data.data(),flag); cufftExec(plan->plan_fast,d_data.data(),d_data.data(),flag);
#else #else
typename FFT_AT::t_FFT_DATA_1d d_tmp = typename FFT_AT::t_FFT_DATA_1d d_tmp =
typename FFT_AT::t_FFT_DATA_1d(Kokkos::view_alloc("fft_3d:tmp",Kokkos::WithoutInitializing),d_in.extent(0)); typename FFT_AT::t_FFT_DATA_1d(Kokkos::view_alloc("fft_3d:tmp",Kokkos::WithoutInitializing),d_data.extent(0));
kiss_fft_functor<DeviceType> f; kiss_fft_functor<DeviceType> f;
if (flag == -1) if (flag == -1)
f = kiss_fft_functor<DeviceType>(d_data,d_tmp,plan->cfg_fast_forward,length); f = kiss_fft_functor<DeviceType>(d_data,d_tmp,plan->cfg_fast_forward,length);
@ -236,7 +236,6 @@ void FFT3dKokkos<DeviceType>::fft_3d_kokkos(typename FFT_AT::t_FFT_DATA_1d d_in,
f = kiss_fft_functor<DeviceType>(d_data,d_tmp,plan->cfg_fast_backward,length); f = kiss_fft_functor<DeviceType>(d_data,d_tmp,plan->cfg_fast_backward,length);
Kokkos::parallel_for(total/length,f); Kokkos::parallel_for(total/length,f);
d_data = d_tmp; d_data = d_tmp;
d_tmp = typename FFT_AT::t_FFT_DATA_1d(Kokkos::view_alloc("fft_3d:tmp",Kokkos::WithoutInitializing),d_in.extent(0));
#endif #endif
@ -273,13 +272,13 @@ void FFT3dKokkos<DeviceType>::fft_3d_kokkos(typename FFT_AT::t_FFT_DATA_1d d_in,
#elif defined(FFT_CUFFT) #elif defined(FFT_CUFFT)
cufftExec(plan->plan_mid,d_data.data(),d_data.data(),flag); cufftExec(plan->plan_mid,d_data.data(),d_data.data(),flag);
#else #else
d_tmp = typename FFT_AT::t_FFT_DATA_1d(Kokkos::view_alloc("fft_3d:tmp",Kokkos::WithoutInitializing),d_data.extent(0));
if (flag == -1) if (flag == -1)
f = kiss_fft_functor<DeviceType>(d_data,d_tmp,plan->cfg_mid_forward,length); f = kiss_fft_functor<DeviceType>(d_data,d_tmp,plan->cfg_mid_forward,length);
else else
f = kiss_fft_functor<DeviceType>(d_data,d_tmp,plan->cfg_mid_backward,length); f = kiss_fft_functor<DeviceType>(d_data,d_tmp,plan->cfg_mid_backward,length);
Kokkos::parallel_for(total/length,f); Kokkos::parallel_for(total/length,f);
d_data = d_tmp; d_data = d_tmp;
d_tmp = typename FFT_AT::t_FFT_DATA_1d(Kokkos::view_alloc("fft_3d:tmp",Kokkos::WithoutInitializing),d_in.extent(0));
#endif #endif
// 2nd mid-remap to prepare for 3rd FFTs // 2nd mid-remap to prepare for 3rd FFTs
@ -315,6 +314,7 @@ void FFT3dKokkos<DeviceType>::fft_3d_kokkos(typename FFT_AT::t_FFT_DATA_1d d_in,
#elif defined(FFT_CUFFT) #elif defined(FFT_CUFFT)
cufftExec(plan->plan_slow,d_data.data(),d_data.data(),flag); cufftExec(plan->plan_slow,d_data.data(),d_data.data(),flag);
#else #else
d_tmp = typename FFT_AT::t_FFT_DATA_1d(Kokkos::view_alloc("fft_3d:tmp",Kokkos::WithoutInitializing),d_data.extent(0));
if (flag == -1) if (flag == -1)
f = kiss_fft_functor<DeviceType>(d_data,d_tmp,plan->cfg_slow_forward,length); f = kiss_fft_functor<DeviceType>(d_data,d_tmp,plan->cfg_slow_forward,length);
else else
@ -866,7 +866,8 @@ void FFT3dKokkos<DeviceType>::fft_3d_1d_only_kokkos(typename FFT_AT::t_FFT_DATA_
cufftExec(plan->plan_slow,d_data.data(),d_data.data(),flag); cufftExec(plan->plan_slow,d_data.data(),d_data.data(),flag);
#else #else
kiss_fft_functor<DeviceType> f; kiss_fft_functor<DeviceType> f;
typename FFT_AT::t_FFT_DATA_1d d_tmp = typename FFT_AT::t_FFT_DATA_1d("fft_3d:tmp",d_data.extent(0)); typename FFT_AT::t_FFT_DATA_1d d_tmp =
typename FFT_AT::t_FFT_DATA_1d(Kokkos::view_alloc("fft_3d:tmp",Kokkos::WithoutInitializing),d_data.extent(0));
if (flag == -1) { if (flag == -1) {
f = kiss_fft_functor<DeviceType>(d_data,d_tmp,plan->cfg_fast_forward,length1); f = kiss_fft_functor<DeviceType>(d_data,d_tmp,plan->cfg_fast_forward,length1);
Kokkos::parallel_for(total1/length1,f); Kokkos::parallel_for(total1/length1,f);

View File

@ -2710,11 +2710,15 @@ double FixRigid::compute_scalar()
void *FixRigid::extract(const char *str, int &dim) void *FixRigid::extract(const char *str, int &dim)
{ {
dim = 0;
if (strcmp(str,"body") == 0) { if (strcmp(str,"body") == 0) {
if (!setupflag) return nullptr;
dim = 1; dim = 1;
return body; return body;
} }
if (strcmp(str,"masstotal") == 0) { if (strcmp(str,"masstotal") == 0) {
if (!setupflag) return nullptr;
dim = 1; dim = 1;
return masstotal; return masstotal;
} }

View File

@ -3424,7 +3424,10 @@ int FixRigidSmall::modify_param(int narg, char **arg)
void *FixRigidSmall::extract(const char *str, int &dim) void *FixRigidSmall::extract(const char *str, int &dim)
{ {
dim = 0;
if (strcmp(str,"body") == 0) { if (strcmp(str,"body") == 0) {
if (!setupflag) return nullptr;
dim = 1; dim = 1;
return atom2body; return atom2body;
} }
@ -3438,6 +3441,7 @@ void *FixRigidSmall::extract(const char *str, int &dim)
// used by granular pair styles, indexed by atom2body // used by granular pair styles, indexed by atom2body
if (strcmp(str,"masstotal") == 0) { if (strcmp(str,"masstotal") == 0) {
if (!setupflag) return nullptr;
dim = 1; dim = 1;
if (nmax_mass < nmax_body) { if (nmax_mass < nmax_body) {

View File

@ -176,6 +176,9 @@ void ComputeSpin::compute_vector()
for (i = 0; i < nlocal; i++) { for (i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) { if (mask[i] & groupbit) {
if (atom->sp_flag) { if (atom->sp_flag) {
// compute first moment
mag[0] += sp[i][0]; mag[0] += sp[i][0];
mag[1] += sp[i][1]; mag[1] += sp[i][1];
mag[2] += sp[i][2]; mag[2] += sp[i][2];
@ -211,11 +214,16 @@ void ComputeSpin::compute_vector()
MPI_Allreduce(&tempdenom,&tempdenomtot,1,MPI_DOUBLE,MPI_SUM,world); MPI_Allreduce(&tempdenom,&tempdenomtot,1,MPI_DOUBLE,MPI_SUM,world);
MPI_Allreduce(&countsp,&countsptot,1,MPI_INT,MPI_SUM,world); MPI_Allreduce(&countsp,&countsptot,1,MPI_INT,MPI_SUM,world);
// compute average magnetization
double scale = 1.0/countsptot; double scale = 1.0/countsptot;
magtot[0] *= scale; magtot[0] *= scale;
magtot[1] *= scale; magtot[1] *= scale;
magtot[2] *= scale; magtot[2] *= scale;
magtot[3] = sqrt((magtot[0]*magtot[0])+(magtot[1]*magtot[1])+(magtot[2]*magtot[2])); magtot[3] = sqrt((magtot[0]*magtot[0])+(magtot[1]*magtot[1])+(magtot[2]*magtot[2]));
// compute spin temperature
spintemperature = hbar*tempnumtot; spintemperature = hbar*tempnumtot;
spintemperature /= (2.0*kb*tempdenomtot); spintemperature /= (2.0*kb*tempdenomtot);
@ -225,7 +233,6 @@ void ComputeSpin::compute_vector()
vector[3] = magtot[3]; vector[3] = magtot[3];
vector[4] = magenergytot; vector[4] = magenergytot;
vector[5] = spintemperature; vector[5] = spintemperature;
} }
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------

View File

@ -45,9 +45,10 @@ PairSpinDipoleCut::PairSpinDipoleCut(LAMMPS *lmp) : PairSpin(lmp)
hbar = force->hplanck/MY_2PI; // eV/(rad.THz) hbar = force->hplanck/MY_2PI; // eV/(rad.THz)
mub = 9.274e-4; // in A.Ang^2 mub = 9.274e-4; // in A.Ang^2
mu_0 = 785.15; // in eV/Ang/A^2 // mu_0 = 785.15; // in eV/Ang/A^2
mu_0 = 784.15; // in eV/Ang/A^2
mub2mu0 = mub * mub * mu_0 / (4.0*MY_PI); // in eV.Ang^3 mub2mu0 = mub * mub * mu_0 / (4.0*MY_PI); // in eV.Ang^3
//mub2mu0 = mub * mub * mu_0 / (4.0*MY_PI); // in eV // mub2mu0 = mub * mub * mu_0 / (4.0*MY_PI); // in eV
mub2mu0hbinv = mub2mu0 / hbar; // in rad.THz mub2mu0hbinv = mub2mu0 / hbar; // in rad.THz
} }
@ -232,22 +233,16 @@ void PairSpinDipoleCut::compute(int eflag, int vflag)
local_cut2 = cut_spin_long[itype][jtype]*cut_spin_long[itype][jtype]; local_cut2 = cut_spin_long[itype][jtype]*cut_spin_long[itype][jtype];
// compute dipolar interaction
if (rsq < local_cut2) { if (rsq < local_cut2) {
r2inv = 1.0/rsq; r2inv = 1.0/rsq;
r3inv = r2inv*rinv; r3inv = r2inv*rinv;
compute_dipolar(i,j,eij,fmi,spi,spj,r3inv); compute_dipolar(i,j,eij,fmi,spi,spj,r3inv);
if (lattice_flag) compute_dipolar_mech(i,j,eij,fi,spi,spj,r2inv);
}
// force accumulation if (lattice_flag)
compute_dipolar_mech(i,j,eij,fi,spi,spj,r2inv);
f[i][0] += fi[0];
f[i][1] += fi[1];
f[i][2] += fi[2];
fm[i][0] += fmi[0];
fm[i][1] += fmi[1];
fm[i][2] += fmi[2];
if (eflag) { if (eflag) {
if (rsq <= local_cut2) { if (rsq <= local_cut2) {
@ -257,11 +252,25 @@ void PairSpinDipoleCut::compute(int eflag, int vflag)
} }
} else evdwl = 0.0; } else evdwl = 0.0;
f[i][0] += fi[0];
f[i][1] += fi[1];
f[i][2] += fi[2];
if (newton_pair || j < nlocal) {
f[j][0] -= fi[0];
f[j][1] -= fi[1];
f[j][2] -= fi[2];
}
fm[i][0] += fmi[0];
fm[i][1] += fmi[1];
fm[i][2] += fmi[2];
if (evflag) ev_tally_xyz(i,j,nlocal,newton_pair, if (evflag) ev_tally_xyz(i,j,nlocal,newton_pair,
evdwl,ecoul,fi[0],fi[1],fi[2],rij[0],rij[1],rij[2]); evdwl,ecoul,fi[0],fi[1],fi[2],rij[0],rij[1],rij[2]);
}
}
}
} if (vflag_fdotr) virial_fdotr_compute();
}
} }
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
@ -390,7 +399,7 @@ void PairSpinDipoleCut::compute_dipolar_mech(int /* i */, int /* j */, double ei
sjeij = spj[0]*eij[0] + spj[1]*eij[1] + spj[2]*eij[2]; sjeij = spj[0]*eij[0] + spj[1]*eij[1] + spj[2]*eij[2];
bij = sisj - 5.0*sieij*sjeij; bij = sisj - 5.0*sieij*sjeij;
pre = 3.0*mub2mu0*gigjri4; pre = 0.5*3.0*mub2mu0*gigjri4;
fi[0] -= pre * (eij[0] * bij + (sjeij*spi[0] + sieij*spj[0])); fi[0] -= pre * (eij[0] * bij + (sjeij*spi[0] + sieij*spj[0]));
fi[1] -= pre * (eij[1] * bij + (sjeij*spi[1] + sieij*spj[1])); fi[1] -= pre * (eij[1] * bij + (sjeij*spi[1] + sieij*spj[1]));

View File

@ -49,7 +49,7 @@ PairSpinDipoleLong::PairSpinDipoleLong(LAMMPS *lmp) : PairSpin(lmp)
hbar = force->hplanck/MY_2PI; // eV/(rad.THz) hbar = force->hplanck/MY_2PI; // eV/(rad.THz)
mub = 9.274e-4; // in A.Ang^2 mub = 9.274e-4; // in A.Ang^2
mu_0 = 785.15; // in eV/Ang/A^2 mu_0 = 784.15; // in eV/Ang/A^2
mub2mu0 = mub * mub * mu_0 / (4.0*MY_PI); // in eV.Ang^3 mub2mu0 = mub * mub * mu_0 / (4.0*MY_PI); // in eV.Ang^3
//mub2mu0 = mub * mub * mu_0 / (4.0*MY_PI); // in eV //mub2mu0 = mub * mub * mu_0 / (4.0*MY_PI); // in eV
mub2mu0hbinv = mub2mu0 / hbar; // in rad.THz mub2mu0hbinv = mub2mu0 / hbar; // in rad.THz
@ -281,17 +281,8 @@ void PairSpinDipoleLong::compute(int eflag, int vflag)
bij[3] = (5.0*bij[2] + pre3*expm2) * r2inv; bij[3] = (5.0*bij[2] + pre3*expm2) * r2inv;
compute_long(i,j,eij,bij,fmi,spi,spj); compute_long(i,j,eij,bij,fmi,spi,spj);
if (lattice_flag)
compute_long_mech(i,j,eij,bij,fmi,spi,spj); compute_long_mech(i,j,eij,bij,fmi,spi,spj);
}
// force accumulation
f[i][0] += fi[0];
f[i][1] += fi[1];
f[i][2] += fi[2];
fm[i][0] += fmi[0];
fm[i][1] += fmi[1];
fm[i][2] += fmi[2];
if (eflag) { if (eflag) {
if (rsq <= local_cut2) { if (rsq <= local_cut2) {
@ -301,12 +292,26 @@ void PairSpinDipoleLong::compute(int eflag, int vflag)
} }
} else evdwl = 0.0; } else evdwl = 0.0;
f[i][0] += fi[0];
f[i][1] += fi[1];
f[i][2] += fi[2];
if (newton_pair || j < nlocal) {
f[j][0] -= fi[0];
f[j][1] -= fi[1];
f[j][2] -= fi[2];
}
fm[i][0] += fmi[0];
fm[i][1] += fmi[1];
fm[i][2] += fmi[2];
if (evflag) ev_tally_xyz(i,j,nlocal,newton_pair, if (evflag) ev_tally_xyz(i,j,nlocal,newton_pair,
evdwl,ecoul,fi[0],fi[1],fi[2],rij[0],rij[1],rij[2]); evdwl,ecoul,fi[0],fi[1],fi[2],rij[0],rij[1],rij[2]);
} }
} }
}
if (vflag_fdotr) virial_fdotr_compute();
} }
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
@ -373,7 +378,6 @@ void PairSpinDipoleLong::compute_single_pair(int ii, double fmi[3])
spi[3] = sp[ii][3]; spi[3] = sp[ii][3];
jlist = firstneigh[ii]; jlist = firstneigh[ii];
jnum = numneigh[ii]; jnum = numneigh[ii];
//itype = type[i];
for (jj = 0; jj < jnum; jj++) { for (jj = 0; jj < jnum; jj++) {
j = jlist[jj]; j = jlist[jj];
@ -459,7 +463,7 @@ void PairSpinDipoleLong::compute_long_mech(int /* i */, int /* j */, double eij[
double g1,g2,g1b2_g2b3,gigj,pre; double g1,g2,g1b2_g2b3,gigj,pre;
gigj = spi[3] * spj[3]; gigj = spi[3] * spj[3];
pre = gigj*mub2mu0; pre = 0.5 * gigj*mub2mu0;
sisj = spi[0]*spj[0] + spi[1]*spj[1] + spi[2]*spj[2]; sisj = spi[0]*spj[0] + spi[1]*spj[1] + spi[2]*spj[2];
sieij = spi[0]*eij[0] + spi[1]*eij[1] + spi[2]*eij[2]; sieij = spi[0]*eij[0] + spi[1]*eij[1] + spi[2]*eij[2];
sjeij = spj[0]*eij[0] + spj[1]*eij[1] + spj[2]*eij[2]; sjeij = spj[0]*eij[0] + spj[1]*eij[1] + spj[2]*eij[2];

View File

@ -244,17 +244,9 @@ void PairSpinDmi::compute(int eflag, int vflag)
if (rsq <= local_cut2) { if (rsq <= local_cut2) {
compute_dmi(i,j,eij,fmi,spj); compute_dmi(i,j,eij,fmi,spj);
if (lattice_flag) {
compute_dmi_mech(i,j,rsq,eij,fi,spi,spj);
}
}
f[i][0] += fi[0]; if (lattice_flag)
f[i][1] += fi[1]; compute_dmi_mech(i,j,rsq,eij,fi,spi,spj);
f[i][2] += fi[2];
fm[i][0] += fmi[0];
fm[i][1] += fmi[1];
fm[i][2] += fmi[2];
if (eflag) { if (eflag) {
evdwl -= (spi[0]*fmi[0] + spi[1]*fmi[1] + spi[2]*fmi[2]); evdwl -= (spi[0]*fmi[0] + spi[1]*fmi[1] + spi[2]*fmi[2]);
@ -262,13 +254,25 @@ void PairSpinDmi::compute(int eflag, int vflag)
emag[i] += evdwl; emag[i] += evdwl;
} else evdwl = 0.0; } else evdwl = 0.0;
f[i][0] += fi[0];
f[i][1] += fi[1];
f[i][2] += fi[2];
if (newton_pair || j < nlocal) {
f[j][0] -= fi[0];
f[j][1] -= fi[1];
f[j][2] -= fi[2];
}
fm[i][0] += fmi[0];
fm[i][1] += fmi[1];
fm[i][2] += fmi[2];
if (evflag) ev_tally_xyz(i,j,nlocal,newton_pair, if (evflag) ev_tally_xyz(i,j,nlocal,newton_pair,
evdwl,ecoul,fi[0],fi[1],fi[2],delx,dely,delz); evdwl,ecoul,fi[0],fi[1],fi[2],delx,dely,delz);
} }
} }
}
if (vflag_fdotr) virial_fdotr_compute(); if (vflag_fdotr) virial_fdotr_compute();
} }
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
@ -405,9 +409,9 @@ void PairSpinDmi::compute_dmi_mech(int i, int j, double rsq, double /*eij*/[3],
cdmy = (dmiz*csx - dmix*csz); cdmy = (dmiz*csx - dmix*csz);
cdmz = (dmix*csy - dmiy*csz); cdmz = (dmix*csy - dmiy*csz);
fi[0] += irij*cdmx; fi[0] += 0.5*irij*cdmx;
fi[1] += irij*cdmy; fi[1] += 0.5*irij*cdmy;
fi[2] += irij*cdmz; fi[2] += 0.5*irij*cdmz;
} }
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------

View File

@ -37,6 +37,14 @@ using namespace LAMMPS_NS;
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
PairSpinExchange::PairSpinExchange(LAMMPS *lmp) :
PairSpin(lmp)
{
e_offset = 0;
}
/* ---------------------------------------------------------------------- */
PairSpinExchange::~PairSpinExchange() PairSpinExchange::~PairSpinExchange()
{ {
if (allocated) { if (allocated) {
@ -59,6 +67,8 @@ void PairSpinExchange::settings(int narg, char **arg)
{ {
PairSpin::settings(narg,arg); PairSpin::settings(narg,arg);
if (narg != 1) error->all(FLERR,"Illegal pair_style command");
cut_spin_exchange_global = utils::numeric(FLERR,arg[0],false,lmp); cut_spin_exchange_global = utils::numeric(FLERR,arg[0],false,lmp);
// reset cutoffs that have been explicitly set // reset cutoffs that have been explicitly set
@ -84,9 +94,9 @@ void PairSpinExchange::coeff(int narg, char **arg)
// check if args correct // check if args correct
if (strcmp(arg[2],"exchange") != 0) if (strcmp(arg[2],"exchange") != 0)
error->all(FLERR,"Incorrect args in pair_style command"); error->all(FLERR,"Incorrect args for pair coefficients");
if (narg != 7) if ((narg != 7) && (narg != 9))
error->all(FLERR,"Incorrect args in pair_style command"); error->all(FLERR,"Incorrect args for pair coefficients");
int ilo,ihi,jlo,jhi; int ilo,ihi,jlo,jhi;
utils::bounds(FLERR,arg[0],1,atom->ntypes,ilo,ihi,error); utils::bounds(FLERR,arg[0],1,atom->ntypes,ilo,ihi,error);
@ -94,11 +104,25 @@ void PairSpinExchange::coeff(int narg, char **arg)
// get exchange arguments from input command // get exchange arguments from input command
int iarg = 7;
const double rc = utils::numeric(FLERR,arg[3],false,lmp); const double rc = utils::numeric(FLERR,arg[3],false,lmp);
const double j1 = utils::numeric(FLERR,arg[4],false,lmp); const double j1 = utils::numeric(FLERR,arg[4],false,lmp);
const double j2 = utils::numeric(FLERR,arg[5],false,lmp); const double j2 = utils::numeric(FLERR,arg[5],false,lmp);
const double j3 = utils::numeric(FLERR,arg[6],false,lmp); const double j3 = utils::numeric(FLERR,arg[6],false,lmp);
// read energy offset flag if specified
while (iarg < narg) {
if (strcmp(arg[7],"offset") == 0) {
if (strcmp(arg[8],"yes") == 0) {
e_offset = 1;
} else if (strcmp(arg[8],"no") == 0) {
e_offset = 0;
} else error->all(FLERR,"Incorrect args for pair coefficients");
iarg += 2;
} else error->all(FLERR,"Incorrect args for pair coefficients");
}
int count = 0; int count = 0;
for (int i = ilo; i <= ihi; i++) { for (int i = ilo; i <= ihi; i++) {
for (int j = MAX(jlo,i); j <= jhi; j++) { for (int j = MAX(jlo,i); j <= jhi; j++) {
@ -228,31 +252,34 @@ void PairSpinExchange::compute(int eflag, int vflag)
if (rsq <= local_cut2) { if (rsq <= local_cut2) {
compute_exchange(i,j,rsq,fmi,spj); compute_exchange(i,j,rsq,fmi,spj);
if (lattice_flag) {
if (lattice_flag)
compute_exchange_mech(i,j,rsq,eij,fi,spi,spj); compute_exchange_mech(i,j,rsq,eij,fi,spi,spj);
}
} if (eflag) {
evdwl -= compute_energy(i,j,rsq,spi,spj);
emag[i] += evdwl;
} else evdwl = 0.0;
f[i][0] += fi[0]; f[i][0] += fi[0];
f[i][1] += fi[1]; f[i][1] += fi[1];
f[i][2] += fi[2]; f[i][2] += fi[2];
if (newton_pair || j < nlocal) {
f[j][0] -= fi[0];
f[j][1] -= fi[1];
f[j][2] -= fi[2];
}
fm[i][0] += fmi[0]; fm[i][0] += fmi[0];
fm[i][1] += fmi[1]; fm[i][1] += fmi[1];
fm[i][2] += fmi[2]; fm[i][2] += fmi[2];
if (eflag) {
evdwl -= (spi[0]*fmi[0] + spi[1]*fmi[1] + spi[2]*fmi[2]);
evdwl *= 0.5*hbar;
emag[i] += evdwl;
} else evdwl = 0.0;
if (evflag) ev_tally_xyz(i,j,nlocal,newton_pair, if (evflag) ev_tally_xyz(i,j,nlocal,newton_pair,
evdwl,ecoul,fi[0],fi[1],fi[2],delx,dely,delz); evdwl,ecoul,fi[0],fi[1],fi[2],delx,dely,delz);
} }
} }
}
if (vflag_fdotr) virial_fdotr_compute(); if (vflag_fdotr) virial_fdotr_compute();
} }
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
@ -361,12 +388,14 @@ void PairSpinExchange::compute_exchange(int i, int j, double rsq, double fmi[3],
compute the mechanical force due to the exchange interaction between atom i and atom j compute the mechanical force due to the exchange interaction between atom i and atom j
------------------------------------------------------------------------- */ ------------------------------------------------------------------------- */
void PairSpinExchange::compute_exchange_mech(int i, int j, double rsq, double eij[3], void PairSpinExchange::compute_exchange_mech(int i, int j, double rsq,
double fi[3], double spi[3], double spj[3]) double eij[3], double fi[3], double spi[3], double spj[3])
{ {
int *type = atom->type; int *type = atom->type;
int itype, jtype; int itype, jtype;
double Jex, Jex_mech, ra, rr, iJ3; double Jex, Jex_mech, ra, sdots;
double rr, iJ3;
double fx, fy, fz;
itype = type[i]; itype = type[i];
jtype = type[j]; jtype = type[j];
@ -378,35 +407,58 @@ void PairSpinExchange::compute_exchange_mech(int i, int j, double rsq, double ei
Jex_mech = 1.0-ra-J2[itype][jtype]*ra*(2.0-ra); Jex_mech = 1.0-ra-J2[itype][jtype]*ra*(2.0-ra);
Jex_mech *= 8.0*Jex*rr*exp(-ra); Jex_mech *= 8.0*Jex*rr*exp(-ra);
Jex_mech *= (spi[0]*spj[0]+spi[1]*spj[1]+spi[2]*spj[2]);
fi[0] -= Jex_mech*eij[0]; sdots = (spi[0]*spj[0]+spi[1]*spj[1]+spi[2]*spj[2]);
fi[1] -= Jex_mech*eij[1];
fi[2] -= Jex_mech*eij[2]; // apply or not energy and force offset
fx = fy = fz = 0.0;
if (e_offset == 1) { // set offset
fx = Jex_mech*(sdots-1.0)*eij[0];
fy = Jex_mech*(sdots-1.0)*eij[1];
fz = Jex_mech*(sdots-1.0)*eij[2];
} else if (e_offset == 0) { // no offset ("normal" calculation)
fx = Jex_mech*sdots*eij[0];
fy = Jex_mech*sdots*eij[1];
fz = Jex_mech*sdots*eij[2];
} else error->all(FLERR,"Illegal option in pair exchange/biquadratic command");
fi[0] -= 0.5*fx;
fi[1] -= 0.5*fy;
fi[2] -= 0.5*fz;
} }
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
compute energy of spin pair i and j compute energy of spin pair i and j
------------------------------------------------------------------------- */ ------------------------------------------------------------------------- */
// double PairSpinExchange::compute_energy(int i, int j, double rsq, double spi[3], double spj[3]) double PairSpinExchange::compute_energy(int i, int j, double rsq, double spi[3], double spj[3])
// { {
// int *type = atom->type; int *type = atom->type;
// int itype, jtype; int itype, jtype;
// double Jex, ra; double Jex, ra, sdots;
// double energy = 0.0; double energy = 0.0;
// itype = type[i]; itype = type[i];
// jtype = type[j]; jtype = type[j];
//
// Jex = J1_mech[itype][jtype]; Jex = J1_mech[itype][jtype];
// ra = rsq/J3[itype][jtype]/J3[itype][jtype]; ra = rsq/J3[itype][jtype]/J3[itype][jtype];
// Jex = 4.0*Jex*ra; Jex = 4.0*Jex*ra;
// Jex *= (1.0-J2[itype][jtype]*ra); Jex *= (1.0-J2[itype][jtype]*ra);
// Jex *= exp(-ra); Jex *= exp(-ra);
//
// energy = Jex*(spi[0]*spj[0]+spi[1]*spj[1]+spi[2]*spj[2]); sdots = (spi[0]*spj[0]+spi[1]*spj[1]+spi[2]*spj[2]);
// return energy;
// } // apply or not energy and force offset
if (e_offset == 1) { // set offset
energy = 0.5*Jex*(sdots-1.0);
} else if (e_offset == 0) { // no offset ("normal" calculation)
energy = 0.5*Jex*sdots;
} else error->all(FLERR,"Illegal option in pair exchange/biquadratic command");
return energy;
}
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
allocate all arrays allocate all arrays
@ -495,6 +547,7 @@ void PairSpinExchange::read_restart(FILE *fp)
void PairSpinExchange::write_restart_settings(FILE *fp) void PairSpinExchange::write_restart_settings(FILE *fp)
{ {
fwrite(&cut_spin_exchange_global,sizeof(double),1,fp); fwrite(&cut_spin_exchange_global,sizeof(double),1,fp);
fwrite(&e_offset,sizeof(int),1,fp);
fwrite(&offset_flag,sizeof(int),1,fp); fwrite(&offset_flag,sizeof(int),1,fp);
fwrite(&mix_flag,sizeof(int),1,fp); fwrite(&mix_flag,sizeof(int),1,fp);
} }
@ -507,10 +560,12 @@ void PairSpinExchange::read_restart_settings(FILE *fp)
{ {
if (comm->me == 0) { if (comm->me == 0) {
utils::sfread(FLERR,&cut_spin_exchange_global,sizeof(double),1,fp,nullptr,error); utils::sfread(FLERR,&cut_spin_exchange_global,sizeof(double),1,fp,nullptr,error);
utils::sfread(FLERR,&e_offset,sizeof(int),1,fp,nullptr,error);
utils::sfread(FLERR,&offset_flag,sizeof(int),1,fp,nullptr,error); utils::sfread(FLERR,&offset_flag,sizeof(int),1,fp,nullptr,error);
utils::sfread(FLERR,&mix_flag,sizeof(int),1,fp,nullptr,error); utils::sfread(FLERR,&mix_flag,sizeof(int),1,fp,nullptr,error);
} }
MPI_Bcast(&cut_spin_exchange_global,1,MPI_DOUBLE,0,world); MPI_Bcast(&cut_spin_exchange_global,1,MPI_DOUBLE,0,world);
MPI_Bcast(&e_offset,1,MPI_INT,0,world);
MPI_Bcast(&offset_flag,1,MPI_INT,0,world); MPI_Bcast(&offset_flag,1,MPI_INT,0,world);
MPI_Bcast(&mix_flag,1,MPI_INT,0,world); MPI_Bcast(&mix_flag,1,MPI_INT,0,world);
} }

View File

@ -26,7 +26,7 @@ namespace LAMMPS_NS {
class PairSpinExchange : public PairSpin { class PairSpinExchange : public PairSpin {
public: public:
PairSpinExchange(LAMMPS *lmp) : PairSpin(lmp) {} PairSpinExchange(class LAMMPS *);
virtual ~PairSpinExchange(); virtual ~PairSpinExchange();
void settings(int, char **); void settings(int, char **);
void coeff(int, char **); void coeff(int, char **);
@ -38,8 +38,7 @@ class PairSpinExchange : public PairSpin {
void compute_exchange(int, int, double, double *, double *); void compute_exchange(int, int, double, double *, double *);
void compute_exchange_mech(int, int, double, double *, double *, double *, double *); void compute_exchange_mech(int, int, double, double *, double *, double *, double *);
double compute_energy(int , int , double , double *, double *);
// double compute_energy(int , int , double , double *, double *);
void write_restart(FILE *); void write_restart(FILE *);
void read_restart(FILE *); void read_restart(FILE *);
@ -49,6 +48,7 @@ class PairSpinExchange : public PairSpin {
double cut_spin_exchange_global; // global exchange cutoff distance double cut_spin_exchange_global; // global exchange cutoff distance
protected: protected:
int e_offset; // apply energy offset
double **J1_mag; // exchange coeffs in eV double **J1_mag; // exchange coeffs in eV
double **J1_mech; // mech exchange coeffs in double **J1_mech; // mech exchange coeffs in
double **J2, **J3; // J1 in eV, J2 adim, J3 in Ang double **J2, **J3; // J1 in eV, J2 adim, J3 in Ang

View File

@ -0,0 +1,635 @@
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
/* ------------------------------------------------------------------------
Contributing authors: Julien Tranchida (SNL)
Aidan Thompson (SNL)
Please cite the related publication:
Tranchida, J., Plimpton, S. J., Thibaudeau, P., & Thompson, A. P. (2018).
Massively parallel symplectic algorithm for coupled magnetic spin dynamics
and molecular dynamics. Journal of Computational Physics.
------------------------------------------------------------------------- */
#include "pair_spin_exchange_biquadratic.h"
#include "atom.h"
#include "comm.h"
#include "error.h"
#include "force.h"
#include "memory.h"
#include "neigh_list.h"
#include <cmath>
#include <cstring>
using namespace LAMMPS_NS;
/* ---------------------------------------------------------------------- */
PairSpinExchangeBiquadratic::PairSpinExchangeBiquadratic(LAMMPS *lmp) :
PairSpin(lmp)
{
e_offset = 0;
}
/* ---------------------------------------------------------------------- */
PairSpinExchangeBiquadratic::~PairSpinExchangeBiquadratic()
{
if (allocated) {
memory->destroy(setflag);
memory->destroy(cut_spin_exchange);
memory->destroy(J1_mag);
memory->destroy(J1_mech);
memory->destroy(J2);
memory->destroy(J3);
memory->destroy(K1_mag);
memory->destroy(K1_mech);
memory->destroy(K2);
memory->destroy(K3);
memory->destroy(cutsq); // to be implemented
memory->destroy(emag);
}
}
/* ----------------------------------------------------------------------
global settings
------------------------------------------------------------------------- */
void PairSpinExchangeBiquadratic::settings(int narg, char **arg)
{
PairSpin::settings(narg,arg);
if (narg != 1) error->all(FLERR,"Illegal pair_style command");
cut_spin_exchange_global = utils::numeric(FLERR,arg[0],false,lmp);
// reset cutoffs that have been explicitly set
if (allocated) {
int i,j;
for (i = 1; i <= atom->ntypes; i++)
for (j = i+1; j <= atom->ntypes; j++)
if (setflag[i][j]) {
cut_spin_exchange[i][j] = cut_spin_exchange_global;
}
}
}
/* ----------------------------------------------------------------------
set coeffs for one or more type spin pairs
------------------------------------------------------------------------- */
void PairSpinExchangeBiquadratic::coeff(int narg, char **arg)
{
if (!allocated) allocate();
// check if args correct
if (strcmp(arg[2],"biquadratic") != 0)
error->all(FLERR,"Incorrect args for pair coefficients");
if ((narg != 10) && (narg != 12))
error->all(FLERR,"Incorrect args for pair coefficients");
int ilo,ihi,jlo,jhi;
utils::bounds(FLERR,arg[0],1,atom->ntypes,ilo,ihi,error);
utils::bounds(FLERR,arg[1],1,atom->ntypes,jlo,jhi,error);
// get exchange arguments from input command
int iarg = 10;
const double rc = utils::numeric(FLERR,arg[3],false,lmp);
const double j1 = utils::numeric(FLERR,arg[4],false,lmp);
const double j2 = utils::numeric(FLERR,arg[5],false,lmp);
const double j3 = utils::numeric(FLERR,arg[6],false,lmp);
const double k1 = utils::numeric(FLERR,arg[7],false,lmp);
const double k2 = utils::numeric(FLERR,arg[8],false,lmp);
const double k3 = utils::numeric(FLERR,arg[9],false,lmp);
// read energy offset flag if specified
while (iarg < narg) {
if (strcmp(arg[10],"offset") == 0) {
if (strcmp(arg[11],"yes") == 0) {
e_offset = 1;
} else if (strcmp(arg[11],"no") == 0) {
e_offset = 0;
} else error->all(FLERR,"Incorrect args for pair coefficients");
iarg += 2;
} else error->all(FLERR,"Incorrect args for pair coefficients");
}
int count = 0;
for (int i = ilo; i <= ihi; i++) {
for (int j = MAX(jlo,i); j <= jhi; j++) {
cut_spin_exchange[i][j] = rc;
J1_mag[i][j] = j1/hbar;
J1_mech[i][j] = j1;
J2[i][j] = j2;
J3[i][j] = j3;
K1_mag[i][j] = k1/hbar;
K1_mech[i][j] = k1;
K2[i][j] = k2;
K3[i][j] = k3;
setflag[i][j] = 1;
count++;
}
}
if (count == 0) error->all(FLERR,"Incorrect args in pair_style command");
}
/* ----------------------------------------------------------------------
init for one type pair i,j and corresponding j,i
------------------------------------------------------------------------- */
double PairSpinExchangeBiquadratic::init_one(int i, int j)
{
if (setflag[i][j] == 0) error->all(FLERR,"All pair coeffs are not set");
J1_mag[j][i] = J1_mag[i][j];
J1_mech[j][i] = J1_mech[i][j];
J2[j][i] = J2[i][j];
J3[j][i] = J3[i][j];
K1_mag[j][i] = K1_mag[i][j];
K1_mech[j][i] = K1_mech[i][j];
K2[j][i] = K2[i][j];
K3[j][i] = K3[i][j];
cut_spin_exchange[j][i] = cut_spin_exchange[i][j];
return cut_spin_exchange_global;
}
/* ----------------------------------------------------------------------
extract the larger cutoff
------------------------------------------------------------------------- */
void *PairSpinExchangeBiquadratic::extract(const char *str, int &dim)
{
dim = 0;
if (strcmp(str,"cut") == 0) return (void *) &cut_spin_exchange_global;
return nullptr;
}
/* ---------------------------------------------------------------------- */
void PairSpinExchangeBiquadratic::compute(int eflag, int vflag)
{
int i,j,ii,jj,inum,jnum,itype,jtype;
double evdwl, ecoul;
double xi[3], eij[3];
double delx,dely,delz;
double spi[3], spj[3];
double fi[3], fmi[3];
double local_cut2;
double rsq, inorm;
int *ilist,*jlist,*numneigh,**firstneigh;
evdwl = ecoul = 0.0;
ev_init(eflag,vflag);
double **x = atom->x;
double **f = atom->f;
double **fm = atom->fm;
double **sp = atom->sp;
int *type = atom->type;
int nlocal = atom->nlocal;
int newton_pair = force->newton_pair;
inum = list->inum;
ilist = list->ilist;
numneigh = list->numneigh;
firstneigh = list->firstneigh;
// checking size of emag
if (nlocal_max < nlocal) { // grow emag lists if necessary
nlocal_max = nlocal;
memory->grow(emag,nlocal_max,"pair/spin:emag");
}
// computation of the exchange interaction
// loop over atoms and their neighbors
for (ii = 0; ii < inum; ii++) {
i = ilist[ii];
itype = type[i];
jlist = firstneigh[i];
jnum = numneigh[i];
xi[0] = x[i][0];
xi[1] = x[i][1];
xi[2] = x[i][2];
spi[0] = sp[i][0];
spi[1] = sp[i][1];
spi[2] = sp[i][2];
emag[i] = 0.0;
// loop on neighbors
for (jj = 0; jj < jnum; jj++) {
j = jlist[jj];
j &= NEIGHMASK;
jtype = type[j];
spj[0] = sp[j][0];
spj[1] = sp[j][1];
spj[2] = sp[j][2];
evdwl = 0.0;
fi[0] = fi[1] = fi[2] = 0.0;
fmi[0] = fmi[1] = fmi[2] = 0.0;
delx = xi[0] - x[j][0];
dely = xi[1] - x[j][1];
delz = xi[2] - x[j][2];
rsq = delx*delx + dely*dely + delz*delz;
inorm = 1.0/sqrt(rsq);
eij[0] = -inorm*delx;
eij[1] = -inorm*dely;
eij[2] = -inorm*delz;
local_cut2 = cut_spin_exchange[itype][jtype]*cut_spin_exchange[itype][jtype];
// compute exchange interaction
if (rsq <= local_cut2) {
compute_exchange(i,j,rsq,fmi,spi,spj);
if (lattice_flag)
compute_exchange_mech(i,j,rsq,eij,fi,spi,spj);
if (eflag) {
evdwl -= compute_energy(i,j,rsq,spi,spj);
emag[i] += evdwl;
} else evdwl = 0.0;
f[i][0] += fi[0];
f[i][1] += fi[1];
f[i][2] += fi[2];
if (newton_pair || j < nlocal) {
f[j][0] -= fi[0];
f[j][1] -= fi[1];
f[j][2] -= fi[2];
}
fm[i][0] += fmi[0];
fm[i][1] += fmi[1];
fm[i][2] += fmi[2];
if (evflag) ev_tally_xyz(i,j,nlocal,newton_pair,
evdwl,ecoul,fi[0],fi[1],fi[2],delx,dely,delz);
}
}
}
if (vflag_fdotr) virial_fdotr_compute();
}
/* ----------------------------------------------------------------------
update the pair interactions fmi acting on the spin ii
------------------------------------------------------------------------- */
void PairSpinExchangeBiquadratic::compute_single_pair(int ii, double fmi[3])
{
int *type = atom->type;
double **x = atom->x;
double **sp = atom->sp;
double local_cut2;
double xi[3];
double delx,dely,delz;
double spi[3],spj[3];
int j,jnum,itype,jtype,ntypes;
int k,locflag;
int *jlist,*numneigh,**firstneigh;
double rsq;
numneigh = list->numneigh;
firstneigh = list->firstneigh;
// check if interaction applies to type of ii
itype = type[ii];
ntypes = atom->ntypes;
locflag = 0;
k = 1;
while (k <= ntypes) {
if (k <= itype) {
if (setflag[k][itype] == 1) {
locflag =1;
break;
}
k++;
} else if (k > itype) {
if (setflag[itype][k] == 1) {
locflag =1;
break;
}
k++;
} else error->all(FLERR,"Wrong type number");
}
// if interaction applies to type ii,
// locflag = 1 and compute pair interaction
if (locflag == 1) {
xi[0] = x[ii][0];
xi[1] = x[ii][1];
xi[2] = x[ii][2];
spi[0] = sp[ii][0];
spi[1] = sp[ii][1];
spi[2] = sp[ii][2];
jlist = firstneigh[ii];
jnum = numneigh[ii];
for (int jj = 0; jj < jnum; jj++) {
j = jlist[jj];
j &= NEIGHMASK;
jtype = type[j];
local_cut2 = cut_spin_exchange[itype][jtype]*cut_spin_exchange[itype][jtype];
spj[0] = sp[j][0];
spj[1] = sp[j][1];
spj[2] = sp[j][2];
delx = xi[0] - x[j][0];
dely = xi[1] - x[j][1];
delz = xi[2] - x[j][2];
rsq = delx*delx + dely*dely + delz*delz;
if (rsq <= local_cut2) {
compute_exchange(ii,j,rsq,fmi,spi,spj);
}
}
}
}
/* ----------------------------------------------------------------------
compute exchange interaction between spins i and j
------------------------------------------------------------------------- */
void PairSpinExchangeBiquadratic::compute_exchange(int i, int j, double rsq,
double fmi[3], double spi[3], double spj[3])
{
int *type = atom->type;
int itype,jtype;
double Jex,Kex,r2j,r2k,sdots;
itype = type[i];
jtype = type[j];
r2j = rsq/J3[itype][jtype]/J3[itype][jtype];
r2k = rsq/J3[itype][jtype]/J3[itype][jtype];
Jex = 4.0*J1_mag[itype][jtype]*r2j;
Jex *= (1.0-J2[itype][jtype]*r2j);
Jex *= exp(-r2j);
Kex = 4.0*K1_mag[itype][jtype]*r2k;
Kex *= (1.0-K2[itype][jtype]*r2k);
Kex *= exp(-r2k);
sdots = (spi[0]*spj[0]+spi[1]*spj[1]+spi[2]*spj[2]);
fmi[0] += (Jex*spj[0] + 2.0*Kex*spj[0]*sdots);
fmi[1] += (Jex*spj[1] + 2.0*Kex*spj[1]*sdots);
fmi[2] += (Jex*spj[2] + 2.0*Kex*spj[2]*sdots);
}
/* ----------------------------------------------------------------------
compute the mechanical force due to the exchange interaction between atom i and atom j
------------------------------------------------------------------------- */
void PairSpinExchangeBiquadratic::compute_exchange_mech(int i, int j,
double rsq, double eij[3], double fi[3], double spi[3], double spj[3])
{
int *type = atom->type;
int itype,jtype;
double Jex,Jex_mech,Kex,Kex_mech,sdots;
double rja,rka,rjr,rkr,iJ3,iK3;
double fx, fy, fz;
itype = type[i];
jtype = type[j];
Jex = J1_mech[itype][jtype];
iJ3 = 1.0/(J3[itype][jtype]*J3[itype][jtype]);
Kex = K1_mech[itype][jtype];
iK3 = 1.0/(K3[itype][jtype]*K3[itype][jtype]);
rja = rsq*iJ3;
rjr = sqrt(rsq)*iJ3;
rka = rsq*iK3;
rkr = sqrt(rsq)*iK3;
Jex_mech = 1.0-rja-J2[itype][jtype]*rja*(2.0-rja);
Jex_mech *= 8.0*Jex*rjr*exp(-rja);
Kex_mech = 1.0-rka-K2[itype][jtype]*rka*(2.0-rka);
Kex_mech *= 8.0*Kex*rkr*exp(-rka);
sdots = (spi[0]*spj[0]+spi[1]*spj[1]+spi[2]*spj[2]);
// apply or not energy and force offset
fx = fy = fz = 0.0;
if (e_offset == 1) { // set offset
fx = (Jex_mech*(sdots-1.0) + Kex_mech*(sdots*sdots-1.0))*eij[0];
fy = (Jex_mech*(sdots-1.0) + Kex_mech*(sdots*sdots-1.0))*eij[1];
fz = (Jex_mech*(sdots-1.0) + Kex_mech*(sdots*sdots-1.0))*eij[2];
} else if (e_offset == 0) { // no offset ("normal" calculation)
fx = (Jex_mech*sdots + Kex_mech*sdots*sdots)*eij[0];
fy = (Jex_mech*sdots + Kex_mech*sdots*sdots)*eij[1];
fz = (Jex_mech*sdots + Kex_mech*sdots*sdots)*eij[2];
} else error->all(FLERR,"Illegal option in pair exchange/biquadratic command");
fi[0] -= 0.5*fx;
fi[1] -= 0.5*fy;
fi[2] -= 0.5*fz;
// fi[0] -= fx;
// fi[1] -= fy;
// fi[2] -= fz;
}
/* ----------------------------------------------------------------------
compute energy of spin pair i and j
------------------------------------------------------------------------- */
double PairSpinExchangeBiquadratic::compute_energy(int i, int j, double rsq,
double spi[3], double spj[3])
{
int *type = atom->type;
int itype,jtype;
double Jex,Kex,ra,sdots;
double rj,rk,r2j,r2k,ir3j,ir3k;
double energy = 0.0;
itype = type[i];
jtype = type[j];
ra = sqrt(rsq);
rj = ra/J3[itype][jtype];
r2j = rsq/J3[itype][jtype]/J3[itype][jtype];
ir3j = 1.0/(rj*rj*rj);
rk = ra/K3[itype][jtype];
r2k = rsq/K3[itype][jtype]/K3[itype][jtype];
ir3k = 1.0/(rk*rk*rk);
Jex = 4.0*J1_mech[itype][jtype]*r2j;
Jex *= (1.0-J2[itype][jtype]*r2j);
Jex *= exp(-r2j);
Kex = 4.0*K1_mech[itype][jtype]*r2k;
Kex *= (1.0-K2[itype][jtype]*r2k);
Kex *= exp(-r2k);
sdots = (spi[0]*spj[0]+spi[1]*spj[1]+spi[2]*spj[2]);
// apply or not energy and force offset
if (e_offset == 1) { // set offset
energy = 0.5*(Jex*(sdots-1.0) + Kex*(sdots*sdots-1.0));
} else if (e_offset == 0) { // no offset ("normal" calculation)
energy = 0.5*(Jex*sdots + Kex*sdots*sdots);
} else error->all(FLERR,"Illegal option in pair exchange/biquadratic command");
return energy;
}
/* ----------------------------------------------------------------------
allocate all arrays
------------------------------------------------------------------------- */
void PairSpinExchangeBiquadratic::allocate()
{
allocated = 1;
int n = atom->ntypes;
memory->create(setflag,n+1,n+1,"pair:setflag");
for (int i = 1; i <= n; i++)
for (int j = i; j <= n; j++)
setflag[i][j] = 0;
memory->create(cut_spin_exchange,n+1,n+1,"pair/spin/exchange:cut_spin_exchange");
memory->create(J1_mag,n+1,n+1,"pair/spin/exchange:J1_mag");
memory->create(J1_mech,n+1,n+1,"pair/spin/exchange:J1_mech");
memory->create(J2,n+1,n+1,"pair/spin/exchange:J2");
memory->create(J3,n+1,n+1,"pair/spin/exchange:J3");
memory->create(K1_mag,n+1,n+1,"pair/spin/exchange:J1_mag");
memory->create(K1_mech,n+1,n+1,"pair/spin/exchange:J1_mech");
memory->create(K2,n+1,n+1,"pair/spin/exchange:J2");
memory->create(K3,n+1,n+1,"pair/spin/exchange:J3");
memory->create(cutsq,n+1,n+1,"pair:cutsq");
}
/* ----------------------------------------------------------------------
proc 0 writes to restart file
------------------------------------------------------------------------- */
void PairSpinExchangeBiquadratic::write_restart(FILE *fp)
{
write_restart_settings(fp);
int i,j;
for (i = 1; i <= atom->ntypes; i++) {
for (j = i; j <= atom->ntypes; j++) {
fwrite(&setflag[i][j],sizeof(int),1,fp);
if (setflag[i][j]) {
fwrite(&J1_mag[i][j],sizeof(double),1,fp);
fwrite(&J1_mech[i][j],sizeof(double),1,fp);
fwrite(&J2[i][j],sizeof(double),1,fp);
fwrite(&J3[i][j],sizeof(double),1,fp);
fwrite(&K1_mag[i][j],sizeof(double),1,fp);
fwrite(&K1_mech[i][j],sizeof(double),1,fp);
fwrite(&K2[i][j],sizeof(double),1,fp);
fwrite(&K3[i][j],sizeof(double),1,fp);
fwrite(&cut_spin_exchange[i][j],sizeof(double),1,fp);
}
}
}
}
/* ----------------------------------------------------------------------
proc 0 reads from restart file, bcasts
------------------------------------------------------------------------- */
void PairSpinExchangeBiquadratic::read_restart(FILE *fp)
{
read_restart_settings(fp);
allocate();
int i,j;
int me = comm->me;
for (i = 1; i <= atom->ntypes; i++) {
for (j = i; j <= atom->ntypes; j++) {
if (me == 0) utils::sfread(FLERR,&setflag[i][j],sizeof(int),1,fp,nullptr,error);
MPI_Bcast(&setflag[i][j],1,MPI_INT,0,world);
if (setflag[i][j]) {
if (me == 0) {
utils::sfread(FLERR,&J1_mag[i][j],sizeof(double),1,fp,nullptr,error);
utils::sfread(FLERR,&J1_mech[i][j],sizeof(double),1,fp,nullptr,error);
utils::sfread(FLERR,&J2[i][j],sizeof(double),1,fp,nullptr,error);
utils::sfread(FLERR,&J3[i][j],sizeof(double),1,fp,nullptr,error);
utils::sfread(FLERR,&K1_mag[i][j],sizeof(double),1,fp,nullptr,error);
utils::sfread(FLERR,&K1_mech[i][j],sizeof(double),1,fp,nullptr,error);
utils::sfread(FLERR,&K2[i][j],sizeof(double),1,fp,nullptr,error);
utils::sfread(FLERR,&K3[i][j],sizeof(double),1,fp,nullptr,error);
utils::sfread(FLERR,&cut_spin_exchange[i][j],sizeof(double),1,fp,nullptr,error);
}
MPI_Bcast(&J1_mag[i][j],1,MPI_DOUBLE,0,world);
MPI_Bcast(&J1_mech[i][j],1,MPI_DOUBLE,0,world);
MPI_Bcast(&J2[i][j],1,MPI_DOUBLE,0,world);
MPI_Bcast(&J3[i][j],1,MPI_DOUBLE,0,world);
MPI_Bcast(&K1_mag[i][j],1,MPI_DOUBLE,0,world);
MPI_Bcast(&K1_mech[i][j],1,MPI_DOUBLE,0,world);
MPI_Bcast(&K2[i][j],1,MPI_DOUBLE,0,world);
MPI_Bcast(&K3[i][j],1,MPI_DOUBLE,0,world);
MPI_Bcast(&cut_spin_exchange[i][j],1,MPI_DOUBLE,0,world);
}
}
}
}
/* ----------------------------------------------------------------------
proc 0 writes to restart file
------------------------------------------------------------------------- */
void PairSpinExchangeBiquadratic::write_restart_settings(FILE *fp)
{
fwrite(&cut_spin_exchange_global,sizeof(double),1,fp);
fwrite(&e_offset,sizeof(int),1,fp);
fwrite(&offset_flag,sizeof(int),1,fp);
fwrite(&mix_flag,sizeof(int),1,fp);
}
/* ----------------------------------------------------------------------
proc 0 reads from restart file, bcasts
------------------------------------------------------------------------- */
void PairSpinExchangeBiquadratic::read_restart_settings(FILE *fp)
{
if (comm->me == 0) {
utils::sfread(FLERR,&cut_spin_exchange_global,sizeof(double),1,fp,nullptr,error);
utils::sfread(FLERR,&e_offset,sizeof(int),1,fp,nullptr,error);
utils::sfread(FLERR,&offset_flag,sizeof(int),1,fp,nullptr,error);
utils::sfread(FLERR,&mix_flag,sizeof(int),1,fp,nullptr,error);
}
MPI_Bcast(&cut_spin_exchange_global,1,MPI_DOUBLE,0,world);
MPI_Bcast(&e_offset,1,MPI_INT,0,world);
MPI_Bcast(&offset_flag,1,MPI_INT,0,world);
MPI_Bcast(&mix_flag,1,MPI_INT,0,world);
}

View File

@ -0,0 +1,87 @@
/* -*- c++ -*- ----------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#ifdef PAIR_CLASS
PairStyle(spin/exchange/biquadratic,PairSpinExchangeBiquadratic)
#else
#ifndef LMP_PAIR_SPIN_EXCHANGE_BIQUADRATIC_H
#define LMP_PAIR_SPIN_EXCHANGE_BIQUADRATIC_H
#include "pair_spin.h"
namespace LAMMPS_NS {
class PairSpinExchangeBiquadratic : public PairSpin {
public:
PairSpinExchangeBiquadratic(class LAMMPS *);
virtual ~PairSpinExchangeBiquadratic();
void settings(int, char **);
void coeff(int, char **);
double init_one(int, int);
void *extract(const char *, int &);
void compute(int, int);
void compute_single_pair(int, double *);
void compute_exchange(int, int, double, double *, double *, double *);
void compute_exchange_mech(int, int, double, double *, double *, double *, double *);
double compute_energy(int , int , double , double *, double *);
void write_restart(FILE *);
void read_restart(FILE *);
void write_restart_settings(FILE *);
void read_restart_settings(FILE *);
double cut_spin_exchange_global; // global exchange cutoff distance
protected:
int e_offset; // apply energy offset
double **J1_mag; // H exchange coeffs in eV
double **J1_mech; // mech exchange coeffs in
double **J2, **J3; // J1 in eV, J2 in Ang-1, J3 in Ang
double **K1_mag; // Bi exchange coeffs in eV
double **K1_mech; // mech exchange coeffs in
double **K2, **K3; // K1 in eV, K2 Ang-1, K3 in Ang
double **cut_spin_exchange; // cutoff distance exchange
void allocate();
};
}
#endif
#endif
/* ERROR/WARNING messages:
E: Incorrect args in pair_spin command
Self-explanatory.
E: Spin simulations require metal unit style
Self-explanatory.
E: Incorrect args for pair coefficients
Self-explanatory. Check the input script or data file.
E: Pair spin requires atom attribute spin
The atom style defined does not have these attributes.
*/

View File

@ -237,17 +237,9 @@ void PairSpinMagelec::compute(int eflag, int vflag)
if (rsq <= local_cut2) { if (rsq <= local_cut2) {
compute_magelec(i,j,eij,fmi,spj); compute_magelec(i,j,eij,fmi,spj);
if (lattice_flag) {
compute_magelec_mech(i,j,fi,spi,spj);
}
}
f[i][0] += fi[0]; if (lattice_flag)
f[i][1] += fi[1]; compute_magelec_mech(i,j,fi,spi,spj);
f[i][2] += fi[2];
fm[i][0] += fmi[0];
fm[i][1] += fmi[1];
fm[i][2] += fmi[2];
if (eflag) { if (eflag) {
evdwl -= (spi[0]*fmi[0] + spi[1]*fmi[1] + spi[2]*fmi[2]); evdwl -= (spi[0]*fmi[0] + spi[1]*fmi[1] + spi[2]*fmi[2]);
@ -255,13 +247,25 @@ void PairSpinMagelec::compute(int eflag, int vflag)
emag[i] += evdwl; emag[i] += evdwl;
} else evdwl = 0.0; } else evdwl = 0.0;
f[i][0] += fi[0];
f[i][1] += fi[1];
f[i][2] += fi[2];
if (newton_pair || j < nlocal) {
f[j][0] -= fi[0];
f[j][1] -= fi[1];
f[j][2] -= fi[2];
}
fm[i][0] += fmi[0];
fm[i][1] += fmi[1];
fm[i][2] += fmi[2];
if (evflag) ev_tally_xyz(i,j,nlocal,newton_pair, if (evflag) ev_tally_xyz(i,j,nlocal,newton_pair,
evdwl,ecoul,fi[0],fi[1],fi[2],delx,dely,delz); evdwl,ecoul,fi[0],fi[1],fi[2],delx,dely,delz);
} }
} }
}
if (vflag_fdotr) virial_fdotr_compute(); if (vflag_fdotr) virial_fdotr_compute();
} }
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
@ -400,9 +404,9 @@ void PairSpinMagelec::compute_magelec_mech(int i, int j, double fi[3], double sp
meiy *= ME_mech[itype][jtype]; meiy *= ME_mech[itype][jtype];
meiz *= ME_mech[itype][jtype]; meiz *= ME_mech[itype][jtype];
fi[0] += (meiy*vz - meiz*vy); fi[0] += 0.5*(meiy*vz - meiz*vy);
fi[1] += (meiz*vx - meix*vz); fi[1] += 0.5*(meiz*vx - meix*vz);
fi[2] += (meix*vy - meiy*vx); fi[2] += 0.5*(meix*vy - meiy*vx);
} }

View File

@ -246,21 +246,23 @@ void PairSpinNeel::compute(int eflag, int vflag)
if (rsq <= local_cut2) { if (rsq <= local_cut2) {
compute_neel(i,j,rsq,eij,fmi,spi,spj); compute_neel(i,j,rsq,eij,fmi,spi,spj);
if (lattice_flag) { if (lattice_flag)
compute_neel_mech(i,j,rsq,eij,fi,spi,spj); compute_neel_mech(i,j,rsq,eij,fi,spi,spj);
}
}
f[i][0] += fi[0]; f[i][0] += fi[0];
f[i][1] += fi[1]; f[i][1] += fi[1];
f[i][2] += fi[2]; f[i][2] += fi[2];
if (newton_pair || j < nlocal) {
f[j][0] -= fi[0];
f[j][1] -= fi[1];
f[j][2] -= fi[2];
}
fm[i][0] += fmi[0]; fm[i][0] += fmi[0];
fm[i][1] += fmi[1]; fm[i][1] += fmi[1];
fm[i][2] += fmi[2]; fm[i][2] += fmi[2];
if (eflag) { if (eflag) {
evdwl = compute_neel_energy(i,j,rsq,eij,spi,spj); evdwl -= compute_neel_energy(i,j,rsq,eij,spi,spj);
evdwl *= 0.5*hbar;
emag[i] += evdwl; emag[i] += evdwl;
} else evdwl = 0.0; } else evdwl = 0.0;
@ -268,9 +270,9 @@ void PairSpinNeel::compute(int eflag, int vflag)
evdwl,ecoul,fi[0],fi[1],fi[2],rij[0],rij[1],rij[2]); evdwl,ecoul,fi[0],fi[1],fi[2],rij[0],rij[1],rij[2]);
} }
} }
}
if (vflag_fdotr) virial_fdotr_compute(); if (vflag_fdotr) virial_fdotr_compute();
} }
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
@ -563,9 +565,9 @@ void PairSpinNeel::compute_neel_mech(int i, int j, double rsq, double eij[3], do
// adding three contributions // adding three contributions
fi[0] = pdx + pq1x + pq2x; fi[0] = 0.5*(pdx + pq1x + pq2x);
fi[1] = pdy + pq1y + pq2y; fi[1] = 0.5*(pdy + pq1y + pq2y);
fi[2] = pdz + pq1z + pq2z; fi[2] = 0.5*(pdz + pq1z + pq2z);
} }
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
@ -585,12 +587,12 @@ double PairSpinNeel::compute_neel_energy(int i, int j, double rsq, double eij[3]
// compute Neel's functions // compute Neel's functions
ra = rsq/g3[itype][jtype]/g3[itype][jtype]; ra = rsq/g3[itype][jtype]/g3[itype][jtype];
gr = 4.0*g1[itype][jtype]*ra; gr = 4.0*g1_mech[itype][jtype]*ra;
gr *= (1.0-g2[itype][jtype]*ra); gr *= (1.0-g2[itype][jtype]*ra);
gr *= exp(-ra); gr *= exp(-ra);
ra = rsq/q3[itype][jtype]/q3[itype][jtype]; ra = rsq/q3[itype][jtype]/q3[itype][jtype];
qr = 4.0*q1[itype][jtype]*ra; qr = 4.0*q1_mech[itype][jtype]*ra;
qr *= (1.0-q2[itype][jtype]*ra); qr *= (1.0-q2[itype][jtype]*ra);
qr *= exp(-ra); qr *= exp(-ra);
@ -609,7 +611,7 @@ double PairSpinNeel::compute_neel_energy(int i, int j, double rsq, double eij[3]
eij_sj_3 = eij_sj*eij_sj_2; eij_sj_3 = eij_sj*eij_sj_2;
epq2 = q2r*(eij_si*eij_sj_3+eij_sj*eij_si_3); epq2 = q2r*(eij_si*eij_sj_3+eij_sj*eij_si_3);
return (epd+epq1+epq2); return 0.5*(epd+epq1+epq2);
} }
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------

View File

@ -9,7 +9,6 @@
the GNU General Public License. the GNU General Public License.
See the README file in the top-level LAMMPS directory. See the README file in the top-level LAMMPS directory.
Contributing author: Maxim Shugaev (UVA), mvs9t@virginia.edu Contributing author: Maxim Shugaev (UVA), mvs9t@virginia.edu
------------------------------------------------------------------------- */ ------------------------------------------------------------------------- */
@ -29,6 +28,7 @@
#include <cstring> #include <cstring>
#include <vector> #include <vector>
#include <cmath> #include <cmath>
#include <array>
#include <fstream> #include <fstream>
#include <sstream> #include <sstream>
@ -36,28 +36,17 @@
using namespace LAMMPS_NS; using namespace LAMMPS_NS;
//since LAMMPS is compiled with C++ 2003, define a substitution for std::array
template<typename T, int N>
class array2003{
public:
T& operator[] (int idx){ return data[idx];};
const T& operator[] (int idx) const{ return data[idx];};
private:
T data[N];
};
class MESONTList { class MESONTList {
public: public:
MESONTList(const Atom* atom, const NeighList* nblist, double rc2); MESONTList(const Atom* atom, const NeighList* nblist);
~MESONTList() {}; ~MESONTList() {};
//list of segments //list of segments
const std::vector<array2003<int,2> >& get_segments() const; const std::vector<std::array<int,2>>& get_segments() const;
//list of triplets //list of triplets
const std::vector<array2003<int,3> >& get_triplets() const; const std::vector<std::array<int,3>>& get_triplets() const;
//list of neighbor chains [start,end] for segments //list of neighbor chains [start,end] for segments
//(use idx() to get real indexes) //(use idx() to get real indexes)
const std::vector<std::vector<array2003<int,2> > >& get_nbs() const; const std::vector<std::vector<std::array<int,2>>>& get_nbs() const;
//convert idx from sorted representation to real idx //convert idx from sorted representation to real idx
int get_idx(int idx) const; int get_idx(int idx) const;
//return list of indexes for conversion from sorted representation //return list of indexes for conversion from sorted representation
@ -69,22 +58,22 @@ public:
//check if the node is the end of the tube //check if the node is the end of the tube
bool is_end(int idx) const; bool is_end(int idx) const;
array2003<int, 2> get_segment(int idx) const; std::array<int,2> get_segment(int idx) const;
array2003<int, 3> get_triplet(int idx) const; std::array<int,3> get_triplet(int idx) const;
static const int cnt_end = -1; static const int cnt_end = -1;
static const int domain_end = -2; static const int domain_end = -2;
static const int not_cnt = -3; static const int not_cnt = -3;
private: private:
std::vector<array2003<int, 2> > chain_list, segments; std::vector<std::array<int,2>> chain_list, segments;
std::vector<array2003<int, 3> > triplets; std::vector<std::array<int,3>> triplets;
std::vector<std::vector<array2003<int, 2> > > nb_chains; std::vector<std::vector<std::array<int,2>>> nb_chains;
std::vector<int> index_list, index_list_b; std::vector<int> index_list, index_list_b;
}; };
//============================================================================= //=============================================================================
inline const std::vector<std::vector<array2003<int, 2> > > & inline const std::vector<std::vector<std::array<int,2>>> &
MESONTList::get_nbs() const { MESONTList::get_nbs() const {
return nb_chains; return nb_chains;
} }
@ -106,25 +95,25 @@ inline const std::vector<int>& MESONTList::get_idxb_list() const {
return index_list_b; return index_list_b;
}; };
inline const std::vector<array2003<int, 2> > & MESONTList::get_segments() inline const std::vector<std::array<int,2>> & MESONTList::get_segments()
const { const {
return segments; return segments;
} }
inline const std::vector<array2003<int, 3> > & MESONTList::get_triplets() inline const std::vector<std::array<int,3>> & MESONTList::get_triplets()
const { const {
return triplets; return triplets;
} }
inline array2003<int, 2> MESONTList::get_segment(int idx) const { inline std::array<int,2> MESONTList::get_segment(int idx) const {
array2003<int, 2> result; std::array<int,2> result;
result[0] = chain_list[idx][0]; result[0] = chain_list[idx][0];
result[1] = idx; result[1] = idx;
return result; return result;
} }
inline array2003<int, 3> MESONTList::get_triplet(int idx) const { inline std::array<int,3> MESONTList::get_triplet(int idx) const {
array2003<int, 3> result; std::array<int,3> result;
result[0] = chain_list[idx][0]; result[0] = chain_list[idx][0];
result[1] = idx; result[1] = idx;
result[2] = chain_list[idx][1]; result[2] = chain_list[idx][1];
@ -165,12 +154,13 @@ void vector_union(std::vector<T>& v1, std::vector<T>& v2,
} }
} }
MESONTList::MESONTList(const Atom* atom, const NeighList* nblist, double /* rc2 */){ MESONTList::MESONTList(const Atom* atom, const NeighList* nblist) {
if (atom == nullptr || nblist == nullptr) return; if (atom == nullptr || nblist == nullptr) return;
//number of local atoms at the node //number of local atoms at the node
int nlocal = atom->nlocal; int nlocal = atom->nlocal;
//total number of atoms in the node and ghost shell //total number of atoms in the node and ghost shell treated as NTs
int nall = nblist->inum + nblist->gnum; int nall = nblist->inum + nblist->gnum;
//total number of atoms in the node and ghost shell
int ntot = atom->nlocal + atom->nghost; int ntot = atom->nlocal + atom->nghost;
tagint* const g_id = atom->tag; tagint* const g_id = atom->tag;
tagint** const bonds = atom->bond_nt; tagint** const bonds = atom->bond_nt;
@ -178,9 +168,7 @@ MESONTList::MESONTList(const Atom* atom, const NeighList* nblist, double /* rc2
int* ilist = nblist->ilist; int* ilist = nblist->ilist;
//convert bonds to local id representation //convert bonds to local id representation
array2003<int, 2> tmp_arr; chain_list.resize(ntot, {not_cnt,not_cnt});
tmp_arr[0] = not_cnt; tmp_arr[1] = not_cnt;
chain_list.resize(ntot, tmp_arr);
for (int ii = 0; ii < nall; ii++) { for (int ii = 0; ii < nall; ii++) {
int i = ilist[ii]; int i = ilist[ii];
chain_list[i][0] = domain_end; chain_list[i][0] = domain_end;
@ -193,7 +181,7 @@ MESONTList::MESONTList(const Atom* atom, const NeighList* nblist, double /* rc2
if (bonds[i][m] == cnt_end) chain_list[i][m] = cnt_end; if (bonds[i][m] == cnt_end) chain_list[i][m] = cnt_end;
for (int j = 0; j < nnb; j++) { for (int j = 0; j < nnb; j++) {
int nb = nblist->firstneigh[i][j]; int nb = nblist->firstneigh[i][j];
if (bonds[i][0] == g_id[nb]){ if (bonds[i][0] == g_id[nb]) {
chain_list[i][0] = nb; chain_list[i][0] = nb;
chain_list[nb][1] = i; chain_list[nb][1] = i;
break; break;
@ -223,22 +211,16 @@ MESONTList::MESONTList(const Atom* atom, const NeighList* nblist, double /* rc2
for (int i = 0; i < nlocal; i++) { for (int i = 0; i < nlocal; i++) {
if (chain_list[i][0] == not_cnt) continue; if (chain_list[i][0] == not_cnt) continue;
if (chain_list[i][0] != cnt_end && chain_list[i][0] != domain_end && if (chain_list[i][0] != cnt_end && chain_list[i][0] != domain_end &&
g_id[i] < g_id[chain_list[i][0]]){ g_id[i] < g_id[chain_list[i][0]])
array2003<int, 2> tmp_c; segments.push_back({i,chain_list[i][0]});
tmp_c[0] = i; tmp_c[1] = chain_list[i][0];
segments.push_back(tmp_c);
}
if (chain_list[i][1] != cnt_end && chain_list[i][1] != domain_end && if (chain_list[i][1] != cnt_end && chain_list[i][1] != domain_end &&
g_id[i] < g_id[chain_list[i][1]]){ g_id[i] < g_id[chain_list[i][1]])
array2003<int, 2> tmp_c; segments.push_back({i,chain_list[i][1]});
tmp_c[0] = i; tmp_c[1] = chain_list[i][1];
segments.push_back(tmp_c);
}
} }
int nbonds = segments.size(); int nbonds = segments.size();
//triplets //triplets
for (int i = 0; i < nlocal; i++){ for (int i = 0; i < nlocal; i++) {
if (chain_list[i][0] == not_cnt) continue; if (chain_list[i][0] == not_cnt) continue;
if (chain_list[i][0] != cnt_end && chain_list[i][0] != domain_end && if (chain_list[i][0] != cnt_end && chain_list[i][0] != domain_end &&
chain_list[i][1] != cnt_end && chain_list[i][1] != domain_end) chain_list[i][1] != cnt_end && chain_list[i][1] != domain_end)
@ -274,7 +256,7 @@ MESONTList::MESONTList(const Atom* atom, const NeighList* nblist, double /* rc2
int idx_next = chain_list[index_list[nb_list[j]]][1]; int idx_next = chain_list[index_list[nb_list[j]]][1];
if ((j == nnb - 1) || (nb_list[j] + 1 != nb_list[j+1]) || if ((j == nnb - 1) || (nb_list[j] + 1 != nb_list[j+1]) ||
(idx_next == cnt_end) || (idx_next == domain_end)) { (idx_next == cnt_end) || (idx_next == domain_end)) {
array2003<int, 2> chain; std::array<int,2> chain;
chain[0] = idx_s; chain[0] = idx_s;
chain[1] = nb_list[j]; chain[1] = nb_list[j];
//make sure that segments having at least one node //make sure that segments having at least one node
@ -285,7 +267,7 @@ MESONTList::MESONTList(const Atom* atom, const NeighList* nblist, double /* rc2
chain_list[idx0][0] != domain_end) chain[0] -= 1; chain_list[idx0][0] != domain_end) chain[0] -= 1;
if (chain_list[idx1][1] != cnt_end && if (chain_list[idx1][1] != cnt_end &&
chain_list[idx1][1] != domain_end) chain[1] += 1; chain_list[idx1][1] != domain_end) chain[1] += 1;
if(chain[0] != chain[1]) nb_chains[i].push_back(chain); if (chain[0] != chain[1]) nb_chains[i].push_back(chain);
idx_s = (j == nnb - 1) ? -1 : nb_list[j + 1]; idx_s = (j == nnb - 1) ? -1 : nb_list[j + 1];
} }
} }
@ -311,8 +293,9 @@ PairMESONTTPM::PairMESONTTPM(LAMMPS *lmp) : Pair(lmp) {
eatom_s = nullptr; eatom_s = nullptr;
eatom_b = nullptr; eatom_b = nullptr;
eatom_t = nullptr; eatom_t = nullptr;
nmax = 0;
instance_count++; instance_count++;
if(instance_count > 1) error->all(FLERR, if (instance_count > 1) error->all(FLERR,
"only a single instance of mesont/tpm pair style can be created"); "only a single instance of mesont/tpm pair style can be created");
} }
@ -335,13 +318,25 @@ PairMESONTTPM::~PairMESONTTPM()
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
void PairMESONTTPM::compute(int eflag, int vflag){ void PairMESONTTPM::compute(int eflag, int vflag) {
// set per atom values and accumulators
// reallocate per-atom arrays if necessary
ev_init(eflag,vflag); ev_init(eflag,vflag);
//total number of atoms in the node and ghost shell if (atom->nmax > nmax && eflag_atom) {
memory->destroy(eatom_s);
memory->create(eatom_s,comm->nthreads*maxeatom,"pair:eatom_s");
memory->destroy(eatom_b);
memory->create(eatom_b,comm->nthreads*maxeatom,"pair:eatom_b");
memory->destroy(eatom_t);
memory->create(eatom_t,comm->nthreads*maxeatom,"pair:eatom_t");
nmax = atom->nmax;
}
//total number of atoms in the node and ghost shell treated as NTs
int nall = list->inum + list->gnum; int nall = list->inum + list->gnum;
//total number of atoms in the node and ghost shell
int ntot = atom->nlocal + atom->nghost; int ntot = atom->nlocal + atom->nghost;
int newton_pair = force->newton_pair; int newton_pair = force->newton_pair;
if(!newton_pair) if (!newton_pair)
error->all(FLERR,"Pair style mesont/tpm requires newton pair on"); error->all(FLERR,"Pair style mesont/tpm requires newton pair on");
double **x = atom->x; double **x = atom->x;
@ -360,7 +355,7 @@ void PairMESONTTPM::compute(int eflag, int vflag){
} }
double Rcut_min = std::max(2.0*Lmax, std::sqrt(0.5*Lmax*Lmax + double Rcut_min = std::max(2.0*Lmax, std::sqrt(0.5*Lmax*Lmax +
std::pow((2.0*RT + TPBRcutoff),2))); std::pow((2.0*RT + TPBRcutoff),2)));
if (cut_global < Rcut_min){ if (cut_global < Rcut_min) {
std::stringstream err; std::stringstream err;
err << "The selected cutoff is too small for the current system : " << err << "The selected cutoff is too small for the current system : " <<
"L_max = " << Lmax << ", R_max = " << RT << ", Rc = " << cut_global << "L_max = " << Lmax << ", R_max = " << RT << ", Rc = " << cut_global <<
@ -369,14 +364,14 @@ void PairMESONTTPM::compute(int eflag, int vflag){
} }
//generate bonds and chain nblist //generate bonds and chain nblist
MESONTList ntlist(atom, list, cut_global*cut_global); MESONTList ntlist(atom, list);
//reorder data to make it contiguous within tubes //reorder data to make it contiguous within tubes
//and compatible with Fortran functions //and compatible with Fortran functions
std::vector<double> x_sort(3*nall), f_sort(3*nall), s_sort(9*nall); std::vector<double> x_sort(3*nall), f_sort(3*nall), s_sort(9*nall);
std::vector<double> u_ts_sort(nall), u_tb_sort(nall), u_tt_sort(nall); std::vector<double> u_ts_sort(nall), u_tb_sort(nall), u_tt_sort(nall);
std::vector<int> b_sort(nall); std::vector<int> b_sort(nall);
for (int i = 0; i < nall; i++){ for (int i = 0; i < nall; i++) {
int idx = ntlist.get_idx(i); int idx = ntlist.get_idx(i);
for (int j = 0; j < 3; j++) x_sort[3*i+j] = x[idx][j]; for (int j = 0; j < 3; j++) x_sort[3*i+j] = x[idx][j];
b_sort[i] = buckling[idx]; b_sort[i] = buckling[idx];
@ -385,7 +380,7 @@ void PairMESONTTPM::compute(int eflag, int vflag){
//bending potential //bending potential
int n_triplets = ntlist.get_triplets().size(); int n_triplets = ntlist.get_triplets().size();
for (int i = 0; i < n_triplets; i++) { for (int i = 0; i < n_triplets; i++) {
const array2003<int,3>& t = ntlist.get_triplets()[i]; const std::array<int,3>& t = ntlist.get_triplets()[i];
//idx of nodes of a triplet in sorted representation //idx of nodes of a triplet in sorted representation
int idx_s0 = ntlist.get_idxb(t[0]); int idx_s0 = ntlist.get_idxb(t[0]);
int idx_s1 = ntlist.get_idxb(t[1]); int idx_s1 = ntlist.get_idxb(t[1]);
@ -412,13 +407,13 @@ void PairMESONTTPM::compute(int eflag, int vflag){
} }
//share new values of buckling //share new values of buckling
if (BendingMode == 1){ if (BendingMode == 1) {
for (int i = 0; i < nall; i++){ for (int i = 0; i < nall; i++) {
int idx = ntlist.get_idx(i); int idx = ntlist.get_idx(i);
buckling[idx] = b_sort[i]; buckling[idx] = b_sort[i];
} }
comm->forward_comm_pair(this); comm->forward_comm_pair(this);
for (int i = 0; i < nall; i++){ for (int i = 0; i < nall; i++) {
int idx = ntlist.get_idx(i); int idx = ntlist.get_idx(i);
b_sort[i] = buckling[idx]; b_sort[i] = buckling[idx];
} }
@ -429,7 +424,7 @@ void PairMESONTTPM::compute(int eflag, int vflag){
double Rmax = 0.0; double Rmax = 0.0;
Lmax = 0.0; Lmax = 0.0;
for (int i = 0; i < n_segments; i++) { for (int i = 0; i < n_segments; i++) {
const array2003<int,2>& s = ntlist.get_segments()[i]; const std::array<int,2>& s = ntlist.get_segments()[i];
//idx of a segment end 1 in sorted representation //idx of a segment end 1 in sorted representation
int idx_s0 = ntlist.get_idxb(s[0]); int idx_s0 = ntlist.get_idxb(s[0]);
//idx of a segment end 2 in sorted representation //idx of a segment end 2 in sorted representation
@ -456,9 +451,9 @@ void PairMESONTTPM::compute(int eflag, int vflag){
mesont_lib_TubeStretchingForceField(U1s, U2s, F1, F2, S1, S2, X1, X2, mesont_lib_TubeStretchingForceField(U1s, U2s, F1, F2, S1, S2, X1, X2,
R12, L12); R12, L12);
for (int nc = 0; nc < (int)ntlist.get_nbs()[i].size(); nc++){ for (int nc = 0; nc < (int)ntlist.get_nbs()[i].size(); nc++) {
//id of the beginning and end of the chain in the sorted representation //id of the beginning and end of the chain in the sorted representation
const array2003<int,2>& chain = ntlist.get_nbs()[i][nc]; const std::array<int,2>& chain = ntlist.get_nbs()[i][nc];
int N = chain[1] - chain[0] + 1; //number of elements in the chain int N = chain[1] - chain[0] + 1; //number of elements in the chain
int end1 = ntlist.get_idx(chain[0]); //chain ends (real representation) int end1 = ntlist.get_idx(chain[0]); //chain ends (real representation)
int end2 = ntlist.get_idx(chain[1]); int end2 = ntlist.get_idx(chain[1]);
@ -475,7 +470,7 @@ void PairMESONTTPM::compute(int eflag, int vflag){
double* Xe = X; double* Fe = F; double* Se = S; double* Xe = X; double* Fe = F; double* Se = S;
if (!E1 && ntlist.get_triplet(end1)[0] != MESONTList::domain_end && if (!E1 && ntlist.get_triplet(end1)[0] != MESONTList::domain_end &&
ntlist.get_triplet(ntlist.get_triplet(end1)[0])[0] == ntlist.get_triplet(ntlist.get_triplet(end1)[0])[0] ==
MESONTList::cnt_end){ MESONTList::cnt_end) {
Ee = 1; Ee = 1;
int idx = ntlist.get_idxb(ntlist.get_triplet(end1)[0]); int idx = ntlist.get_idxb(ntlist.get_triplet(end1)[0]);
Xe = &(x_sort[3*idx]); Xe = &(x_sort[3*idx]);
@ -484,7 +479,7 @@ void PairMESONTTPM::compute(int eflag, int vflag){
} }
else if (!E2 && ntlist.get_triplet(end2)[2] != MESONTList::domain_end && else if (!E2 && ntlist.get_triplet(end2)[2] != MESONTList::domain_end &&
ntlist.get_triplet(ntlist.get_triplet(end2)[2])[2] == ntlist.get_triplet(ntlist.get_triplet(end2)[2])[2] ==
MESONTList::cnt_end){ MESONTList::cnt_end) {
Ee = 2; Ee = 2;
int idx = ntlist.get_idxb(ntlist.get_triplet(end2)[2]); int idx = ntlist.get_idxb(ntlist.get_triplet(end2)[2]);
Xe = &(x_sort[3*idx]); Xe = &(x_sort[3*idx]);
@ -500,7 +495,7 @@ void PairMESONTTPM::compute(int eflag, int vflag){
//check if cutoff is chosen correctly //check if cutoff is chosen correctly
Rcut_min = std::max(2.0*Lmax, std::sqrt(0.5*Lmax*Lmax + Rcut_min = std::max(2.0*Lmax, std::sqrt(0.5*Lmax*Lmax +
std::pow((2.0*Rmax + TPBRcutoff),2))); std::pow((2.0*Rmax + TPBRcutoff),2)));
if (cut_global < Rcut_min){ if (cut_global < Rcut_min) {
std::stringstream err; std::stringstream err;
err << "The selected cutoff is too small for the current system : " << err << "The selected cutoff is too small for the current system : " <<
"L_max = " << Lmax << ", R_max = " << RT << ", Rc = " << cut_global << "L_max = " << Lmax << ", R_max = " << RT << ", Rc = " << cut_global <<
@ -508,65 +503,69 @@ void PairMESONTTPM::compute(int eflag, int vflag){
error->all(FLERR, err.str().c_str()); error->all(FLERR, err.str().c_str());
} }
// set per atom values and accumulators //convert from sorted representation
// reallocate per-atom arrays if necessary for (int i = 0; i < nall; i++) {
if (atom->nmax > maxeatom) { int idx = ntlist.get_idx(i);
maxeatom = atom->nmax; for (int j = 0; j < 3; j++) f[idx][j] += f_sort[3*i+j];
memory->destroy(eatom); buckling[idx] = b_sort[i];
memory->create(eatom,comm->nthreads*maxeatom,"pair:eatom");
memory->destroy(eatom_s);
memory->create(eatom_s,comm->nthreads*maxeatom,"pair:eatom_s");
memory->destroy(eatom_b);
memory->create(eatom_b,comm->nthreads*maxeatom,"pair:eatom_b");
memory->destroy(eatom_t);
memory->create(eatom_t,comm->nthreads*maxeatom,"pair:eatom_t");
} }
if (eflag_global) {
if (atom->nmax > maxvatom) {
maxvatom = atom->nmax;
memory->destroy(vatom);
memory->create(vatom,comm->nthreads*maxvatom,6,"pair:vatom");
}
// zero accumulators
eng_vdwl = 0.0; energy_s = 0.0; eng_vdwl = 0.0; energy_s = 0.0;
energy_b = 0.0; energy_t = 0.0; energy_b = 0.0; energy_t = 0.0;
for (int i = 0; i < 6; i++) virial[i] = 0.0; for (int i = 0; i < nall; i++) {
for (int i = 0; i < ntot; i++){ int idx = ntlist.get_idx(i);
energy_s += u_ts_sort[i];
energy_b += u_tb_sort[i];
energy_t += u_tt_sort[i];
}
eng_vdwl = energy_s + energy_b + energy_t;
}
if (eflag_atom) {
for (int i = 0; i < ntot; i++) {
eatom[i] = 0.0; eatom_s[i] = 0.0; eatom[i] = 0.0; eatom_s[i] = 0.0;
eatom_b[i] = 0.0; eatom_t[i] = 0.0; eatom_b[i] = 0.0; eatom_t[i] = 0.0;
} }
for (int i = 0; i < ntot; i++) for (int i = 0; i < nall; i++) {
for (int j = 0; j < 6; j++) vatom[i][j] = 0.0;
//convert from sorted representation
for (int i = 0; i < nall; i++){
int idx = ntlist.get_idx(i); int idx = ntlist.get_idx(i);
for (int j = 0; j < 3; j++) f[idx][j] += f_sort[3*i+j];
eatom_s[idx] = u_ts_sort[i]; eatom_s[idx] = u_ts_sort[i];
eatom_b[idx] = u_tb_sort[i]; eatom_b[idx] = u_tb_sort[i];
eatom_t[idx] = u_tt_sort[i]; eatom_t[idx] = u_tt_sort[i];
eatom[idx] = u_ts_sort[i] + u_tb_sort[i] + u_tt_sort[i]; eatom[idx] = u_ts_sort[i] + u_tb_sort[i] + u_tt_sort[i];
energy_s += u_ts_sort[i]; }
energy_b += u_tb_sort[i]; }
energy_t += u_tt_sort[i]; if (vflag_global) {
for (int i = 0; i < 6; i++) virial[i] = 0.0;
for (int i = 0; i < nall; i++) {
int idx = ntlist.get_idx(i);
virial[0] += s_sort[9*i+0]; //xx
virial[1] += s_sort[9*i+4]; //yy
virial[2] += s_sort[9*i+8]; //zz
virial[3] += s_sort[9*i+1]; //xy
virial[4] += s_sort[9*i+2]; //xz
virial[5] += s_sort[9*i+5]; //yz
}
}
if (vflag_atom) {
for (int i = 0; i < ntot; i++)
for (int j = 0; j < 6; j++) vatom[i][j] = 0.0;
for (int i = 0; i < nall; i++) {
int idx = ntlist.get_idx(i);
vatom[idx][0] = s_sort[9*i+0]; //xx vatom[idx][0] = s_sort[9*i+0]; //xx
vatom[idx][1] = s_sort[9*i+4]; //yy vatom[idx][1] = s_sort[9*i+4]; //yy
vatom[idx][2] = s_sort[9*i+8]; //zz vatom[idx][2] = s_sort[9*i+8]; //zz
vatom[idx][3] = s_sort[9*i+1]; //xy vatom[idx][3] = s_sort[9*i+1]; //xy
vatom[idx][4] = s_sort[9*i+2]; //xz vatom[idx][4] = s_sort[9*i+2]; //xz
vatom[idx][5] = s_sort[9*i+5]; //yz vatom[idx][5] = s_sort[9*i+5]; //yz
for (int j = 0; j < 6; j++) virial[j] += vatom[idx][j];
buckling[idx] = b_sort[i];
} }
eng_vdwl = energy_s + energy_b + energy_t; }
} }
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
allocate all arrays allocate all arrays
------------------------------------------------------------------------- */ ------------------------------------------------------------------------- */
void PairMESONTTPM::allocate(){ void PairMESONTTPM::allocate() {
allocated = 1; allocated = 1;
int n = atom->ntypes; int n = atom->ntypes;
@ -583,7 +582,7 @@ void PairMESONTTPM::allocate(){
global settings global settings
------------------------------------------------------------------------- */ ------------------------------------------------------------------------- */
void PairMESONTTPM::settings(int narg, char **arg){ void PairMESONTTPM::settings(int narg, char **arg) {
if ((narg == 0) || (narg > 4)) if ((narg == 0) || (narg > 4))
error->all(FLERR,"Illegal pair_style command"); error->all(FLERR,"Illegal pair_style command");
cut_global = utils::numeric(FLERR,arg[0],false,lmp); cut_global = utils::numeric(FLERR,arg[0],false,lmp);
@ -632,7 +631,7 @@ void PairMESONTTPM::settings(int narg, char **arg){
set coeffs for one or more type pairs set coeffs for one or more type pairs
------------------------------------------------------------------------- */ ------------------------------------------------------------------------- */
void PairMESONTTPM::coeff(int narg, char **arg){ void PairMESONTTPM::coeff(int narg, char **arg) {
if ((narg < 2) || (narg > 3)) if ((narg < 2) || (narg > 3))
error->all(FLERR,"Incorrect args for pair coefficients"); error->all(FLERR,"Incorrect args for pair coefficients");
@ -661,7 +660,7 @@ void PairMESONTTPM::coeff(int narg, char **arg){
init for one type pair i,j and corresponding j,i init for one type pair i,j and corresponding j,i
------------------------------------------------------------------------- */ ------------------------------------------------------------------------- */
double PairMESONTTPM::init_one(int i, int j){ double PairMESONTTPM::init_one(int i, int j) {
if (setflag[i][j] == 0) { if (setflag[i][j] == 0) {
cut[i][j] = mix_distance(cut[i][i],cut[j][j]); cut[i][j] = mix_distance(cut[i][i],cut[j][j]);
} }
@ -673,7 +672,7 @@ double PairMESONTTPM::init_one(int i, int j){
proc 0 writes to restart file proc 0 writes to restart file
------------------------------------------------------------------------- */ ------------------------------------------------------------------------- */
void PairMESONTTPM::write_restart(FILE *fp){ void PairMESONTTPM::write_restart(FILE *fp) {
write_restart_settings(fp); write_restart_settings(fp);
int i,j; int i,j;
@ -690,7 +689,7 @@ void PairMESONTTPM::write_restart(FILE *fp){
proc 0 reads from restart file, bcasts proc 0 reads from restart file, bcasts
------------------------------------------------------------------------- */ ------------------------------------------------------------------------- */
void PairMESONTTPM::read_restart(FILE *fp){ void PairMESONTTPM::read_restart(FILE *fp) {
read_restart_settings(fp); read_restart_settings(fp);
allocate(); allocate();
@ -713,7 +712,7 @@ void PairMESONTTPM::read_restart(FILE *fp){
proc 0 writes to restart file proc 0 writes to restart file
------------------------------------------------------------------------- */ ------------------------------------------------------------------------- */
void PairMESONTTPM::write_restart_settings(FILE *fp){ void PairMESONTTPM::write_restart_settings(FILE *fp) {
fwrite(&BendingMode,sizeof(int),1,fp); fwrite(&BendingMode,sizeof(int),1,fp);
fwrite(&TPMType,sizeof(int),1,fp); fwrite(&TPMType,sizeof(int),1,fp);
fwrite(&cut_global,sizeof(double),1,fp); fwrite(&cut_global,sizeof(double),1,fp);
@ -725,7 +724,7 @@ void PairMESONTTPM::write_restart_settings(FILE *fp){
proc 0 reads from restart file, bcasts proc 0 reads from restart file, bcasts
------------------------------------------------------------------------- */ ------------------------------------------------------------------------- */
void PairMESONTTPM::read_restart_settings(FILE *fp){ void PairMESONTTPM::read_restart_settings(FILE *fp) {
int me = comm->me; int me = comm->me;
if (me == 0) { if (me == 0) {
fread(&BendingMode,sizeof(int),1,fp); fread(&BendingMode,sizeof(int),1,fp);
@ -761,7 +760,7 @@ void PairMESONTTPM::read_restart_settings(FILE *fp){
proc 0 writes to data file proc 0 writes to data file
------------------------------------------------------------------------- */ ------------------------------------------------------------------------- */
void PairMESONTTPM::write_data(FILE *fp){ void PairMESONTTPM::write_data(FILE *fp) {
for (int i = 1; i <= atom->ntypes; i++) for (int i = 1; i <= atom->ntypes; i++)
fprintf(fp,"%d\n",i); fprintf(fp,"%d\n",i);
} }
@ -770,7 +769,7 @@ void PairMESONTTPM::write_data(FILE *fp){
proc 0 writes all pairs to data file proc 0 writes all pairs to data file
------------------------------------------------------------------------- */ ------------------------------------------------------------------------- */
void PairMESONTTPM::write_data_all(FILE *fp){ void PairMESONTTPM::write_data_all(FILE *fp) {
for (int i = 1; i <= atom->ntypes; i++) for (int i = 1; i <= atom->ntypes; i++)
for (int j = i; j <= atom->ntypes; j++) for (int j = i; j <= atom->ntypes; j++)
fprintf(fp,"%d %d %g\n",i,j,cut[i][j]); fprintf(fp,"%d %d %g\n",i,j,cut[i][j]);
@ -778,7 +777,7 @@ void PairMESONTTPM::write_data_all(FILE *fp){
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
void PairMESONTTPM::init_style(){ void PairMESONTTPM::init_style() {
//make sure that a full list is created (including ghost nodes) //make sure that a full list is created (including ghost nodes)
int r = neighbor->request(this,instance_me); int r = neighbor->request(this,instance_me);
neighbor->requests[r]->half = false; neighbor->requests[r]->half = false;
@ -786,7 +785,7 @@ void PairMESONTTPM::init_style(){
neighbor->requests[r]->ghost = true; neighbor->requests[r]->ghost = true;
} }
void* PairMESONTTPM::extract(const char *str, int &){ void* PairMESONTTPM::extract(const char *str, int &) {
if (strcmp(str,"mesonttpm_Es_tot") == 0) return &energy_s; if (strcmp(str,"mesonttpm_Es_tot") == 0) return &energy_s;
else if (strcmp(str,"mesonttpm_Eb_tot") == 0) return &energy_b; else if (strcmp(str,"mesonttpm_Eb_tot") == 0) return &energy_b;
else if (strcmp(str,"mesonttpm_Et_tot") == 0) return &energy_t; else if (strcmp(str,"mesonttpm_Et_tot") == 0) return &energy_t;

View File

@ -54,6 +54,7 @@ class PairMESONTTPM : public Pair {
double cut_global; double cut_global;
double **cut; double **cut;
static int instance_count; static int instance_count;
int nmax;
virtual void allocate(); virtual void allocate();
virtual void *extract(const char *, int &); virtual void *extract(const char *, int &);

View File

@ -183,7 +183,7 @@ DumpNetCDF::DumpNetCDF(LAMMPS *lmp, int narg, char **arg) :
for (int j = 0; j < DUMP_NC_MAX_DIMS; j++) { for (int j = 0; j < DUMP_NC_MAX_DIMS; j++) {
perat[inc].field[j] = -1; perat[inc].field[j] = -1;
} }
strcpy(perat[inc].name, mangled); strncpy(perat[inc].name, mangled, NC_FIELD_NAME_MAX);
n_perat++; n_perat++;
} }

View File

@ -30,7 +30,7 @@ even if PLUMED is not in the path if as long as the input does not contain a fix
plumed command. plumed command.
If you wish to statically link PLUMED you must download PLUMED to the /lib/plumed directory before compiling LAMMPS. You can If you wish to statically link PLUMED you must download PLUMED to the /lib/plumed directory before compiling LAMMPS. You can
download a tar ball into that directory or you can clone the plumed2 repository from github there. Once you have created a download a tar ball into that directory or you can clone the plumed2 repository from GitHub there. Once you have created a
directory containing a distribution of PLUMED within /lib/plumed you then must build PLUMED within that directory by issuing directory containing a distribution of PLUMED within /lib/plumed you then must build PLUMED within that directory by issuing
the usual commands. It is worth noting that we have provided a script that will download and build PLUMED for you with the usual commands. It is worth noting that we have provided a script that will download and build PLUMED for you with
a minimal set of options. To run this script you need to issue the following command: a minimal set of options. To run this script you need to issue the following command:

View File

@ -42,6 +42,12 @@
#define ENERGY_MASK 0x00010000 #define ENERGY_MASK 0x00010000
#define VIRIAL_MASK 0x00020000 #define VIRIAL_MASK 0x00020000
// SPIN
#define SP_MASK 0x00000001
#define FM_MASK 0x00000002
#define FML_MASK 0x00000004
// DPD // DPD
#define DPDRHO_MASK 0x00040000 #define DPDRHO_MASK 0x00040000

View File

@ -533,8 +533,7 @@ void AtomVecEllipsoid::write_data_bonus(FILE *fp, int n, double *buf, int /*flag
this may create or delete entry in bonus data this may create or delete entry in bonus data
------------------------------------------------------------------------- */ ------------------------------------------------------------------------- */
void AtomVecEllipsoid:: void AtomVecEllipsoid::set_shape(int i, double shapex, double shapey, double shapez)
set_shape(int i, double shapex, double shapey, double shapez)
{ {
if (ellipsoid[i] < 0) { if (ellipsoid[i] < 0) {
if (shapex == 0.0 && shapey == 0.0 && shapez == 0.0) return; if (shapex == 0.0 && shapey == 0.0 && shapez == 0.0) return;

View File

@ -1873,6 +1873,7 @@ void Domain::set_boundary(int narg, char **arg, int flag)
else zperiodic = 0; else zperiodic = 0;
// record if we changed a periodic boundary to a non-periodic one // record if we changed a periodic boundary to a non-periodic one
int pflag=0; int pflag=0;
if ((periodicity[0] && !xperiodic) if ((periodicity[0] && !xperiodic)
|| (periodicity[1] && !yperiodic) || (periodicity[1] && !yperiodic)
@ -1889,23 +1890,27 @@ void Domain::set_boundary(int narg, char **arg, int flag)
boundary[1][0] >= 2 || boundary[1][1] >= 2 || boundary[1][0] >= 2 || boundary[1][1] >= 2 ||
boundary[2][0] >= 2 || boundary[2][1] >= 2) nonperiodic = 2; boundary[2][0] >= 2 || boundary[2][1] >= 2) nonperiodic = 2;
} }
// force non-zero image flags to zero for non-periodic dimensions
// keep track if a change was made, so we can print a warning message
if (pflag) { if (pflag) {
pflag = 0; pflag = 0;
for (int i=0; i < atom->nlocal; ++i) { for (int i=0; i < atom->nlocal; ++i) {
int xbox = (atom->image[i] & IMGMASK) - IMGMAX; int xbox = (atom->image[i] & IMGMASK) - IMGMAX;
int ybox = (atom->image[i] >> IMGBITS & IMGMASK) - IMGMAX; int ybox = (atom->image[i] >> IMGBITS & IMGMASK) - IMGMAX;
int zbox = (atom->image[i] >> IMG2BITS) - IMGMAX; int zbox = (atom->image[i] >> IMG2BITS) - IMGMAX;
if (!xperiodic) { xbox = 0; pflag = 1; } if ((!xperiodic) && (xbox != 0)) { xbox = 0; pflag = 1; }
if (!yperiodic) { ybox = 0; pflag = 1; } if ((!yperiodic) && (ybox != 0)) { ybox = 0; pflag = 1; }
if (!zperiodic) { zbox = 0; pflag = 1; } if ((!zperiodic) && (zbox != 0)) { zbox = 0; pflag = 1; }
atom->image[i] = ((imageint) (xbox + IMGMAX) & IMGMASK) | atom->image[i] = ((imageint) (xbox + IMGMAX) & IMGMASK) |
(((imageint) (ybox + IMGMAX) & IMGMASK) << IMGBITS) | (((imageint) (ybox + IMGMAX) & IMGMASK) << IMGBITS) |
(((imageint) (zbox + IMGMAX) & IMGMASK) << IMG2BITS); (((imageint) (zbox + IMGMAX) & IMGMASK) << IMG2BITS);
} }
int flag_all; int flag_all;
MPI_Allreduce(&flag,&flag_all, 1, MPI_INT, MPI_SUM, world); MPI_Allreduce(&pflag,&flag_all, 1, MPI_INT, MPI_SUM, world);
if ((flag_all > 0) && (comm->me == 0)) if ((flag_all > 0) && (comm->me == 0))
error->warning(FLERR,"Reset image flags for non-periodic boundary"); error->warning(FLERR,"Resetting image flags for non-periodic dimensions");
} }
} }

View File

@ -72,43 +72,27 @@ FMT_CONSTEXPR To lossless_integral_conversion(const From from, int& ec) {
static_assert(F::is_integer, "From must be integral"); static_assert(F::is_integer, "From must be integral");
static_assert(T::is_integer, "To must be integral"); static_assert(T::is_integer, "To must be integral");
if (F::is_signed && !T::is_signed) { if (detail::const_check(F::is_signed && !T::is_signed)) {
// From may be negative, not allowed! // From may be negative, not allowed!
if (fmt::detail::is_negative(from)) { if (fmt::detail::is_negative(from)) {
ec = 1; ec = 1;
return {}; return {};
} }
// From is positive. Can it always fit in To? // From is positive. Can it always fit in To?
if (F::digits <= T::digits) { if (F::digits > T::digits &&
// yes, From always fits in To. from > static_cast<From>(detail::max_value<To>())) {
} else {
// from may not fit in To, we have to do a dynamic check
if (from > static_cast<From>((T::max)())) {
ec = 1; ec = 1;
return {}; return {};
} }
} }
}
if (!F::is_signed && T::is_signed) { if (!F::is_signed && T::is_signed && F::digits >= T::digits &&
// can from be held in To? from > static_cast<From>(detail::max_value<To>())) {
if (F::digits < T::digits) {
// yes, From always fits in To.
} else {
// from may not fit in To, we have to do a dynamic check
if (from > static_cast<From>((T::max)())) {
// outside range.
ec = 1; ec = 1;
return {}; return {};
} }
} return static_cast<To>(from); // Lossless conversion.
} }
// reaching here means all is ok for lossless conversion.
return static_cast<To>(from);
} // function
template <typename To, typename From, template <typename To, typename From,
FMT_ENABLE_IF(std::is_same<From, To>::value)> FMT_ENABLE_IF(std::is_same<From, To>::value)>
@ -190,11 +174,9 @@ To safe_duration_cast(std::chrono::duration<FromRep, FromPeriod> from,
// safe conversion to IntermediateRep // safe conversion to IntermediateRep
IntermediateRep count = IntermediateRep count =
lossless_integral_conversion<IntermediateRep>(from.count(), ec); lossless_integral_conversion<IntermediateRep>(from.count(), ec);
if (ec) { if (ec) return {};
return {};
}
// multiply with Factor::num without overflow or underflow // multiply with Factor::num without overflow or underflow
if (Factor::num != 1) { if (detail::const_check(Factor::num != 1)) {
const auto max1 = detail::max_value<IntermediateRep>() / Factor::num; const auto max1 = detail::max_value<IntermediateRep>() / Factor::num;
if (count > max1) { if (count > max1) {
ec = 1; ec = 1;
@ -209,17 +191,9 @@ To safe_duration_cast(std::chrono::duration<FromRep, FromPeriod> from,
count *= Factor::num; count *= Factor::num;
} }
// this can't go wrong, right? den>0 is checked earlier. if (detail::const_check(Factor::den != 1)) count /= Factor::den;
if (Factor::den != 1) { auto tocount = lossless_integral_conversion<typename To::rep>(count, ec);
count /= Factor::den; return ec ? To() : To(tocount);
}
// convert to the to type, safely
using ToRep = typename To::rep;
const ToRep tocount = lossless_integral_conversion<ToRep>(count, ec);
if (ec) {
return {};
}
return To{tocount};
} }
/** /**
@ -351,6 +325,11 @@ inline std::tm localtime(std::time_t time) {
return lt.tm_; return lt.tm_;
} }
inline std::tm localtime(
std::chrono::time_point<std::chrono::system_clock> time_point) {
return localtime(std::chrono::system_clock::to_time_t(time_point));
}
// Thread-safe replacement for std::gmtime // Thread-safe replacement for std::gmtime
inline std::tm gmtime(std::time_t time) { inline std::tm gmtime(std::time_t time) {
struct dispatcher { struct dispatcher {
@ -387,6 +366,11 @@ inline std::tm gmtime(std::time_t time) {
return gt.tm_; return gt.tm_;
} }
inline std::tm gmtime(
std::chrono::time_point<std::chrono::system_clock> time_point) {
return gmtime(std::chrono::system_clock::to_time_t(time_point));
}
namespace detail { namespace detail {
inline size_t strftime(char* str, size_t count, const char* format, inline size_t strftime(char* str, size_t count, const char* format,
const std::tm* time) { const std::tm* time) {
@ -399,6 +383,17 @@ inline size_t strftime(wchar_t* str, size_t count, const wchar_t* format,
} }
} // namespace detail } // namespace detail
template <typename Char>
struct formatter<std::chrono::time_point<std::chrono::system_clock>, Char>
: formatter<std::tm, Char> {
template <typename FormatContext>
auto format(std::chrono::time_point<std::chrono::system_clock> val,
FormatContext& ctx) -> decltype(ctx.out()) {
std::tm time = localtime(val);
return formatter<std::tm, Char>::format(time, ctx);
}
};
template <typename Char> struct formatter<std::tm, Char> { template <typename Char> struct formatter<std::tm, Char> {
template <typename ParseContext> template <typename ParseContext>
auto parse(ParseContext& ctx) -> decltype(ctx.begin()) { auto parse(ParseContext& ctx) -> decltype(ctx.begin()) {

View File

@ -463,16 +463,16 @@ template <> inline void reset_color<wchar_t>(FILE* stream) FMT_NOEXCEPT {
} }
template <typename Char> template <typename Char>
inline void reset_color(basic_memory_buffer<Char>& buffer) FMT_NOEXCEPT { inline void reset_color(buffer<Char>& buffer) FMT_NOEXCEPT {
const char* begin = data::reset_color; const char* begin = data::reset_color;
const char* end = begin + sizeof(data::reset_color) - 1; const char* end = begin + sizeof(data::reset_color) - 1;
buffer.append(begin, end); buffer.append(begin, end);
} }
template <typename Char> template <typename Char>
void vformat_to(basic_memory_buffer<Char>& buf, const text_style& ts, void vformat_to(buffer<Char>& buf, const text_style& ts,
basic_string_view<Char> format_str, basic_string_view<Char> format_str,
basic_format_args<buffer_context<Char>> args) { basic_format_args<buffer_context<type_identity_t<Char>>> args) {
bool has_style = false; bool has_style = false;
if (ts.has_emphasis()) { if (ts.has_emphasis()) {
has_style = true; has_style = true;
@ -496,7 +496,7 @@ void vformat_to(basic_memory_buffer<Char>& buf, const text_style& ts,
template <typename S, typename Char = char_t<S>> template <typename S, typename Char = char_t<S>>
void vprint(std::FILE* f, const text_style& ts, const S& format, void vprint(std::FILE* f, const text_style& ts, const S& format,
basic_format_args<buffer_context<Char>> args) { basic_format_args<buffer_context<type_identity_t<Char>>> args) {
basic_memory_buffer<Char> buf; basic_memory_buffer<Char> buf;
detail::vformat_to(buf, ts, to_string_view(format), args); detail::vformat_to(buf, ts, to_string_view(format), args);
buf.push_back(Char(0)); buf.push_back(Char(0));
@ -504,20 +504,22 @@ void vprint(std::FILE* f, const text_style& ts, const S& format,
} }
/** /**
\rst
Formats a string and prints it to the specified file stream using ANSI Formats a string and prints it to the specified file stream using ANSI
escape sequences to specify text formatting. escape sequences to specify text formatting.
Example:
**Example**::
fmt::print(fmt::emphasis::bold | fg(fmt::color::red), fmt::print(fmt::emphasis::bold | fg(fmt::color::red),
"Elapsed time: {0:.2f} seconds", 1.23); "Elapsed time: {0:.2f} seconds", 1.23);
\endrst
*/ */
template <typename S, typename... Args, template <typename S, typename... Args,
FMT_ENABLE_IF(detail::is_string<S>::value)> FMT_ENABLE_IF(detail::is_string<S>::value)>
void print(std::FILE* f, const text_style& ts, const S& format_str, void print(std::FILE* f, const text_style& ts, const S& format_str,
const Args&... args) { const Args&... args) {
detail::check_format_string<Args...>(format_str); vprint(f, ts, format_str,
using context = buffer_context<char_t<S>>; fmt::make_args_checked<Args...>(format_str, args...));
format_arg_store<context, Args...> as{args...};
vprint(f, ts, format_str, basic_format_args<context>(as));
} }
/** /**
@ -558,7 +560,42 @@ template <typename S, typename... Args, typename Char = char_t<S>>
inline std::basic_string<Char> format(const text_style& ts, const S& format_str, inline std::basic_string<Char> format(const text_style& ts, const S& format_str,
const Args&... args) { const Args&... args) {
return vformat(ts, to_string_view(format_str), return vformat(ts, to_string_view(format_str),
detail::make_args_checked<Args...>(format_str, args...)); fmt::make_args_checked<Args...>(format_str, args...));
}
/**
Formats a string with the given text_style and writes the output to ``out``.
*/
template <typename OutputIt, typename Char,
FMT_ENABLE_IF(detail::is_output_iterator<OutputIt, Char>::value)>
OutputIt vformat_to(
OutputIt out, const text_style& ts, basic_string_view<Char> format_str,
basic_format_args<buffer_context<type_identity_t<Char>>> args) {
decltype(detail::get_buffer<Char>(out)) buf(detail::get_buffer_init(out));
detail::vformat_to(buf, ts, format_str, args);
return detail::get_iterator(buf);
}
/**
\rst
Formats arguments with the given text_style, writes the result to the output
iterator ``out`` and returns the iterator past the end of the output range.
**Example**::
std::vector<char> out;
fmt::format_to(std::back_inserter(out),
fmt::emphasis::bold | fg(fmt::color::red), "{}", 42);
\endrst
*/
template <typename OutputIt, typename S, typename... Args,
bool enable = detail::is_output_iterator<OutputIt, char_t<S>>::value&&
detail::is_string<S>::value>
inline auto format_to(OutputIt out, const text_style& ts, const S& format_str,
Args&&... args) ->
typename std::enable_if<enable, OutputIt>::type {
return vformat_to(out, ts, to_string_view(format_str),
fmt::make_args_checked<Args...>(format_str, args...));
} }
FMT_END_NAMESPACE FMT_END_NAMESPACE

View File

@ -368,7 +368,8 @@ template <typename... Args> struct type_list {};
// Returns a reference to the argument at index N from [first, rest...]. // Returns a reference to the argument at index N from [first, rest...].
template <int N, typename T, typename... Args> template <int N, typename T, typename... Args>
constexpr const auto& get(const T& first, const Args&... rest) { constexpr const auto& get([[maybe_unused]] const T& first,
[[maybe_unused]] const Args&... rest) {
static_assert(N < 1 + sizeof...(Args), "index is out of bounds"); static_assert(N < 1 + sizeof...(Args), "index is out of bounds");
if constexpr (N == 0) if constexpr (N == 0)
return first; return first;
@ -406,6 +407,19 @@ constexpr text<Char> make_text(basic_string_view<Char> s, size_t pos,
return {{&s[pos], size}}; return {{&s[pos], size}};
} }
template <typename Char> struct code_unit {
Char value;
using char_type = Char;
template <typename OutputIt, typename... Args>
OutputIt format(OutputIt out, const Args&...) const {
return write<Char>(out, value);
}
};
template <typename Char>
struct is_compiled_format<code_unit<Char>> : std::true_type {};
// A replacement field that refers to argument N. // A replacement field that refers to argument N.
template <typename Char, typename T, int N> struct field { template <typename Char, typename T, int N> struct field {
using char_type = Char; using char_type = Char;
@ -430,7 +444,9 @@ template <typename Char, typename T, int N> struct spec_field {
OutputIt format(OutputIt out, const Args&... args) const { OutputIt format(OutputIt out, const Args&... args) const {
// This ensures that the argument type is convertile to `const T&`. // This ensures that the argument type is convertile to `const T&`.
const T& arg = get<N>(args...); const T& arg = get<N>(args...);
basic_format_context<OutputIt, Char> ctx(out, {}); const auto& vargs =
make_format_args<basic_format_context<OutputIt, Char>>(args...);
basic_format_context<OutputIt, Char> ctx(out, vargs);
return fmt.format(arg, ctx); return fmt.format(arg, ctx);
} }
}; };
@ -489,16 +505,17 @@ constexpr auto parse_tail(T head, S format_str) {
template <typename T, typename Char> struct parse_specs_result { template <typename T, typename Char> struct parse_specs_result {
formatter<T, Char> fmt; formatter<T, Char> fmt;
size_t end; size_t end;
int next_arg_id;
}; };
template <typename T, typename Char> template <typename T, typename Char>
constexpr parse_specs_result<T, Char> parse_specs(basic_string_view<Char> str, constexpr parse_specs_result<T, Char> parse_specs(basic_string_view<Char> str,
size_t pos) { size_t pos, int arg_id) {
str.remove_prefix(pos); str.remove_prefix(pos);
auto ctx = basic_format_parse_context<Char>(str); auto ctx = basic_format_parse_context<Char>(str, {}, arg_id + 1);
auto f = formatter<T, Char>(); auto f = formatter<T, Char>();
auto end = f.parse(ctx); auto end = f.parse(ctx);
return {f, pos + (end - str.data()) + 1}; return {f, pos + (end - str.data()) + 1, ctx.next_arg_id()};
} }
// Compiles a non-empty format string and returns the compiled representation // Compiles a non-empty format string and returns the compiled representation
@ -518,8 +535,8 @@ constexpr auto compile_format_string(S format_str) {
format_str); format_str);
} else if constexpr (str[POS + 1] == ':') { } else if constexpr (str[POS + 1] == ':') {
using type = get_type<ID, Args>; using type = get_type<ID, Args>;
constexpr auto result = parse_specs<type>(str, POS + 2); constexpr auto result = parse_specs<type>(str, POS + 2, ID);
return parse_tail<Args, result.end, ID + 1>( return parse_tail<Args, result.end, result.next_arg_id>(
spec_field<char_type, type, ID>{result.fmt}, format_str); spec_field<char_type, type, ID>{result.fmt}, format_str);
} else { } else {
return unknown_format(); return unknown_format();
@ -530,8 +547,13 @@ constexpr auto compile_format_string(S format_str) {
return parse_tail<Args, POS + 2, ID>(make_text(str, POS, 1), format_str); return parse_tail<Args, POS + 2, ID>(make_text(str, POS, 1), format_str);
} else { } else {
constexpr auto end = parse_text(str, POS + 1); constexpr auto end = parse_text(str, POS + 1);
if constexpr (end - POS > 1) {
return parse_tail<Args, end, ID>(make_text(str, POS, end - POS), return parse_tail<Args, end, ID>(make_text(str, POS, end - POS),
format_str); format_str);
} else {
return parse_tail<Args, end, ID>(code_unit<char_type>{str[POS]},
format_str);
}
} }
} }
@ -587,8 +609,7 @@ template <typename CompiledFormat, typename... Args,
FMT_INLINE std::basic_string<Char> format(const CompiledFormat& cf, FMT_INLINE std::basic_string<Char> format(const CompiledFormat& cf,
const Args&... args) { const Args&... args) {
basic_memory_buffer<Char> buffer; basic_memory_buffer<Char> buffer;
detail::buffer<Char>& base = buffer; cf.format(detail::buffer_appender<Char>(buffer), args...);
cf.format(std::back_inserter(base), args...);
return to_string(buffer); return to_string(buffer);
} }
@ -608,8 +629,7 @@ template <typename CompiledFormat, typename... Args,
std::basic_string<Char> format(const CompiledFormat& cf, const Args&... args) { std::basic_string<Char> format(const CompiledFormat& cf, const Args&... args) {
basic_memory_buffer<Char> buffer; basic_memory_buffer<Char> buffer;
using context = buffer_context<Char>; using context = buffer_context<Char>;
detail::buffer<Char>& base = buffer; detail::cf::vformat_to<context>(detail::buffer_appender<Char>(buffer), cf,
detail::cf::vformat_to<context>(std::back_inserter(base), cf,
make_format_args<context>(args...)); make_format_args<context>(args...));
return to_string(buffer); return to_string(buffer);
} }
@ -618,9 +638,13 @@ template <typename S, typename... Args,
FMT_ENABLE_IF(detail::is_compiled_string<S>::value)> FMT_ENABLE_IF(detail::is_compiled_string<S>::value)>
FMT_INLINE std::basic_string<typename S::char_type> format(const S&, FMT_INLINE std::basic_string<typename S::char_type> format(const S&,
Args&&... args) { Args&&... args) {
#ifdef __cpp_if_constexpr
if constexpr (std::is_same<typename S::char_type, char>::value) {
constexpr basic_string_view<typename S::char_type> str = S(); constexpr basic_string_view<typename S::char_type> str = S();
if (str.size() == 2 && str[0] == '{' && str[1] == '}') if (str.size() == 2 && str[0] == '{' && str[1] == '}')
return fmt::to_string(detail::first(args...)); return fmt::to_string(detail::first(args...));
}
#endif
constexpr auto compiled = detail::compile<Args...>(S()); constexpr auto compiled = detail::compile<Args...>(S());
return format(compiled, std::forward<Args>(args)...); return format(compiled, std::forward<Args>(args)...);
} }
@ -643,18 +667,30 @@ OutputIt format_to(OutputIt out, const S&, const Args&... args) {
return format_to(out, compiled, args...); return format_to(out, compiled, args...);
} }
template < template <typename OutputIt, typename CompiledFormat, typename... Args>
typename OutputIt, typename CompiledFormat, typename... Args, auto format_to_n(OutputIt out, size_t n, const CompiledFormat& cf,
FMT_ENABLE_IF(detail::is_output_iterator<OutputIt>::value&& std::is_base_of< const Args&... args) ->
detail::basic_compiled_format, CompiledFormat>::value)> typename std::enable_if<
format_to_n_result<OutputIt> format_to_n(OutputIt out, size_t n, detail::is_output_iterator<OutputIt,
const CompiledFormat& cf, typename CompiledFormat::char_type>::value &&
const Args&... args) { std::is_base_of<detail::basic_compiled_format,
CompiledFormat>::value,
format_to_n_result<OutputIt>>::type {
auto it = auto it =
format_to(detail::truncating_iterator<OutputIt>(out, n), cf, args...); format_to(detail::truncating_iterator<OutputIt>(out, n), cf, args...);
return {it.base(), it.count()}; return {it.base(), it.count()};
} }
template <typename OutputIt, typename S, typename... Args,
FMT_ENABLE_IF(detail::is_compiled_string<S>::value)>
format_to_n_result<OutputIt> format_to_n(OutputIt out, size_t n, const S&,
const Args&... args) {
constexpr auto compiled = detail::compile<Args...>(S());
auto it = format_to(detail::truncating_iterator<OutputIt>(out, n), compiled,
args...);
return {it.base(), it.count()};
}
template <typename CompiledFormat, typename... Args> template <typename CompiledFormat, typename... Args>
size_t formatted_size(const CompiledFormat& cf, const Args&... args) { size_t formatted_size(const CompiledFormat& cf, const Args&... args) {
return format_to(detail::counting_iterator(), cf, args...).count(); return format_to(detail::counting_iterator(), cf, args...).count();

View File

@ -18,7 +18,7 @@
#include <vector> #include <vector>
// The fmt library version in the form major * 10000 + minor * 100 + patch. // The fmt library version in the form major * 10000 + minor * 100 + patch.
#define FMT_VERSION 70003 #define FMT_VERSION 70102
#ifdef __clang__ #ifdef __clang__
# define FMT_CLANG_VERSION (__clang_major__ * 100 + __clang_minor__) # define FMT_CLANG_VERSION (__clang_major__ * 100 + __clang_minor__)
@ -57,6 +57,7 @@
# define FMT_MSC_VER 0 # define FMT_MSC_VER 0
# define FMT_SUPPRESS_MSC_WARNING(n) # define FMT_SUPPRESS_MSC_WARNING(n)
#endif #endif
#ifdef __has_feature #ifdef __has_feature
# define FMT_HAS_FEATURE(x) __has_feature(x) # define FMT_HAS_FEATURE(x) __has_feature(x)
#else #else
@ -64,7 +65,7 @@
#endif #endif
#if defined(__has_include) && !defined(__INTELLISENSE__) && \ #if defined(__has_include) && !defined(__INTELLISENSE__) && \
!(FMT_ICC_VERSION && FMT_ICC_VERSION < 1600) (!FMT_ICC_VERSION || FMT_ICC_VERSION >= 1600)
# define FMT_HAS_INCLUDE(x) __has_include(x) # define FMT_HAS_INCLUDE(x) __has_include(x)
#else #else
# define FMT_HAS_INCLUDE(x) 0 # define FMT_HAS_INCLUDE(x) 0
@ -99,7 +100,7 @@
#endif #endif
#ifndef FMT_OVERRIDE #ifndef FMT_OVERRIDE
# if FMT_HAS_FEATURE(cxx_override) || \ # if FMT_HAS_FEATURE(cxx_override_control) || \
(FMT_GCC_VERSION >= 408 && FMT_HAS_GXX_CXX11) || FMT_MSC_VER >= 1900 (FMT_GCC_VERSION >= 408 && FMT_HAS_GXX_CXX11) || FMT_MSC_VER >= 1900
# define FMT_OVERRIDE override # define FMT_OVERRIDE override
# else # else
@ -152,7 +153,7 @@
# if FMT_HAS_CPP14_ATTRIBUTE(deprecated) || FMT_MSC_VER >= 1900 # if FMT_HAS_CPP14_ATTRIBUTE(deprecated) || FMT_MSC_VER >= 1900
# define FMT_DEPRECATED [[deprecated]] # define FMT_DEPRECATED [[deprecated]]
# else # else
# if defined(__GNUC__) || defined(__clang__) # if (defined(__GNUC__) && !defined(__LCC__)) || defined(__clang__)
# define FMT_DEPRECATED __attribute__((deprecated)) # define FMT_DEPRECATED __attribute__((deprecated))
# elif FMT_MSC_VER # elif FMT_MSC_VER
# define FMT_DEPRECATED __declspec(deprecated) # define FMT_DEPRECATED __declspec(deprecated)
@ -177,6 +178,15 @@
# endif # endif
#endif #endif
#ifndef FMT_USE_INLINE_NAMESPACES
# if FMT_HAS_FEATURE(cxx_inline_namespaces) || FMT_GCC_VERSION >= 404 || \
(FMT_MSC_VER >= 1900 && !_MANAGED)
# define FMT_USE_INLINE_NAMESPACES 1
# else
# define FMT_USE_INLINE_NAMESPACES 0
# endif
#endif
// LAMMPS customization // LAMMPS customization
// use 'v7_lmp' namespace instead of 'v7' so that our // use 'v7_lmp' namespace instead of 'v7' so that our
// bundled copy does not collide with linking other code // bundled copy does not collide with linking other code
@ -184,8 +194,7 @@
// a different version. // a different version.
#ifndef FMT_BEGIN_NAMESPACE #ifndef FMT_BEGIN_NAMESPACE
# if FMT_HAS_FEATURE(cxx_inline_namespaces) || FMT_GCC_VERSION >= 404 || \ # if FMT_USE_INLINE_NAMESPACES
FMT_MSC_VER >= 1900
# define FMT_INLINE_NAMESPACE inline namespace # define FMT_INLINE_NAMESPACE inline namespace
# define FMT_END_NAMESPACE \ # define FMT_END_NAMESPACE \
} \ } \
@ -275,8 +284,7 @@ struct monostate {};
namespace detail { namespace detail {
// A helper function to suppress bogus "conditional expression is constant" // A helper function to suppress "conditional expression is constant" warnings.
// warnings.
template <typename T> constexpr T const_check(T value) { return value; } template <typename T> constexpr T const_check(T value) { return value; }
FMT_NORETURN FMT_API void assert_fail(const char* file, int line, FMT_NORETURN FMT_API void assert_fail(const char* file, int line,
@ -305,7 +313,8 @@ template <typename T> struct std_string_view {};
#ifdef FMT_USE_INT128 #ifdef FMT_USE_INT128
// Do nothing. // Do nothing.
#elif defined(__SIZEOF_INT128__) && !FMT_NVCC && !(FMT_CLANG_VERSION && FMT_MSC_VER) #elif defined(__SIZEOF_INT128__) && !FMT_NVCC && \
!(FMT_CLANG_VERSION && FMT_MSC_VER)
# define FMT_USE_INT128 1 # define FMT_USE_INT128 1
using int128_t = __int128_t; using int128_t = __int128_t;
using uint128_t = __uint128_t; using uint128_t = __uint128_t;
@ -514,6 +523,18 @@ template <typename S> struct char_t_impl<S, enable_if_t<is_string<S>::value>> {
using type = typename result::value_type; using type = typename result::value_type;
}; };
// Reports a compile-time error if S is not a valid format string.
template <typename..., typename S, FMT_ENABLE_IF(!is_compile_string<S>::value)>
FMT_INLINE void check_format_string(const S&) {
#ifdef FMT_ENFORCE_COMPILE_STRING
static_assert(is_compile_string<S>::value,
"FMT_ENFORCE_COMPILE_STRING requires all format strings to use "
"FMT_STRING.");
#endif
}
template <typename..., typename S, FMT_ENABLE_IF(is_compile_string<S>::value)>
void check_format_string(S);
struct error_handler { struct error_handler {
constexpr error_handler() = default; constexpr error_handler() = default;
constexpr error_handler(const error_handler&) = default; constexpr error_handler(const error_handler&) = default;
@ -553,8 +574,9 @@ class basic_format_parse_context : private ErrorHandler {
using iterator = typename basic_string_view<Char>::iterator; using iterator = typename basic_string_view<Char>::iterator;
explicit constexpr basic_format_parse_context( explicit constexpr basic_format_parse_context(
basic_string_view<Char> format_str, ErrorHandler eh = {}) basic_string_view<Char> format_str, ErrorHandler eh = {},
: ErrorHandler(eh), format_str_(format_str), next_arg_id_(0) {} int next_arg_id = 0)
: ErrorHandler(eh), format_str_(format_str), next_arg_id_(next_arg_id) {}
/** /**
Returns an iterator to the beginning of the format string range being Returns an iterator to the beginning of the format string range being
@ -624,8 +646,24 @@ template <typename T, typename Context>
using has_formatter = using has_formatter =
std::is_constructible<typename Context::template formatter_type<T>>; std::is_constructible<typename Context::template formatter_type<T>>;
// Checks whether T is a container with contiguous storage.
template <typename T> struct is_contiguous : std::false_type {};
template <typename Char>
struct is_contiguous<std::basic_string<Char>> : std::true_type {};
namespace detail { namespace detail {
// Extracts a reference to the container from back_insert_iterator.
template <typename Container>
inline Container& get_container(std::back_insert_iterator<Container> it) {
using bi_iterator = std::back_insert_iterator<Container>;
struct accessor : bi_iterator {
accessor(bi_iterator iter) : bi_iterator(iter) {}
using bi_iterator::container;
};
return *accessor(it).container;
}
/** /**
\rst \rst
A contiguous memory buffer with an optional growing ability. It is an internal A contiguous memory buffer with an optional growing ability. It is an internal
@ -648,6 +686,8 @@ template <typename T> class buffer {
size_(sz), size_(sz),
capacity_(cap) {} capacity_(cap) {}
~buffer() = default;
/** Sets the buffer data and capacity. */ /** Sets the buffer data and capacity. */
void set(T* buf_data, size_t buf_capacity) FMT_NOEXCEPT { void set(T* buf_data, size_t buf_capacity) FMT_NOEXCEPT {
ptr_ = buf_data; ptr_ = buf_data;
@ -663,7 +703,6 @@ template <typename T> class buffer {
buffer(const buffer&) = delete; buffer(const buffer&) = delete;
void operator=(const buffer&) = delete; void operator=(const buffer&) = delete;
virtual ~buffer() = default;
T* begin() FMT_NOEXCEPT { return ptr_; } T* begin() FMT_NOEXCEPT { return ptr_; }
T* end() FMT_NOEXCEPT { return ptr_ + size_; } T* end() FMT_NOEXCEPT { return ptr_ + size_; }
@ -683,24 +722,26 @@ template <typename T> class buffer {
/** Returns a pointer to the buffer data. */ /** Returns a pointer to the buffer data. */
const T* data() const FMT_NOEXCEPT { return ptr_; } const T* data() const FMT_NOEXCEPT { return ptr_; }
/**
Resizes the buffer. If T is a POD type new elements may not be initialized.
*/
void resize(size_t new_size) {
reserve(new_size);
size_ = new_size;
}
/** Clears this buffer. */ /** Clears this buffer. */
void clear() { size_ = 0; } void clear() { size_ = 0; }
/** Reserves space to store at least *capacity* elements. */ // Tries resizing the buffer to contain *count* elements. If T is a POD type
void reserve(size_t new_capacity) { // the new elements may not be initialized.
void try_resize(size_t count) {
try_reserve(count);
size_ = count <= capacity_ ? count : capacity_;
}
// Tries increasing the buffer capacity to *new_capacity*. It can increase the
// capacity by a smaller amount than requested but guarantees there is space
// for at least one additional element either by increasing the capacity or by
// flushing the buffer if it is full.
void try_reserve(size_t new_capacity) {
if (new_capacity > capacity_) grow(new_capacity); if (new_capacity > capacity_) grow(new_capacity);
} }
void push_back(const T& value) { void push_back(const T& value) {
reserve(size_ + 1); try_reserve(size_ + 1);
ptr_[size_++] = value; ptr_[size_++] = value;
} }
@ -713,32 +754,150 @@ template <typename T> class buffer {
} }
}; };
// A container-backed buffer. struct buffer_traits {
explicit buffer_traits(size_t) {}
size_t count() const { return 0; }
size_t limit(size_t size) { return size; }
};
class fixed_buffer_traits {
private:
size_t count_ = 0;
size_t limit_;
public:
explicit fixed_buffer_traits(size_t limit) : limit_(limit) {}
size_t count() const { return count_; }
size_t limit(size_t size) {
size_t n = limit_ - count_;
count_ += size;
return size < n ? size : n;
}
};
// A buffer that writes to an output iterator when flushed.
template <typename OutputIt, typename T, typename Traits = buffer_traits>
class iterator_buffer final : public Traits, public buffer<T> {
private:
OutputIt out_;
enum { buffer_size = 256 };
T data_[buffer_size];
protected:
void grow(size_t) final FMT_OVERRIDE {
if (this->size() == buffer_size) flush();
}
void flush();
public:
explicit iterator_buffer(OutputIt out, size_t n = buffer_size)
: Traits(n),
buffer<T>(data_, 0, n < size_t(buffer_size) ? n : size_t(buffer_size)),
out_(out) {}
~iterator_buffer() { flush(); }
OutputIt out() {
flush();
return out_;
}
size_t count() const { return Traits::count() + this->size(); }
};
template <typename T> class iterator_buffer<T*, T> final : public buffer<T> {
protected:
void grow(size_t) final FMT_OVERRIDE {}
public:
explicit iterator_buffer(T* out, size_t = 0) : buffer<T>(out, 0, ~size_t()) {}
T* out() { return &*this->end(); }
};
// A buffer that writes to a container with the contiguous storage.
template <typename Container> template <typename Container>
class container_buffer : public buffer<typename Container::value_type> { class iterator_buffer<std::back_insert_iterator<Container>,
enable_if_t<is_contiguous<Container>::value,
typename Container::value_type>>
final : public buffer<typename Container::value_type> {
private: private:
Container& container_; Container& container_;
protected: protected:
void grow(size_t capacity) FMT_OVERRIDE { void grow(size_t capacity) final FMT_OVERRIDE {
container_.resize(capacity); container_.resize(capacity);
this->set(&container_[0], capacity); this->set(&container_[0], capacity);
} }
public: public:
explicit container_buffer(Container& c) explicit iterator_buffer(Container& c)
: buffer<typename Container::value_type>(c.size()), container_(c) {} : buffer<typename Container::value_type>(c.size()), container_(c) {}
explicit iterator_buffer(std::back_insert_iterator<Container> out, size_t = 0)
: iterator_buffer(get_container(out)) {}
std::back_insert_iterator<Container> out() {
return std::back_inserter(container_);
}
}; };
// Extracts a reference to the container from back_insert_iterator. // A buffer that counts the number of code units written discarding the output.
template <typename Container> template <typename T = char> class counting_buffer final : public buffer<T> {
inline Container& get_container(std::back_insert_iterator<Container> it) { private:
using bi_iterator = std::back_insert_iterator<Container>; enum { buffer_size = 256 };
struct accessor : bi_iterator { T data_[buffer_size];
accessor(bi_iterator iter) : bi_iterator(iter) {} size_t count_ = 0;
using bi_iterator::container;
}; protected:
return *accessor(it).container; void grow(size_t) final FMT_OVERRIDE {
if (this->size() != buffer_size) return;
count_ += this->size();
this->clear();
}
public:
counting_buffer() : buffer<T>(data_, 0, buffer_size) {}
size_t count() { return count_ + this->size(); }
};
// An output iterator that appends to the buffer.
// It is used to reduce symbol sizes for the common case.
template <typename T>
class buffer_appender : public std::back_insert_iterator<buffer<T>> {
using base = std::back_insert_iterator<buffer<T>>;
public:
explicit buffer_appender(buffer<T>& buf) : base(buf) {}
buffer_appender(base it) : base(it) {}
buffer_appender& operator++() {
base::operator++();
return *this;
}
buffer_appender operator++(int) {
buffer_appender tmp = *this;
++*this;
return tmp;
}
};
// Maps an output iterator into a buffer.
template <typename T, typename OutputIt>
iterator_buffer<OutputIt, T> get_buffer(OutputIt);
template <typename T> buffer<T>& get_buffer(buffer_appender<T>);
template <typename OutputIt> OutputIt get_buffer_init(OutputIt out) {
return out;
}
template <typename T> buffer<T>& get_buffer_init(buffer_appender<T> out) {
return get_container(out);
}
template <typename Buffer>
auto get_iterator(Buffer& buf) -> decltype(buf.out()) {
return buf.out();
}
template <typename T> buffer_appender<T> get_iterator(buffer<T>& buf) {
return buffer_appender<T>(buf);
} }
template <typename T, typename Char = char, typename Enable = void> template <typename T, typename Char = char, typename Enable = void>
@ -767,7 +926,8 @@ template <typename Char> struct named_arg_info {
template <typename T, typename Char, size_t NUM_ARGS, size_t NUM_NAMED_ARGS> template <typename T, typename Char, size_t NUM_ARGS, size_t NUM_NAMED_ARGS>
struct arg_data { struct arg_data {
// args_[0].named_args points to named_args_ to avoid bloating format_args. // args_[0].named_args points to named_args_ to avoid bloating format_args.
T args_[1 + (NUM_ARGS != 0 ? NUM_ARGS : 1)]; // +1 to workaround a bug in gcc 7.5 that causes duplicated-branches warning.
T args_[1 + (NUM_ARGS != 0 ? NUM_ARGS : +1)];
named_arg_info<Char> named_args_[NUM_NAMED_ARGS]; named_arg_info<Char> named_args_[NUM_NAMED_ARGS];
template <typename... U> template <typename... U>
@ -779,7 +939,8 @@ struct arg_data {
template <typename T, typename Char, size_t NUM_ARGS> template <typename T, typename Char, size_t NUM_ARGS>
struct arg_data<T, Char, NUM_ARGS, 0> { struct arg_data<T, Char, NUM_ARGS, 0> {
T args_[NUM_ARGS != 0 ? NUM_ARGS : 1]; // +1 to workaround a bug in gcc 7.5 that causes duplicated-branches warning.
T args_[NUM_ARGS != 0 ? NUM_ARGS : +1];
template <typename... U> template <typename... U>
FMT_INLINE arg_data(const U&... init) : args_{init...} {} FMT_INLINE arg_data(const U&... init) : args_{init...} {}
@ -967,6 +1128,8 @@ enum { long_short = sizeof(long) == sizeof(int) };
using long_type = conditional_t<long_short, int, long long>; using long_type = conditional_t<long_short, int, long long>;
using ulong_type = conditional_t<long_short, unsigned, unsigned long long>; using ulong_type = conditional_t<long_short, unsigned, unsigned long long>;
struct unformattable {};
// Maps formatting arguments to core types. // Maps formatting arguments to core types.
template <typename Context> struct arg_mapper { template <typename Context> struct arg_mapper {
using char_type = typename Context::char_type; using char_type = typename Context::char_type;
@ -1075,15 +1238,7 @@ template <typename Context> struct arg_mapper {
return map(val.value); return map(val.value);
} }
int map(...) { unformattable map(...) { return {}; }
constexpr bool formattable = sizeof(Context) == 0;
static_assert(
formattable,
"Cannot format argument. To make type T formattable provide a "
"formatter<T> specialization: "
"https://fmt.dev/latest/api.html#formatting-user-defined-types");
return 0;
}
}; };
// A type constant after applying arg_mapper<Context>. // A type constant after applying arg_mapper<Context>.
@ -1207,15 +1362,25 @@ FMT_CONSTEXPR_DECL FMT_INLINE auto visit_format_arg(
return vis(monostate()); return vis(monostate());
} }
// Checks whether T is a container with contiguous storage. template <typename T> struct formattable : std::false_type {};
template <typename T> struct is_contiguous : std::false_type {};
template <typename Char>
struct is_contiguous<std::basic_string<Char>> : std::true_type {};
template <typename Char>
struct is_contiguous<detail::buffer<Char>> : std::true_type {};
namespace detail { namespace detail {
// A workaround for gcc 4.8 to make void_t work in a SFINAE context.
template <typename... Ts> struct void_t_impl { using type = void; };
template <typename... Ts>
using void_t = typename detail::void_t_impl<Ts...>::type;
template <typename It, typename T, typename Enable = void>
struct is_output_iterator : std::false_type {};
template <typename It, typename T>
struct is_output_iterator<
It, T,
void_t<typename std::iterator_traits<It>::iterator_category,
decltype(*std::declval<It>() = std::declval<T>())>>
: std::true_type {};
template <typename OutputIt> template <typename OutputIt>
struct is_back_insert_iterator : std::false_type {}; struct is_back_insert_iterator : std::false_type {};
template <typename Container> template <typename Container>
@ -1227,6 +1392,9 @@ struct is_contiguous_back_insert_iterator : std::false_type {};
template <typename Container> template <typename Container>
struct is_contiguous_back_insert_iterator<std::back_insert_iterator<Container>> struct is_contiguous_back_insert_iterator<std::back_insert_iterator<Container>>
: is_contiguous<Container> {}; : is_contiguous<Container> {};
template <typename Char>
struct is_contiguous_back_insert_iterator<buffer_appender<Char>>
: std::true_type {};
// A type-erased reference to an std::locale to avoid heavy <locale> include. // A type-erased reference to an std::locale to avoid heavy <locale> include.
class locale_ref { class locale_ref {
@ -1258,13 +1426,24 @@ FMT_CONSTEXPR basic_format_arg<Context> make_arg(const T& value) {
return arg; return arg;
} }
template <typename T> int check(unformattable) {
static_assert(
formattable<T>(),
"Cannot format an argument. To make type T formattable provide a "
"formatter<T> specialization: https://fmt.dev/latest/api.html#udt");
return 0;
}
template <typename T, typename U> inline const U& check(const U& val) {
return val;
}
// The type template parameter is there to avoid an ODR violation when using // The type template parameter is there to avoid an ODR violation when using
// a fallback formatter in one translation unit and an implicit conversion in // a fallback formatter in one translation unit and an implicit conversion in
// another (not recommended). // another (not recommended).
template <bool IS_PACKED, typename Context, type, typename T, template <bool IS_PACKED, typename Context, type, typename T,
FMT_ENABLE_IF(IS_PACKED)> FMT_ENABLE_IF(IS_PACKED)>
inline value<Context> make_arg(const T& val) { inline value<Context> make_arg(const T& val) {
return arg_mapper<Context>().map(val); return check<T>(arg_mapper<Context>().map(val));
} }
template <bool IS_PACKED, typename Context, type, typename T, template <bool IS_PACKED, typename Context, type, typename T,
@ -1364,13 +1543,13 @@ template <typename OutputIt, typename Char> class basic_format_context {
template <typename Char> template <typename Char>
using buffer_context = using buffer_context =
basic_format_context<std::back_insert_iterator<detail::buffer<Char>>, Char>; basic_format_context<detail::buffer_appender<Char>, Char>;
using format_context = buffer_context<char>; using format_context = buffer_context<char>;
using wformat_context = buffer_context<wchar_t>; using wformat_context = buffer_context<wchar_t>;
// Workaround a bug in gcc: https://stackoverflow.com/q/62767544/471164. // Workaround an alias issue: https://stackoverflow.com/q/62767544/471164.
#define FMT_BUFFER_CONTEXT(Char) \ #define FMT_BUFFER_CONTEXT(Char) \
basic_format_context<std::back_insert_iterator<detail::buffer<Char>>, Char> basic_format_context<detail::buffer_appender<Char>, Char>
/** /**
\rst \rst
@ -1422,7 +1601,7 @@ class format_arg_store
/** /**
\rst \rst
Constructs an `~fmt::format_arg_store` object that contains references to Constructs a `~fmt::format_arg_store` object that contains references to
arguments and can be implicitly converted to `~fmt::format_args`. `Context` arguments and can be implicitly converted to `~fmt::format_args`. `Context`
can be omitted in which case it defaults to `~fmt::context`. can be omitted in which case it defaults to `~fmt::context`.
See `~fmt::arg` for lifetime considerations. See `~fmt::arg` for lifetime considerations.
@ -1434,6 +1613,27 @@ inline format_arg_store<Context, Args...> make_format_args(
return {args...}; return {args...};
} }
/**
\rst
Constructs a `~fmt::format_arg_store` object that contains references
to arguments and can be implicitly converted to `~fmt::format_args`.
If ``format_str`` is a compile-time string then `make_args_checked` checks
its validity at compile time.
\endrst
*/
template <typename... Args, typename S, typename Char = char_t<S>>
inline auto make_args_checked(const S& format_str,
const remove_reference_t<Args>&... args)
-> format_arg_store<buffer_context<Char>, remove_reference_t<Args>...> {
static_assert(
detail::count<(
std::is_base_of<detail::view, remove_reference_t<Args>>::value &&
std::is_reference<Args>::value)...>() == 0,
"passing views as lvalues is disallowed");
detail::check_format_string<Args...>(format_str);
return {args...};
}
/** /**
\rst \rst
Returns a named argument to be used in a formatting function. It should only Returns a named argument to be used in a formatting function. It should only
@ -1737,7 +1937,14 @@ template <typename Context> class basic_format_args {
} }
}; };
/** An alias to ``basic_format_args<context>``. */ #ifdef FMT_ARM_ABI_COMPATIBILITY
/** An alias to ``basic_format_args<format_context>``. */
// Separate types would result in shorter symbols but break ABI compatibility
// between clang and gcc on ARM (#1919).
using format_args = basic_format_args<format_context>;
using wformat_args = basic_format_args<wformat_context>;
#else
// DEPRECATED! These are kept for ABI compatibility.
// It is a separate type rather than an alias to make symbols readable. // It is a separate type rather than an alias to make symbols readable.
struct format_args : basic_format_args<format_context> { struct format_args : basic_format_args<format_context> {
template <typename... Args> template <typename... Args>
@ -1746,32 +1953,10 @@ struct format_args : basic_format_args<format_context> {
struct wformat_args : basic_format_args<wformat_context> { struct wformat_args : basic_format_args<wformat_context> {
using basic_format_args::basic_format_args; using basic_format_args::basic_format_args;
}; };
#endif
namespace detail { namespace detail {
// Reports a compile-time error if S is not a valid format string.
template <typename..., typename S, FMT_ENABLE_IF(!is_compile_string<S>::value)>
FMT_INLINE void check_format_string(const S&) {
#ifdef FMT_ENFORCE_COMPILE_STRING
static_assert(is_compile_string<S>::value,
"FMT_ENFORCE_COMPILE_STRING requires all format strings to use "
"FMT_STRING.");
#endif
}
template <typename..., typename S, FMT_ENABLE_IF(is_compile_string<S>::value)>
void check_format_string(S);
template <typename... Args, typename S, typename Char = char_t<S>>
inline format_arg_store<buffer_context<Char>, remove_reference_t<Args>...>
make_args_checked(const S& format_str,
const remove_reference_t<Args>&... args) {
static_assert(count<(std::is_base_of<view, remove_reference_t<Args>>::value &&
std::is_reference<Args>::value)...>() == 0,
"passing views as lvalues is disallowed");
check_format_string<Args...>(format_str);
return {args...};
}
template <typename Char, FMT_ENABLE_IF(!std::is_same<Char, char>::value)> template <typename Char, FMT_ENABLE_IF(!std::is_same<Char, char>::value)>
std::basic_string<Char> vformat( std::basic_string<Char> vformat(
basic_string_view<Char> format_str, basic_string_view<Char> format_str,
@ -1780,9 +1965,10 @@ std::basic_string<Char> vformat(
FMT_API std::string vformat(string_view format_str, format_args args); FMT_API std::string vformat(string_view format_str, format_args args);
template <typename Char> template <typename Char>
typename FMT_BUFFER_CONTEXT(Char)::iterator vformat_to( void vformat_to(
buffer<Char>& buf, basic_string_view<Char> format_str, buffer<Char>& buf, basic_string_view<Char> format_str,
basic_format_args<FMT_BUFFER_CONTEXT(type_identity_t<Char>)> args); basic_format_args<FMT_BUFFER_CONTEXT(type_identity_t<Char>)> args,
detail::locale_ref loc = {});
template <typename Char, typename Args, template <typename Char, typename Args,
FMT_ENABLE_IF(!std::is_same<Char, char>::value)> FMT_ENABLE_IF(!std::is_same<Char, char>::value)>
@ -1797,26 +1983,80 @@ inline void vprint_mojibake(std::FILE*, string_view, format_args) {}
/** Formats a string and writes the output to ``out``. */ /** Formats a string and writes the output to ``out``. */
// GCC 8 and earlier cannot handle std::back_insert_iterator<Container> with // GCC 8 and earlier cannot handle std::back_insert_iterator<Container> with
// vformat_to<ArgFormatter>(...) overload, so SFINAE on iterator type instead. // vformat_to<ArgFormatter>(...) overload, so SFINAE on iterator type instead.
template < template <typename OutputIt, typename S, typename Char = char_t<S>,
typename OutputIt, typename S, typename Char = char_t<S>, bool enable = detail::is_output_iterator<OutputIt, Char>::value>
FMT_ENABLE_IF(detail::is_contiguous_back_insert_iterator<OutputIt>::value)> auto vformat_to(OutputIt out, const S& format_str,
OutputIt vformat_to( basic_format_args<buffer_context<type_identity_t<Char>>> args)
OutputIt out, const S& format_str, -> typename std::enable_if<enable, OutputIt>::type {
basic_format_args<buffer_context<type_identity_t<Char>>> args) { decltype(detail::get_buffer<Char>(out)) buf(detail::get_buffer_init(out));
auto& c = detail::get_container(out);
detail::container_buffer<remove_reference_t<decltype(c)>> buf(c);
detail::vformat_to(buf, to_string_view(format_str), args); detail::vformat_to(buf, to_string_view(format_str), args);
return out; return detail::get_iterator(buf);
} }
template <typename Container, typename S, typename... Args, /**
FMT_ENABLE_IF( \rst
is_contiguous<Container>::value&& detail::is_string<S>::value)> Formats arguments, writes the result to the output iterator ``out`` and returns
inline std::back_insert_iterator<Container> format_to( the iterator past the end of the output range.
std::back_insert_iterator<Container> out, const S& format_str,
Args&&... args) { **Example**::
return vformat_to(out, to_string_view(format_str),
detail::make_args_checked<Args...>(format_str, args...)); std::vector<char> out;
fmt::format_to(std::back_inserter(out), "{}", 42);
\endrst
*/
// We cannot use FMT_ENABLE_IF because of a bug in gcc 8.3.
template <typename OutputIt, typename S, typename... Args,
bool enable = detail::is_output_iterator<OutputIt, char_t<S>>::value>
inline auto format_to(OutputIt out, const S& format_str, Args&&... args) ->
typename std::enable_if<enable, OutputIt>::type {
const auto& vargs = fmt::make_args_checked<Args...>(format_str, args...);
return vformat_to(out, to_string_view(format_str), vargs);
}
template <typename OutputIt> struct format_to_n_result {
/** Iterator past the end of the output range. */
OutputIt out;
/** Total (not truncated) output size. */
size_t size;
};
template <typename OutputIt, typename Char, typename... Args,
FMT_ENABLE_IF(detail::is_output_iterator<OutputIt, Char>::value)>
inline format_to_n_result<OutputIt> vformat_to_n(
OutputIt out, size_t n, basic_string_view<Char> format_str,
basic_format_args<buffer_context<type_identity_t<Char>>> args) {
detail::iterator_buffer<OutputIt, Char, detail::fixed_buffer_traits> buf(out,
n);
detail::vformat_to(buf, format_str, args);
return {buf.out(), buf.count()};
}
/**
\rst
Formats arguments, writes up to ``n`` characters of the result to the output
iterator ``out`` and returns the total output size and the iterator past the
end of the output range.
\endrst
*/
template <typename OutputIt, typename S, typename... Args,
bool enable = detail::is_output_iterator<OutputIt, char_t<S>>::value>
inline auto format_to_n(OutputIt out, size_t n, const S& format_str,
const Args&... args) ->
typename std::enable_if<enable, format_to_n_result<OutputIt>>::type {
const auto& vargs = fmt::make_args_checked<Args...>(format_str, args...);
return vformat_to_n(out, n, to_string_view(format_str), vargs);
}
/**
Returns the number of characters in the output of
``format(format_str, args...)``.
*/
template <typename... Args>
inline size_t formatted_size(string_view format_str, Args&&... args) {
const auto& vargs = fmt::make_args_checked<Args...>(format_str, args...);
detail::counting_buffer<> buf;
detail::vformat_to(buf, format_str, vargs);
return buf.count();
} }
template <typename S, typename Char = char_t<S>> template <typename S, typename Char = char_t<S>>
@ -1840,7 +2080,7 @@ FMT_INLINE std::basic_string<Char> vformat(
// std::basic_string<char_t<S>> to reduce the symbol size. // std::basic_string<char_t<S>> to reduce the symbol size.
template <typename S, typename... Args, typename Char = char_t<S>> template <typename S, typename... Args, typename Char = char_t<S>>
FMT_INLINE std::basic_string<Char> format(const S& format_str, Args&&... args) { FMT_INLINE std::basic_string<Char> format(const S& format_str, Args&&... args) {
const auto& vargs = detail::make_args_checked<Args...>(format_str, args...); const auto& vargs = fmt::make_args_checked<Args...>(format_str, args...);
return detail::vformat(to_string_view(format_str), vargs); return detail::vformat(to_string_view(format_str), vargs);
} }
@ -1860,7 +2100,7 @@ FMT_API void vprint(std::FILE*, string_view, format_args);
*/ */
template <typename S, typename... Args, typename Char = char_t<S>> template <typename S, typename... Args, typename Char = char_t<S>>
inline void print(std::FILE* f, const S& format_str, Args&&... args) { inline void print(std::FILE* f, const S& format_str, Args&&... args) {
const auto& vargs = detail::make_args_checked<Args...>(format_str, args...); const auto& vargs = fmt::make_args_checked<Args...>(format_str, args...);
return detail::is_unicode<Char>() return detail::is_unicode<Char>()
? vprint(f, to_string_view(format_str), vargs) ? vprint(f, to_string_view(format_str), vargs)
: detail::vprint_mojibake(f, to_string_view(format_str), vargs); : detail::vprint_mojibake(f, to_string_view(format_str), vargs);
@ -1879,7 +2119,7 @@ inline void print(std::FILE* f, const S& format_str, Args&&... args) {
*/ */
template <typename S, typename... Args, typename Char = char_t<S>> template <typename S, typename... Args, typename Char = char_t<S>>
inline void print(const S& format_str, Args&&... args) { inline void print(const S& format_str, Args&&... args) {
const auto& vargs = detail::make_args_checked<Args...>(format_str, args...); const auto& vargs = fmt::make_args_checked<Args...>(format_str, args...);
return detail::is_unicode<Char>() return detail::is_unicode<Char>()
? vprint(to_string_view(format_str), vargs) ? vprint(to_string_view(format_str), vargs)
: detail::vprint_mojibake(stdout, to_string_view(format_str), : detail::vprint_mojibake(stdout, to_string_view(format_str),

File diff suppressed because it is too large Load Diff

File diff suppressed because it is too large Load Diff

View File

@ -15,22 +15,12 @@
FMT_BEGIN_NAMESPACE FMT_BEGIN_NAMESPACE
namespace detail { namespace detail {
template <typename Char>
typename buffer_context<Char>::iterator vformat_to(
const std::locale& loc, buffer<Char>& buf,
basic_string_view<Char> format_str,
basic_format_args<buffer_context<type_identity_t<Char>>> args) {
using af = arg_formatter<typename buffer_context<Char>::iterator, Char>;
return vformat_to<af>(std::back_inserter(buf), to_string_view(format_str),
args, detail::locale_ref(loc));
}
template <typename Char> template <typename Char>
std::basic_string<Char> vformat( std::basic_string<Char> vformat(
const std::locale& loc, basic_string_view<Char> format_str, const std::locale& loc, basic_string_view<Char> format_str,
basic_format_args<buffer_context<type_identity_t<Char>>> args) { basic_format_args<buffer_context<type_identity_t<Char>>> args) {
basic_memory_buffer<Char> buffer; basic_memory_buffer<Char> buffer;
detail::vformat_to(loc, buffer, format_str, args); detail::vformat_to(buffer, format_str, args, detail::locale_ref(loc));
return fmt::to_string(buffer); return fmt::to_string(buffer);
} }
} // namespace detail } // namespace detail
@ -45,32 +35,28 @@ inline std::basic_string<Char> vformat(
template <typename S, typename... Args, typename Char = char_t<S>> template <typename S, typename... Args, typename Char = char_t<S>>
inline std::basic_string<Char> format(const std::locale& loc, inline std::basic_string<Char> format(const std::locale& loc,
const S& format_str, Args&&... args) { const S& format_str, Args&&... args) {
return detail::vformat( return detail::vformat(loc, to_string_view(format_str),
loc, to_string_view(format_str), fmt::make_args_checked<Args...>(format_str, args...));
detail::make_args_checked<Args...>(format_str, args...));
} }
template <typename S, typename OutputIt, typename... Args, template <typename S, typename OutputIt, typename... Args,
typename Char = enable_if_t< typename Char = char_t<S>,
detail::is_output_iterator<OutputIt>::value, char_t<S>>> FMT_ENABLE_IF(detail::is_output_iterator<OutputIt, Char>::value)>
inline OutputIt vformat_to( inline OutputIt vformat_to(
OutputIt out, const std::locale& loc, const S& format_str, OutputIt out, const std::locale& loc, const S& format_str,
format_args_t<type_identity_t<OutputIt>, Char> args) { basic_format_args<buffer_context<type_identity_t<Char>>> args) {
using af = detail::arg_formatter<OutputIt, Char>; decltype(detail::get_buffer<Char>(out)) buf(detail::get_buffer_init(out));
return vformat_to<af>(out, to_string_view(format_str), args, vformat_to(buf, to_string_view(format_str), args, detail::locale_ref(loc));
detail::locale_ref(loc)); return detail::get_iterator(buf);
} }
template <typename OutputIt, typename S, typename... Args, template <typename OutputIt, typename S, typename... Args,
FMT_ENABLE_IF(detail::is_output_iterator<OutputIt>::value&& bool enable = detail::is_output_iterator<OutputIt, char_t<S>>::value>
detail::is_string<S>::value)> inline auto format_to(OutputIt out, const std::locale& loc,
inline OutputIt format_to(OutputIt out, const std::locale& loc, const S& format_str, Args&&... args) ->
const S& format_str, Args&&... args) { typename std::enable_if<enable, OutputIt>::type {
detail::check_format_string<Args...>(format_str); const auto& vargs = fmt::make_args_checked<Args...>(format_str, args...);
using context = format_context_t<OutputIt, char_t<S>>; return vformat_to(out, loc, to_string_view(format_str), vargs);
format_arg_store<context, Args...> as{args...};
return vformat_to(out, loc, to_string_view(format_str),
basic_format_args<context>(as));
} }
FMT_END_NAMESPACE FMT_END_NAMESPACE

View File

@ -29,7 +29,8 @@
#if FMT_HAS_INCLUDE("winapifamily.h") #if FMT_HAS_INCLUDE("winapifamily.h")
# include <winapifamily.h> # include <winapifamily.h>
#endif #endif
#if FMT_HAS_INCLUDE("fcntl.h") && \ #if (FMT_HAS_INCLUDE(<fcntl.h>) || defined(__APPLE__) || \
defined(__linux__)) && \
(!defined(WINAPI_FAMILY) || (WINAPI_FAMILY == WINAPI_FAMILY_DESKTOP_APP)) (!defined(WINAPI_FAMILY) || (WINAPI_FAMILY == WINAPI_FAMILY_DESKTOP_APP))
# include <fcntl.h> // for O_RDONLY # include <fcntl.h> // for O_RDONLY
# define FMT_USE_FCNTL 1 # define FMT_USE_FCNTL 1
@ -278,7 +279,8 @@ class file {
RDONLY = FMT_POSIX(O_RDONLY), // Open for reading only. RDONLY = FMT_POSIX(O_RDONLY), // Open for reading only.
WRONLY = FMT_POSIX(O_WRONLY), // Open for writing only. WRONLY = FMT_POSIX(O_WRONLY), // Open for writing only.
RDWR = FMT_POSIX(O_RDWR), // Open for reading and writing. RDWR = FMT_POSIX(O_RDWR), // Open for reading and writing.
CREATE = FMT_POSIX(O_CREAT) // Create if the file doesn't exist. CREATE = FMT_POSIX(O_CREAT), // Create if the file doesn't exist.
APPEND = FMT_POSIX(O_APPEND) // Open in append mode.
}; };
// Constructs a file object which doesn't represent any file. // Constructs a file object which doesn't represent any file.
@ -343,36 +345,69 @@ class file {
// Returns the memory page size. // Returns the memory page size.
long getpagesize(); long getpagesize();
class direct_buffered_file; namespace detail {
template <typename S, typename... Args> struct buffer_size {
void print(direct_buffered_file& f, const S& format_str, size_t value = 0;
const Args&... args); buffer_size operator=(size_t val) const {
auto bs = buffer_size();
bs.value = val;
return bs;
}
};
// A buffered file with a direct buffer access and no synchronization. struct ostream_params {
class direct_buffered_file { int oflag = file::WRONLY | file::CREATE;
size_t buffer_size = BUFSIZ > 32768 ? BUFSIZ : 32768;
ostream_params() {}
template <typename... T>
ostream_params(T... params, int oflag) : ostream_params(params...) {
this->oflag = oflag;
}
template <typename... T>
ostream_params(T... params, detail::buffer_size bs)
: ostream_params(params...) {
this->buffer_size = bs.value;
}
};
} // namespace detail
static constexpr detail::buffer_size buffer_size;
// A fast output stream which is not thread-safe.
class ostream final : private detail::buffer<char> {
private: private:
file file_; file file_;
enum { buffer_size = 4096 };
char buffer_[buffer_size];
int pos_;
void flush() { void flush() {
if (pos_ == 0) return; if (size() == 0) return;
file_.write(buffer_, pos_); file_.write(data(), size());
pos_ = 0; clear();
} }
int free_capacity() const { return buffer_size - pos_; } void grow(size_t) final;
ostream(cstring_view path, const detail::ostream_params& params)
: file_(path, params.oflag) {
set(new char[params.buffer_size], params.buffer_size);
}
public: public:
direct_buffered_file(cstring_view path, int oflag) ostream(ostream&& other)
: file_(path, oflag), pos_(0) {} : detail::buffer<char>(other.data(), other.size(), other.capacity()),
file_(std::move(other.file_)) {
~direct_buffered_file() { other.set(nullptr, 0);
flush();
} }
~ostream() {
flush();
delete[] data();
}
template <typename... T>
friend ostream output_file(cstring_view path, T... params);
void close() { void close() {
flush(); flush();
@ -380,25 +415,20 @@ class direct_buffered_file {
} }
template <typename S, typename... Args> template <typename S, typename... Args>
friend void print(direct_buffered_file& f, const S& format_str, void print(const S& format_str, const Args&... args) {
const Args&... args) { format_to(detail::buffer_appender<char>(*this), format_str, args...);
// We could avoid double buffering.
auto buf = fmt::memory_buffer();
fmt::format_to(std::back_inserter(buf), format_str, args...);
auto remaining_pos = 0;
auto remaining_size = buf.size();
while (remaining_size > detail::to_unsigned(f.free_capacity())) {
auto size = f.free_capacity();
memcpy(f.buffer_ + f.pos_, buf.data() + remaining_pos, size);
f.pos_ += size;
f.flush();
remaining_pos += size;
remaining_size -= size;
}
memcpy(f.buffer_ + f.pos_, buf.data() + remaining_pos, remaining_size);
f.pos_ += static_cast<int>(remaining_size);
} }
}; };
/**
Opens a file for writing. Supported parameters passed in `params`:
* ``<integer>``: Output flags (``file::WRONLY | file::CREATE`` by default)
* ``buffer_size=<integer>``: Output buffer size
*/
template <typename... T>
inline ostream output_file(cstring_view path, T... params) {
return {path, detail::ostream_params(params...)};
}
#endif // FMT_USE_FCNTL #endif // FMT_USE_FCNTL
#ifdef FMT_LOCALE #ifdef FMT_LOCALE

View File

@ -49,17 +49,27 @@ template <class Char> class formatbuf : public std::basic_streambuf<Char> {
} }
}; };
struct converter {
template <typename T, FMT_ENABLE_IF(is_integral<T>::value)> converter(T);
};
template <typename Char> struct test_stream : std::basic_ostream<Char> { template <typename Char> struct test_stream : std::basic_ostream<Char> {
private: private:
// Hide all operator<< from std::basic_ostream<Char>. void_t<> operator<<(converter);
void_t<> operator<<(null<>);
void_t<> operator<<(const Char*);
template <typename T, FMT_ENABLE_IF(std::is_convertible<T, int>::value &&
!std::is_enum<T>::value)>
void_t<> operator<<(T);
}; };
// Hide insertion operators for built-in types.
template <typename Char, typename Traits>
void_t<> operator<<(std::basic_ostream<Char, Traits>&, Char);
template <typename Char, typename Traits>
void_t<> operator<<(std::basic_ostream<Char, Traits>&, char);
template <typename Traits>
void_t<> operator<<(std::basic_ostream<char, Traits>&, char);
template <typename Traits>
void_t<> operator<<(std::basic_ostream<char, Traits>&, signed char);
template <typename Traits>
void_t<> operator<<(std::basic_ostream<char, Traits>&, unsigned char);
// Checks if T has a user-defined operator<< (e.g. not a member of // Checks if T has a user-defined operator<< (e.g. not a member of
// std::ostream). // std::ostream).
template <typename T, typename Char> class is_streamable { template <typename T, typename Char> class is_streamable {
@ -103,7 +113,7 @@ void format_value(buffer<Char>& buf, const T& value,
#endif #endif
output << value; output << value;
output.exceptions(std::ios_base::failbit | std::ios_base::badbit); output.exceptions(std::ios_base::failbit | std::ios_base::badbit);
buf.resize(buf.size()); buf.try_resize(buf.size());
} }
// Formats an object of type T that has an overloaded ostream operator<<. // Formats an object of type T that has an overloaded ostream operator<<.
@ -160,7 +170,7 @@ template <typename S, typename... Args,
typename Char = enable_if_t<detail::is_string<S>::value, char_t<S>>> typename Char = enable_if_t<detail::is_string<S>::value, char_t<S>>>
void print(std::basic_ostream<Char>& os, const S& format_str, Args&&... args) { void print(std::basic_ostream<Char>& os, const S& format_str, Args&&... args) {
vprint(os, to_string_view(format_str), vprint(os, to_string_view(format_str),
detail::make_args_checked<Args...>(format_str, args...)); fmt::make_args_checked<Args...>(format_str, args...));
} }
FMT_END_NAMESPACE FMT_END_NAMESPACE

View File

@ -181,7 +181,7 @@ template <typename Char> class printf_width_handler {
template <typename Char, typename Context> template <typename Char, typename Context>
void vprintf(buffer<Char>& buf, basic_string_view<Char> format, void vprintf(buffer<Char>& buf, basic_string_view<Char> format,
basic_format_args<Context> args) { basic_format_args<Context> args) {
Context(std::back_inserter(buf), format, args).format(); Context(buffer_appender<Char>(buf), format, args).format();
} }
} // namespace detail } // namespace detail
@ -598,7 +598,7 @@ OutputIt basic_printf_context<OutputIt, Char>::format() {
template <typename Char> template <typename Char>
using basic_printf_context_t = using basic_printf_context_t =
basic_printf_context<std::back_insert_iterator<detail::buffer<Char>>, Char>; basic_printf_context<detail::buffer_appender<Char>, Char>;
using printf_context = basic_printf_context_t<char>; using printf_context = basic_printf_context_t<char>;
using wprintf_context = basic_printf_context_t<wchar_t>; using wprintf_context = basic_printf_context_t<wchar_t>;

View File

@ -157,6 +157,9 @@ template <class Tuple, class F> void for_each(Tuple&& tup, F&& f) {
for_each(indexes, std::forward<Tuple>(tup), std::forward<F>(f)); for_each(indexes, std::forward<Tuple>(tup), std::forward<F>(f));
} }
template <typename Range>
using value_type = remove_cvref_t<decltype(*std::declval<Range>().begin())>;
template <typename Arg, FMT_ENABLE_IF(!is_like_std_string< template <typename Arg, FMT_ENABLE_IF(!is_like_std_string<
typename std::decay<Arg>::type>::value)> typename std::decay<Arg>::type>::value)>
FMT_CONSTEXPR const char* format_str_quoted(bool add_space, const Arg&) { FMT_CONSTEXPR const char* format_str_quoted(bool add_space, const Arg&) {
@ -182,7 +185,6 @@ FMT_CONSTEXPR const char* format_str_quoted(bool add_space, const char) {
FMT_CONSTEXPR const wchar_t* format_str_quoted(bool add_space, const wchar_t) { FMT_CONSTEXPR const wchar_t* format_str_quoted(bool add_space, const wchar_t) {
return add_space ? L" '{}'" : L"'{}'"; return add_space ? L" '{}'" : L"'{}'";
} }
} // namespace detail } // namespace detail
template <typename T> struct is_tuple_like { template <typename T> struct is_tuple_like {
@ -246,9 +248,15 @@ template <typename T, typename Char> struct is_range {
!std::is_constructible<detail::std_string_view<Char>, T>::value; !std::is_constructible<detail::std_string_view<Char>, T>::value;
}; };
template <typename RangeT, typename Char> template <typename T, typename Char>
struct formatter<RangeT, Char, struct formatter<
enable_if_t<fmt::is_range<RangeT, Char>::value>> { T, Char,
enable_if_t<fmt::is_range<T, Char>::value
// Workaround a bug in MSVC 2017 and earlier.
#if !FMT_MSC_VER || FMT_MSC_VER >= 1927
&& has_formatter<detail::value_type<T>, format_context>::value
#endif
>> {
formatting_range<Char> formatting; formatting_range<Char> formatting;
template <typename ParseContext> template <typename ParseContext>
@ -257,8 +265,7 @@ struct formatter<RangeT, Char,
} }
template <typename FormatContext> template <typename FormatContext>
typename FormatContext::iterator format(const RangeT& values, typename FormatContext::iterator format(const T& values, FormatContext& ctx) {
FormatContext& ctx) {
auto out = detail::copy(formatting.prefix, ctx.out()); auto out = detail::copy(formatting.prefix, ctx.out());
size_t i = 0; size_t i = 0;
auto it = values.begin(); auto it = values.begin();

View File

@ -23,6 +23,36 @@ int format_float(char* buf, std::size_t size, const char* format, int precision,
return precision < 0 ? snprintf_ptr(buf, size, format, value) return precision < 0 ? snprintf_ptr(buf, size, format, value)
: snprintf_ptr(buf, size, format, precision, value); : snprintf_ptr(buf, size, format, precision, value);
} }
template dragonbox::decimal_fp<float> dragonbox::to_decimal(float x)
FMT_NOEXCEPT;
template dragonbox::decimal_fp<double> dragonbox::to_decimal(double x)
FMT_NOEXCEPT;
// DEPRECATED! This function exists for ABI compatibility.
template <typename Char>
typename basic_format_context<std::back_insert_iterator<buffer<Char>>,
Char>::iterator
vformat_to(buffer<Char>& buf, basic_string_view<Char> format_str,
basic_format_args<basic_format_context<
std::back_insert_iterator<buffer<type_identity_t<Char>>>,
type_identity_t<Char>>>
args) {
using iterator = std::back_insert_iterator<buffer<char>>;
using context = basic_format_context<
std::back_insert_iterator<buffer<type_identity_t<Char>>>,
type_identity_t<Char>>;
auto out = iterator(buf);
format_handler<iterator, Char, context> h(out, format_str, args, {});
parse_format_string<false>(format_str, h);
return out;
}
template basic_format_context<std::back_insert_iterator<buffer<char>>,
char>::iterator
vformat_to(buffer<char>&, string_view,
basic_format_args<basic_format_context<
std::back_insert_iterator<buffer<type_identity_t<char>>>,
type_identity_t<char>>>);
} // namespace detail } // namespace detail
template struct FMT_INSTANTIATION_DEF_API detail::basic_data<void>; template struct FMT_INSTANTIATION_DEF_API detail::basic_data<void>;
@ -44,9 +74,9 @@ template FMT_API char detail::decimal_point_impl(locale_ref);
template FMT_API void detail::buffer<char>::append(const char*, const char*); template FMT_API void detail::buffer<char>::append(const char*, const char*);
template FMT_API FMT_BUFFER_CONTEXT(char)::iterator detail::vformat_to( template FMT_API void detail::vformat_to(
detail::buffer<char>&, string_view, detail::buffer<char>&, string_view,
basic_format_args<FMT_BUFFER_CONTEXT(char)>); basic_format_args<FMT_BUFFER_CONTEXT(char)>, detail::locale_ref);
template FMT_API int detail::snprintf_float(double, int, detail::float_specs, template FMT_API int detail::snprintf_float(double, int, detail::float_specs,
detail::buffer<char>&); detail::buffer<char>&);

View File

@ -62,7 +62,7 @@ using RWResult = int;
inline unsigned convert_rwcount(std::size_t count) { inline unsigned convert_rwcount(std::size_t count) {
return count <= UINT_MAX ? static_cast<unsigned>(count) : UINT_MAX; return count <= UINT_MAX ? static_cast<unsigned>(count) : UINT_MAX;
} }
#else #elif FMT_USE_FCNTL
// Return type of read and write functions. // Return type of read and write functions.
using RWResult = ssize_t; using RWResult = ssize_t;
@ -124,7 +124,8 @@ void detail::format_windows_error(detail::buffer<char>& out, int error_code,
if (result != 0) { if (result != 0) {
utf16_to_utf8 utf8_message; utf16_to_utf8 utf8_message;
if (utf8_message.convert(system_message) == ERROR_SUCCESS) { if (utf8_message.convert(system_message) == ERROR_SUCCESS) {
format_to(std::back_inserter(out), "{}: {}", message, utf8_message); format_to(buffer_appender<char>(out), "{}: {}", message,
utf8_message);
return; return;
} }
break; break;
@ -288,12 +289,12 @@ void file::pipe(file& read_end, file& write_end) {
} }
buffered_file file::fdopen(const char* mode) { buffered_file file::fdopen(const char* mode) {
// Don't retry as fdopen doesn't return EINTR. // Don't retry as fdopen doesn't return EINTR.
#if defined(__MINGW32__) && defined(_POSIX_) # if defined(__MINGW32__) && defined(_POSIX_)
FILE* f = ::fdopen(fd_, mode); FILE* f = ::fdopen(fd_, mode);
#else # else
FILE* f = FMT_POSIX_CALL(fdopen(fd_, mode)); FILE* f = FMT_POSIX_CALL(fdopen(fd_, mode));
#endif # endif
if (!f) if (!f)
FMT_THROW( FMT_THROW(
system_error(errno, "cannot associate stream with file descriptor")); system_error(errno, "cannot associate stream with file descriptor"));
@ -313,5 +314,9 @@ long getpagesize() {
return size; return size;
# endif # endif
} }
void ostream::grow(size_t) {
if (this->size() == this->capacity()) flush();
}
#endif // FMT_USE_FCNTL #endif // FMT_USE_FCNTL
FMT_END_NAMESPACE FMT_END_NAMESPACE

View File

@ -23,7 +23,7 @@ using namespace LAMMPS_NS;
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
Integrate::Integrate(LAMMPS *lmp, int /*narg*/, char **/*arg*/) : Pointers(lmp) Integrate::Integrate(LAMMPS *lmp, int /*narg*/, char ** /*arg*/) : Pointers(lmp)
{ {
elist_global = elist_atom = nullptr; elist_global = elist_atom = nullptr;
vlist_global = vlist_atom = cvlist_atom = nullptr; vlist_global = vlist_atom = cvlist_atom = nullptr;

View File

@ -1976,7 +1976,7 @@ void lammps_gather_atoms(void *handle, char *name, int type, int count, void *da
// use atom ID to insert each atom's values into copy // use atom ID to insert each atom's values into copy
// MPI_Allreduce with MPI_SUM to merge into data, ordered by atom ID // MPI_Allreduce with MPI_SUM to merge into data, ordered by atom ID
if (type == LAMMPS_INT) { if (type == 0) {
int *vector = nullptr; int *vector = nullptr;
int **array = nullptr; int **array = nullptr;
const int imgunpack = (count == 3) && (strcmp(name,"image") == 0); const int imgunpack = (count == 3) && (strcmp(name,"image") == 0);
@ -2015,7 +2015,7 @@ void lammps_gather_atoms(void *handle, char *name, int type, int count, void *da
MPI_Allreduce(copy,data,count*natoms,MPI_INT,MPI_SUM,lmp->world); MPI_Allreduce(copy,data,count*natoms,MPI_INT,MPI_SUM,lmp->world);
lmp->memory->destroy(copy); lmp->memory->destroy(copy);
} else if (type == LAMMPS_DOUBLE) { } else if (type == 1) {
double *vector = nullptr; double *vector = nullptr;
double **array = nullptr; double **array = nullptr;
if (count == 1) vector = (double *) vptr; if (count == 1) vector = (double *) vptr;

View File

@ -998,6 +998,23 @@ void Modify::replace_fix(const char *replaceID,
add_fix(narg,arg,trysuffix); add_fix(narg,arg,trysuffix);
} }
/* ----------------------------------------------------------------------
convenience function to allow replacing a fix from a single string
------------------------------------------------------------------------- */
void Modify::replace_fix(const std::string &oldfix,
const std::string &fixcmd, int trysuffix)
{
auto args = utils::split_words(fixcmd);
char **newarg = new char*[args.size()];
int i=0;
for (const auto &arg : args) {
newarg[i++] = (char *)arg.c_str();
}
replace_fix(oldfix.c_str(),args.size(),newarg,trysuffix);
delete[] newarg;
}
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
one instance per fix in style_fix.h one instance per fix in style_fix.h
------------------------------------------------------------------------- */ ------------------------------------------------------------------------- */

View File

@ -101,6 +101,7 @@ class Modify : protected Pointers {
void add_fix(int, char **, int trysuffix=1); void add_fix(int, char **, int trysuffix=1);
void add_fix(const std::string &, int trysuffix=1); void add_fix(const std::string &, int trysuffix=1);
void replace_fix(const char *, int, char **, int trysuffix=1); void replace_fix(const char *, int, char **, int trysuffix=1);
void replace_fix(const std::string &, const std::string &, int trysuffix=1);
void modify_fix(int, char **); void modify_fix(int, char **);
void delete_fix(const std::string &); void delete_fix(const std::string &);
void delete_fix(int); void delete_fix(int);

View File

@ -323,7 +323,12 @@ void Update::create_integrate(int narg, char **arg, int trysuffix)
delete integrate; delete integrate;
int sflag; int sflag;
if(narg-1 > 0) {
new_integrate(arg[0],narg-1,&arg[1],trysuffix,sflag); new_integrate(arg[0],narg-1,&arg[1],trysuffix,sflag);
} else {
new_integrate(arg[0],0,nullptr,trysuffix,sflag);
}
std::string estyle = arg[0]; std::string estyle = arg[0];
if (sflag) { if (sflag) {

View File

@ -862,6 +862,7 @@ std::string utils::get_potential_date(const std::string &path, const std::string
reader.ignore_comments = false; reader.ignore_comments = false;
char *line = reader.next_line(); char *line = reader.next_line();
if (line == nullptr) return "";
Tokenizer words(line); Tokenizer words(line);
while (words.has_next()) { while (words.has_next()) {
if (words.next() == "DATE:") { if (words.next() == "DATE:") {
@ -881,6 +882,7 @@ std::string utils::get_potential_units(const std::string &path, const std::strin
reader.ignore_comments = false; reader.ignore_comments = false;
char *line = reader.next_line(); char *line = reader.next_line();
if (line == nullptr) return "";
Tokenizer words(line); Tokenizer words(line);
while (words.has_next()) { while (words.has_next()) {
if (words.next() == "UNITS:") { if (words.next() == "UNITS:") {

View File

@ -1,8 +1,8 @@
find_package(YAML) find_package(YAML)
if(NOT YAML_FOUND) if(NOT YAML_FOUND)
message(STATUS "Skipping tests because libyaml is not found") # download and build a local copy of libyaml
return() include(YAML)
endif() endif()
if(CMAKE_VERSION VERSION_LESS 3.12) if(CMAKE_VERSION VERSION_LESS 3.12)

View File

@ -26,7 +26,7 @@ TEST(ErrorStats, test)
std::stringstream out; std::stringstream out;
out << stats; out << stats;
ASSERT_EQ(out.str(), "Average: 5.800e-01 StdDev: 7.305e-01 MaxErr: 2.000e+00 @ item: 3.0"); ASSERT_EQ(out.str(), "Average: 5.800e-01 StdDev: 7.305e-01 MaxErr: 2.000e+00 @ item: 3");
stats.reset(); stats.reset();
ASSERT_EQ(stats.has_data(), false); ASSERT_EQ(stats.has_data(), false);