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

Merge branch 'master' into reset-invalid-image-flags

Browse Source
This commit is contained in:
Axel Kohlmeyer
2020-11-18 10:36:54 -05:00
parent 20159fff23 64910d636c
commit 898ffc7d80
46 changed files with 1335 additions and 382 deletions
Download Patch File Download Diff File Expand all files Collapse all files

2
.github/CONTRIBUTING.md vendored
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.
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.
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.

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)

2
doc/github-development-workflow.md
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
rationale behind choices made. Exceptions to this policy are technical
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

6
doc/src/Build_development.rst
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
of the LAMMPS source code distribution itself. It can be enabled by
setting ``-D ENABLE_TESTING=on`` during the CMake configuration step.
It requires the `PyYAML <http://pyyaml.org/>`_ library and development
headers to compile and will download and compile a recent version of the
It requires the `YAML <http://pyyaml.org/>`_ library and development
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
for implementing the tests.

1
doc/src/Commands_pair.rst
View File

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

1
doc/src/Errors_warnings.rst
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>`.
*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
comprising a bond, angle, dihedral or improper must be present on any
sub-domain that "owns" the atom with the information, either as a

4
doc/src/Howto_github.rst
View File

@ -72,7 +72,7 @@ explained in more detail here: `feature branch workflow <https://www.atlassian.c
**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:
.. code-block:: bash
@ -155,7 +155,7 @@ useful message that explains the change.
.. 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:

2
doc/src/Howto_walls.rst
View File

@ -67,5 +67,5 @@ rotate.
The only frictional idealized walls currently in LAMMPS are flat or
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.

2
doc/src/Packages_details.rst
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/dmi <pair_spin_dmi>`
* :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/neel <pair_spin_neel>`
* :doc:`fix nve/spin <fix_nve_spin>`
* :doc:`fix langevin/spin <fix_langevin_spin>`
* :doc:`fix precession/spin <fix_precession_spin>`
* :doc:`compute spin <compute_spin>`
* :doc:`neb/spin <neb_spin>`

2
doc/src/atc_output.rst
View File

@ -14,7 +14,7 @@ Syntax
* AtC fixID = ID of :doc:`fix atc <fix_atc>` instance
* *output* or *output index* = name of the AtC sub-command
* 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*:
- text = creates text output of index, step and nodal variable values for unique nodes

2
doc/src/fix_filter_corotate.rst
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
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.
----------

2
doc/src/fix_precession_spin.rst
View File

@ -62,7 +62,7 @@ with:
The field value in Tesla is multiplied by the gyromagnetic
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
single spin (of norm :math:`\mu_i = 1.0`) submitted to an external

4
doc/src/fix_rx.rst
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.
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
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
and run-time per ODE as well as the average/RMS/min/max per process. If the

161
doc/src/pair_spin_exchange.rst
View File

@ -1,14 +1,19 @@
.. index:: pair_style spin/exchange
.. index:: pair_style spin/exchange/biquadratic
pair_style spin/exchange command
================================
pair_style spin/exchange/biquadratic command
============================================
Syntax
""""""
.. code-block:: LAMMPS
pair_style spin/exchange cutoff
pair_style spin/exchange/biquadratic cutoff
* cutoff = global cutoff pair (distance in metal units)
@ -19,7 +24,11 @@ Examples
pair_style spin/exchange 4.0
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
"""""""""""
@ -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
where :math:`\vec{s}_i` and :math:`\vec{s}_j` are two neighboring magnetic spins of two particles,
:math:`r_{ij} = \vert \vec{r}_i - \vec{r}_j \vert` is the inter-atomic distance between the two
particles. The summation is over pairs of nearest neighbors.
:math:`J(r_{ij})` is a function defining the intensity and the sign of the exchange
interaction for different neighboring shells. This function is defined as:
where :math:`\vec{s}_i` and :math:`\vec{s}_j` are two unit vectors representing
the magnetic spins of two particles (usually atoms), and
:math:`r_{ij} = \vert \vec{r}_i - \vec{r}_j \vert` is the inter-atomic distance
between those two particles. The summation is over pairs of nearest neighbors.
: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::
{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
"pair_coeff" command, and :math:`R_c` is the radius cutoff associated to
the pair interaction (see below for more explanations).
where :math:`\vec{s}_i`, :math:`\vec{s}_j`, :math:`r_{ij}` and
:math:`J(r_{ij})` have the same definitions as above, and :math:`K(r_{ij})` is
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 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>`.
The interatomic dependence of :math:`J(r_{ij})` and :math:`K(r_{ij})` in both
interactions above is defined by the following function:
.. 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
to a magnetic torque :math:`\vec{\omega}`, and its associated atom can be submitted to a
force :math:`\vec{F}` for spin-lattice calculations (see :doc:`fix nve/spin <fix_nve_spin>`),
such as:
to a magnetic torque :math:`\vec{\omega}_{i}`, and its associated atom can be
submitted to a force :math:`\vec{F}_{i}` for spin-lattice calculations (see
:doc:`fix nve/spin <fix_nve_spin>`), such as:
.. math::
\vec{\omega}_{i} = \frac{1}{\hbar} \sum_{j}^{Neighb} {J}
\left(r_{ij} \right)\,\vec{s}_{j}
~~{\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`.
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
:ref:`(Tranchida) <Tranchida3>`.
For the *spin/exchange* pair style, the following coefficients must be defined
for each pair of atoms types via the :doc:`pair_coeff <pair_coeff>` command as in
the examples above, or in the data file or restart files read by the
:doc:`read_data <read_data>` or :doc:`read_restart <read_restart>` commands, and
set in the following order:
For the *spin/exchange* and *spin/exchange/biquadratic* pair styles, the
following coefficients must be defined for each pair of atoms types via the
:doc:`pair_coeff <pair_coeff>` command as in the examples above, or in the data
file or restart files read by the :doc:`read_data <read_data>` or
:doc:`read_restart <read_restart>` commands, and set in the following order:
* :math:`R_c` (distance units)
* :math:`a` (energy units)
* :math:`b` (adim parameter)
* :math:`d` (distance units)
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.
for the *spin/exchange* pair style, and:
* :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
*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
""""""""""""
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
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
""""""""""""""""
@ -105,7 +208,7 @@ Default
"""""""
none
The default *offset* keyword value is *no*.
----------

1
doc/src/pair_style.rst
View File

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

2
doc/utils/sphinx-config/false_positives.txt
View File

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

2
examples/SPIN/cobalt_hcp/in.spin.cobalt_hcp
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
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
neighbor 0.1 bin

2
examples/SPIN/iron/in.spin.iron
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_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
neigh_modify every 10 check yes delay 20

4
examples/SPIN/iron/in.spin.iron_cubic
View File

@ -21,9 +21,9 @@ mass 1 55.845
set group all spin 2.2 -1.0 0.0 0.0
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 * * 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
neigh_modify every 10 check yes delay 20

20
examples/SPIN/test_problems/validation_damped_exchange/llg_exchange.py
View File

@ -6,9 +6,17 @@ import matplotlib.pyplot as plt
import mpmath as mp
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])
S2 = np.array([0.0, 1.0, 0.0])
alpha=0.01 # damping coefficient
pi=math.pi
@ -30,6 +38,14 @@ def rotation_matrix(axis, theta):
[2 * (bc - ad), aa + cc - bb - dd, 2 * (cd + ab)],
[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
def calc_rot_vector(Sr,Sf):
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
Sm = (S1+S2)*0.5
# 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(t*dt/1000.0,Sm[0],Sm[1],Sm[2],en)

4
examples/SPIN/test_problems/validation_damped_exchange/run-test-exchange.sh
View File

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

23
examples/SPIN/test_problems/validation_damped_exchange/test-spin-precession.in
View File

@ -5,22 +5,24 @@ atom_style spin
atom_modify map array
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_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 RUN equal 30000
fix 1 all nve/spin lattice no
fix 2 all precession/spin zeeman ${H} 0.0 0.0 1.0 anisotropy ${Kan} 0.0 0.0 1.0
fix_modify 2 energy yes
fix 3 all langevin/spin ${Temperature} 0.01 12345
fix 1 all nve/spin lattice frozen
fix 2 all langevin/spin ${Temperature} 0.01 12345
compute out_mag all spin
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 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
run ${RUN}

22
examples/SPIN/test_problems/validation_damped_exchange/two_spins.data
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

2
examples/SPIN/test_problems/validation_nve/run-test-nve.sh
View File

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

2
examples/SPIN/test_problems/validation_nvt/in.spin.nvt_lattice
View File

@ -30,7 +30,7 @@ neighbor 0.1 bin
neigh_modify every 10 check yes delay 20
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 4 all nve/spin lattice moving
timestep 0.001

2
examples/SPIN/test_problems/validation_nvt/in.spin.nvt_spin
View File

@ -29,7 +29,7 @@ neighbor 0.1 bin
neigh_modify every 10 check yes delay 20
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
timestep 0.001

2
examples/SPIN/test_problems/validation_nvt/plot_nvt.py
View File

@ -39,5 +39,5 @@ plt.xlabel('Time (in ps)')
plt.legend()
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)

4
lib/kokkos/README.md
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.
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
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.
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
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.
Which means contributing to Kokkos allows anyone else to use your contributions

2
lib/quip/README
View File

@ -17,7 +17,7 @@ Building LAMMPS with QUIP support:
1) Building 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

2
lib/scafacos/README
View File

@ -3,7 +3,7 @@ is required to use the KSPACE scafacos and its kspace_style
scafacos command in a LAMMPS input script.
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
Germany (Bonn, Chemnitz, Stuttgart, Wuppertal) and
the Research Centre Juelich (Juelich Supercomputing Centre).

9
src/KOKKOS/fft3d_kokkos.cpp
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);
#else
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;
if (flag == -1)
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);
Kokkos::parallel_for(total/length,f);
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
@ -273,13 +272,13 @@ void FFT3dKokkos<DeviceType>::fft_3d_kokkos(typename FFT_AT::t_FFT_DATA_1d d_in,
#elif defined(FFT_CUFFT)
cufftExec(plan->plan_mid,d_data.data(),d_data.data(),flag);
#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)
f = kiss_fft_functor<DeviceType>(d_data,d_tmp,plan->cfg_mid_forward,length);
else
f = kiss_fft_functor<DeviceType>(d_data,d_tmp,plan->cfg_mid_backward,length);
Kokkos::parallel_for(total/length,f);
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
// 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)
cufftExec(plan->plan_slow,d_data.data(),d_data.data(),flag);
#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)
f = kiss_fft_functor<DeviceType>(d_data,d_tmp,plan->cfg_slow_forward,length);
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);
#else
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) {
f = kiss_fft_functor<DeviceType>(d_data,d_tmp,plan->cfg_fast_forward,length1);
Kokkos::parallel_for(total1/length1,f);

9
src/SPIN/compute_spin.cpp
View File

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

51
src/SPIN/pair_spin_dipole_cut.cpp
View File

@ -45,9 +45,10 @@ PairSpinDipoleCut::PairSpinDipoleCut(LAMMPS *lmp) : PairSpin(lmp)
hbar = force->hplanck/MY_2PI; // eV/(rad.THz)
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
// mub2mu0 = mub * mub * mu_0 / (4.0*MY_PI); // in eV
mub2mu0hbinv = mub2mu0 / hbar; // in rad.THz
}
@ -232,36 +233,44 @@ void PairSpinDipoleCut::compute(int eflag, int vflag)
local_cut2 = cut_spin_long[itype][jtype]*cut_spin_long[itype][jtype];
// compute dipolar interaction
if (rsq < local_cut2) {
r2inv = 1.0/rsq;
r3inv = r2inv*rinv;
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 (rsq <= local_cut2) {
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 (eflag) {
if (rsq <= local_cut2) {
evdwl -= (spi[0]*fmi[0] + spi[1]*fmi[1] + spi[2]*fmi[2]);
evdwl *= 0.5*hbar;
emag[i] += evdwl;
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];
}
} else evdwl = 0.0;
if (evflag) ev_tally_xyz(i,j,nlocal,newton_pair,
evdwl,ecoul,fi[0],fi[1],fi[2],rij[0],rij[1],rij[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],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];
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[1] -= pre * (eij[1] * bij + (sjeij*spi[1] + sieij*spj[1]));

48
src/SPIN/pair_spin_dipole_long.cpp
View File

@ -49,7 +49,7 @@ PairSpinDipoleLong::PairSpinDipoleLong(LAMMPS *lmp) : PairSpin(lmp)
hbar = force->hplanck/MY_2PI; // eV/(rad.THz)
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
mub2mu0hbinv = mub2mu0 / hbar; // in rad.THz
@ -281,32 +281,37 @@ void PairSpinDipoleLong::compute(int eflag, int vflag)
bij[3] = (5.0*bij[2] + pre3*expm2) * r2inv;
compute_long(i,j,eij,bij,fmi,spi,spj);
compute_long_mech(i,j,eij,bij,fmi,spi,spj);
}
if (lattice_flag)
compute_long_mech(i,j,eij,bij,fmi,spi,spj);
// force accumulation
if (eflag) {
if (rsq <= local_cut2) {
evdwl -= (spi[0]*fmi[0] + spi[1]*fmi[1] + spi[2]*fmi[2]);
evdwl *= 0.5*hbar;
emag[i] += evdwl;
}
} else evdwl = 0.0;
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 (rsq <= local_cut2) {
evdwl -= (spi[0]*fmi[0] + spi[1]*fmi[1] + spi[2]*fmi[2]);
evdwl *= 0.5*hbar;
emag[i] += evdwl;
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];
}
} else evdwl = 0.0;
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],rij[0],rij[1],rij[2]);
if (evflag) ev_tally_xyz(i,j,nlocal,newton_pair,
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];
jlist = firstneigh[ii];
jnum = numneigh[ii];
//itype = type[i];
for (jj = 0; jj < jnum; 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;
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];
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];

46
src/SPIN/pair_spin_dmi.cpp
View File

@ -244,31 +244,35 @@ void PairSpinDmi::compute(int eflag, int vflag)
if (rsq <= local_cut2) {
compute_dmi(i,j,eij,fmi,spj);
if (lattice_flag) {
if (lattice_flag)
compute_dmi_mech(i,j,rsq,eij,fi,spi,spj);
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;
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);
}
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) {
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,
evdwl,ecoul,fi[0],fi[1],fi[2],delx,dely,delz);
}
}
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);
cdmz = (dmix*csy - dmiy*csz);
fi[0] += irij*cdmx;
fi[1] += irij*cdmy;
fi[2] += irij*cdmz;
fi[0] += 0.5*irij*cdmx;
fi[1] += 0.5*irij*cdmy;
fi[2] += 0.5*irij*cdmz;
}
/* ----------------------------------------------------------------------

147
src/SPIN/pair_spin_exchange.cpp
View File

@ -37,6 +37,14 @@ using namespace LAMMPS_NS;
/* ---------------------------------------------------------------------- */
PairSpinExchange::PairSpinExchange(LAMMPS *lmp) :
PairSpin(lmp)
{
e_offset = 0;
}
/* ---------------------------------------------------------------------- */
PairSpinExchange::~PairSpinExchange()
{
if (allocated) {
@ -59,6 +67,8 @@ void PairSpinExchange::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
@ -84,9 +94,9 @@ void PairSpinExchange::coeff(int narg, char **arg)
// check if args correct
if (strcmp(arg[2],"exchange") != 0)
error->all(FLERR,"Incorrect args in pair_style command");
if (narg != 7)
error->all(FLERR,"Incorrect args in pair_style command");
error->all(FLERR,"Incorrect args for pair coefficients");
if ((narg != 7) && (narg != 9))
error->all(FLERR,"Incorrect args for pair coefficients");
int ilo,ihi,jlo,jhi;
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
int iarg = 7;
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);
// 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;
for (int i = ilo; i <= ihi; i++) {
for (int j = MAX(jlo,i); j <= jhi; j++) {
@ -228,31 +252,34 @@ void PairSpinExchange::compute(int eflag, int vflag)
if (rsq <= local_cut2) {
compute_exchange(i,j,rsq,fmi,spj);
if (lattice_flag) {
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);
}
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) {
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,
evdwl,ecoul,fi[0],fi[1],fi[2],delx,dely,delz);
}
}
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
------------------------------------------------------------------------- */
void PairSpinExchange::compute_exchange_mech(int i, int j, double rsq, double eij[3],
double fi[3], double spi[3], double spj[3])
void PairSpinExchange::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, ra, rr, iJ3;
double Jex, Jex_mech, ra, sdots;
double rr, iJ3;
double fx, fy, fz;
itype = type[i];
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 *= 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];
fi[1] -= Jex_mech*eij[1];
fi[2] -= Jex_mech*eij[2];
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)*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
------------------------------------------------------------------------- */
// double PairSpinExchange::compute_energy(int i, int j, double rsq, double spi[3], double spj[3])
// {
// int *type = atom->type;
// int itype, jtype;
// double Jex, ra;
// double energy = 0.0;
// itype = type[i];
// jtype = type[j];
//
// Jex = J1_mech[itype][jtype];
// ra = rsq/J3[itype][jtype]/J3[itype][jtype];
// Jex = 4.0*Jex*ra;
// Jex *= (1.0-J2[itype][jtype]*ra);
// Jex *= exp(-ra);
//
// energy = Jex*(spi[0]*spj[0]+spi[1]*spj[1]+spi[2]*spj[2]);
// return energy;
// }
double PairSpinExchange::compute_energy(int i, int j, double rsq, double spi[3], double spj[3])
{
int *type = atom->type;
int itype, jtype;
double Jex, ra, sdots;
double energy = 0.0;
itype = type[i];
jtype = type[j];
Jex = J1_mech[itype][jtype];
ra = rsq/J3[itype][jtype]/J3[itype][jtype];
Jex = 4.0*Jex*ra;
Jex *= (1.0-J2[itype][jtype]*ra);
Jex *= exp(-ra);
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);
} 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
@ -495,6 +547,7 @@ void PairSpinExchange::read_restart(FILE *fp)
void PairSpinExchange::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);
}
@ -507,10 +560,12 @@ void PairSpinExchange::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);
}

6
src/SPIN/pair_spin_exchange.h
View File

@ -26,7 +26,7 @@ namespace LAMMPS_NS {
class PairSpinExchange : public PairSpin {
public:
PairSpinExchange(LAMMPS *lmp) : PairSpin(lmp) {}
PairSpinExchange(class LAMMPS *);
virtual ~PairSpinExchange();
void settings(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_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 read_restart(FILE *);
@ -49,6 +48,7 @@ class PairSpinExchange : public PairSpin {
double cut_spin_exchange_global; // global exchange cutoff distance
protected:
int e_offset; // apply energy offset
double **J1_mag; // exchange coeffs in eV
double **J1_mech; // mech exchange coeffs in
double **J2, **J3; // J1 in eV, J2 adim, J3 in Ang

635
src/SPIN/pair_spin_exchange_biquadratic.cpp Normal file
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);
}

87
src/SPIN/pair_spin_exchange_biquadratic.h Normal file
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.
*/

46
src/SPIN/pair_spin_magelec.cpp
View File

@ -237,31 +237,35 @@ void PairSpinMagelec::compute(int eflag, int vflag)
if (rsq <= local_cut2) {
compute_magelec(i,j,eij,fmi,spj);
if (lattice_flag) {
if (lattice_flag)
compute_magelec_mech(i,j,fi,spi,spj);
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;
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);
}
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) {
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,
evdwl,ecoul,fi[0],fi[1],fi[2],delx,dely,delz);
}
}
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];
meiz *= ME_mech[itype][jtype];
fi[0] += (meiy*vz - meiz*vy);
fi[1] += (meiz*vx - meix*vz);
fi[2] += (meix*vy - meiy*vx);
fi[0] += 0.5*(meiy*vz - meiz*vy);
fi[1] += 0.5*(meiz*vx - meix*vz);
fi[2] += 0.5*(meix*vy - meiy*vx);
}

50
src/SPIN/pair_spin_neel.cpp
View File

@ -246,31 +246,33 @@ void PairSpinNeel::compute(int eflag, int vflag)
if (rsq <= local_cut2) {
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);
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 (eflag) {
evdwl -= compute_neel_energy(i,j,rsq,eij,spi,spj);
emag[i] += evdwl;
} else evdwl = 0.0;
if (evflag) ev_tally_xyz(i,j,nlocal,newton_pair,
evdwl,ecoul,fi[0],fi[1],fi[2],rij[0],rij[1],rij[2]);
}
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) {
evdwl = compute_neel_energy(i,j,rsq,eij,spi,spj);
evdwl *= 0.5*hbar;
emag[i] += evdwl;
} else evdwl = 0.0;
if (evflag) ev_tally_xyz(i,j,nlocal,newton_pair,
evdwl,ecoul,fi[0],fi[1],fi[2],rij[0],rij[1],rij[2]);
}
}
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
fi[0] = pdx + pq1x + pq2x;
fi[1] = pdy + pq1y + pq2y;
fi[2] = pdz + pq1z + pq2z;
fi[0] = 0.5*(pdx + pq1x + pq2x);
fi[1] = 0.5*(pdy + pq1y + pq2y);
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
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 *= exp(-ra);
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 *= 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;
epq2 = q2r*(eij_si*eij_sj_3+eij_sj*eij_si_3);
return (epd+epq1+epq2);
return 0.5*(epd+epq1+epq2);
}
/* ----------------------------------------------------------------------

251
src/USER-MESONT/pair_mesont_tpm.cpp
View File

@ -9,7 +9,6 @@
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
Contributing author: Maxim Shugaev (UVA), mvs9t@virginia.edu
------------------------------------------------------------------------- */
@ -29,6 +28,7 @@
#include <cstring>
#include <vector>
#include <cmath>
#include <array>
#include <fstream>
#include <sstream>
@ -36,28 +36,17 @@
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 {
public:
MESONTList(const Atom* atom, const NeighList* nblist, double rc2);
MESONTList(const Atom* atom, const NeighList* nblist);
~MESONTList() {};
//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
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
//(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
int get_idx(int idx) const;
//return list of indexes for conversion from sorted representation
@ -69,22 +58,22 @@ public:
//check if the node is the end of the tube
bool is_end(int idx) const;
array2003<int, 2> get_segment(int idx) const;
array2003<int, 3> get_triplet(int idx) const;
std::array<int,2> get_segment(int idx) const;
std::array<int,3> get_triplet(int idx) const;
static const int cnt_end = -1;
static const int domain_end = -2;
static const int not_cnt = -3;
private:
std::vector<array2003<int, 2> > chain_list, segments;
std::vector<array2003<int, 3> > triplets;
std::vector<std::vector<array2003<int, 2> > > nb_chains;
std::vector<std::array<int,2>> chain_list, segments;
std::vector<std::array<int,3>> triplets;
std::vector<std::vector<std::array<int,2>>> nb_chains;
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 {
return nb_chains;
}
@ -106,25 +95,25 @@ inline const std::vector<int>& MESONTList::get_idxb_list() const {
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 {
return segments;
}
inline const std::vector<array2003<int, 3> > & MESONTList::get_triplets()
inline const std::vector<std::array<int,3>> & MESONTList::get_triplets()
const {
return triplets;
}
inline array2003<int, 2> MESONTList::get_segment(int idx) const {
array2003<int, 2> result;
inline std::array<int,2> MESONTList::get_segment(int idx) const {
std::array<int,2> result;
result[0] = chain_list[idx][0];
result[1] = idx;
return result;
}
inline array2003<int, 3> MESONTList::get_triplet(int idx) const {
array2003<int, 3> result;
inline std::array<int,3> MESONTList::get_triplet(int idx) const {
std::array<int,3> result;
result[0] = chain_list[idx][0];
result[1] = idx;
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;
//number of local atoms at the node
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;
//total number of atoms in the node and ghost shell
int ntot = atom->nlocal + atom->nghost;
tagint* const g_id = atom->tag;
tagint** const bonds = atom->bond_nt;
@ -178,9 +168,7 @@ MESONTList::MESONTList(const Atom* atom, const NeighList* nblist, double /* rc2
int* ilist = nblist->ilist;
//convert bonds to local id representation
array2003<int, 2> tmp_arr;
tmp_arr[0] = not_cnt; tmp_arr[1] = not_cnt;
chain_list.resize(ntot, tmp_arr);
chain_list.resize(ntot, {not_cnt,not_cnt});
for (int ii = 0; ii < nall; ii++) {
int i = ilist[ii];
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;
for (int j = 0; j < nnb; 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[nb][1] = i;
break;
@ -223,22 +211,16 @@ MESONTList::MESONTList(const Atom* atom, const NeighList* nblist, double /* rc2
for (int i = 0; i < nlocal; i++) {
if (chain_list[i][0] == not_cnt) continue;
if (chain_list[i][0] != cnt_end && chain_list[i][0] != domain_end &&
g_id[i] < g_id[chain_list[i][0]]){
array2003<int, 2> tmp_c;
tmp_c[0] = i; tmp_c[1] = chain_list[i][0];
segments.push_back(tmp_c);
}
g_id[i] < g_id[chain_list[i][0]])
segments.push_back({i,chain_list[i][0]});
if (chain_list[i][1] != cnt_end && chain_list[i][1] != domain_end &&
g_id[i] < g_id[chain_list[i][1]]){
array2003<int, 2> tmp_c;
tmp_c[0] = i; tmp_c[1] = chain_list[i][1];
segments.push_back(tmp_c);
}
g_id[i] < g_id[chain_list[i][1]])
segments.push_back({i,chain_list[i][1]});
}
int nbonds = segments.size();
//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] != cnt_end && chain_list[i][0] != 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];
if ((j == nnb - 1) || (nb_list[j] + 1 != nb_list[j+1]) ||
(idx_next == cnt_end) || (idx_next == domain_end)) {
array2003<int, 2> chain;
std::array<int,2> chain;
chain[0] = idx_s;
chain[1] = nb_list[j];
//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;
if (chain_list[idx1][1] != cnt_end &&
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];
}
}
@ -311,8 +293,9 @@ PairMESONTTPM::PairMESONTTPM(LAMMPS *lmp) : Pair(lmp) {
eatom_s = nullptr;
eatom_b = nullptr;
eatom_t = nullptr;
nmax = 0;
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");
}
@ -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);
//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;
//total number of atoms in the node and ghost shell
int ntot = atom->nlocal + atom->nghost;
int newton_pair = force->newton_pair;
if(!newton_pair)
if (!newton_pair)
error->all(FLERR,"Pair style mesont/tpm requires newton pair on");
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 +
std::pow((2.0*RT + TPBRcutoff),2)));
if (cut_global < Rcut_min){
if (cut_global < Rcut_min) {
std::stringstream err;
err << "The selected cutoff is too small for the current system : " <<
"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
MESONTList ntlist(atom, list, cut_global*cut_global);
MESONTList ntlist(atom, list);
//reorder data to make it contiguous within tubes
//and compatible with Fortran functions
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<int> b_sort(nall);
for (int i = 0; i < nall; i++){
for (int i = 0; i < nall; i++) {
int idx = ntlist.get_idx(i);
for (int j = 0; j < 3; j++) x_sort[3*i+j] = x[idx][j];
b_sort[i] = buckling[idx];
@ -385,7 +380,7 @@ void PairMESONTTPM::compute(int eflag, int vflag){
//bending potential
int n_triplets = ntlist.get_triplets().size();
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
int idx_s0 = ntlist.get_idxb(t[0]);
int idx_s1 = ntlist.get_idxb(t[1]);
@ -412,13 +407,13 @@ void PairMESONTTPM::compute(int eflag, int vflag){
}
//share new values of buckling
if (BendingMode == 1){
for (int i = 0; i < nall; i++){
if (BendingMode == 1) {
for (int i = 0; i < nall; i++) {
int idx = ntlist.get_idx(i);
buckling[idx] = b_sort[i];
}
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);
b_sort[i] = buckling[idx];
}
@ -429,7 +424,7 @@ void PairMESONTTPM::compute(int eflag, int vflag){
double Rmax = 0.0;
Lmax = 0.0;
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
int idx_s0 = ntlist.get_idxb(s[0]);
//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,
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
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 end1 = ntlist.get_idx(chain[0]); //chain ends (real representation)
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;
if (!E1 && ntlist.get_triplet(end1)[0] != MESONTList::domain_end &&
ntlist.get_triplet(ntlist.get_triplet(end1)[0])[0] ==
MESONTList::cnt_end){
MESONTList::cnt_end) {
Ee = 1;
int idx = ntlist.get_idxb(ntlist.get_triplet(end1)[0]);
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 &&
ntlist.get_triplet(ntlist.get_triplet(end2)[2])[2] ==
MESONTList::cnt_end){
MESONTList::cnt_end) {
Ee = 2;
int idx = ntlist.get_idxb(ntlist.get_triplet(end2)[2]);
Xe = &(x_sort[3*idx]);
@ -500,7 +495,7 @@ void PairMESONTTPM::compute(int eflag, int vflag){
//check if cutoff is chosen correctly
Rcut_min = std::max(2.0*Lmax, std::sqrt(0.5*Lmax*Lmax +
std::pow((2.0*Rmax + TPBRcutoff),2)));
if (cut_global < Rcut_min){
if (cut_global < Rcut_min) {
std::stringstream err;
err << "The selected cutoff is too small for the current system : " <<
"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());
}
// set per atom values and accumulators
// reallocate per-atom arrays if necessary
if (atom->nmax > maxeatom) {
maxeatom = atom->nmax;
memory->destroy(eatom);
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 (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;
energy_b = 0.0; energy_t = 0.0;
for (int i = 0; i < 6; i++) virial[i] = 0.0;
for (int i = 0; i < ntot; i++){
eatom[i] = 0.0; eatom_s[i] = 0.0;
eatom_b[i] = 0.0; eatom_t[i] = 0.0;
}
for (int i = 0; i < ntot; 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);
for (int j = 0; j < 3; j++) f[idx][j] += f_sort[3*i+j];
eatom_s[idx] = u_ts_sort[i];
eatom_b[idx] = u_tb_sort[i];
eatom_t[idx] = 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];
vatom[idx][0] = s_sort[9*i+0]; //xx
vatom[idx][1] = s_sort[9*i+4]; //yy
vatom[idx][2] = s_sort[9*i+8]; //zz
vatom[idx][3] = s_sort[9*i+1]; //xy
vatom[idx][4] = s_sort[9*i+2]; //xz
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];
for (int i = 0; i < nall; i++) {
int idx = ntlist.get_idx(i);
for (int j = 0; j < 3; j++) f[idx][j] += f_sort[3*i+j];
buckling[idx] = b_sort[i];
}
eng_vdwl = energy_s + energy_b + energy_t;
if (eflag_global) {
eng_vdwl = 0.0; energy_s = 0.0;
energy_b = 0.0; energy_t = 0.0;
for (int i = 0; i < nall; 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_b[i] = 0.0; eatom_t[i] = 0.0;
}
for (int i = 0; i < nall; i++) {
int idx = ntlist.get_idx(i);
eatom_s[idx] = u_ts_sort[i];
eatom_b[idx] = u_tb_sort[i];
eatom_t[idx] = u_tt_sort[i];
eatom[idx] = u_ts_sort[i] + u_tb_sort[i] + 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][1] = s_sort[9*i+4]; //yy
vatom[idx][2] = s_sort[9*i+8]; //zz
vatom[idx][3] = s_sort[9*i+1]; //xy
vatom[idx][4] = s_sort[9*i+2]; //xz
vatom[idx][5] = s_sort[9*i+5]; //yz
}
}
}
/* ----------------------------------------------------------------------
allocate all arrays
------------------------------------------------------------------------- */
void PairMESONTTPM::allocate(){
void PairMESONTTPM::allocate() {
allocated = 1;
int n = atom->ntypes;
@ -583,7 +582,7 @@ void PairMESONTTPM::allocate(){
global settings
------------------------------------------------------------------------- */
void PairMESONTTPM::settings(int narg, char **arg){
void PairMESONTTPM::settings(int narg, char **arg) {
if ((narg == 0) || (narg > 4))
error->all(FLERR,"Illegal pair_style command");
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
------------------------------------------------------------------------- */
void PairMESONTTPM::coeff(int narg, char **arg){
void PairMESONTTPM::coeff(int narg, char **arg) {
if ((narg < 2) || (narg > 3))
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
------------------------------------------------------------------------- */
double PairMESONTTPM::init_one(int i, int j){
double PairMESONTTPM::init_one(int i, int j) {
if (setflag[i][j] == 0) {
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
------------------------------------------------------------------------- */
void PairMESONTTPM::write_restart(FILE *fp){
void PairMESONTTPM::write_restart(FILE *fp) {
write_restart_settings(fp);
int i,j;
@ -690,7 +689,7 @@ void PairMESONTTPM::write_restart(FILE *fp){
proc 0 reads from restart file, bcasts
------------------------------------------------------------------------- */
void PairMESONTTPM::read_restart(FILE *fp){
void PairMESONTTPM::read_restart(FILE *fp) {
read_restart_settings(fp);
allocate();
@ -713,7 +712,7 @@ void PairMESONTTPM::read_restart(FILE *fp){
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(&TPMType,sizeof(int),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
------------------------------------------------------------------------- */
void PairMESONTTPM::read_restart_settings(FILE *fp){
void PairMESONTTPM::read_restart_settings(FILE *fp) {
int me = comm->me;
if (me == 0) {
fread(&BendingMode,sizeof(int),1,fp);
@ -761,7 +760,7 @@ void PairMESONTTPM::read_restart_settings(FILE *fp){
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++)
fprintf(fp,"%d\n",i);
}
@ -770,7 +769,7 @@ void PairMESONTTPM::write_data(FILE *fp){
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 j = i; j <= atom->ntypes; 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)
int r = neighbor->request(this,instance_me);
neighbor->requests[r]->half = false;
@ -786,7 +785,7 @@ void PairMESONTTPM::init_style(){
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;
else if (strcmp(str,"mesonttpm_Eb_tot") == 0) return &energy_b;
else if (strcmp(str,"mesonttpm_Et_tot") == 0) return &energy_t;

1
src/USER-MESONT/pair_mesont_tpm.h
View File

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

2
src/USER-PLUMED/README
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.
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
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:

6
src/atom_masks.h
View File

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

4
unittest/force-styles/CMakeLists.txt
View File

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