From f48b71f46ba66655f1b9a20753737391ba3e42c5 Mon Sep 17 00:00:00 2001
From: Steve Plimpton <sjplimp@sandia.gov>
Date: Fri, 4 Nov 2016 10:56:04 -0600
Subject: [PATCH] added examples/threebody, fix reaxc/speceies/kk
---
doc/src/compute_temp_profile.txt | 4 +-
doc/src/fix_reax_bonds.txt | 26 +++
doc/src/fix_reaxc_species.txt | 26 +++
doc/src/suffix.txt | 7 +
examples/threebody/BNC.tersoff | 74 +++++++
examples/threebody/CdTeZnSeHgS0.sw | 233 ++++++++++++++++++++
examples/threebody/InP.vashishta | 38 ++++
examples/threebody/Si.sw | 18 ++
examples/threebody/in.threebody | 116 ++++++++++
potentials/CdTe.sw | 8 +-
potentials/CdTeZnSeHgS0.sw | 14 +-
src/KOKKOS/fix_qeq_reax_kokkos.cpp | 76 +++++--
src/KOKKOS/fix_qeq_reax_kokkos.h | 2 +-
src/KOKKOS/fix_reaxc_species_kokkos.cpp | 20 +-
src/KOKKOS/fix_wall_reflect_kokkos.h | 2 -
src/KOKKOS/neigh_full_kokkos.h | 1 +
src/KOKKOS/neigh_list_kokkos.cpp | 5 +-
src/KOKKOS/neigh_list_kokkos.h | 3 +-
src/KOKKOS/pair_reax_c_kokkos.h | 2 -
src/KOKKOS/pair_sw_kokkos.cpp | 68 ++++--
src/KOKKOS/pair_sw_kokkos.h | 8 +
src/KOKKOS/pair_tersoff_kokkos.cpp | 270 +++++++++++++++---------
src/KOKKOS/pair_tersoff_kokkos.h | 10 +-
src/KOKKOS/pair_tersoff_mod_kokkos.cpp | 237 +++++++++++++--------
src/KOKKOS/pair_tersoff_mod_kokkos.h | 10 +-
src/KOKKOS/pair_tersoff_zbl_kokkos.cpp | 235 +++++++++++++--------
src/KOKKOS/pair_tersoff_zbl_kokkos.h | 9 +-
src/pair_morse.cpp | 3 +-
28 files changed, 1184 insertions(+), 341 deletions(-)
create mode 100644 examples/threebody/BNC.tersoff
create mode 100644 examples/threebody/CdTeZnSeHgS0.sw
create mode 100644 examples/threebody/InP.vashishta
create mode 100644 examples/threebody/Si.sw
create mode 100644 examples/threebody/in.threebody
diff --git a/doc/src/compute_temp_profile.txt b/doc/src/compute_temp_profile.txt
index fedbd5cb32..54eebd6d8f 100644
--- a/doc/src/compute_temp_profile.txt
+++ b/doc/src/compute_temp_profile.txt
@@ -69,8 +69,8 @@ velocity for each atom. Note that if there is only one atom in the
bin, its thermal velocity will thus be 0.0.
After the spatially-averaged velocity field has been subtracted from
-each atom, the temperature is calculated by the formula KE = (dim/2 N
-- dim*Nx*Ny*Nz) k T, where KE = total kinetic energy of the group of
+each atom, the temperature is calculated by the formula KE = (dim*N
+- dim*Nx*Ny*Nz) k T/2, where KE = total kinetic energy of the group of
atoms (sum of 1/2 m v^2), dim = 2 or 3 = dimensionality of the
simulation, N = number of atoms in the group, k = Boltzmann constant,
and T = temperature. The dim*Nx*Ny*Nz term are degrees of freedom
diff --git a/doc/src/fix_reax_bonds.txt b/doc/src/fix_reax_bonds.txt
index a0396ce210..a85e140b62 100644
--- a/doc/src/fix_reax_bonds.txt
+++ b/doc/src/fix_reax_bonds.txt
@@ -8,6 +8,7 @@
fix reax/bonds command :h3
fix reax/c/bonds command :h3
+fix reax/c/bonds/kk command :h3
[Syntax:]
@@ -47,6 +48,31 @@ commands"_Section_howto.html#howto_15. No parameter of this fix can
be used with the {start/stop} keywords of the "run"_run.html command.
This fix is not invoked during "energy minimization"_minimize.html.
+:line
+
+Styles with a {gpu}, {intel}, {kk}, {omp}, or {opt} suffix are
+functionally the same as the corresponding style without the suffix.
+They have been optimized to run faster, depending on your available
+hardware, as discussed in "Section_accelerate"_Section_accelerate.html
+of the manual. The accelerated styles take the same arguments and
+should produce the same results, except for round-off and precision
+issues.
+
+These accelerated styles are part of the GPU, USER-INTEL, KOKKOS,
+USER-OMP and OPT packages, respectively. They are only enabled if
+LAMMPS was built with those packages. See the "Making
+LAMMPS"_Section_start.html#start_3 section for more info.
+
+You can specify the accelerated styles explicitly in your input script
+by including their suffix, or you can use the "-suffix command-line
+switch"_Section_start.html#start_7 when you invoke LAMMPS, or you can
+use the "suffix"_suffix.html command in your input script.
+
+See "Section_accelerate"_Section_accelerate.html of the manual for
+more instructions on how to use the accelerated styles effectively.
+
+:line
+
[Restrictions:]
The fix reax/bonds command requires that the "pair_style
diff --git a/doc/src/fix_reaxc_species.txt b/doc/src/fix_reaxc_species.txt
index 630c802a83..00db91900e 100644
--- a/doc/src/fix_reaxc_species.txt
+++ b/doc/src/fix_reaxc_species.txt
@@ -7,6 +7,7 @@
:line
fix reax/c/species command :h3
+fix reax/c/species/kk command :h3
[Syntax:]
@@ -129,6 +130,31 @@ No parameter of this fix can be used with the {start/stop} keywords of
the "run"_run.html command. This fix is not invoked during "energy
minimization"_minimize.html.
+:line
+
+Styles with a {gpu}, {intel}, {kk}, {omp}, or {opt} suffix are
+functionally the same as the corresponding style without the suffix.
+They have been optimized to run faster, depending on your available
+hardware, as discussed in "Section_accelerate"_Section_accelerate.html
+of the manual. The accelerated styles take the same arguments and
+should produce the same results, except for round-off and precision
+issues.
+
+These accelerated styles are part of the GPU, USER-INTEL, KOKKOS,
+USER-OMP and OPT packages, respectively. They are only enabled if
+LAMMPS was built with those packages. See the "Making
+LAMMPS"_Section_start.html#start_3 section for more info.
+
+You can specify the accelerated styles explicitly in your input script
+by including their suffix, or you can use the "-suffix command-line
+switch"_Section_start.html#start_7 when you invoke LAMMPS, or you can
+use the "suffix"_suffix.html command in your input script.
+
+See "Section_accelerate"_Section_accelerate.html of the manual for
+more instructions on how to use the accelerated styles effectively.
+
+:line
+
[Restrictions:]
The fix species currently only works with
diff --git a/doc/src/suffix.txt b/doc/src/suffix.txt
index 669d483f00..70c6a3031b 100644
--- a/doc/src/suffix.txt
+++ b/doc/src/suffix.txt
@@ -94,6 +94,13 @@ Of course this is also possible by not using any suffix commands, and
explictly appending or not appending the suffix to the relevant
commands in your input script.
+NOTE: The default "run_style"_run_style.html verlet is invoked prior to
+reading the input script and is therefore not affected by a suffix command
+in the input script. The KOKKOS package requires "run_style verlet/kk",
+so when using the KOKKOS package it is necessary to either use the command
+line "-sf kk" command or add an explicit "run_style verlet" command to the
+input script.
+
[Restrictions:] none
[Related commands:]
diff --git a/examples/threebody/BNC.tersoff b/examples/threebody/BNC.tersoff
new file mode 100644
index 0000000000..1630062524
--- /dev/null
+++ b/examples/threebody/BNC.tersoff
@@ -0,0 +1,74 @@
+# DATE: 2013-03-21 CONTRIBUTOR: Cem Sevik CITATION: Kinaci, Haskins, Sevik and Cagin, Phys Rev B, 86, 115410 (2012)
+# Tersoff parameters for B, C, and BN-C hybrid based graphene like nano structures
+# multiple entries can be added to this file, LAMMPS reads the ones it needs
+
+# these entries are in LAMMPS "metal" units:
+# A,B = eV; lambda1,lambda2,lambda3 = 1/Angstroms; R,D = Angstroms
+# other quantities are unitless
+
+# Cem Sevik (csevik at anadolu.edu.tr) takes full blame for this
+# file. It specifies B-N, B-C, and N-C interaction parameters
+# generated and published by the reseacrh group of Prof. Tahir Cagin.
+
+# 1. Physical Review B 84, 085409 2011
+# Characterization of thermal transport in low-dimensional boron nitride nanostructures,
+#
+
+# 2. Physical Review B 86, 075403 2012
+# Influence of disorder on thermal transport properties of boron nitride nanostructures
+#
+
+# 3. Physical Review B 86, 075403 2012, Please see for further information about B-C and N-C parameters
+# Thermal conductivity of BN-C nanostructures
+#
+
+# The file also specifies C-C, interaction parameters
+# generated and published by the reseacrh group of Dr. D. A. Broido
+# Physical Review B 81, 205441 2010
+# Optimized Tersoff and Brenner empirical potential parameters for
+# lattice dynamics and phonon thermal transport in carbon nanotubes and graphene
+
+# Users in referring the full parameters can cite the full parameter paper (3) as:
+# A. Kinaci, J. B. Haskins, C. Sevik, T. Cagin, Physical Review B 86, 115410 (2012)
+# Thermal conductivity of BN-C nanostructures
+#
+
+# format of a single entry (one or more lines):
+# element 1, element 2, element 3,
+# m, gamma, lambda3, c, d, costheta0, n, beta, lambda2, B, R, D, lambda1, A
+
+N B B 3.0 1.0 0.0 25000 4.3484 -0.89000 0.72751 1.25724e-7 2.199 340.00 1.95 0.05 3.568 1380.0
+N B N 3.0 1.0 0.0 25000 4.3484 -0.89000 0.72751 1.25724e-7 2.199 340.00 1.95 0.05 3.568 1380.0
+N B C 3.0 1.0 0.0 25000 4.3484 -0.89000 0.72751 1.25724e-7 2.199 340.00 1.95 0.05 3.568 1380.0
+
+B N B 3.0 1.0 0.0 25000 4.3484 -0.89000 0.72751 1.25724e-7 2.199 340.00 1.95 0.05 3.568 1380.0
+B N N 3.0 1.0 0.0 25000 4.3484 -0.89000 0.72751 1.25724e-7 2.199 340.00 1.95 0.05 3.568 1380.0
+B N C 3.0 1.0 0.0 25000 4.3484 -0.89000 0.72751 1.25724e-7 2.199 340.00 1.95 0.05 3.568 1380.0
+
+N N B 3.0 1.0 0.0 17.7959 5.9484 0.00000 0.6184432 0.019251 2.6272721 138.77866 2.0 0.1 2.8293093 128.86866
+N N N 3.0 1.0 0.0 17.7959 5.9484 0.00000 0.6184432 0.019251 2.6272721 138.77866 2.0 0.1 2.8293093 128.86866
+N N C 3.0 1.0 0.0 17.7959 5.9484 0.00000 0.6184432 0.019251 2.6272721 138.77866 2.0 0.1 2.8293093 128.86866
+
+B B B 3.0 1.0 0.0 0.52629 0.001587 0.5 3.9929061 1.6e-6 2.0774982 43.132016 2.0 0.1 2.2372578 40.0520156
+B B N 3.0 1.0 0.0 0.52629 0.001587 0.5 3.9929061 1.6e-6 2.0774982 43.132016 2.0 0.1 2.2372578 40.0520156
+B B C 3.0 1.0 0.0 0.52629 0.001587 0.5 3.9929061 1.6e-6 2.0774982 43.132016 2.0 0.1 2.2372578 40.0520156
+
+C C C 3.0 1.0 0.0 3.8049e4 4.3484 -0.93000 0.72751 1.5724e-7 2.2119 430.00 1.95 0.15 3.4879 1393.6
+C C B 3.0 1.0 0.0 3.8049e4 4.3484 -0.93000 0.72751 1.5724e-7 2.2119 430.00 1.95 0.15 3.4879 1393.6
+C C N 3.0 1.0 0.0 3.8049e4 4.3484 -0.93000 0.72751 1.5724e-7 2.2119 430.00 1.95 0.15 3.4879 1393.6
+
+C B B 3.0 1.0 0.0 3.8049e4 4.3484 -0.93000 0.72751 1.5724e-7 2.2054 339.068910 1.95 0.10 3.5279 1386.78
+C B N 3.0 1.0 0.0 3.8049e4 4.3484 -0.93000 0.72751 1.5724e-7 2.2054 339.068910 1.95 0.10 3.5279 1386.78
+C B C 3.0 1.0 0.0 3.8049e4 4.3484 -0.93000 0.72751 1.5724e-7 2.2054 339.068910 1.95 0.10 3.5279 1386.78
+
+C N B 3.0 1.0 0.0 3.8049e4 4.3484 -0.93000 0.72751 1.5724e-7 2.2054 387.575152 1.95 0.10 3.5279 1386.78
+C N N 3.0 1.0 0.0 3.8049e4 4.3484 -0.93000 0.72751 1.5724e-7 2.2054 387.575152 1.95 0.10 3.5279 1386.78
+C N C 3.0 1.0 0.0 3.8049e4 4.3484 -0.93000 0.72751 1.5724e-7 2.2054 387.575152 1.95 0.10 3.5279 1386.78
+
+B C C 3.0 1.0 0.0 25000 4.3484 -0.89000 0.72751 1.25724e-7 2.2054 339.068910 1.95 0.10 3.5279 1386.78
+B C B 3.0 1.0 0.0 25000 4.3484 -0.89000 0.72751 1.25724e-7 2.2054 339.068910 1.95 0.10 3.5279 1386.78
+B C N 3.0 1.0 0.0 25000 4.3484 -0.89000 0.72751 1.25724e-7 2.2054 339.068910 1.95 0.10 3.5279 1386.78
+
+N C C 3.0 1.0 0.0 25000 4.3484 -0.89000 0.72751 1.25724e-7 2.2054 387.575152 1.95 0.10 3.5279 1386.78
+N C B 3.0 1.0 0.0 25000 4.3484 -0.89000 0.72751 1.25724e-7 2.2054 387.575152 1.95 0.10 3.5279 1386.78
+N C N 3.0 1.0 0.0 25000 4.3484 -0.89000 0.72751 1.25724e-7 2.2054 387.575152 1.95 0.10 3.5279 1386.78
diff --git a/examples/threebody/CdTeZnSeHgS0.sw b/examples/threebody/CdTeZnSeHgS0.sw
new file mode 100644
index 0000000000..d6f05d41df
--- /dev/null
+++ b/examples/threebody/CdTeZnSeHgS0.sw
@@ -0,0 +1,233 @@
+### DATE: 2013-08-09 CONTRIBUTOR: X. W. Zhou, xzhou@sandia.gov, CITATION: Zhou, Ward, Martin, van Swol, Cruz-Campa, and D. Zubia, Phys. Rev. B, 88, 085309 (2013).
+#
+# Note that the way the parameters can be entered is not unique.
+# As one way, we assume that eps_ijk is equal to eps_ik and
+# lambda_ijk is equal to sqrt(lambda_ij*eps_ij*lambda_ik*eps_ik)/eps_ik,
+# and all other parameters in the ijk line are for ik.
+#
+# The twobody ik pair parameters are entered on the i*k lines, where *
+# can be any species. This is consistent with the LAMMPS requirement
+# that twobody ik parameters be defined on the ikk line. Entries on all
+# the other i*k lines are ignored by LAMMPS
+#
+# These entries are in LAMMPS "metal" units: epsilon = eV;
+# sigma = Angstroms; other quantities are unitless
+#
+# cutoff distance = 4.632
+# eps sigma a lambda gamma cos(theta) A B p q tol
+Cd Cd Cd 1.182358e+00 2.663951e+00 1.527956e+00 3.250000e+01 1.200000e+00 -3.333333e-01 7.917000e+00 7.674460e-01 4.000000e+00 0.000000e+00 0.000000e+00
+Cd Cd Te 1.385284e+00 2.352141e+00 1.810919e+00 3.002537e+01 1.200000e+00 -3.333333e-01 7.049600e+00 8.861252e-01 4.000000e+00 0.000000e+00 0.000000e+00
+Cd Cd Zn 6.908179e-01 2.238699e+00 1.812616e+00 4.251831e+01 1.200000e+00 -3.333333e-01 7.049600e+00 1.010632e+00 4.000000e+00 0.000000e+00 0.000000e+00
+Cd Cd Se 1.352371e+00 2.045165e+00 1.953387e+00 3.038855e+01 1.200000e+00 -3.333333e-01 7.049600e+00 1.116149e+00 4.000000e+00 0.000000e+00 0.000000e+00
+Cd Cd Hg 4.881231e-01 2.432694e+00 1.677987e+00 5.058167e+01 1.200000e+00 -3.333333e-01 7.049600e+00 6.250999e-01 4.000000e+00 0.000000e+00 0.000000e+00
+Cd Cd S 1.300376e+00 1.804151e+00 2.124568e+00 3.099013e+01 1.200000e+00 -3.333333e-01 7.049600e+00 1.540087e+00 4.000000e+00 0.000000e+00 0.000000e+00
+Cd Te Cd 1.182358e+00 2.663951e+00 1.527956e+00 3.517858e+01 1.200000e+00 -3.333333e-01 7.917000e+00 7.674460e-01 4.000000e+00 0.000000e+00 0.000000e+00
+Cd Te Te 1.385284e+00 2.352141e+00 1.810919e+00 3.250000e+01 1.200000e+00 -3.333333e-01 7.049600e+00 8.861252e-01 4.000000e+00 0.000000e+00 0.000000e+00
+Cd Te Zn 6.908179e-01 2.238699e+00 1.812616e+00 4.602259e+01 1.200000e+00 -3.333333e-01 7.049600e+00 1.010632e+00 4.000000e+00 0.000000e+00 0.000000e+00
+Cd Te Se 1.352371e+00 2.045165e+00 1.953387e+00 3.289311e+01 1.200000e+00 -3.333333e-01 7.049600e+00 1.116149e+00 4.000000e+00 0.000000e+00 0.000000e+00
+Cd Te Hg 4.881231e-01 2.432694e+00 1.677987e+00 5.475051e+01 1.200000e+00 -3.333333e-01 7.049600e+00 6.250999e-01 4.000000e+00 0.000000e+00 0.000000e+00
+Cd Te S 1.300376e+00 1.804151e+00 2.124568e+00 3.354428e+01 1.200000e+00 -3.333333e-01 7.049600e+00 1.540087e+00 4.000000e+00 0.000000e+00 0.000000e+00
+Cd Zn Cd 1.182358e+00 2.663951e+00 1.527956e+00 2.484224e+01 1.200000e+00 -3.333333e-01 7.917000e+00 7.674460e-01 4.000000e+00 0.000000e+00 0.000000e+00
+Cd Zn Te 1.385284e+00 2.352141e+00 1.810919e+00 2.295069e+01 1.200000e+00 -3.333333e-01 7.049600e+00 8.861252e-01 4.000000e+00 0.000000e+00 0.000000e+00
+Cd Zn Zn 6.908179e-01 2.238699e+00 1.812616e+00 3.250000e+01 1.200000e+00 -3.333333e-01 7.049600e+00 1.010632e+00 4.000000e+00 0.000000e+00 0.000000e+00
+Cd Zn Se 1.352371e+00 2.045165e+00 1.953387e+00 2.322829e+01 1.200000e+00 -3.333333e-01 7.049600e+00 1.116149e+00 4.000000e+00 0.000000e+00 0.000000e+00
+Cd Zn Hg 4.881231e-01 2.432694e+00 1.677987e+00 3.866344e+01 1.200000e+00 -3.333333e-01 7.049600e+00 6.250999e-01 4.000000e+00 0.000000e+00 0.000000e+00
+Cd Zn S 1.300376e+00 1.804151e+00 2.124568e+00 2.368813e+01 1.200000e+00 -3.333333e-01 7.049600e+00 1.540087e+00 4.000000e+00 0.000000e+00 0.000000e+00
+Cd Se Cd 1.182358e+00 2.663951e+00 1.527956e+00 3.475816e+01 1.200000e+00 -3.333333e-01 7.917000e+00 7.674460e-01 4.000000e+00 0.000000e+00 0.000000e+00
+Cd Se Te 1.385284e+00 2.352141e+00 1.810919e+00 3.211159e+01 1.200000e+00 -3.333333e-01 7.049600e+00 8.861252e-01 4.000000e+00 0.000000e+00 0.000000e+00
+Cd Se Zn 6.908179e-01 2.238699e+00 1.812616e+00 4.547256e+01 1.200000e+00 -3.333333e-01 7.049600e+00 1.010632e+00 4.000000e+00 0.000000e+00 0.000000e+00
+Cd Se Se 1.352371e+00 2.045165e+00 1.953387e+00 3.250000e+01 1.200000e+00 -3.333333e-01 7.049600e+00 1.116149e+00 4.000000e+00 0.000000e+00 0.000000e+00
+Cd Se Hg 4.881231e-01 2.432694e+00 1.677987e+00 5.409618e+01 1.200000e+00 -3.333333e-01 7.049600e+00 6.250999e-01 4.000000e+00 0.000000e+00 0.000000e+00
+Cd Se S 1.300376e+00 1.804151e+00 2.124568e+00 3.314338e+01 1.200000e+00 -3.333333e-01 7.049600e+00 1.540087e+00 4.000000e+00 0.000000e+00 0.000000e+00
+Cd Hg Cd 1.182358e+00 2.663951e+00 1.527956e+00 2.088207e+01 1.200000e+00 -3.333333e-01 7.917000e+00 7.674460e-01 4.000000e+00 0.000000e+00 0.000000e+00
+Cd Hg Te 1.385284e+00 2.352141e+00 1.810919e+00 1.929206e+01 1.200000e+00 -3.333333e-01 7.049600e+00 8.861252e-01 4.000000e+00 0.000000e+00 0.000000e+00
+Cd Hg Zn 6.908179e-01 2.238699e+00 1.812616e+00 2.731909e+01 1.200000e+00 -3.333333e-01 7.049600e+00 1.010632e+00 4.000000e+00 0.000000e+00 0.000000e+00
+Cd Hg Se 1.352371e+00 2.045165e+00 1.953387e+00 1.952541e+01 1.200000e+00 -3.333333e-01 7.049600e+00 1.116149e+00 4.000000e+00 0.000000e+00 0.000000e+00
+Cd Hg Hg 4.881231e-01 2.432694e+00 1.677987e+00 3.250000e+01 1.200000e+00 -3.333333e-01 7.049600e+00 6.250999e-01 4.000000e+00 0.000000e+00 0.000000e+00
+Cd Hg S 1.300376e+00 1.804151e+00 2.124568e+00 1.991194e+01 1.200000e+00 -3.333333e-01 7.049600e+00 1.540087e+00 4.000000e+00 0.000000e+00 0.000000e+00
+Cd S Cd 1.182358e+00 2.663951e+00 1.527956e+00 3.408343e+01 1.200000e+00 -3.333333e-01 7.917000e+00 7.674460e-01 4.000000e+00 0.000000e+00 0.000000e+00
+Cd S Te 1.385284e+00 2.352141e+00 1.810919e+00 3.148823e+01 1.200000e+00 -3.333333e-01 7.049600e+00 8.861252e-01 4.000000e+00 0.000000e+00 0.000000e+00
+Cd S Zn 6.908179e-01 2.238699e+00 1.812616e+00 4.458985e+01 1.200000e+00 -3.333333e-01 7.049600e+00 1.010632e+00 4.000000e+00 0.000000e+00 0.000000e+00
+Cd S Se 1.352371e+00 2.045165e+00 1.953387e+00 3.186911e+01 1.200000e+00 -3.333333e-01 7.049600e+00 1.116149e+00 4.000000e+00 0.000000e+00 0.000000e+00
+Cd S Hg 4.881231e-01 2.432694e+00 1.677987e+00 5.304605e+01 1.200000e+00 -3.333333e-01 7.049600e+00 6.250999e-01 4.000000e+00 0.000000e+00 0.000000e+00
+Cd S S 1.300376e+00 1.804151e+00 2.124568e+00 3.250000e+01 1.200000e+00 -3.333333e-01 7.049600e+00 1.540087e+00 4.000000e+00 0.000000e+00 0.000000e+00
+Te Cd Cd 1.385284e+00 2.352141e+00 1.810919e+00 3.250000e+01 1.200000e+00 -3.333333e-01 7.049600e+00 8.861252e-01 4.000000e+00 0.000000e+00 0.000000e+00
+Te Cd Te 1.849775e+00 2.905254e+00 1.594353e+00 2.812506e+01 1.200000e+00 -3.333333e-01 7.917000e+00 7.307283e-01 4.000000e+00 0.000000e+00 0.000000e+00
+Te Cd Zn 1.546239e+00 2.056363e+00 1.907922e+00 3.076200e+01 1.200000e+00 -3.333333e-01 7.049600e+00 1.255846e+00 4.000000e+00 0.000000e+00 0.000000e+00
+Te Cd Se 1.295053e+00 2.231716e+00 1.809645e+00 3.361313e+01 1.200000e+00 -3.333333e-01 7.049600e+00 1.005396e+00 4.000000e+00 0.000000e+00 0.000000e+00
+Te Cd Hg 1.204715e+00 2.135591e+00 1.892491e+00 3.485063e+01 1.200000e+00 -3.333333e-01 7.049600e+00 1.445180e+00 4.000000e+00 0.000000e+00 0.000000e+00
+Te Cd S 1.450015e+00 2.297301e+00 1.726905e+00 3.176630e+01 1.200000e+00 -3.333333e-01 7.049600e+00 7.794685e-01 4.000000e+00 0.000000e+00 0.000000e+00
+Te Te Cd 1.385284e+00 2.352141e+00 1.810919e+00 3.755548e+01 1.200000e+00 -3.333333e-01 7.049600e+00 8.861252e-01 4.000000e+00 0.000000e+00 0.000000e+00
+Te Te Te 1.849775e+00 2.905254e+00 1.594353e+00 3.250000e+01 1.200000e+00 -3.333333e-01 7.917000e+00 7.307283e-01 4.000000e+00 0.000000e+00 0.000000e+00
+Te Te Zn 1.546239e+00 2.056363e+00 1.907922e+00 3.554713e+01 1.200000e+00 -3.333333e-01 7.049600e+00 1.255846e+00 4.000000e+00 0.000000e+00 0.000000e+00
+Te Te Se 1.295053e+00 2.231716e+00 1.809645e+00 3.884177e+01 1.200000e+00 -3.333333e-01 7.049600e+00 1.005396e+00 4.000000e+00 0.000000e+00 0.000000e+00
+Te Te Hg 1.204715e+00 2.135591e+00 1.892491e+00 4.027176e+01 1.200000e+00 -3.333333e-01 7.049600e+00 1.445180e+00 4.000000e+00 0.000000e+00 0.000000e+00
+Te Te S 1.450015e+00 2.297301e+00 1.726905e+00 3.670765e+01 1.200000e+00 -3.333333e-01 7.049600e+00 7.794685e-01 4.000000e+00 0.000000e+00 0.000000e+00
+Te Zn Cd 1.385284e+00 2.352141e+00 1.810919e+00 3.433620e+01 1.200000e+00 -3.333333e-01 7.049600e+00 8.861252e-01 4.000000e+00 0.000000e+00 0.000000e+00
+Te Zn Te 1.849775e+00 2.905254e+00 1.594353e+00 2.971408e+01 1.200000e+00 -3.333333e-01 7.917000e+00 7.307283e-01 4.000000e+00 0.000000e+00 0.000000e+00
+Te Zn Zn 1.546239e+00 2.056363e+00 1.907922e+00 3.250000e+01 1.200000e+00 -3.333333e-01 7.049600e+00 1.255846e+00 4.000000e+00 0.000000e+00 0.000000e+00
+Te Zn Se 1.295053e+00 2.231716e+00 1.809645e+00 3.551222e+01 1.200000e+00 -3.333333e-01 7.049600e+00 1.005396e+00 4.000000e+00 0.000000e+00 0.000000e+00
+Te Zn Hg 1.204715e+00 2.135591e+00 1.892491e+00 3.681964e+01 1.200000e+00 -3.333333e-01 7.049600e+00 1.445180e+00 4.000000e+00 0.000000e+00 0.000000e+00
+Te Zn S 1.450015e+00 2.297301e+00 1.726905e+00 3.356105e+01 1.200000e+00 -3.333333e-01 7.049600e+00 7.794685e-01 4.000000e+00 0.000000e+00 0.000000e+00
+Te Se Cd 1.385284e+00 2.352141e+00 1.810919e+00 3.142373e+01 1.200000e+00 -3.333333e-01 7.049600e+00 8.861252e-01 4.000000e+00 0.000000e+00 0.000000e+00
+Te Se Te 1.849775e+00 2.905254e+00 1.594353e+00 2.719366e+01 1.200000e+00 -3.333333e-01 7.917000e+00 7.307283e-01 4.000000e+00 0.000000e+00 0.000000e+00
+Te Se Zn 1.546239e+00 2.056363e+00 1.907922e+00 2.974328e+01 1.200000e+00 -3.333333e-01 7.049600e+00 1.255846e+00 4.000000e+00 0.000000e+00 0.000000e+00
+Te Se Se 1.295053e+00 2.231716e+00 1.809645e+00 3.250000e+01 1.200000e+00 -3.333333e-01 7.049600e+00 1.005396e+00 4.000000e+00 0.000000e+00 0.000000e+00
+Te Se Hg 1.204715e+00 2.135591e+00 1.892491e+00 3.369652e+01 1.200000e+00 -3.333333e-01 7.049600e+00 1.445180e+00 4.000000e+00 0.000000e+00 0.000000e+00
+Te Se S 1.450015e+00 2.297301e+00 1.726905e+00 3.071433e+01 1.200000e+00 -3.333333e-01 7.049600e+00 7.794685e-01 4.000000e+00 0.000000e+00 0.000000e+00
+Te Hg Cd 1.385284e+00 2.352141e+00 1.810919e+00 3.030791e+01 1.200000e+00 -3.333333e-01 7.049600e+00 8.861252e-01 4.000000e+00 0.000000e+00 0.000000e+00
+Te Hg Te 1.849775e+00 2.905254e+00 1.594353e+00 2.622805e+01 1.200000e+00 -3.333333e-01 7.917000e+00 7.307283e-01 4.000000e+00 0.000000e+00 0.000000e+00
+Te Hg Zn 1.546239e+00 2.056363e+00 1.907922e+00 2.868714e+01 1.200000e+00 -3.333333e-01 7.049600e+00 1.255846e+00 4.000000e+00 0.000000e+00 0.000000e+00
+Te Hg Se 1.295053e+00 2.231716e+00 1.809645e+00 3.134597e+01 1.200000e+00 -3.333333e-01 7.049600e+00 1.005396e+00 4.000000e+00 0.000000e+00 0.000000e+00
+Te Hg Hg 1.204715e+00 2.135591e+00 1.892491e+00 3.250000e+01 1.200000e+00 -3.333333e-01 7.049600e+00 1.445180e+00 4.000000e+00 0.000000e+00 0.000000e+00
+Te Hg S 1.450015e+00 2.297301e+00 1.726905e+00 2.962370e+01 1.200000e+00 -3.333333e-01 7.049600e+00 7.794685e-01 4.000000e+00 0.000000e+00 0.000000e+00
+Te S Cd 1.385284e+00 2.352141e+00 1.810919e+00 3.325065e+01 1.200000e+00 -3.333333e-01 7.049600e+00 8.861252e-01 4.000000e+00 0.000000e+00 0.000000e+00
+Te S Te 1.849775e+00 2.905254e+00 1.594353e+00 2.877465e+01 1.200000e+00 -3.333333e-01 7.917000e+00 7.307283e-01 4.000000e+00 0.000000e+00 0.000000e+00
+Te S Zn 1.546239e+00 2.056363e+00 1.907922e+00 3.147250e+01 1.200000e+00 -3.333333e-01 7.049600e+00 1.255846e+00 4.000000e+00 0.000000e+00 0.000000e+00
+Te S Se 1.295053e+00 2.231716e+00 1.809645e+00 3.438949e+01 1.200000e+00 -3.333333e-01 7.049600e+00 1.005396e+00 4.000000e+00 0.000000e+00 0.000000e+00
+Te S Hg 1.204715e+00 2.135591e+00 1.892491e+00 3.565557e+01 1.200000e+00 -3.333333e-01 7.049600e+00 1.445180e+00 4.000000e+00 0.000000e+00 0.000000e+00
+Te S S 1.450015e+00 2.297301e+00 1.726905e+00 3.250000e+01 1.200000e+00 -3.333333e-01 7.049600e+00 7.794685e-01 4.000000e+00 0.000000e+00 0.000000e+00
+Zn Cd Cd 6.908179e-01 2.238699e+00 1.812616e+00 3.250000e+01 1.200000e+00 -3.333333e-01 7.049600e+00 1.010632e+00 4.000000e+00 0.000000e+00 0.000000e+00
+Zn Cd Te 1.546239e+00 2.056363e+00 1.907922e+00 2.172335e+01 1.200000e+00 -3.333333e-01 7.049600e+00 1.255846e+00 4.000000e+00 0.000000e+00 0.000000e+00
+Zn Cd Zn 1.392961e+00 2.367650e+00 1.525521e+00 2.288736e+01 1.200000e+00 -3.333333e-01 7.917000e+00 7.676279e-01 4.000000e+00 0.000000e+00 0.000000e+00
+Zn Cd Se 1.691181e+00 2.028827e+00 1.836907e+00 2.077161e+01 1.200000e+00 -3.333333e-01 7.049600e+00 9.510930e-01 4.000000e+00 0.000000e+00 0.000000e+00
+Zn Cd Hg 4.951616e-01 2.239186e+00 1.761363e+00 3.838766e+01 1.200000e+00 -3.333333e-01 7.049600e+00 7.461167e-01 4.000000e+00 0.000000e+00 0.000000e+00
+Zn Cd S 2.208390e+00 2.323783e+00 1.589241e+00 1.817721e+01 1.200000e+00 -3.333333e-01 7.049600e+00 4.643181e-01 4.000000e+00 0.000000e+00 0.000000e+00
+Zn Te Cd 6.908179e-01 2.238699e+00 1.812616e+00 4.862279e+01 1.200000e+00 -3.333333e-01 7.049600e+00 1.010632e+00 4.000000e+00 0.000000e+00 0.000000e+00
+Zn Te Te 1.546239e+00 2.056363e+00 1.907922e+00 3.250000e+01 1.200000e+00 -3.333333e-01 7.049600e+00 1.255846e+00 4.000000e+00 0.000000e+00 0.000000e+00
+Zn Te Zn 1.392961e+00 2.367650e+00 1.525521e+00 3.424146e+01 1.200000e+00 -3.333333e-01 7.917000e+00 7.676279e-01 4.000000e+00 0.000000e+00 0.000000e+00
+Zn Te Se 1.691181e+00 2.028827e+00 1.836907e+00 3.107611e+01 1.200000e+00 -3.333333e-01 7.049600e+00 9.510930e-01 4.000000e+00 0.000000e+00 0.000000e+00
+Zn Te Hg 4.951616e-01 2.239186e+00 1.761363e+00 5.743124e+01 1.200000e+00 -3.333333e-01 7.049600e+00 7.461167e-01 4.000000e+00 0.000000e+00 0.000000e+00
+Zn Te S 2.208390e+00 2.323783e+00 1.589241e+00 2.719467e+01 1.200000e+00 -3.333333e-01 7.049600e+00 4.643181e-01 4.000000e+00 0.000000e+00 0.000000e+00
+Zn Zn Cd 6.908179e-01 2.238699e+00 1.812616e+00 4.614993e+01 1.200000e+00 -3.333333e-01 7.049600e+00 1.010632e+00 4.000000e+00 0.000000e+00 0.000000e+00
+Zn Zn Te 1.546239e+00 2.056363e+00 1.907922e+00 3.084711e+01 1.200000e+00 -3.333333e-01 7.049600e+00 1.255846e+00 4.000000e+00 0.000000e+00 0.000000e+00
+Zn Zn Zn 1.392961e+00 2.367650e+00 1.525521e+00 3.250000e+01 1.200000e+00 -3.333333e-01 7.917000e+00 7.676279e-01 4.000000e+00 0.000000e+00 0.000000e+00
+Zn Zn Se 1.691181e+00 2.028827e+00 1.836907e+00 2.949563e+01 1.200000e+00 -3.333333e-01 7.049600e+00 9.510930e-01 4.000000e+00 0.000000e+00 0.000000e+00
+Zn Zn Hg 4.951616e-01 2.239186e+00 1.761363e+00 5.451040e+01 1.200000e+00 -3.333333e-01 7.049600e+00 7.461167e-01 4.000000e+00 0.000000e+00 0.000000e+00
+Zn Zn S 2.208390e+00 2.323783e+00 1.589241e+00 2.581160e+01 1.200000e+00 -3.333333e-01 7.049600e+00 4.643181e-01 4.000000e+00 0.000000e+00 0.000000e+00
+Zn Se Cd 6.908179e-01 2.238699e+00 1.812616e+00 5.085067e+01 1.200000e+00 -3.333333e-01 7.049600e+00 1.010632e+00 4.000000e+00 0.000000e+00 0.000000e+00
+Zn Se Te 1.546239e+00 2.056363e+00 1.907922e+00 3.398914e+01 1.200000e+00 -3.333333e-01 7.049600e+00 1.255846e+00 4.000000e+00 0.000000e+00 0.000000e+00
+Zn Se Zn 1.392961e+00 2.367650e+00 1.525521e+00 3.581039e+01 1.200000e+00 -3.333333e-01 7.917000e+00 7.676279e-01 4.000000e+00 0.000000e+00 0.000000e+00
+Zn Se Se 1.691181e+00 2.028827e+00 1.836907e+00 3.250000e+01 1.200000e+00 -3.333333e-01 7.049600e+00 9.510930e-01 4.000000e+00 0.000000e+00 0.000000e+00
+Zn Se Hg 4.951616e-01 2.239186e+00 1.761363e+00 6.006272e+01 1.200000e+00 -3.333333e-01 7.049600e+00 7.461167e-01 4.000000e+00 0.000000e+00 0.000000e+00
+Zn Se S 2.208390e+00 2.323783e+00 1.589241e+00 2.844072e+01 1.200000e+00 -3.333333e-01 7.049600e+00 4.643181e-01 4.000000e+00 0.000000e+00 0.000000e+00
+Zn Hg Cd 6.908179e-01 2.238699e+00 1.812616e+00 2.751535e+01 1.200000e+00 -3.333333e-01 7.049600e+00 1.010632e+00 4.000000e+00 0.000000e+00 0.000000e+00
+Zn Hg Te 1.546239e+00 2.056363e+00 1.907922e+00 1.839156e+01 1.200000e+00 -3.333333e-01 7.049600e+00 1.255846e+00 4.000000e+00 0.000000e+00 0.000000e+00
+Zn Hg Zn 1.392961e+00 2.367650e+00 1.525521e+00 1.937704e+01 1.200000e+00 -3.333333e-01 7.917000e+00 7.676279e-01 4.000000e+00 0.000000e+00 0.000000e+00
+Zn Hg Se 1.691181e+00 2.028827e+00 1.836907e+00 1.758578e+01 1.200000e+00 -3.333333e-01 7.049600e+00 9.510930e-01 4.000000e+00 0.000000e+00 0.000000e+00
+Zn Hg Hg 4.951616e-01 2.239186e+00 1.761363e+00 3.250000e+01 1.200000e+00 -3.333333e-01 7.049600e+00 7.461167e-01 4.000000e+00 0.000000e+00 0.000000e+00
+Zn Hg S 2.208390e+00 2.323783e+00 1.589241e+00 1.538930e+01 1.200000e+00 -3.333333e-01 7.049600e+00 4.643181e-01 4.000000e+00 0.000000e+00 0.000000e+00
+Zn S Cd 6.908179e-01 2.238699e+00 1.812616e+00 5.810847e+01 1.200000e+00 -3.333333e-01 7.049600e+00 1.010632e+00 4.000000e+00 0.000000e+00 0.000000e+00
+Zn S Te 1.546239e+00 2.056363e+00 1.907922e+00 3.884033e+01 1.200000e+00 -3.333333e-01 7.049600e+00 1.255846e+00 4.000000e+00 0.000000e+00 0.000000e+00
+Zn S Zn 1.392961e+00 2.367650e+00 1.525521e+00 4.092153e+01 1.200000e+00 -3.333333e-01 7.917000e+00 7.676279e-01 4.000000e+00 0.000000e+00 0.000000e+00
+Zn S Se 1.691181e+00 2.028827e+00 1.836907e+00 3.713865e+01 1.200000e+00 -3.333333e-01 7.049600e+00 9.510930e-01 4.000000e+00 0.000000e+00 0.000000e+00
+Zn S Hg 4.951616e-01 2.239186e+00 1.761363e+00 6.863534e+01 1.200000e+00 -3.333333e-01 7.049600e+00 7.461167e-01 4.000000e+00 0.000000e+00 0.000000e+00
+Zn S S 2.208390e+00 2.323783e+00 1.589241e+00 3.250000e+01 1.200000e+00 -3.333333e-01 7.049600e+00 4.643181e-01 4.000000e+00 0.000000e+00 0.000000e+00
+Se Cd Cd 1.352371e+00 2.045165e+00 1.953387e+00 3.250000e+01 1.200000e+00 -3.333333e-01 7.049600e+00 1.116149e+00 4.000000e+00 0.000000e+00 0.000000e+00
+Se Cd Te 1.295053e+00 2.231716e+00 1.809645e+00 3.321142e+01 1.200000e+00 -3.333333e-01 7.049600e+00 1.005396e+00 4.000000e+00 0.000000e+00 0.000000e+00
+Se Cd Zn 1.691181e+00 2.028827e+00 1.836907e+00 2.906271e+01 1.200000e+00 -3.333333e-01 7.049600e+00 9.510930e-01 4.000000e+00 0.000000e+00 0.000000e+00
+Se Cd Se 2.400781e+00 2.789002e+00 1.544925e+00 2.439242e+01 1.200000e+00 -3.333333e-01 7.917000e+00 7.672131e-01 4.000000e+00 0.000000e+00 0.000000e+00
+Se Cd Hg 1.299758e+00 2.113406e+00 1.831821e+00 3.315126e+01 1.200000e+00 -3.333333e-01 7.049600e+00 1.150200e+00 4.000000e+00 0.000000e+00 0.000000e+00
+Se Cd S 1.307592e+00 2.229392e+00 1.747782e+00 3.305180e+01 1.200000e+00 -3.333333e-01 7.049600e+00 6.932325e-01 4.000000e+00 0.000000e+00 0.000000e+00
+Se Te Cd 1.352371e+00 2.045165e+00 1.953387e+00 3.180382e+01 1.200000e+00 -3.333333e-01 7.049600e+00 1.116149e+00 4.000000e+00 0.000000e+00 0.000000e+00
+Se Te Te 1.295053e+00 2.231716e+00 1.809645e+00 3.250000e+01 1.200000e+00 -3.333333e-01 7.049600e+00 1.005396e+00 4.000000e+00 0.000000e+00 0.000000e+00
+Se Te Zn 1.691181e+00 2.028827e+00 1.836907e+00 2.844016e+01 1.200000e+00 -3.333333e-01 7.049600e+00 9.510930e-01 4.000000e+00 0.000000e+00 0.000000e+00
+Se Te Se 2.400781e+00 2.789002e+00 1.544925e+00 2.386992e+01 1.200000e+00 -3.333333e-01 7.917000e+00 7.672131e-01 4.000000e+00 0.000000e+00 0.000000e+00
+Se Te Hg 1.299758e+00 2.113406e+00 1.831821e+00 3.244113e+01 1.200000e+00 -3.333333e-01 7.049600e+00 1.150200e+00 4.000000e+00 0.000000e+00 0.000000e+00
+Se Te S 1.307592e+00 2.229392e+00 1.747782e+00 3.234380e+01 1.200000e+00 -3.333333e-01 7.049600e+00 6.932325e-01 4.000000e+00 0.000000e+00 0.000000e+00
+Se Zn Cd 1.352371e+00 2.045165e+00 1.953387e+00 3.634382e+01 1.200000e+00 -3.333333e-01 7.049600e+00 1.116149e+00 4.000000e+00 0.000000e+00 0.000000e+00
+Se Zn Te 1.295053e+00 2.231716e+00 1.809645e+00 3.713938e+01 1.200000e+00 -3.333333e-01 7.049600e+00 1.005396e+00 4.000000e+00 0.000000e+00 0.000000e+00
+Se Zn Zn 1.691181e+00 2.028827e+00 1.836907e+00 3.250000e+01 1.200000e+00 -3.333333e-01 7.049600e+00 9.510930e-01 4.000000e+00 0.000000e+00 0.000000e+00
+Se Zn Se 2.400781e+00 2.789002e+00 1.544925e+00 2.727735e+01 1.200000e+00 -3.333333e-01 7.917000e+00 7.672131e-01 4.000000e+00 0.000000e+00 0.000000e+00
+Se Zn Hg 1.299758e+00 2.113406e+00 1.831821e+00 3.707211e+01 1.200000e+00 -3.333333e-01 7.049600e+00 1.150200e+00 4.000000e+00 0.000000e+00 0.000000e+00
+Se Zn S 1.307592e+00 2.229392e+00 1.747782e+00 3.696088e+01 1.200000e+00 -3.333333e-01 7.049600e+00 6.932325e-01 4.000000e+00 0.000000e+00 0.000000e+00
+Se Se Cd 1.352371e+00 2.045165e+00 1.953387e+00 4.330238e+01 1.200000e+00 -3.333333e-01 7.049600e+00 1.116149e+00 4.000000e+00 0.000000e+00 0.000000e+00
+Se Se Te 1.295053e+00 2.231716e+00 1.809645e+00 4.425026e+01 1.200000e+00 -3.333333e-01 7.049600e+00 1.005396e+00 4.000000e+00 0.000000e+00 0.000000e+00
+Se Se Zn 1.691181e+00 2.028827e+00 1.836907e+00 3.872260e+01 1.200000e+00 -3.333333e-01 7.049600e+00 9.510930e-01 4.000000e+00 0.000000e+00 0.000000e+00
+Se Se Se 2.400781e+00 2.789002e+00 1.544925e+00 3.250000e+01 1.200000e+00 -3.333333e-01 7.917000e+00 7.672131e-01 4.000000e+00 0.000000e+00 0.000000e+00
+Se Se Hg 1.299758e+00 2.113406e+00 1.831821e+00 4.417011e+01 1.200000e+00 -3.333333e-01 7.049600e+00 1.150200e+00 4.000000e+00 0.000000e+00 0.000000e+00
+Se Se S 1.307592e+00 2.229392e+00 1.747782e+00 4.403758e+01 1.200000e+00 -3.333333e-01 7.049600e+00 6.932325e-01 4.000000e+00 0.000000e+00 0.000000e+00
+Se Hg Cd 1.352371e+00 2.045165e+00 1.953387e+00 3.186153e+01 1.200000e+00 -3.333333e-01 7.049600e+00 1.116149e+00 4.000000e+00 0.000000e+00 0.000000e+00
+Se Hg Te 1.295053e+00 2.231716e+00 1.809645e+00 3.255898e+01 1.200000e+00 -3.333333e-01 7.049600e+00 1.005396e+00 4.000000e+00 0.000000e+00 0.000000e+00
+Se Hg Zn 1.691181e+00 2.028827e+00 1.836907e+00 2.849177e+01 1.200000e+00 -3.333333e-01 7.049600e+00 9.510930e-01 4.000000e+00 0.000000e+00 0.000000e+00
+Se Hg Se 2.400781e+00 2.789002e+00 1.544925e+00 2.391323e+01 1.200000e+00 -3.333333e-01 7.917000e+00 7.672131e-01 4.000000e+00 0.000000e+00 0.000000e+00
+Se Hg Hg 1.299758e+00 2.113406e+00 1.831821e+00 3.250000e+01 1.200000e+00 -3.333333e-01 7.049600e+00 1.150200e+00 4.000000e+00 0.000000e+00 0.000000e+00
+Se Hg S 1.307592e+00 2.229392e+00 1.747782e+00 3.240249e+01 1.200000e+00 -3.333333e-01 7.049600e+00 6.932325e-01 4.000000e+00 0.000000e+00 0.000000e+00
+Se S Cd 1.352371e+00 2.045165e+00 1.953387e+00 3.195742e+01 1.200000e+00 -3.333333e-01 7.049600e+00 1.116149e+00 4.000000e+00 0.000000e+00 0.000000e+00
+Se S Te 1.295053e+00 2.231716e+00 1.809645e+00 3.265696e+01 1.200000e+00 -3.333333e-01 7.049600e+00 1.005396e+00 4.000000e+00 0.000000e+00 0.000000e+00
+Se S Zn 1.691181e+00 2.028827e+00 1.836907e+00 2.857751e+01 1.200000e+00 -3.333333e-01 7.049600e+00 9.510930e-01 4.000000e+00 0.000000e+00 0.000000e+00
+Se S Se 2.400781e+00 2.789002e+00 1.544925e+00 2.398520e+01 1.200000e+00 -3.333333e-01 7.917000e+00 7.672131e-01 4.000000e+00 0.000000e+00 0.000000e+00
+Se S Hg 1.299758e+00 2.113406e+00 1.831821e+00 3.259780e+01 1.200000e+00 -3.333333e-01 7.049600e+00 1.150200e+00 4.000000e+00 0.000000e+00 0.000000e+00
+Se S S 1.307592e+00 2.229392e+00 1.747782e+00 3.250000e+01 1.200000e+00 -3.333333e-01 7.049600e+00 6.932325e-01 4.000000e+00 0.000000e+00 0.000000e+00
+Hg Cd Cd 4.881231e-01 2.432694e+00 1.677987e+00 3.250000e+01 1.200000e+00 -3.333333e-01 7.049600e+00 6.250999e-01 4.000000e+00 0.000000e+00 0.000000e+00
+Hg Cd Te 1.204715e+00 2.135591e+00 1.892491e+00 2.068740e+01 1.200000e+00 -3.333333e-01 7.049600e+00 1.445180e+00 4.000000e+00 0.000000e+00 0.000000e+00
+Hg Cd Zn 4.951616e-01 2.239186e+00 1.761363e+00 3.226819e+01 1.200000e+00 -3.333333e-01 7.049600e+00 7.461167e-01 4.000000e+00 0.000000e+00 0.000000e+00
+Hg Cd Se 1.299758e+00 2.113406e+00 1.831821e+00 1.991668e+01 1.200000e+00 -3.333333e-01 7.049600e+00 1.150200e+00 4.000000e+00 0.000000e+00 0.000000e+00
+Hg Cd Hg 1.272807e+00 2.699097e+00 1.498503e+00 2.012643e+01 1.200000e+00 -3.333333e-01 7.917000e+00 1.211532e+00 4.000000e+00 0.000000e+00 0.000000e+00
+Hg Cd S 1.531211e+00 2.025045e+00 1.833708e+00 1.834976e+01 1.200000e+00 -3.333333e-01 7.049600e+00 1.184541e+00 4.000000e+00 0.000000e+00 0.000000e+00
+Hg Te Cd 4.881231e-01 2.432694e+00 1.677987e+00 5.105765e+01 1.200000e+00 -3.333333e-01 7.049600e+00 6.250999e-01 4.000000e+00 0.000000e+00 0.000000e+00
+Hg Te Te 1.204715e+00 2.135591e+00 1.892491e+00 3.250000e+01 1.200000e+00 -3.333333e-01 7.049600e+00 1.445180e+00 4.000000e+00 0.000000e+00 0.000000e+00
+Hg Te Zn 4.951616e-01 2.239186e+00 1.761363e+00 5.069347e+01 1.200000e+00 -3.333333e-01 7.049600e+00 7.461167e-01 4.000000e+00 0.000000e+00 0.000000e+00
+Hg Te Se 1.299758e+00 2.113406e+00 1.831821e+00 3.128919e+01 1.200000e+00 -3.333333e-01 7.049600e+00 1.150200e+00 4.000000e+00 0.000000e+00 0.000000e+00
+Hg Te Hg 1.272807e+00 2.699097e+00 1.498503e+00 3.161872e+01 1.200000e+00 -3.333333e-01 7.917000e+00 1.211532e+00 4.000000e+00 0.000000e+00 0.000000e+00
+Hg Te S 1.531211e+00 2.025045e+00 1.833708e+00 2.882756e+01 1.200000e+00 -3.333333e-01 7.049600e+00 1.184541e+00 4.000000e+00 0.000000e+00 0.000000e+00
+Hg Zn Cd 4.881231e-01 2.432694e+00 1.677987e+00 3.273348e+01 1.200000e+00 -3.333333e-01 7.049600e+00 6.250999e-01 4.000000e+00 0.000000e+00 0.000000e+00
+Hg Zn Te 1.204715e+00 2.135591e+00 1.892491e+00 2.083602e+01 1.200000e+00 -3.333333e-01 7.049600e+00 1.445180e+00 4.000000e+00 0.000000e+00 0.000000e+00
+Hg Zn Zn 4.951616e-01 2.239186e+00 1.761363e+00 3.250000e+01 1.200000e+00 -3.333333e-01 7.049600e+00 7.461167e-01 4.000000e+00 0.000000e+00 0.000000e+00
+Hg Zn Se 1.299758e+00 2.113406e+00 1.831821e+00 2.005976e+01 1.200000e+00 -3.333333e-01 7.049600e+00 1.150200e+00 4.000000e+00 0.000000e+00 0.000000e+00
+Hg Zn Hg 1.272807e+00 2.699097e+00 1.498503e+00 2.027102e+01 1.200000e+00 -3.333333e-01 7.917000e+00 1.211532e+00 4.000000e+00 0.000000e+00 0.000000e+00
+Hg Zn S 1.531211e+00 2.025045e+00 1.833708e+00 1.848159e+01 1.200000e+00 -3.333333e-01 7.049600e+00 1.184541e+00 4.000000e+00 0.000000e+00 0.000000e+00
+Hg Se Cd 4.881231e-01 2.432694e+00 1.677987e+00 5.303345e+01 1.200000e+00 -3.333333e-01 7.049600e+00 6.250999e-01 4.000000e+00 0.000000e+00 0.000000e+00
+Hg Se Te 1.204715e+00 2.135591e+00 1.892491e+00 3.375766e+01 1.200000e+00 -3.333333e-01 7.049600e+00 1.445180e+00 4.000000e+00 0.000000e+00 0.000000e+00
+Hg Se Zn 4.951616e-01 2.239186e+00 1.761363e+00 5.265518e+01 1.200000e+00 -3.333333e-01 7.049600e+00 7.461167e-01 4.000000e+00 0.000000e+00 0.000000e+00
+Hg Se Se 1.299758e+00 2.113406e+00 1.831821e+00 3.250000e+01 1.200000e+00 -3.333333e-01 7.049600e+00 1.150200e+00 4.000000e+00 0.000000e+00 0.000000e+00
+Hg Se Hg 1.272807e+00 2.699097e+00 1.498503e+00 3.284228e+01 1.200000e+00 -3.333333e-01 7.917000e+00 1.211532e+00 4.000000e+00 0.000000e+00 0.000000e+00
+Hg Se S 1.531211e+00 2.025045e+00 1.833708e+00 2.994311e+01 1.200000e+00 -3.333333e-01 7.049600e+00 1.184541e+00 4.000000e+00 0.000000e+00 0.000000e+00
+Hg Hg Cd 4.881231e-01 2.432694e+00 1.677987e+00 5.248074e+01 1.200000e+00 -3.333333e-01 7.049600e+00 6.250999e-01 4.000000e+00 0.000000e+00 0.000000e+00
+Hg Hg Te 1.204715e+00 2.135591e+00 1.892491e+00 3.340584e+01 1.200000e+00 -3.333333e-01 7.049600e+00 1.445180e+00 4.000000e+00 0.000000e+00 0.000000e+00
+Hg Hg Zn 4.951616e-01 2.239186e+00 1.761363e+00 5.210641e+01 1.200000e+00 -3.333333e-01 7.049600e+00 7.461167e-01 4.000000e+00 0.000000e+00 0.000000e+00
+Hg Hg Se 1.299758e+00 2.113406e+00 1.831821e+00 3.216129e+01 1.200000e+00 -3.333333e-01 7.049600e+00 1.150200e+00 4.000000e+00 0.000000e+00 0.000000e+00
+Hg Hg Hg 1.272807e+00 2.699097e+00 1.498503e+00 3.250000e+01 1.200000e+00 -3.333333e-01 7.917000e+00 1.211532e+00 4.000000e+00 0.000000e+00 0.000000e+00
+Hg Hg S 1.531211e+00 2.025045e+00 1.833708e+00 2.963105e+01 1.200000e+00 -3.333333e-01 7.049600e+00 1.184541e+00 4.000000e+00 0.000000e+00 0.000000e+00
+Hg S Cd 4.881231e-01 2.432694e+00 1.677987e+00 5.756205e+01 1.200000e+00 -3.333333e-01 7.049600e+00 6.250999e-01 4.000000e+00 0.000000e+00 0.000000e+00
+Hg S Te 1.204715e+00 2.135591e+00 1.892491e+00 3.664028e+01 1.200000e+00 -3.333333e-01 7.049600e+00 1.445180e+00 4.000000e+00 0.000000e+00 0.000000e+00
+Hg S Zn 4.951616e-01 2.239186e+00 1.761363e+00 5.715148e+01 1.200000e+00 -3.333333e-01 7.049600e+00 7.461167e-01 4.000000e+00 0.000000e+00 0.000000e+00
+Hg S Se 1.299758e+00 2.113406e+00 1.831821e+00 3.527522e+01 1.200000e+00 -3.333333e-01 7.049600e+00 1.150200e+00 4.000000e+00 0.000000e+00 0.000000e+00
+Hg S Hg 1.272807e+00 2.699097e+00 1.498503e+00 3.564673e+01 1.200000e+00 -3.333333e-01 7.917000e+00 1.211532e+00 4.000000e+00 0.000000e+00 0.000000e+00
+Hg S S 1.531211e+00 2.025045e+00 1.833708e+00 3.250000e+01 1.200000e+00 -3.333333e-01 7.049600e+00 1.184541e+00 4.000000e+00 0.000000e+00 0.000000e+00
+S Cd Cd 1.300376e+00 1.804151e+00 2.124568e+00 3.250000e+01 1.200000e+00 -3.333333e-01 7.049600e+00 1.540087e+00 4.000000e+00 0.000000e+00 0.000000e+00
+S Cd Te 1.450015e+00 2.297301e+00 1.726905e+00 3.077737e+01 1.200000e+00 -3.333333e-01 7.049600e+00 7.794685e-01 4.000000e+00 0.000000e+00 0.000000e+00
+S Cd Zn 2.208390e+00 2.323783e+00 1.589241e+00 2.493905e+01 1.200000e+00 -3.333333e-01 7.049600e+00 4.643181e-01 4.000000e+00 0.000000e+00 0.000000e+00
+S Cd Se 1.307592e+00 2.229392e+00 1.747782e+00 3.241019e+01 1.200000e+00 -3.333333e-01 7.049600e+00 6.932325e-01 4.000000e+00 0.000000e+00 0.000000e+00
+S Cd Hg 1.531211e+00 2.025045e+00 1.833708e+00 2.995023e+01 1.200000e+00 -3.333333e-01 7.049600e+00 1.184541e+00 4.000000e+00 0.000000e+00 0.000000e+00
+S Cd S 2.434871e+00 2.423171e+00 1.711097e+00 2.375088e+01 1.200000e+00 -3.333333e-01 7.917000e+00 1.049688e+00 4.000000e+00 0.000000e+00 0.000000e+00
+S Te Cd 1.300376e+00 1.804151e+00 2.124568e+00 3.431904e+01 1.200000e+00 -3.333333e-01 7.049600e+00 1.540087e+00 4.000000e+00 0.000000e+00 0.000000e+00
+S Te Te 1.450015e+00 2.297301e+00 1.726905e+00 3.250000e+01 1.200000e+00 -3.333333e-01 7.049600e+00 7.794685e-01 4.000000e+00 0.000000e+00 0.000000e+00
+S Te Zn 2.208390e+00 2.323783e+00 1.589241e+00 2.633490e+01 1.200000e+00 -3.333333e-01 7.049600e+00 4.643181e-01 4.000000e+00 0.000000e+00 0.000000e+00
+S Te Se 1.307592e+00 2.229392e+00 1.747782e+00 3.422421e+01 1.200000e+00 -3.333333e-01 7.049600e+00 6.932325e-01 4.000000e+00 0.000000e+00 0.000000e+00
+S Te Hg 1.531211e+00 2.025045e+00 1.833708e+00 3.162656e+01 1.200000e+00 -3.333333e-01 7.049600e+00 1.184541e+00 4.000000e+00 0.000000e+00 0.000000e+00
+S Te S 2.434871e+00 2.423171e+00 1.711097e+00 2.508023e+01 1.200000e+00 -3.333333e-01 7.917000e+00 1.049688e+00 4.000000e+00 0.000000e+00 0.000000e+00
+S Zn Cd 1.300376e+00 1.804151e+00 2.124568e+00 4.235326e+01 1.200000e+00 -3.333333e-01 7.049600e+00 1.540087e+00 4.000000e+00 0.000000e+00 0.000000e+00
+S Zn Te 1.450015e+00 2.297301e+00 1.726905e+00 4.010837e+01 1.200000e+00 -3.333333e-01 7.049600e+00 7.794685e-01 4.000000e+00 0.000000e+00 0.000000e+00
+S Zn Zn 2.208390e+00 2.323783e+00 1.589241e+00 3.250000e+01 1.200000e+00 -3.333333e-01 7.049600e+00 4.643181e-01 4.000000e+00 0.000000e+00 0.000000e+00
+S Zn Se 1.307592e+00 2.229392e+00 1.747782e+00 4.223622e+01 1.200000e+00 -3.333333e-01 7.049600e+00 6.932325e-01 4.000000e+00 0.000000e+00 0.000000e+00
+S Zn Hg 1.531211e+00 2.025045e+00 1.833708e+00 3.903046e+01 1.200000e+00 -3.333333e-01 7.049600e+00 1.184541e+00 4.000000e+00 0.000000e+00 0.000000e+00
+S Zn S 2.434871e+00 2.423171e+00 1.711097e+00 3.095161e+01 1.200000e+00 -3.333333e-01 7.917000e+00 1.049688e+00 4.000000e+00 0.000000e+00 0.000000e+00
+S Se Cd 1.300376e+00 1.804151e+00 2.124568e+00 3.259006e+01 1.200000e+00 -3.333333e-01 7.049600e+00 1.540087e+00 4.000000e+00 0.000000e+00 0.000000e+00
+S Se Te 1.450015e+00 2.297301e+00 1.726905e+00 3.086266e+01 1.200000e+00 -3.333333e-01 7.049600e+00 7.794685e-01 4.000000e+00 0.000000e+00 0.000000e+00
+S Se Zn 2.208390e+00 2.323783e+00 1.589241e+00 2.500815e+01 1.200000e+00 -3.333333e-01 7.049600e+00 4.643181e-01 4.000000e+00 0.000000e+00 0.000000e+00
+S Se Se 1.307592e+00 2.229392e+00 1.747782e+00 3.250000e+01 1.200000e+00 -3.333333e-01 7.049600e+00 6.932325e-01 4.000000e+00 0.000000e+00 0.000000e+00
+S Se Hg 1.531211e+00 2.025045e+00 1.833708e+00 3.003322e+01 1.200000e+00 -3.333333e-01 7.049600e+00 1.184541e+00 4.000000e+00 0.000000e+00 0.000000e+00
+S Se S 2.434871e+00 2.423171e+00 1.711097e+00 2.381670e+01 1.200000e+00 -3.333333e-01 7.917000e+00 1.049688e+00 4.000000e+00 0.000000e+00 0.000000e+00
+S Hg Cd 1.300376e+00 1.804151e+00 2.124568e+00 3.526684e+01 1.200000e+00 -3.333333e-01 7.049600e+00 1.540087e+00 4.000000e+00 0.000000e+00 0.000000e+00
+S Hg Te 1.450015e+00 2.297301e+00 1.726905e+00 3.339756e+01 1.200000e+00 -3.333333e-01 7.049600e+00 7.794685e-01 4.000000e+00 0.000000e+00 0.000000e+00
+S Hg Zn 2.208390e+00 2.323783e+00 1.589241e+00 2.706220e+01 1.200000e+00 -3.333333e-01 7.049600e+00 4.643181e-01 4.000000e+00 0.000000e+00 0.000000e+00
+S Hg Se 1.307592e+00 2.229392e+00 1.747782e+00 3.516939e+01 1.200000e+00 -3.333333e-01 7.049600e+00 6.932325e-01 4.000000e+00 0.000000e+00 0.000000e+00
+S Hg Hg 1.531211e+00 2.025045e+00 1.833708e+00 3.250000e+01 1.200000e+00 -3.333333e-01 7.049600e+00 1.184541e+00 4.000000e+00 0.000000e+00 0.000000e+00
+S Hg S 2.434871e+00 2.423171e+00 1.711097e+00 2.577288e+01 1.200000e+00 -3.333333e-01 7.917000e+00 1.049688e+00 4.000000e+00 0.000000e+00 0.000000e+00
+S S Cd 1.300376e+00 1.804151e+00 2.124568e+00 4.447203e+01 1.200000e+00 -3.333333e-01 7.049600e+00 1.540087e+00 4.000000e+00 0.000000e+00 0.000000e+00
+S S Te 1.450015e+00 2.297301e+00 1.726905e+00 4.211484e+01 1.200000e+00 -3.333333e-01 7.049600e+00 7.794685e-01 4.000000e+00 0.000000e+00 0.000000e+00
+S S Zn 2.208390e+00 2.323783e+00 1.589241e+00 3.412585e+01 1.200000e+00 -3.333333e-01 7.049600e+00 4.643181e-01 4.000000e+00 0.000000e+00 0.000000e+00
+S S Se 1.307592e+00 2.229392e+00 1.747782e+00 4.434914e+01 1.200000e+00 -3.333333e-01 7.049600e+00 6.932325e-01 4.000000e+00 0.000000e+00 0.000000e+00
+S S Hg 1.531211e+00 2.025045e+00 1.833708e+00 4.098300e+01 1.200000e+00 -3.333333e-01 7.049600e+00 1.184541e+00 4.000000e+00 0.000000e+00 0.000000e+00
+S S S 2.434871e+00 2.423171e+00 1.711097e+00 3.250000e+01 1.200000e+00 -3.333333e-01 7.917000e+00 1.049688e+00 4.000000e+00 0.000000e+00 0.000000e+00
diff --git a/examples/threebody/InP.vashishta b/examples/threebody/InP.vashishta
new file mode 100644
index 0000000000..9fefd4ef19
--- /dev/null
+++ b/examples/threebody/InP.vashishta
@@ -0,0 +1,38 @@
+# DATE: 2015-10-14 CONTRIBUTOR: Aidan Thompson, athomps@sandia.gov CITATION: Branicio, Rino, Gan and Tsuzuki, J. Phys Condensed Matter 21 (2009) 095002
+#
+# Vashishta potential file for InP, Branicio, Rino, Gan and Tsuzuki,
+# J. Phys Condensed Matter 21 (2009) 095002
+#
+# These entries are in LAMMPS "metal" units:
+# H = eV*Angstroms^eta; Zi, Zj = |e| (e = electronic charge);
+# lambda1, lambda4, rc, r0, gamma = Angstroms;
+# D = eV*Angstroms^4; W = eV*Angstroms^6; B = eV;
+# other quantities are unitless
+
+# element1 element2 element3
+# H eta Zi Zj lambda1 D lambda4
+# W rc B gamma r0 C cos(theta)
+
+In In In 273.584 7 -1.21 -1.21 4.5 0.0 2.75
+ 0.0 6.0 0.0 0.0 0.0 0.0 0.0
+
+P P P 1813.06 7 1.21 1.21 4.5 52.7067 2.75
+ 0.0 6.0 0.0 0.0 0.0 0.0 -0.333333333333
+
+In P P 4847.09 9 1.21 -1.21 4.5 26.3533 2.75
+ 270.105 6.0 4.34967 1.0 3.55 7.0 -0.333333333333
+
+P In In 4847.09 9 1.21 -1.21 4.5 26.3533 2.75
+ 270.105 6.0 4.34967 1.0 3.55 7.0 -0.333333333333
+
+In In P 0.0 0.0 0.0 0.0 0.0 0.0 0.0
+ 0.0 0.0 0.0 0.0 0.0 0.0 0.0
+
+In P In 0.0 0.0 0.0 0.0 0.0 0.0 0.0
+ 0.0 0.0 0.0 0.0 0.0 0.0 0.0
+
+P In P 0.0 0.0 0.0 0.0 0.0 0.0 0.0
+ 0.0 0.0 0.0 0.0 0.0 0.0 0.0
+
+P P In 0.0 0.0 0.0 0.0 0.0 0.0 0.0
+ 0.0 0.0 0.0 0.0 0.0 0.0 0.0
diff --git a/examples/threebody/Si.sw b/examples/threebody/Si.sw
new file mode 100644
index 0000000000..db4be100ef
--- /dev/null
+++ b/examples/threebody/Si.sw
@@ -0,0 +1,18 @@
+# DATE: 2007-06-11 CONTRIBUTOR: Aidan Thompson, athomps@sandia.gov CITATION: Stillinger and Weber, Phys Rev B, 31, 5262, (1985)
+# Stillinger-Weber parameters for various elements and mixtures
+# multiple entries can be added to this file, LAMMPS reads the ones it needs
+# these entries are in LAMMPS "metal" units:
+# epsilon = eV; sigma = Angstroms
+# other quantities are unitless
+
+# format of a single entry (one or more lines):
+# element 1, element 2, element 3,
+# epsilon, sigma, a, lambda, gamma, costheta0, A, B, p, q, tol
+
+# Here are the original parameters in metal units, for Silicon from:
+#
+# Stillinger and Weber, Phys. Rev. B, v. 31, p. 5262, (1985)
+#
+
+Si Si Si 2.1683 2.0951 1.80 21.0 1.20 -0.333333333333
+ 7.049556277 0.6022245584 4.0 0.0 0.0
diff --git a/examples/threebody/in.threebody b/examples/threebody/in.threebody
new file mode 100644
index 0000000000..6eff77d6d6
--- /dev/null
+++ b/examples/threebody/in.threebody
@@ -0,0 +1,116 @@
+# Simple regression tests for threebody potentials
+
+# NOTE: These are not intended to represent real materials
+
+units metal
+
+atom_style atomic
+atom_modify map array
+boundary p p p
+atom_modify sort 0 0.0
+
+# temperature
+
+variable t equal 1800.0
+
+# cubic diamond unit cell
+
+variable a equal 5.431
+lattice custom $a &
+ a1 1.0 0.0 0.0 &
+ a2 0.0 1.0 0.0 &
+ a3 0.0 0.0 1.0 &
+ basis 0.0 0.0 0.0 &
+ basis 0.0 0.5 0.5 &
+ basis 0.5 0.0 0.5 &
+ basis 0.5 0.5 0.0 &
+ basis 0.25 0.25 0.25 &
+ basis 0.25 0.75 0.75 &
+ basis 0.75 0.25 0.75 &
+ basis 0.75 0.75 0.25
+
+region myreg block 0 4 &
+ 0 4 &
+ 0 4
+
+create_box 8 myreg
+create_atoms 1 region myreg &
+ basis 1 1 &
+ basis 2 2 &
+ basis 3 3 &
+ basis 4 4 &
+ basis 5 5 &
+ basis 6 6 &
+ basis 7 7 &
+ basis 8 8
+
+mass * 28.06
+
+velocity all create $t 5287287 mom yes rot yes dist gaussian
+
+# Equilibrate using Stillinger-Weber model for silicon
+
+pair_style sw
+pair_coeff * * Si.sw Si Si Si Si Si Si Si Si
+
+thermo 10
+fix 1 all nvt temp $t $t 0.1
+fix_modify 1 energy yes
+timestep 1.0e-3
+neighbor 1.0 bin
+neigh_modify every 1 delay 10 check yes
+run 100
+
+write_restart restart.equil
+
+# Test Stillinger-Weber model for Cd/Te/Zn/Se/Hg/S
+
+clear
+read_restart restart.equil
+
+pair_style sw
+pair_coeff * * CdTeZnSeHgS0.sw Cd Zn Hg Cd Te S Se Te
+
+thermo 10
+fix 1 all nvt temp $t $t 0.1
+fix_modify 1 energy yes
+timestep 1.0e-3
+neighbor 1.0 bin
+neigh_modify every 1 delay 10 check yes
+run 100
+
+# Test Vashishta model for In/P
+
+clear
+read_restart restart.equil
+
+pair_style vashishta
+pair_coeff * * InP.vashishta In In In In P P P P
+
+thermo 10
+fix 1 all nvt temp $t $t 0.1
+fix_modify 1 energy yes
+timestep 1.0e-3
+neighbor 1.0 bin
+neigh_modify every 1 delay 10 check yes
+run 100
+
+# Test Tersoff model for B/N/C
+
+clear
+read_restart restart.equil
+
+variable fac equal 0.6
+change_box all x scale ${fac} y scale ${fac} z scale ${fac} remap
+
+pair_style tersoff
+pair_coeff * * BNC.tersoff N N N C B B C B
+
+thermo 10
+fix 1 all nvt temp $t $t 0.1
+fix_modify 1 energy yes
+timestep 1.0e-3
+neighbor 1.0 bin
+neigh_modify every 1 delay 10 check yes
+run 100
+
diff --git a/potentials/CdTe.sw b/potentials/CdTe.sw
index ed79f7d384..3259657d06 100644
--- a/potentials/CdTe.sw
+++ b/potentials/CdTe.sw
@@ -18,9 +18,11 @@
# sqrt(lambda_ij*epsilon_ij*lambda_ik*epsilon_ik)/lambda_ik, and the
# results are directly entered in this table. Obviously, this
# conversion does not change the two-body parameters epsilon_ijj.
-# All other ik pair parameters are entered on the i*k line, where *
-# can be any species. This is consistent with the requirement of
-# the ik parameter being on the ikk line.
+
+# The twobody ik pair parameters are entered on the i*k lines, where *
+# can be any species. This is consistent with the LAMMPS requirement
+# that twobody ik parameters be defined on the ikk line. Entries on all
+# the other i*k lines are ignored by LAMMPS
# These entries are in LAMMPS "metal" units: epsilon = eV;
# sigma = Angstroms; other quantities are unitless
diff --git a/potentials/CdTeZnSeHgS0.sw b/potentials/CdTeZnSeHgS0.sw
index cae54097dd..d6f05d41df 100644
--- a/potentials/CdTeZnSeHgS0.sw
+++ b/potentials/CdTeZnSeHgS0.sw
@@ -1,9 +1,17 @@
### DATE: 2013-08-09 CONTRIBUTOR: X. W. Zhou, xzhou@sandia.gov, CITATION: Zhou, Ward, Martin, van Swol, Cruz-Campa, and D. Zubia, Phys. Rev. B, 88, 085309 (2013).
#
-# Note that the way the parameters can be entered is not unique. As one way, we assume that eps_ijk is equal to eps_ik and lambda_ijk is equal to
-# sqrt(lambda_ij*eps_ij*lambda_ik*eps_ik)/eps_ik, and all other parameters in the ijk line are for ik.
+# Note that the way the parameters can be entered is not unique.
+# As one way, we assume that eps_ijk is equal to eps_ik and
+# lambda_ijk is equal to sqrt(lambda_ij*eps_ij*lambda_ik*eps_ik)/eps_ik,
+# and all other parameters in the ijk line are for ik.
+#
+# The twobody ik pair parameters are entered on the i*k lines, where *
+# can be any species. This is consistent with the LAMMPS requirement
+# that twobody ik parameters be defined on the ikk line. Entries on all
+# the other i*k lines are ignored by LAMMPS
#
-# These entries are in LAMMPS "metal" units: epsilon = eV; sigma = Angstroms; other quantities are unitless;
+# These entries are in LAMMPS "metal" units: epsilon = eV;
+# sigma = Angstroms; other quantities are unitless
#
# cutoff distance = 4.632
# eps sigma a lambda gamma cos(theta) A B p q tol
diff --git a/src/KOKKOS/fix_qeq_reax_kokkos.cpp b/src/KOKKOS/fix_qeq_reax_kokkos.cpp
index d39d5cc84a..0c0039a18a 100644
--- a/src/KOKKOS/fix_qeq_reax_kokkos.cpp
+++ b/src/KOKKOS/fix_qeq_reax_kokkos.cpp
@@ -37,7 +37,7 @@
#include "math_const.h"
#include "memory.h"
#include "error.h"
-#include "pair_kokkos.h"
+#include "pair_reax_c_kokkos.h"
using namespace LAMMPS_NS;
using namespace FixConst;
@@ -62,6 +62,11 @@ FixQEqReaxKokkos<DeviceType>::FixQEqReaxKokkos(LAMMPS *lmp, int narg, char **arg
nmax = nmax = m_cap = 0;
allocated_flag = 0;
+
+ reaxc = (PairReaxC *) force->pair_match("reax/c/kk",1);
+ use_pair_list = 0;
+ if (reaxc->execution_space == this->execution_space)
+ use_pair_list = 1;
}
/* ---------------------------------------------------------------------- */
@@ -80,29 +85,54 @@ void FixQEqReaxKokkos<DeviceType>::init()
atomKK->k_q.modify<LMPHostType>();
atomKK->k_q.sync<LMPDeviceType>();
- FixQEqReax::init();
+ //FixQEqReax::init();
+ {
+ if (!atom->q_flag) error->all(FLERR,"Fix qeq/reax requires atom attribute q");
+
+ ngroup = group->count(igroup);
+ if (ngroup == 0) error->all(FLERR,"Fix qeq/reax group has no atoms");
+
+ // need a half neighbor list w/ Newton off and ghost neighbors
+ // built whenever re-neighboring occurs
+
+ if (!use_pair_list) {
+ int irequest = neighbor->request(this,instance_me);
+ neighbor->requests[irequest]->pair = 0;
+ neighbor->requests[irequest]->fix = 1;
+ neighbor->requests[irequest]->newton = 2;
+ neighbor->requests[irequest]->ghost = 1;
+ }
+
+ init_shielding();
+ init_taper();
+
+ if (strstr(update->integrate_style,"respa"))
+ nlevels_respa = ((Respa *) update->integrate)->nlevels;
+ }
neighflag = lmp->kokkos->neighflag;
- int irequest = neighbor->nrequest - 1;
-
- neighbor->requests[irequest]->
- kokkos_host = Kokkos::Impl::is_same<DeviceType,LMPHostType>::value &&
- !Kokkos::Impl::is_same<DeviceType,LMPDeviceType>::value;
- neighbor->requests[irequest]->
- kokkos_device = Kokkos::Impl::is_same<DeviceType,LMPDeviceType>::value;
-
- if (neighflag == FULL) {
- neighbor->requests[irequest]->fix = 1;
- neighbor->requests[irequest]->pair = 0;
- neighbor->requests[irequest]->full = 1;
- neighbor->requests[irequest]->half = 0;
- neighbor->requests[irequest]->full_cluster = 0;
- } else { //if (neighflag == HALF || neighflag == HALFTHREAD)
- neighbor->requests[irequest]->fix = 1;
- neighbor->requests[irequest]->full = 0;
- neighbor->requests[irequest]->half = 1;
- neighbor->requests[irequest]->full_cluster = 0;
- neighbor->requests[irequest]->ghost = 1;
+
+ if (!use_pair_list) {
+ int irequest = neighbor->nrequest - 1;
+ neighbor->requests[irequest]->
+ kokkos_host = Kokkos::Impl::is_same<DeviceType,LMPHostType>::value &&
+ !Kokkos::Impl::is_same<DeviceType,LMPDeviceType>::value;
+ neighbor->requests[irequest]->
+ kokkos_device = Kokkos::Impl::is_same<DeviceType,LMPDeviceType>::value;
+
+ if (neighflag == FULL) {
+ neighbor->requests[irequest]->fix = 1;
+ neighbor->requests[irequest]->pair = 0;
+ neighbor->requests[irequest]->full = 1;
+ neighbor->requests[irequest]->half = 0;
+ neighbor->requests[irequest]->full_cluster = 0;
+ } else { //if (neighflag == HALF || neighflag == HALFTHREAD)
+ neighbor->requests[irequest]->fix = 1;
+ neighbor->requests[irequest]->full = 0;
+ neighbor->requests[irequest]->half = 1;
+ neighbor->requests[irequest]->full_cluster = 0;
+ neighbor->requests[irequest]->ghost = 1;
+ }
}
int ntypes = atom->ntypes;
@@ -213,6 +243,8 @@ void FixQEqReaxKokkos<DeviceType>::pre_force(int vflag)
k_shield.template sync<DeviceType>();
k_tap.template sync<DeviceType>();
+ if (use_pair_list)
+ list = reaxc->list;
NeighListKokkos<DeviceType>* k_list = static_cast<NeighListKokkos<DeviceType>*>(list);
d_numneigh = k_list->d_numneigh;
d_neighbors = k_list->d_neighbors;
diff --git a/src/KOKKOS/fix_qeq_reax_kokkos.h b/src/KOKKOS/fix_qeq_reax_kokkos.h
index eca0d761b7..fcfc28fa74 100644
--- a/src/KOKKOS/fix_qeq_reax_kokkos.h
+++ b/src/KOKKOS/fix_qeq_reax_kokkos.h
@@ -146,7 +146,7 @@ class FixQEqReaxKokkos : public FixQEqReax {
double memory_usage();
protected:
- int inum;
+ int inum,use_pair_list;
int allocated_flag;
typedef Kokkos::DualView<int***,DeviceType> tdual_int_1d;
diff --git a/src/KOKKOS/fix_reaxc_species_kokkos.cpp b/src/KOKKOS/fix_reaxc_species_kokkos.cpp
index d3de0c998c..17b42174c6 100644
--- a/src/KOKKOS/fix_reaxc_species_kokkos.cpp
+++ b/src/KOKKOS/fix_reaxc_species_kokkos.cpp
@@ -71,12 +71,6 @@ void FixReaxCSpeciesKokkos::init()
"pair_style reax/c/kk");
FixReaxCSpecies::init();
-
- int irequest = neighbor->request(this,instance_me);
- neighbor->requests[irequest]->pair = 0;
- neighbor->requests[irequest]->fix = 1;
- neighbor->requests[irequest]->newton = 2;
- neighbor->requests[irequest]->ghost = 1;
}
/* ---------------------------------------------------------------------- */
@@ -86,12 +80,20 @@ void FixReaxCSpeciesKokkos::FindMolecule()
int i,j,ii,jj,inum,itype,jtype,loop,looptot;
int change,done,anychange;
int *mask = atom->mask;
- int *ilist;
double bo_tmp,bo_cut;
double **spec_atom = f_SPECBOND->array_atom;
- inum = list->inum;
- ilist = list->ilist;
+ inum = reaxc->list->inum;
+ typename ArrayTypes<LMPHostType>::t_int_1d ilist;
+ if (reaxc->execution_space == Host) {
+ NeighListKokkos<LMPHostType>* k_list = static_cast<NeighListKokkos<LMPHostType>*>(reaxc->list);
+ k_list->k_ilist.sync<LMPHostType>();
+ ilist = k_list->k_ilist.h_view;
+ } else {
+ NeighListKokkos<LMPDeviceType>* k_list = static_cast<NeighListKokkos<LMPDeviceType>*>(reaxc->list);
+ k_list->k_ilist.sync<LMPHostType>();
+ ilist = k_list->k_ilist.h_view;
+ }
for (ii = 0; ii < inum; ii++) {
i = ilist[ii];
diff --git a/src/KOKKOS/fix_wall_reflect_kokkos.h b/src/KOKKOS/fix_wall_reflect_kokkos.h
index f5e28796fd..d59088b592 100644
--- a/src/KOKKOS/fix_wall_reflect_kokkos.h
+++ b/src/KOKKOS/fix_wall_reflect_kokkos.h
@@ -41,8 +41,6 @@ class FixWallReflectKokkos : public FixWallReflect {
void operator()(TagFixWallReflectPostIntegrate, const int&) const;
protected:
- class AtomKokkos *atomKK;
-
typename AT::t_x_array x;
typename AT::t_v_array v;
typename AT::t_int_1d_randomread mask;
diff --git a/src/KOKKOS/neigh_full_kokkos.h b/src/KOKKOS/neigh_full_kokkos.h
index 883b6dc115..9125b5fbe2 100644
--- a/src/KOKKOS/neigh_full_kokkos.h
+++ b/src/KOKKOS/neigh_full_kokkos.h
@@ -164,6 +164,7 @@ if (GHOST) {
list->gnum = 0;
}
+ list->k_ilist.template modify<DeviceType>();
}
/* ---------------------------------------------------------------------- */
diff --git a/src/KOKKOS/neigh_list_kokkos.cpp b/src/KOKKOS/neigh_list_kokkos.cpp
index 5fe796f84d..cbba2120bd 100644
--- a/src/KOKKOS/neigh_list_kokkos.cpp
+++ b/src/KOKKOS/neigh_list_kokkos.cpp
@@ -49,8 +49,9 @@ void NeighListKokkos<Device>::grow(int nmax)
if (nmax <= maxatoms) return;
maxatoms = nmax;
- d_ilist =
- typename ArrayTypes<Device>::t_int_1d("neighlist:ilist",maxatoms);
+ k_ilist =
+ DAT::tdual_int_1d("neighlist:ilist",maxatoms);
+ d_ilist = k_ilist.view<Device>();
d_numneigh =
typename ArrayTypes<Device>::t_int_1d("neighlist:numneigh",maxatoms);
d_neighbors =
diff --git a/src/KOKKOS/neigh_list_kokkos.h b/src/KOKKOS/neigh_list_kokkos.h
index 5200b24595..85f0f38d2c 100644
--- a/src/KOKKOS/neigh_list_kokkos.h
+++ b/src/KOKKOS/neigh_list_kokkos.h
@@ -71,7 +71,8 @@ public:
void clean_copy();
void grow(int nmax);
typename ArrayTypes<Device>::t_neighbors_2d d_neighbors;
- typename ArrayTypes<Device>::t_int_1d d_ilist; // local indices of I atoms
+ typename DAT::tdual_int_1d k_ilist; // local indices of I atoms
+ typename ArrayTypes<Device>::t_int_1d d_ilist;
typename ArrayTypes<Device>::t_int_1d d_numneigh; // # of J neighs for each I
typename ArrayTypes<Device>::t_int_1d d_stencil; // # of J neighs for each I
typename ArrayTypes<LMPHostType>::t_int_1d h_stencil; // # of J neighs per I
diff --git a/src/KOKKOS/pair_reax_c_kokkos.h b/src/KOKKOS/pair_reax_c_kokkos.h
index 204ad7732f..2d746dee0d 100644
--- a/src/KOKKOS/pair_reax_c_kokkos.h
+++ b/src/KOKKOS/pair_reax_c_kokkos.h
@@ -426,8 +426,6 @@ class PairReaxCKokkos : public PairReaxC {
typename AT::t_ffloat_2d_dl d_sum_ovun;
typename AT::t_ffloat_2d_dl d_dBOpx, d_dBOpy, d_dBOpz;
- class AtomKokkos *atomKK;
-
int neighflag,newton_pair, maxnumneigh, maxhb, maxbo;
int nlocal,nall,eflag,vflag;
F_FLOAT cut_nbsq, cut_hbsq, cut_bosq, bo_cut, thb_cut, thb_cutsq;
diff --git a/src/KOKKOS/pair_sw_kokkos.cpp b/src/KOKKOS/pair_sw_kokkos.cpp
index bb81a29c0a..d2cda316be 100644
--- a/src/KOKKOS/pair_sw_kokkos.cpp
+++ b/src/KOKKOS/pair_sw_kokkos.cpp
@@ -121,6 +121,18 @@ void PairSWKokkos<DeviceType>::compute(int eflag_in, int vflag_in)
EV_FLOAT ev;
EV_FLOAT ev_all;
+ // build short neighbor list
+
+ int max_neighs = d_neighbors.dimension_1();
+
+ if ((d_neighbors_short.dimension_1() != max_neighs) ||
+ (d_neighbors_short.dimension_0() != ignum)) {
+ d_neighbors_short = Kokkos::View<int**,DeviceType>("SW::neighbors_short",ignum,max_neighs);
+ }
+ if (d_numneigh_short.dimension_0()!=ignum)
+ d_numneigh_short = Kokkos::View<int*,DeviceType>("SW::numneighs_short",ignum);
+ Kokkos::parallel_for(Kokkos::RangePolicy<DeviceType,TagPairSWComputeShortNeigh>(0,neighflag==FULL?ignum:inum), *this);
+
// loop over neighbor list of my atoms
if (neighflag == HALF) {
@@ -174,6 +186,38 @@ void PairSWKokkos<DeviceType>::compute(int eflag_in, int vflag_in)
copymode = 0;
}
+
+/* ---------------------------------------------------------------------- */
+
+template<class DeviceType>
+KOKKOS_INLINE_FUNCTION
+void PairSWKokkos<DeviceType>::operator()(TagPairSWComputeShortNeigh, const int& ii) const {
+ const int i = d_ilist[ii];
+ const X_FLOAT xtmp = x(i,0);
+ const X_FLOAT ytmp = x(i,1);
+ const X_FLOAT ztmp = x(i,2);
+
+ const int jnum = d_numneigh[i];
+ int inside = 0;
+ for (int jj = 0; jj < jnum; jj++) {
+ int j = d_neighbors(i,jj);
+ j &= NEIGHMASK;
+
+ const X_FLOAT delx = xtmp - x(j,0);
+ const X_FLOAT dely = ytmp - x(j,1);
+ const X_FLOAT delz = ztmp - x(j,2);
+ const F_FLOAT rsq = delx*delx + dely*dely + delz*delz;
+
+ if (rsq < cutmax*cutmax) {
+ d_neighbors_short(i,inside) = j;
+ inside++;
+ }
+ }
+ d_numneigh_short(i) = inside;
+}
+
+/* ---------------------------------------------------------------------- */
+
template<class DeviceType>
template<int NEIGHFLAG, int EVFLAG>
KOKKOS_INLINE_FUNCTION
@@ -196,14 +240,14 @@ void PairSWKokkos<DeviceType>::operator()(TagPairSWComputeHalf<NEIGHFLAG,EVFLAG>
// two-body interactions, skip half of them
- const int jnum = d_numneigh[i];
+ const int jnum = d_numneigh_short[i];
F_FLOAT fxtmpi = 0.0;
F_FLOAT fytmpi = 0.0;
F_FLOAT fztmpi = 0.0;
for (int jj = 0; jj < jnum; jj++) {
- int j = d_neighbors(i,jj);
+ int j = d_neighbors_short(i,jj);
j &= NEIGHMASK;
const tagint jtag = tag[j];
@@ -245,7 +289,7 @@ void PairSWKokkos<DeviceType>::operator()(TagPairSWComputeHalf<NEIGHFLAG,EVFLAG>
const int jnumm1 = jnum - 1;
for (int jj = 0; jj < jnumm1; jj++) {
- int j = d_neighbors(i,jj);
+ int j = d_neighbors_short(i,jj);
j &= NEIGHMASK;
const int jtype = d_map[type[j]];
const int ijparam = d_elem2param(itype,jtype,jtype);
@@ -260,7 +304,7 @@ void PairSWKokkos<DeviceType>::operator()(TagPairSWComputeHalf<NEIGHFLAG,EVFLAG>
F_FLOAT fztmpj = 0.0;
for (int kk = jj+1; kk < jnum; kk++) {
- int k = d_neighbors(i,kk);
+ int k = d_neighbors_short(i,kk);
k &= NEIGHMASK;
const int ktype = d_map[type[k]];
const int ikparam = d_elem2param(itype,ktype,ktype);
@@ -331,14 +375,14 @@ void PairSWKokkos<DeviceType>::operator()(TagPairSWComputeFullA<NEIGHFLAG,EVFLAG
// two-body interactions
- const int jnum = d_numneigh[i];
+ const int jnum = d_numneigh_short[i];
F_FLOAT fxtmpi = 0.0;
F_FLOAT fytmpi = 0.0;
F_FLOAT fztmpi = 0.0;
for (int jj = 0; jj < jnum; jj++) {
- int j = d_neighbors(i,jj);
+ int j = d_neighbors_short(i,jj);
j &= NEIGHMASK;
const tagint jtag = tag[j];
@@ -368,7 +412,7 @@ void PairSWKokkos<DeviceType>::operator()(TagPairSWComputeFullA<NEIGHFLAG,EVFLAG
const int jnumm1 = jnum - 1;
for (int jj = 0; jj < jnumm1; jj++) {
- int j = d_neighbors(i,jj);
+ int j = d_neighbors_short(i,jj);
j &= NEIGHMASK;
const int jtype = d_map[type[j]];
const int ijparam = d_elem2param(itype,jtype,jtype);
@@ -380,7 +424,7 @@ void PairSWKokkos<DeviceType>::operator()(TagPairSWComputeFullA<NEIGHFLAG,EVFLAG
if (rsq1 >= d_params[ijparam].cutsq) continue;
for (int kk = jj+1; kk < jnum; kk++) {
- int k = d_neighbors(i,kk);
+ int k = d_neighbors_short(i,kk);
k &= NEIGHMASK;
const int ktype = d_map[type[k]];
const int ikparam = d_elem2param(itype,ktype,ktype);
@@ -438,14 +482,14 @@ void PairSWKokkos<DeviceType>::operator()(TagPairSWComputeFullB<NEIGHFLAG,EVFLAG
const X_FLOAT ytmpi = x(i,1);
const X_FLOAT ztmpi = x(i,2);
- const int jnum = d_numneigh[i];
+ const int jnum = d_numneigh_short[i];
F_FLOAT fxtmpi = 0.0;
F_FLOAT fytmpi = 0.0;
F_FLOAT fztmpi = 0.0;
for (int jj = 0; jj < jnum; jj++) {
- int j = d_neighbors(i,jj);
+ int j = d_neighbors_short(i,jj);
j &= NEIGHMASK;
if (j >= nlocal) continue;
const int jtype = d_map[type[j]];
@@ -461,10 +505,10 @@ void PairSWKokkos<DeviceType>::operator()(TagPairSWComputeFullB<NEIGHFLAG,EVFLAG
if (rsq1 >= d_params[jiparam].cutsq) continue;
- const int j_jnum = d_numneigh[j];
+ const int j_jnum = d_numneigh_short[j];
for (int kk = 0; kk < j_jnum; kk++) {
- int k = d_neighbors(j,kk);
+ int k = d_neighbors_short(j,kk);
k &= NEIGHMASK;
if (k == i) continue;
const int ktype = d_map[type[k]];
diff --git a/src/KOKKOS/pair_sw_kokkos.h b/src/KOKKOS/pair_sw_kokkos.h
index cec240e064..c722d9d52c 100644
--- a/src/KOKKOS/pair_sw_kokkos.h
+++ b/src/KOKKOS/pair_sw_kokkos.h
@@ -34,6 +34,8 @@ struct TagPairSWComputeFullA{};
template<int NEIGHFLAG, int EVFLAG>
struct TagPairSWComputeFullB{};
+struct TagPairSWComputeShortNeigh{};
+
namespace LAMMPS_NS {
template<class DeviceType>
@@ -75,6 +77,9 @@ class PairSWKokkos : public PairSW {
KOKKOS_INLINE_FUNCTION
void operator()(TagPairSWComputeFullB<NEIGHFLAG,EVFLAG>, const int&) const;
+ KOKKOS_INLINE_FUNCTION
+ void operator()(TagPairSWComputeShortNeigh, const int&) const;
+
template<int NEIGHFLAG>
KOKKOS_INLINE_FUNCTION
void ev_tally(EV_FLOAT &ev, const int &i, const int &j,
@@ -136,6 +141,9 @@ class PairSWKokkos : public PairSW {
int nlocal,nall,eflag,vflag;
int inum;
+ Kokkos::View<int**,DeviceType> d_neighbors_short;
+ Kokkos::View<int*,DeviceType> d_numneigh_short;
+
friend void pair_virial_fdotr_compute<PairSWKokkos>(PairSWKokkos*);
};
diff --git a/src/KOKKOS/pair_tersoff_kokkos.cpp b/src/KOKKOS/pair_tersoff_kokkos.cpp
index cf9b510ed6..162661430b 100644
--- a/src/KOKKOS/pair_tersoff_kokkos.cpp
+++ b/src/KOKKOS/pair_tersoff_kokkos.cpp
@@ -12,7 +12,7 @@
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
- Contributing author: Ray Shan (SNL)
+ Contributing author: Ray Shan (SNL) and Christian Trott (SNL)
------------------------------------------------------------------------- */
#include <math.h>
@@ -191,7 +191,7 @@ void PairTersoffKokkos<DeviceType>::compute(int eflag_in, int vflag_in)
nall = atom->nlocal + atom->nghost;
newton_pair = force->newton_pair;
- const int inum = list->inum;
+ inum = list->inum;
const int ignum = inum + list->gnum;
NeighListKokkos<DeviceType>* k_list = static_cast<NeighListKokkos<DeviceType>*>(list);
d_numneigh = k_list->d_numneigh;
@@ -204,6 +204,18 @@ void PairTersoffKokkos<DeviceType>::compute(int eflag_in, int vflag_in)
EV_FLOAT ev;
EV_FLOAT ev_all;
+ // build short neighbor list
+
+ int max_neighs = d_neighbors.dimension_1();
+
+ if ((d_neighbors_short.dimension_1() != max_neighs) ||
+ (d_neighbors_short.dimension_0() != ignum)) {
+ d_neighbors_short = Kokkos::View<int**,DeviceType>("Tersoff::neighbors_short",ignum,max_neighs);
+ }
+ if (d_numneigh_short.dimension_0()!=ignum)
+ d_numneigh_short = Kokkos::View<int*,DeviceType>("Tersoff::numneighs_short",ignum);
+ Kokkos::parallel_for(Kokkos::RangePolicy<DeviceType,TagPairTersoffComputeShortNeigh>(0,neighflag==FULL?ignum:inum), *this);
+
if (neighflag == HALF) {
if (evflag)
Kokkos::parallel_reduce(Kokkos::RangePolicy<DeviceType, TagPairTersoffComputeHalf<HALF,1> >(0,inum),*this,ev);
@@ -257,6 +269,35 @@ void PairTersoffKokkos<DeviceType>::compute(int eflag_in, int vflag_in)
/* ---------------------------------------------------------------------- */
+template<class DeviceType>
+KOKKOS_INLINE_FUNCTION
+void PairTersoffKokkos<DeviceType>::operator()(TagPairTersoffComputeShortNeigh, const int& ii) const {
+ const int i = d_ilist[ii];
+ const X_FLOAT xtmp = x(i,0);
+ const X_FLOAT ytmp = x(i,1);
+ const X_FLOAT ztmp = x(i,2);
+
+ const int jnum = d_numneigh[i];
+ int inside = 0;
+ for (int jj = 0; jj < jnum; jj++) {
+ int j = d_neighbors(i,jj);
+ j &= NEIGHMASK;
+
+ const X_FLOAT delx = xtmp - x(j,0);
+ const X_FLOAT dely = ytmp - x(j,1);
+ const X_FLOAT delz = ztmp - x(j,2);
+ const F_FLOAT rsq = delx*delx + dely*dely + delz*delz;
+
+ if (rsq < cutmax*cutmax) {
+ d_neighbors_short(i,inside) = j;
+ inside++;
+ }
+ }
+ d_numneigh_short(i) = inside;
+}
+
+/* ---------------------------------------------------------------------- */
+
template<class DeviceType>
template<int NEIGHFLAG, int EVFLAG>
KOKKOS_INLINE_FUNCTION
@@ -273,21 +314,22 @@ void PairTersoffKokkos<DeviceType>::operator()(TagPairTersoffComputeHalf<NEIGHFL
const int itype = type(i);
const int itag = tag(i);
- int j,k,jj,kk,jtag,jtype,ktype;
- F_FLOAT rsq1, cutsq1, rsq2, cutsq2, rij, rik, bo_ij;
F_FLOAT fi[3], fj[3], fk[3];
- X_FLOAT delx1, dely1, delz1, delx2, dely2, delz2;
//const AtomNeighborsConst d_neighbors_i = k_list.get_neighbors_const(i);
- const int jnum = d_numneigh[i];
+ const int jnum = d_numneigh_short[i];
// repulsive
- for (jj = 0; jj < jnum; jj++) {
- j = d_neighbors(i,jj);
+ F_FLOAT f_x = 0.0;
+ F_FLOAT f_y = 0.0;
+ F_FLOAT f_z = 0.0;
+
+ for (int jj = 0; jj < jnum; jj++) {
+ int j = d_neighbors_short(i,jj);
j &= NEIGHMASK;
- jtype = type(j);
- jtag = tag(j);
+ const int jtype = type(j);
+ const int jtag = tag(j);
if (itag > jtag) {
if ((itag+jtag) % 2 == 0) continue;
@@ -315,9 +357,9 @@ void PairTersoffKokkos<DeviceType>::operator()(TagPairTersoffComputeHalf<NEIGHFL
(tmp_fcd - tmp_fce*paramskk(itype,jtype,jtype).lam1) / r;
const F_FLOAT eng = tmp_fce * paramskk(itype,jtype,jtype).biga * tmp_exp;
- a_f(i,0) += delx*frep;
- a_f(i,1) += dely*frep;
- a_f(i,2) += delz*frep;
+ f_x += delx*frep;
+ f_y += dely*frep;
+ f_z += delz*frep;
a_f(j,0) -= delx*frep;
a_f(j,1) -= dely*frep;
a_f(j,2) -= delz*frep;
@@ -330,35 +372,35 @@ void PairTersoffKokkos<DeviceType>::operator()(TagPairTersoffComputeHalf<NEIGHFL
// attractive: bond order
- for (jj = 0; jj < jnum; jj++) {
- j = d_neighbors(i,jj);
+ for (int jj = 0; jj < jnum; jj++) {
+ int j = d_neighbors_short(i,jj);
j &= NEIGHMASK;
- jtype = type(j);
+ const int jtype = type(j);
- delx1 = xtmp - x(j,0);
- dely1 = ytmp - x(j,1);
- delz1 = ztmp - x(j,2);
- rsq1 = delx1*delx1 + dely1*dely1 + delz1*delz1;
- cutsq1 = paramskk(itype,jtype,jtype).cutsq;
+ const F_FLOAT delx1 = xtmp - x(j,0);
+ const F_FLOAT dely1 = ytmp - x(j,1);
+ const F_FLOAT delz1 = ztmp - x(j,2);
+ const F_FLOAT rsq1 = delx1*delx1 + dely1*dely1 + delz1*delz1;
+ const F_FLOAT cutsq1 = paramskk(itype,jtype,jtype).cutsq;
- bo_ij = 0.0;
+ F_FLOAT bo_ij = 0.0;
if (rsq1 > cutsq1) continue;
- rij = sqrt(rsq1);
+ const F_FLOAT rij = sqrt(rsq1);
- for (kk = 0; kk < jnum; kk++) {
+ for (int kk = 0; kk < jnum; kk++) {
if (jj == kk) continue;
- k = d_neighbors(i,kk);
+ int k = d_neighbors_short(i,kk);
k &= NEIGHMASK;
- ktype = type(k);
+ const int ktype = type(k);
- delx2 = xtmp - x(k,0);
- dely2 = ytmp - x(k,1);
- delz2 = ztmp - x(k,2);
- rsq2 = delx2*delx2 + dely2*dely2 + delz2*delz2;
- cutsq2 = paramskk(itype,jtype,ktype).cutsq;
+ const F_FLOAT delx2 = xtmp - x(k,0);
+ const F_FLOAT dely2 = ytmp - x(k,1);
+ const F_FLOAT delz2 = ztmp - x(k,2);
+ const F_FLOAT rsq2 = delx2*delx2 + dely2*dely2 + delz2*delz2;
+ const F_FLOAT cutsq2 = paramskk(itype,jtype,ktype).cutsq;
if (rsq2 > cutsq2) continue;
- rik = sqrt(rsq2);
+ const F_FLOAT rik = sqrt(rsq2);
bo_ij += bondorder(itype,jtype,ktype,rij,delx1,dely1,delz1,rik,delx2,dely2,delz2);
}
@@ -369,58 +411,64 @@ void PairTersoffKokkos<DeviceType>::operator()(TagPairTersoffComputeHalf<NEIGHFL
const F_FLOAT bij = ters_bij_k(itype,jtype,jtype,bo_ij);
const F_FLOAT fatt = -0.5*bij * dfa / rij;
const F_FLOAT prefactor = 0.5*fa * ters_dbij(itype,jtype,jtype,bo_ij);
- const F_FLOAT eng = 0.5*bij * fa;
- a_f(i,0) += delx1*fatt;
- a_f(i,1) += dely1*fatt;
- a_f(i,2) += delz1*fatt;
- a_f(j,0) -= delx1*fatt;
- a_f(j,1) -= dely1*fatt;
- a_f(j,2) -= delz1*fatt;
+ f_x += delx1*fatt;
+ f_y += dely1*fatt;
+ f_z += delz1*fatt;
+ F_FLOAT fj_x = -delx1*fatt;
+ F_FLOAT fj_y = -dely1*fatt;
+ F_FLOAT fj_z = -delz1*fatt;
if (EVFLAG) {
+ const F_FLOAT eng = 0.5*bij * fa;
if (eflag) ev.evdwl += eng;
if (vflag_either || eflag_atom)
- this->template ev_tally<NEIGHFLAG>(ev,i,j,eng,fatt,delx1,dely1,delz1);
+ this->template ev_tally<NEIGHFLAG>(ev,i,j,eng,fatt,delx1,dely1,delz1);
}
// attractive: three-body force
- for (kk = 0; kk < jnum; kk++) {
+ for (int kk = 0; kk < jnum; kk++) {
if (jj == kk) continue;
- k = d_neighbors(i,kk);
+ int k = d_neighbors_short(i,kk);
k &= NEIGHMASK;
- ktype = type(k);
+ const int ktype = type(k);
- delx2 = xtmp - x(k,0);
- dely2 = ytmp - x(k,1);
- delz2 = ztmp - x(k,2);
- rsq2 = delx2*delx2 + dely2*dely2 + delz2*delz2;
- cutsq2 = paramskk(itype,jtype,ktype).cutsq;
+ const F_FLOAT delx2 = xtmp - x(k,0);
+ const F_FLOAT dely2 = ytmp - x(k,1);
+ const F_FLOAT delz2 = ztmp - x(k,2);
+ const F_FLOAT rsq2 = delx2*delx2 + dely2*dely2 + delz2*delz2;
+ const F_FLOAT cutsq2 = paramskk(itype,jtype,ktype).cutsq;
if (rsq2 > cutsq2) continue;
- rik = sqrt(rsq2);
+ const F_FLOAT rik = sqrt(rsq2);
ters_dthb(itype,jtype,ktype,prefactor,rij,delx1,dely1,delz1,
- rik,delx2,dely2,delz2,fi,fj,fk);
-
- a_f(i,0) += fi[0];
- a_f(i,1) += fi[1];
- a_f(i,2) += fi[2];
- a_f(j,0) += fj[0];
- a_f(j,1) += fj[1];
- a_f(j,2) += fj[2];
+ rik,delx2,dely2,delz2,fi,fj,fk);
+
+ f_x += fi[0];
+ f_y += fi[1];
+ f_z += fi[2];
+ fj_x += fj[0];
+ fj_y += fj[1];
+ fj_z += fj[2];
a_f(k,0) += fk[0];
a_f(k,1) += fk[1];
a_f(k,2) += fk[2];
if (vflag_atom) {
- F_FLOAT delrij[3], delrik[3];
- delrij[0] = -delx1; delrij[1] = -dely1; delrij[2] = -delz1;
- delrik[0] = -delx2; delrik[1] = -dely2; delrik[2] = -delz2;
- if (vflag_either) this->template v_tally3<NEIGHFLAG>(ev,i,j,k,fj,fk,delrij,delrik);
+ F_FLOAT delrij[3], delrik[3];
+ delrij[0] = -delx1; delrij[1] = -dely1; delrij[2] = -delz1;
+ delrik[0] = -delx2; delrik[1] = -dely2; delrik[2] = -delz2;
+ if (vflag_either) this->template v_tally3<NEIGHFLAG>(ev,i,j,k,fj,fk,delrij,delrik);
}
}
+ a_f(j,0) += fj_x;
+ a_f(j,1) += fj_y;
+ a_f(j,2) += fj_z;
}
+ a_f(i,0) += f_x;
+ a_f(i,1) += f_y;
+ a_f(i,2) += f_z;
}
template<class DeviceType>
@@ -450,12 +498,15 @@ void PairTersoffKokkos<DeviceType>::operator()(TagPairTersoffComputeFullA<NEIGHF
X_FLOAT delx1, dely1, delz1, delx2, dely2, delz2;
//const AtomNeighborsConst d_neighbors_i = k_list.get_neighbors_const(i);
- const int jnum = d_numneigh[i];
+ const int jnum = d_numneigh_short[i];
// repulsive
+ F_FLOAT f_x = 0.0;
+ F_FLOAT f_y = 0.0;
+ F_FLOAT f_z = 0.0;
for (jj = 0; jj < jnum; jj++) {
- j = d_neighbors(i,jj);
+ j = d_neighbors_short(i,jj);
j &= NEIGHMASK;
const int jtype = type(j);
@@ -475,9 +526,9 @@ void PairTersoffKokkos<DeviceType>::operator()(TagPairTersoffComputeFullA<NEIGHF
(tmp_fcd - tmp_fce*paramskk(itype,jtype,jtype).lam1) / r;
const F_FLOAT eng = tmp_fce * paramskk(itype,jtype,jtype).biga * tmp_exp;
- f(i,0) += delx*frep;
- f(i,1) += dely*frep;
- f(i,2) += delz*frep;
+ f_x += delx*frep;
+ f_y += dely*frep;
+ f_z += delz*frep;
if (EVFLAG) {
if (eflag)
@@ -490,7 +541,7 @@ void PairTersoffKokkos<DeviceType>::operator()(TagPairTersoffComputeFullA<NEIGHF
// attractive: bond order
for (jj = 0; jj < jnum; jj++) {
- j = d_neighbors(i,jj);
+ j = d_neighbors_short(i,jj);
j &= NEIGHMASK;
jtype = type(j);
@@ -506,7 +557,7 @@ void PairTersoffKokkos<DeviceType>::operator()(TagPairTersoffComputeFullA<NEIGHF
for (kk = 0; kk < jnum; kk++) {
if (jj == kk) continue;
- k = d_neighbors(i,kk);
+ k = d_neighbors_short(i,kk);
k &= NEIGHMASK;
ktype = type(k);
@@ -530,9 +581,9 @@ void PairTersoffKokkos<DeviceType>::operator()(TagPairTersoffComputeFullA<NEIGHF
const F_FLOAT prefactor = 0.5*fa * ters_dbij(itype,jtype,jtype,bo_ij);
const F_FLOAT eng = 0.5*bij * fa;
- f(i,0) += delx1*fatt;
- f(i,1) += dely1*fatt;
- f(i,2) += delz1*fatt;
+ f_x += delx1*fatt;
+ f_y += dely1*fatt;
+ f_z += delz1*fatt;
if (EVFLAG) {
if (eflag) ev.evdwl += 0.5*eng;
@@ -544,7 +595,7 @@ void PairTersoffKokkos<DeviceType>::operator()(TagPairTersoffComputeFullA<NEIGHF
for (kk = 0; kk < jnum; kk++) {
if (jj == kk) continue;
- k = d_neighbors(i,kk);
+ k = d_neighbors_short(i,kk);
k &= NEIGHMASK;
ktype = type(k);
@@ -559,9 +610,9 @@ void PairTersoffKokkos<DeviceType>::operator()(TagPairTersoffComputeFullA<NEIGHF
ters_dthb(itype,jtype,ktype,prefactor,rij,delx1,dely1,delz1,
rik,delx2,dely2,delz2,fi,fj,fk);
- f(i,0) += fi[0];
- f(i,1) += fi[1];
- f(i,2) += fi[2];
+ f_x += fi[0];
+ f_y += fi[1];
+ f_z += fi[2];
if (vflag_atom) {
F_FLOAT delrij[3], delrik[3];
@@ -571,6 +622,9 @@ void PairTersoffKokkos<DeviceType>::operator()(TagPairTersoffComputeFullA<NEIGHF
}
}
}
+ f(i,0) += f_x;
+ f(i,1) += f_y;
+ f(i,2) += f_z;
}
template<class DeviceType>
@@ -599,12 +653,16 @@ void PairTersoffKokkos<DeviceType>::operator()(TagPairTersoffComputeFullB<NEIGHF
F_FLOAT fj[3], fk[3];
X_FLOAT delx1, dely1, delz1, delx2, dely2, delz2;
- const int jnum = d_numneigh[i];
+ const int jnum = d_numneigh_short[i];
+
+ F_FLOAT f_x = 0.0;
+ F_FLOAT f_y = 0.0;
+ F_FLOAT f_z = 0.0;
// attractive: bond order
for (jj = 0; jj < jnum; jj++) {
- j = d_neighbors(i,jj);
+ j = d_neighbors_short(i,jj);
j &= NEIGHMASK;
if (j >= nlocal) continue;
jtype = type(j);
@@ -619,10 +677,10 @@ void PairTersoffKokkos<DeviceType>::operator()(TagPairTersoffComputeFullB<NEIGHF
if (rsq1 > cutsq1) continue;
rij = sqrt(rsq1);
- j_jnum = d_numneigh[j];
+ j_jnum = d_numneigh_short[j];
for (kk = 0; kk < j_jnum; kk++) {
- k = d_neighbors(j,kk);
+ k = d_neighbors_short(j,kk);
if (k == i) continue;
k &= NEIGHMASK;
ktype = type(k);
@@ -648,9 +706,9 @@ void PairTersoffKokkos<DeviceType>::operator()(TagPairTersoffComputeFullB<NEIGHF
const F_FLOAT prefactor = 0.5*fa * ters_dbij(jtype,itype,itype,bo_ij);
const F_FLOAT eng = 0.5*bij * fa;
- f(i,0) -= delx1*fatt;
- f(i,1) -= dely1*fatt;
- f(i,2) -= delz1*fatt;
+ f_x -= delx1*fatt;
+ f_y -= dely1*fatt;
+ f_z -= delz1*fatt;
if (EVFLAG) {
if (eflag)
@@ -662,7 +720,7 @@ void PairTersoffKokkos<DeviceType>::operator()(TagPairTersoffComputeFullB<NEIGHF
// attractive: three-body force
for (kk = 0; kk < j_jnum; kk++) {
- k = d_neighbors(j,kk);
+ k = d_neighbors_short(j,kk);
if (k == i) continue;
k &= NEIGHMASK;
ktype = type(k);
@@ -677,9 +735,9 @@ void PairTersoffKokkos<DeviceType>::operator()(TagPairTersoffComputeFullB<NEIGHF
rik = sqrt(rsq2);
ters_dthbj(jtype,itype,ktype,prefactor,rij,delx1,dely1,delz1,
rik,delx2,dely2,delz2,fj,fk);
- f(i,0) += fj[0];
- f(i,1) += fj[1];
- f(i,2) += fj[2];
+ f_x += fj[0];
+ f_y += fj[1];
+ f_z += fj[2];
if (vflag_atom) {
F_FLOAT delrji[3], delrjk[3];
@@ -692,11 +750,14 @@ void PairTersoffKokkos<DeviceType>::operator()(TagPairTersoffComputeFullB<NEIGHF
const F_FLOAT prefactor_jk = 0.5*fa_jk * ters_dbij(jtype,ktype,itype,bo_ij);
ters_dthbk(jtype,ktype,itype,prefactor_jk,rik,delx2,dely2,delz2,
rij,delx1,dely1,delz1,fk);
- f(i,0) += fk[0];
- f(i,1) += fk[1];
- f(i,2) += fk[2];
+ f_x += fk[0];
+ f_y += fk[1];
+ f_z += fk[2];
}
}
+ f(i,0) += f_x;
+ f(i,1) += f_y;
+ f(i,2) += f_z;
}
template<class DeviceType>
@@ -749,8 +810,9 @@ double PairTersoffKokkos<DeviceType>::bondorder(const int &i, const int &j, cons
const F_FLOAT costheta = (dx1*dx2 + dy1*dy2 + dz1*dz2)/(rij*rik);
- if (int(paramskk(i,j,k).powerm) == 3) arg = pow(paramskk(i,j,k).lam3 * (rij-rik),3.0);
- else arg = paramskk(i,j,k).lam3 * (rij-rik);
+ const F_FLOAT param = paramskk(i,j,k).lam3 * (rij-rik);
+ if (int(paramskk(i,j,k).powerm) == 3) arg = param*param*param;//pow(paramskk(i,j,k).lam3 * (rij-rik),3.0);
+ else arg = param;
if (arg > 69.0776) ex_delr = 1.e30;
else if (arg < -69.0776) ex_delr = 0.0;
@@ -780,7 +842,6 @@ KOKKOS_INLINE_FUNCTION
double PairTersoffKokkos<DeviceType>::
ters_dgijk(const int &i, const int &j, const int &k, const F_FLOAT &cos) const
{
-
const F_FLOAT ters_c = paramskk(i,j,k).c * paramskk(i,j,k).c;
const F_FLOAT ters_d = paramskk(i,j,k).d * paramskk(i,j,k).d;
const F_FLOAT hcth = paramskk(i,j,k).h - cos;
@@ -838,9 +899,9 @@ double PairTersoffKokkos<DeviceType>::ters_dbij(const int &i, const int &j,
const int &k, const F_FLOAT &bo) const
{
const F_FLOAT tmp = paramskk(i,j,k).beta * bo;
- if (tmp > paramskk(i,j,k).c1) return paramskk(i,j,k).beta * -0.5*pow(tmp,-1.5);
+ if (tmp > paramskk(i,j,k).c1) return paramskk(i,j,k).beta * -0.5/sqrt(tmp*tmp);//*pow(tmp,-1.5);
if (tmp > paramskk(i,j,k).c2)
- return paramskk(i,j,k).beta * (-0.5*pow(tmp,-1.5) *
+ return paramskk(i,j,k).beta * (-0.5/sqrt(tmp*tmp) * //*pow(tmp,-1.5) *
(1.0 - 0.5*(1.0 + 1.0/(2.0*paramskk(i,j,k).powern)) *
pow(tmp,-paramskk(i,j,k).powern)));
if (tmp < paramskk(i,j,k).c4) return 0.0;
@@ -883,15 +944,16 @@ void PairTersoffKokkos<DeviceType>::ters_dthb(
fc = ters_fc_k(i,j,k,rik);
dfc = ters_dfc(i,j,k,rik);
- if (int(paramskk(i,j,k).powerm) == 3) tmp = pow(paramskk(i,j,k).lam3 * (rij-rik),3.0);
- else tmp = paramskk(i,j,k).lam3 * (rij-rik);
+ const F_FLOAT param = paramskk(i,j,k).lam3 * (rij-rik);
+ if (int(paramskk(i,j,k).powerm) == 3) tmp = param*param*param;//pow(paramskk(i,j,k).lam3 * (rij-rik),3.0);
+ else tmp = param;
if (tmp > 69.0776) ex_delr = 1.e30;
else if (tmp < -69.0776) ex_delr = 0.0;
else ex_delr = exp(tmp);
if (int(paramskk(i,j,k).powerm) == 3)
- dex_delr = 3.0*pow(paramskk(i,j,k).lam3,3.0) * pow(rij-rik,2.0)*ex_delr;
+ dex_delr = 3.0*param*param*paramskk(i,j,k).lam3*ex_delr;//pow(rij-rik,2.0)*ex_delr;
else dex_delr = paramskk(i,j,k).lam3 * ex_delr;
cos = vec3_dot(rij_hat,rik_hat);
@@ -951,15 +1013,16 @@ void PairTersoffKokkos<DeviceType>::ters_dthbj(
fc = ters_fc_k(i,j,k,rik);
dfc = ters_dfc(i,j,k,rik);
- if (int(paramskk(i,j,k).powerm) == 3) tmp = pow(paramskk(i,j,k).lam3 * (rij-rik),3.0);
- else tmp = paramskk(i,j,k).lam3 * (rij-rik);
+ const F_FLOAT param = paramskk(i,j,k).lam3 * (rij-rik);
+ if (int(paramskk(i,j,k).powerm) == 3) tmp = param*param*param;//pow(paramskk(i,j,k).lam3 * (rij-rik),3.0);
+ else tmp = param;//paramskk(i,j,k).lam3 * (rij-rik);
if (tmp > 69.0776) ex_delr = 1.e30;
else if (tmp < -69.0776) ex_delr = 0.0;
else ex_delr = exp(tmp);
if (int(paramskk(i,j,k).powerm) == 3)
- dex_delr = 3.0*pow(paramskk(i,j,k).lam3,3.0) * pow(rij-rik,2.0)*ex_delr;
+ dex_delr = 3.0*param*param*paramskk(i,j,k).lam3*ex_delr;//pow(paramskk(i,j,k).lam3,3.0) * pow(rij-rik,2.0)*ex_delr;
else dex_delr = paramskk(i,j,k).lam3 * ex_delr;
cos = vec3_dot(rij_hat,rik_hat);
@@ -1012,15 +1075,16 @@ void PairTersoffKokkos<DeviceType>::ters_dthbk(
fc = ters_fc_k(i,j,k,rik);
dfc = ters_dfc(i,j,k,rik);
- if (int(paramskk(i,j,k).powerm) == 3) tmp = pow(paramskk(i,j,k).lam3 * (rij-rik),3.0);
- else tmp = paramskk(i,j,k).lam3 * (rij-rik);
+ const F_FLOAT param = paramskk(i,j,k).lam3 * (rij-rik);
+ if (int(paramskk(i,j,k).powerm) == 3) tmp = param*param*param;//pow(paramskk(i,j,k).lam3 * (rij-rik),3.0);
+ else tmp = param;//paramskk(i,j,k).lam3 * (rij-rik);
if (tmp > 69.0776) ex_delr = 1.e30;
else if (tmp < -69.0776) ex_delr = 0.0;
else ex_delr = exp(tmp);
if (int(paramskk(i,j,k).powerm) == 3)
- dex_delr = 3.0*pow(paramskk(i,j,k).lam3,3.0) * pow(rij-rik,2.0)*ex_delr;
+ dex_delr = 3.0*param*param*paramskk(i,j,k).lam3*ex_delr;//pow(paramskk(i,j,k).lam3,3.0) * pow(rij-rik,2.0)*ex_delr;
else dex_delr = paramskk(i,j,k).lam3 * ex_delr;
cos = vec3_dot(rij_hat,rik_hat);
diff --git a/src/KOKKOS/pair_tersoff_kokkos.h b/src/KOKKOS/pair_tersoff_kokkos.h
index c810a052dc..7490d3f45a 100644
--- a/src/KOKKOS/pair_tersoff_kokkos.h
+++ b/src/KOKKOS/pair_tersoff_kokkos.h
@@ -39,6 +39,8 @@ struct TagPairTersoffComputeFullA{};
template<int NEIGHFLAG, int EVFLAG>
struct TagPairTersoffComputeFullB{};
+struct TagPairTersoffComputeShortNeigh{};
+
template<class DeviceType>
class PairTersoffKokkos : public PairTersoff {
public:
@@ -77,6 +79,9 @@ class PairTersoffKokkos : public PairTersoff {
KOKKOS_INLINE_FUNCTION
void operator()(TagPairTersoffComputeFullB<NEIGHFLAG,EVFLAG>, const int&) const;
+ KOKKOS_INLINE_FUNCTION
+ void operator()(TagPairTersoffComputeShortNeigh, const int&) const;
+
KOKKOS_INLINE_FUNCTION
double ters_fc_k(const int &i, const int &j, const int &k, const F_FLOAT &r) const;
@@ -184,6 +189,7 @@ class PairTersoffKokkos : public PairTersoff {
// hardwired to space for 15 atom types
//params_ters m_params[MAX_TYPES_STACKPARAMS+1][MAX_TYPES_STACKPARAMS+1][MAX_TYPES_STACKPARAMS+1];
+ int inum;
typename AT::t_x_array_randomread x;
typename AT::t_f_array f;
typename AT::t_int_1d_randomread type;
@@ -203,10 +209,12 @@ class PairTersoffKokkos : public PairTersoff {
typename ArrayTypes<DeviceType>::t_int_1d_randomread d_numneigh;
//NeighListKokkos<DeviceType> k_list;
- class AtomKokkos *atomKK;
int neighflag,newton_pair;
int nlocal,nall,eflag,vflag;
+ Kokkos::View<int**,DeviceType> d_neighbors_short;
+ Kokkos::View<int*,DeviceType> d_numneigh_short;
+
friend void pair_virial_fdotr_compute<PairTersoffKokkos>(PairTersoffKokkos*);
};
diff --git a/src/KOKKOS/pair_tersoff_mod_kokkos.cpp b/src/KOKKOS/pair_tersoff_mod_kokkos.cpp
index 8cfaa63e7a..0c8f46a30d 100644
--- a/src/KOKKOS/pair_tersoff_mod_kokkos.cpp
+++ b/src/KOKKOS/pair_tersoff_mod_kokkos.cpp
@@ -191,7 +191,7 @@ void PairTersoffMODKokkos<DeviceType>::compute(int eflag_in, int vflag_in)
nall = atom->nlocal + atom->nghost;
newton_pair = force->newton_pair;
- const int inum = list->inum;
+ inum = list->inum;
const int ignum = inum + list->gnum;
NeighListKokkos<DeviceType>* k_list = static_cast<NeighListKokkos<DeviceType>*>(list);
d_numneigh = k_list->d_numneigh;
@@ -204,6 +204,18 @@ void PairTersoffMODKokkos<DeviceType>::compute(int eflag_in, int vflag_in)
EV_FLOAT ev;
EV_FLOAT ev_all;
+ // build short neighbor list
+
+ int max_neighs = d_neighbors.dimension_1();
+
+ if ((d_neighbors_short.dimension_1() != max_neighs) ||
+ (d_neighbors_short.dimension_0() != ignum)) {
+ d_neighbors_short = Kokkos::View<int**,DeviceType>("Tersoff::neighbors_short",ignum,max_neighs);
+ }
+ if (d_numneigh_short.dimension_0()!=ignum)
+ d_numneigh_short = Kokkos::View<int*,DeviceType>("Tersoff::numneighs_short",ignum);
+ Kokkos::parallel_for(Kokkos::RangePolicy<DeviceType,TagPairTersoffMODComputeShortNeigh>(0,neighflag==FULL?ignum:inum), *this);
+
if (neighflag == HALF) {
if (evflag)
Kokkos::parallel_reduce(Kokkos::RangePolicy<DeviceType, TagPairTersoffMODComputeHalf<HALF,1> >(0,inum),*this,ev);
@@ -257,6 +269,35 @@ void PairTersoffMODKokkos<DeviceType>::compute(int eflag_in, int vflag_in)
/* ---------------------------------------------------------------------- */
+template<class DeviceType>
+KOKKOS_INLINE_FUNCTION
+void PairTersoffMODKokkos<DeviceType>::operator()(TagPairTersoffMODComputeShortNeigh, const int& ii) const {
+ const int i = d_ilist[ii];
+ const X_FLOAT xtmp = x(i,0);
+ const X_FLOAT ytmp = x(i,1);
+ const X_FLOAT ztmp = x(i,2);
+
+ const int jnum = d_numneigh[i];
+ int inside = 0;
+ for (int jj = 0; jj < jnum; jj++) {
+ int j = d_neighbors(i,jj);
+ j &= NEIGHMASK;
+
+ const X_FLOAT delx = xtmp - x(j,0);
+ const X_FLOAT dely = ytmp - x(j,1);
+ const X_FLOAT delz = ztmp - x(j,2);
+ const F_FLOAT rsq = delx*delx + dely*dely + delz*delz;
+
+ if (rsq < cutmax*cutmax) {
+ d_neighbors_short(i,inside) = j;
+ inside++;
+ }
+ }
+ d_numneigh_short(i) = inside;
+}
+
+/* ---------------------------------------------------------------------- */
+
template<class DeviceType>
template<int NEIGHFLAG, int EVFLAG>
KOKKOS_INLINE_FUNCTION
@@ -273,21 +314,22 @@ void PairTersoffMODKokkos<DeviceType>::operator()(TagPairTersoffMODComputeHalf<N
const int itype = type(i);
const int itag = tag(i);
- int j,k,jj,kk,jtag,jtype,ktype;
- F_FLOAT rsq1, cutsq1, rsq2, cutsq2, rij, rik, bo_ij;
F_FLOAT fi[3], fj[3], fk[3];
- X_FLOAT delx1, dely1, delz1, delx2, dely2, delz2;
//const AtomNeighborsConst d_neighbors_i = k_list.get_neighbors_const(i);
- const int jnum = d_numneigh[i];
+ const int jnum = d_numneigh_short[i];
// repulsive
- for (jj = 0; jj < jnum; jj++) {
- j = d_neighbors(i,jj);
+ F_FLOAT f_x = 0.0;
+ F_FLOAT f_y = 0.0;
+ F_FLOAT f_z = 0.0;
+
+ for (int jj = 0; jj < jnum; jj++) {
+ int j = d_neighbors_short(i,jj);
j &= NEIGHMASK;
- jtype = type(j);
- jtag = tag(j);
+ const int jtype = type(j);
+ const int jtag = tag(j);
if (itag > jtag) {
if ((itag+jtag) % 2 == 0) continue;
@@ -315,9 +357,9 @@ void PairTersoffMODKokkos<DeviceType>::operator()(TagPairTersoffMODComputeHalf<N
(tmp_fcd - tmp_fce*paramskk(itype,jtype,jtype).lam1) / r;
const F_FLOAT eng = tmp_fce * paramskk(itype,jtype,jtype).biga * tmp_exp;
- a_f(i,0) += delx*frep;
- a_f(i,1) += dely*frep;
- a_f(i,2) += delz*frep;
+ f_x += delx*frep;
+ f_y += dely*frep;
+ f_z += delz*frep;
a_f(j,0) -= delx*frep;
a_f(j,1) -= dely*frep;
a_f(j,2) -= delz*frep;
@@ -330,35 +372,35 @@ void PairTersoffMODKokkos<DeviceType>::operator()(TagPairTersoffMODComputeHalf<N
// attractive: bond order
- for (jj = 0; jj < jnum; jj++) {
- j = d_neighbors(i,jj);
+ for (int jj = 0; jj < jnum; jj++) {
+ int j = d_neighbors_short(i,jj);
j &= NEIGHMASK;
- jtype = type(j);
+ const int jtype = type(j);
- delx1 = xtmp - x(j,0);
- dely1 = ytmp - x(j,1);
- delz1 = ztmp - x(j,2);
- rsq1 = delx1*delx1 + dely1*dely1 + delz1*delz1;
- cutsq1 = paramskk(itype,jtype,jtype).cutsq;
+ const F_FLOAT delx1 = xtmp - x(j,0);
+ const F_FLOAT dely1 = ytmp - x(j,1);
+ const F_FLOAT delz1 = ztmp - x(j,2);
+ const F_FLOAT rsq1 = delx1*delx1 + dely1*dely1 + delz1*delz1;
+ const F_FLOAT cutsq1 = paramskk(itype,jtype,jtype).cutsq;
- bo_ij = 0.0;
+ F_FLOAT bo_ij = 0.0;
if (rsq1 > cutsq1) continue;
- rij = sqrt(rsq1);
+ const F_FLOAT rij = sqrt(rsq1);
- for (kk = 0; kk < jnum; kk++) {
+ for (int kk = 0; kk < jnum; kk++) {
if (jj == kk) continue;
- k = d_neighbors(i,kk);
+ int k = d_neighbors_short(i,kk);
k &= NEIGHMASK;
- ktype = type(k);
+ const int ktype = type(k);
- delx2 = xtmp - x(k,0);
- dely2 = ytmp - x(k,1);
- delz2 = ztmp - x(k,2);
- rsq2 = delx2*delx2 + dely2*dely2 + delz2*delz2;
- cutsq2 = paramskk(itype,jtype,ktype).cutsq;
+ const F_FLOAT delx2 = xtmp - x(k,0);
+ const F_FLOAT dely2 = ytmp - x(k,1);
+ const F_FLOAT delz2 = ztmp - x(k,2);
+ const F_FLOAT rsq2 = delx2*delx2 + dely2*dely2 + delz2*delz2;
+ const F_FLOAT cutsq2 = paramskk(itype,jtype,ktype).cutsq;
if (rsq2 > cutsq2) continue;
- rik = sqrt(rsq2);
+ const F_FLOAT rik = sqrt(rsq2);
bo_ij += bondorder(itype,jtype,ktype,rij,delx1,dely1,delz1,rik,delx2,dely2,delz2);
}
@@ -369,58 +411,64 @@ void PairTersoffMODKokkos<DeviceType>::operator()(TagPairTersoffMODComputeHalf<N
const F_FLOAT bij = ters_bij_k(itype,jtype,jtype,bo_ij);
const F_FLOAT fatt = -0.5*bij * dfa / rij;
const F_FLOAT prefactor = 0.5*fa * ters_dbij(itype,jtype,jtype,bo_ij);
- const F_FLOAT eng = 0.5*bij * fa;
- a_f(i,0) += delx1*fatt;
- a_f(i,1) += dely1*fatt;
- a_f(i,2) += delz1*fatt;
- a_f(j,0) -= delx1*fatt;
- a_f(j,1) -= dely1*fatt;
- a_f(j,2) -= delz1*fatt;
+ f_x += delx1*fatt;
+ f_y += dely1*fatt;
+ f_z += delz1*fatt;
+ F_FLOAT fj_x = -delx1*fatt;
+ F_FLOAT fj_y = -dely1*fatt;
+ F_FLOAT fj_z = -delz1*fatt;
if (EVFLAG) {
+ const F_FLOAT eng = 0.5*bij * fa;
if (eflag) ev.evdwl += eng;
if (vflag_either || eflag_atom)
- this->template ev_tally<NEIGHFLAG>(ev,i,j,eng,fatt,delx1,dely1,delz1);
+ this->template ev_tally<NEIGHFLAG>(ev,i,j,eng,fatt,delx1,dely1,delz1);
}
// attractive: three-body force
- for (kk = 0; kk < jnum; kk++) {
+ for (int kk = 0; kk < jnum; kk++) {
if (jj == kk) continue;
- k = d_neighbors(i,kk);
+ int k = d_neighbors_short(i,kk);
k &= NEIGHMASK;
- ktype = type(k);
+ const int ktype = type(k);
- delx2 = xtmp - x(k,0);
- dely2 = ytmp - x(k,1);
- delz2 = ztmp - x(k,2);
- rsq2 = delx2*delx2 + dely2*dely2 + delz2*delz2;
- cutsq2 = paramskk(itype,jtype,ktype).cutsq;
+ const F_FLOAT delx2 = xtmp - x(k,0);
+ const F_FLOAT dely2 = ytmp - x(k,1);
+ const F_FLOAT delz2 = ztmp - x(k,2);
+ const F_FLOAT rsq2 = delx2*delx2 + dely2*dely2 + delz2*delz2;
+ const F_FLOAT cutsq2 = paramskk(itype,jtype,ktype).cutsq;
if (rsq2 > cutsq2) continue;
- rik = sqrt(rsq2);
+ const F_FLOAT rik = sqrt(rsq2);
ters_dthb(itype,jtype,ktype,prefactor,rij,delx1,dely1,delz1,
- rik,delx2,dely2,delz2,fi,fj,fk);
-
- a_f(i,0) += fi[0];
- a_f(i,1) += fi[1];
- a_f(i,2) += fi[2];
- a_f(j,0) += fj[0];
- a_f(j,1) += fj[1];
- a_f(j,2) += fj[2];
+ rik,delx2,dely2,delz2,fi,fj,fk);
+
+ f_x += fi[0];
+ f_y += fi[1];
+ f_z += fi[2];
+ fj_x += fj[0];
+ fj_y += fj[1];
+ fj_z += fj[2];
a_f(k,0) += fk[0];
a_f(k,1) += fk[1];
a_f(k,2) += fk[2];
if (vflag_atom) {
- F_FLOAT delrij[3], delrik[3];
- delrij[0] = -delx1; delrij[1] = -dely1; delrij[2] = -delz1;
- delrik[0] = -delx2; delrik[1] = -dely2; delrik[2] = -delz2;
- if (vflag_either) this->template v_tally3<NEIGHFLAG>(ev,i,j,k,fj,fk,delrij,delrik);
+ F_FLOAT delrij[3], delrik[3];
+ delrij[0] = -delx1; delrij[1] = -dely1; delrij[2] = -delz1;
+ delrik[0] = -delx2; delrik[1] = -dely2; delrik[2] = -delz2;
+ if (vflag_either) this->template v_tally3<NEIGHFLAG>(ev,i,j,k,fj,fk,delrij,delrik);
}
}
+ a_f(j,0) += fj_x;
+ a_f(j,1) += fj_y;
+ a_f(j,2) += fj_z;
}
+ a_f(i,0) += f_x;
+ a_f(i,1) += f_y;
+ a_f(i,2) += f_z;
}
template<class DeviceType>
@@ -450,12 +498,15 @@ void PairTersoffMODKokkos<DeviceType>::operator()(TagPairTersoffMODComputeFullA<
X_FLOAT delx1, dely1, delz1, delx2, dely2, delz2;
//const AtomNeighborsConst d_neighbors_i = k_list.get_neighbors_const(i);
- const int jnum = d_numneigh[i];
+ const int jnum = d_numneigh_short[i];
// repulsive
+ F_FLOAT f_x = 0.0;
+ F_FLOAT f_y = 0.0;
+ F_FLOAT f_z = 0.0;
for (jj = 0; jj < jnum; jj++) {
- j = d_neighbors(i,jj);
+ j = d_neighbors_short(i,jj);
j &= NEIGHMASK;
const int jtype = type(j);
@@ -475,9 +526,9 @@ void PairTersoffMODKokkos<DeviceType>::operator()(TagPairTersoffMODComputeFullA<
(tmp_fcd - tmp_fce*paramskk(itype,jtype,jtype).lam1) / r;
const F_FLOAT eng = tmp_fce * paramskk(itype,jtype,jtype).biga * tmp_exp;
- f(i,0) += delx*frep;
- f(i,1) += dely*frep;
- f(i,2) += delz*frep;
+ f_x += delx*frep;
+ f_y += dely*frep;
+ f_z += delz*frep;
if (EVFLAG) {
if (eflag)
@@ -490,7 +541,7 @@ void PairTersoffMODKokkos<DeviceType>::operator()(TagPairTersoffMODComputeFullA<
// attractive: bond order
for (jj = 0; jj < jnum; jj++) {
- j = d_neighbors(i,jj);
+ j = d_neighbors_short(i,jj);
j &= NEIGHMASK;
jtype = type(j);
@@ -506,7 +557,7 @@ void PairTersoffMODKokkos<DeviceType>::operator()(TagPairTersoffMODComputeFullA<
for (kk = 0; kk < jnum; kk++) {
if (jj == kk) continue;
- k = d_neighbors(i,kk);
+ k = d_neighbors_short(i,kk);
k &= NEIGHMASK;
ktype = type(k);
@@ -530,9 +581,9 @@ void PairTersoffMODKokkos<DeviceType>::operator()(TagPairTersoffMODComputeFullA<
const F_FLOAT prefactor = 0.5*fa * ters_dbij(itype,jtype,jtype,bo_ij);
const F_FLOAT eng = 0.5*bij * fa;
- f(i,0) += delx1*fatt;
- f(i,1) += dely1*fatt;
- f(i,2) += delz1*fatt;
+ f_x += delx1*fatt;
+ f_y += dely1*fatt;
+ f_z += delz1*fatt;
if (EVFLAG) {
if (eflag) ev.evdwl += 0.5*eng;
@@ -544,7 +595,7 @@ void PairTersoffMODKokkos<DeviceType>::operator()(TagPairTersoffMODComputeFullA<
for (kk = 0; kk < jnum; kk++) {
if (jj == kk) continue;
- k = d_neighbors(i,kk);
+ k = d_neighbors_short(i,kk);
k &= NEIGHMASK;
ktype = type(k);
@@ -559,9 +610,9 @@ void PairTersoffMODKokkos<DeviceType>::operator()(TagPairTersoffMODComputeFullA<
ters_dthb(itype,jtype,ktype,prefactor,rij,delx1,dely1,delz1,
rik,delx2,dely2,delz2,fi,fj,fk);
- f(i,0) += fi[0];
- f(i,1) += fi[1];
- f(i,2) += fi[2];
+ f_x += fi[0];
+ f_y += fi[1];
+ f_z += fi[2];
if (vflag_atom) {
F_FLOAT delrij[3], delrik[3];
@@ -571,6 +622,9 @@ void PairTersoffMODKokkos<DeviceType>::operator()(TagPairTersoffMODComputeFullA<
}
}
}
+ f(i,0) += f_x;
+ f(i,1) += f_y;
+ f(i,2) += f_z;
}
template<class DeviceType>
@@ -599,12 +653,16 @@ void PairTersoffMODKokkos<DeviceType>::operator()(TagPairTersoffMODComputeFullB<
F_FLOAT fj[3], fk[3];
X_FLOAT delx1, dely1, delz1, delx2, dely2, delz2;
- const int jnum = d_numneigh[i];
+ const int jnum = d_numneigh_short[i];
+
+ F_FLOAT f_x = 0.0;
+ F_FLOAT f_y = 0.0;
+ F_FLOAT f_z = 0.0;
// attractive: bond order
for (jj = 0; jj < jnum; jj++) {
- j = d_neighbors(i,jj);
+ j = d_neighbors_short(i,jj);
j &= NEIGHMASK;
if (j >= nlocal) continue;
jtype = type(j);
@@ -619,10 +677,10 @@ void PairTersoffMODKokkos<DeviceType>::operator()(TagPairTersoffMODComputeFullB<
if (rsq1 > cutsq1) continue;
rij = sqrt(rsq1);
- j_jnum = d_numneigh[j];
+ j_jnum = d_numneigh_short[j];
for (kk = 0; kk < j_jnum; kk++) {
- k = d_neighbors(j,kk);
+ k = d_neighbors_short(j,kk);
if (k == i) continue;
k &= NEIGHMASK;
ktype = type(k);
@@ -648,9 +706,9 @@ void PairTersoffMODKokkos<DeviceType>::operator()(TagPairTersoffMODComputeFullB<
const F_FLOAT prefactor = 0.5*fa * ters_dbij(jtype,itype,itype,bo_ij);
const F_FLOAT eng = 0.5*bij * fa;
- f(i,0) -= delx1*fatt;
- f(i,1) -= dely1*fatt;
- f(i,2) -= delz1*fatt;
+ f_x -= delx1*fatt;
+ f_y -= dely1*fatt;
+ f_z -= delz1*fatt;
if (EVFLAG) {
if (eflag)
@@ -662,7 +720,7 @@ void PairTersoffMODKokkos<DeviceType>::operator()(TagPairTersoffMODComputeFullB<
// attractive: three-body force
for (kk = 0; kk < j_jnum; kk++) {
- k = d_neighbors(j,kk);
+ k = d_neighbors_short(j,kk);
if (k == i) continue;
k &= NEIGHMASK;
ktype = type(k);
@@ -677,9 +735,9 @@ void PairTersoffMODKokkos<DeviceType>::operator()(TagPairTersoffMODComputeFullB<
rik = sqrt(rsq2);
ters_dthbj(jtype,itype,ktype,prefactor,rij,delx1,dely1,delz1,
rik,delx2,dely2,delz2,fj,fk);
- f(i,0) += fj[0];
- f(i,1) += fj[1];
- f(i,2) += fj[2];
+ f_x += fj[0];
+ f_y += fj[1];
+ f_z += fj[2];
if (vflag_atom) {
F_FLOAT delrji[3], delrjk[3];
@@ -692,11 +750,14 @@ void PairTersoffMODKokkos<DeviceType>::operator()(TagPairTersoffMODComputeFullB<
const F_FLOAT prefactor_jk = 0.5*fa_jk * ters_dbij(jtype,ktype,itype,bo_ij);
ters_dthbk(jtype,ktype,itype,prefactor_jk,rik,delx2,dely2,delz2,
rij,delx1,dely1,delz1,fk);
- f(i,0) += fk[0];
- f(i,1) += fk[1];
- f(i,2) += fk[2];
+ f_x += fk[0];
+ f_y += fk[1];
+ f_z += fk[2];
}
}
+ f(i,0) += f_x;
+ f(i,1) += f_y;
+ f(i,2) += f_z;
}
template<class DeviceType>
diff --git a/src/KOKKOS/pair_tersoff_mod_kokkos.h b/src/KOKKOS/pair_tersoff_mod_kokkos.h
index 07a76b4e66..5a26fa1557 100644
--- a/src/KOKKOS/pair_tersoff_mod_kokkos.h
+++ b/src/KOKKOS/pair_tersoff_mod_kokkos.h
@@ -39,6 +39,8 @@ struct TagPairTersoffMODComputeFullA{};
template<int NEIGHFLAG, int EVFLAG>
struct TagPairTersoffMODComputeFullB{};
+struct TagPairTersoffMODComputeShortNeigh{};
+
template<class DeviceType>
class PairTersoffMODKokkos : public PairTersoffMOD {
public:
@@ -77,6 +79,9 @@ class PairTersoffMODKokkos : public PairTersoffMOD {
KOKKOS_INLINE_FUNCTION
void operator()(TagPairTersoffMODComputeFullB<NEIGHFLAG,EVFLAG>, const int&) const;
+ KOKKOS_INLINE_FUNCTION
+ void operator()(TagPairTersoffMODComputeShortNeigh, const int&) const;
+
KOKKOS_INLINE_FUNCTION
double ters_fc_k(const int &i, const int &j, const int &k, const F_FLOAT &r) const;
@@ -184,6 +189,7 @@ class PairTersoffMODKokkos : public PairTersoffMOD {
// hardwired to space for 15 atom types
//params_ters m_params[MAX_TYPES_STACKPARAMS+1][MAX_TYPES_STACKPARAMS+1][MAX_TYPES_STACKPARAMS+1];
+ int inum;
typename AT::t_x_array_randomread x;
typename AT::t_f_array f;
typename AT::t_int_1d_randomread type;
@@ -203,10 +209,12 @@ class PairTersoffMODKokkos : public PairTersoffMOD {
typename ArrayTypes<DeviceType>::t_int_1d_randomread d_numneigh;
//NeighListKokkos<DeviceType> k_list;
- class AtomKokkos *atomKK;
int neighflag,newton_pair;
int nlocal,nall,eflag,vflag;
+ Kokkos::View<int**,DeviceType> d_neighbors_short;
+ Kokkos::View<int*,DeviceType> d_numneigh_short;
+
friend void pair_virial_fdotr_compute<PairTersoffMODKokkos>(PairTersoffMODKokkos*);
};
diff --git a/src/KOKKOS/pair_tersoff_zbl_kokkos.cpp b/src/KOKKOS/pair_tersoff_zbl_kokkos.cpp
index 382cff8ef5..6f011aecf4 100644
--- a/src/KOKKOS/pair_tersoff_zbl_kokkos.cpp
+++ b/src/KOKKOS/pair_tersoff_zbl_kokkos.cpp
@@ -205,7 +205,7 @@ void PairTersoffZBLKokkos<DeviceType>::compute(int eflag_in, int vflag_in)
nall = atom->nlocal + atom->nghost;
newton_pair = force->newton_pair;
- const int inum = list->inum;
+ inum = list->inum;
const int ignum = inum + list->gnum;
NeighListKokkos<DeviceType>* k_list = static_cast<NeighListKokkos<DeviceType>*>(list);
d_numneigh = k_list->d_numneigh;
@@ -218,6 +218,18 @@ void PairTersoffZBLKokkos<DeviceType>::compute(int eflag_in, int vflag_in)
EV_FLOAT ev;
EV_FLOAT ev_all;
+ // build short neighbor list
+
+ int max_neighs = d_neighbors.dimension_1();
+
+ if ((d_neighbors_short.dimension_1() != max_neighs) ||
+ (d_neighbors_short.dimension_0() != ignum)) {
+ d_neighbors_short = Kokkos::View<int**,DeviceType>("Tersoff::neighbors_short",ignum,max_neighs);
+ }
+ if (d_numneigh_short.dimension_0()!=ignum)
+ d_numneigh_short = Kokkos::View<int*,DeviceType>("Tersoff::numneighs_short",ignum);
+ Kokkos::parallel_for(Kokkos::RangePolicy<DeviceType,TagPairTersoffZBLComputeShortNeigh>(0,neighflag==FULL?ignum:inum), *this);
+
if (neighflag == HALF) {
if (evflag)
Kokkos::parallel_reduce(Kokkos::RangePolicy<DeviceType, TagPairTersoffZBLComputeHalf<HALF,1> >(0,inum),*this,ev);
@@ -271,6 +283,35 @@ void PairTersoffZBLKokkos<DeviceType>::compute(int eflag_in, int vflag_in)
/* ---------------------------------------------------------------------- */
+template<class DeviceType>
+KOKKOS_INLINE_FUNCTION
+void PairTersoffZBLKokkos<DeviceType>::operator()(TagPairTersoffZBLComputeShortNeigh, const int& ii) const {
+ const int i = d_ilist[ii];
+ const X_FLOAT xtmp = x(i,0);
+ const X_FLOAT ytmp = x(i,1);
+ const X_FLOAT ztmp = x(i,2);
+
+ const int jnum = d_numneigh[i];
+ int inside = 0;
+ for (int jj = 0; jj < jnum; jj++) {
+ int j = d_neighbors(i,jj);
+ j &= NEIGHMASK;
+
+ const X_FLOAT delx = xtmp - x(j,0);
+ const X_FLOAT dely = ytmp - x(j,1);
+ const X_FLOAT delz = ztmp - x(j,2);
+ const F_FLOAT rsq = delx*delx + dely*dely + delz*delz;
+
+ if (rsq < cutmax*cutmax) {
+ d_neighbors_short(i,inside) = j;
+ inside++;
+ }
+ }
+ d_numneigh_short(i) = inside;
+}
+
+/* ---------------------------------------------------------------------- */
+
template<class DeviceType>
template<int NEIGHFLAG, int EVFLAG>
KOKKOS_INLINE_FUNCTION
@@ -287,21 +328,22 @@ void PairTersoffZBLKokkos<DeviceType>::operator()(TagPairTersoffZBLComputeHalf<N
const int itype = type(i);
const int itag = tag(i);
- int j,k,jj,kk,jtag,jtype,ktype;
- F_FLOAT rsq1, cutsq1, rsq2, cutsq2, rij, rik, bo_ij;
F_FLOAT fi[3], fj[3], fk[3];
- X_FLOAT delx1, dely1, delz1, delx2, dely2, delz2;
//const AtomNeighborsConst d_neighbors_i = k_list.get_neighbors_const(i);
- const int jnum = d_numneigh[i];
+ const int jnum = d_numneigh_short[i];
// repulsive
- for (jj = 0; jj < jnum; jj++) {
- j = d_neighbors(i,jj);
+ F_FLOAT f_x = 0.0;
+ F_FLOAT f_y = 0.0;
+ F_FLOAT f_z = 0.0;
+
+ for (int jj = 0; jj < jnum; jj++) {
+ int j = d_neighbors_short(i,jj);
j &= NEIGHMASK;
- jtype = type(j);
- jtag = tag(j);
+ const int jtype = type(j);
+ const int jtag = tag(j);
if (itag > jtag) {
if ((itag+jtag) % 2 == 0) continue;
@@ -359,9 +401,9 @@ void PairTersoffZBLKokkos<DeviceType>::operator()(TagPairTersoffZBLComputeHalf<N
eng = (1.0 - fermi_k(itype,jtype,jtype,r)) * eng_z +
fermi_k(itype,jtype,jtype,r) * eng_t;
- a_f(i,0) += delx*frep;
- a_f(i,1) += dely*frep;
- a_f(i,2) += delz*frep;
+ f_x += delx*frep;
+ f_y += dely*frep;
+ f_z += delz*frep;
a_f(j,0) -= delx*frep;
a_f(j,1) -= dely*frep;
a_f(j,2) -= delz*frep;
@@ -374,35 +416,35 @@ void PairTersoffZBLKokkos<DeviceType>::operator()(TagPairTersoffZBLComputeHalf<N
// attractive: bond order
- for (jj = 0; jj < jnum; jj++) {
- j = d_neighbors(i,jj);
+ for (int jj = 0; jj < jnum; jj++) {
+ int j = d_neighbors_short(i,jj);
j &= NEIGHMASK;
- jtype = type(j);
+ const int jtype = type(j);
- delx1 = xtmp - x(j,0);
- dely1 = ytmp - x(j,1);
- delz1 = ztmp - x(j,2);
- rsq1 = delx1*delx1 + dely1*dely1 + delz1*delz1;
- cutsq1 = paramskk(itype,jtype,jtype).cutsq;
+ const F_FLOAT delx1 = xtmp - x(j,0);
+ const F_FLOAT dely1 = ytmp - x(j,1);
+ const F_FLOAT delz1 = ztmp - x(j,2);
+ const F_FLOAT rsq1 = delx1*delx1 + dely1*dely1 + delz1*delz1;
+ const F_FLOAT cutsq1 = paramskk(itype,jtype,jtype).cutsq;
- bo_ij = 0.0;
+ F_FLOAT bo_ij = 0.0;
if (rsq1 > cutsq1) continue;
- rij = sqrt(rsq1);
+ const F_FLOAT rij = sqrt(rsq1);
- for (kk = 0; kk < jnum; kk++) {
+ for (int kk = 0; kk < jnum; kk++) {
if (jj == kk) continue;
- k = d_neighbors(i,kk);
+ int k = d_neighbors_short(i,kk);
k &= NEIGHMASK;
- ktype = type(k);
+ const int ktype = type(k);
- delx2 = xtmp - x(k,0);
- dely2 = ytmp - x(k,1);
- delz2 = ztmp - x(k,2);
- rsq2 = delx2*delx2 + dely2*dely2 + delz2*delz2;
- cutsq2 = paramskk(itype,jtype,ktype).cutsq;
+ const F_FLOAT delx2 = xtmp - x(k,0);
+ const F_FLOAT dely2 = ytmp - x(k,1);
+ const F_FLOAT delz2 = ztmp - x(k,2);
+ const F_FLOAT rsq2 = delx2*delx2 + dely2*dely2 + delz2*delz2;
+ const F_FLOAT cutsq2 = paramskk(itype,jtype,ktype).cutsq;
if (rsq2 > cutsq2) continue;
- rik = sqrt(rsq2);
+ const F_FLOAT rik = sqrt(rsq2);
bo_ij += bondorder(itype,jtype,ktype,rij,delx1,dely1,delz1,rik,delx2,dely2,delz2);
}
@@ -413,58 +455,64 @@ void PairTersoffZBLKokkos<DeviceType>::operator()(TagPairTersoffZBLComputeHalf<N
const F_FLOAT bij = ters_bij_k(itype,jtype,jtype,bo_ij);
const F_FLOAT fatt = -0.5*bij * dfa / rij;
const F_FLOAT prefactor = 0.5*fa * ters_dbij(itype,jtype,jtype,bo_ij);
- const F_FLOAT eng = 0.5*bij * fa;
- a_f(i,0) += delx1*fatt;
- a_f(i,1) += dely1*fatt;
- a_f(i,2) += delz1*fatt;
- a_f(j,0) -= delx1*fatt;
- a_f(j,1) -= dely1*fatt;
- a_f(j,2) -= delz1*fatt;
+ f_x += delx1*fatt;
+ f_y += dely1*fatt;
+ f_z += delz1*fatt;
+ F_FLOAT fj_x = -delx1*fatt;
+ F_FLOAT fj_y = -dely1*fatt;
+ F_FLOAT fj_z = -delz1*fatt;
if (EVFLAG) {
+ const F_FLOAT eng = 0.5*bij * fa;
if (eflag) ev.evdwl += eng;
if (vflag_either || eflag_atom)
- this->template ev_tally<NEIGHFLAG>(ev,i,j,eng,fatt,delx1,dely1,delz1);
+ this->template ev_tally<NEIGHFLAG>(ev,i,j,eng,fatt,delx1,dely1,delz1);
}
// attractive: three-body force
- for (kk = 0; kk < jnum; kk++) {
+ for (int kk = 0; kk < jnum; kk++) {
if (jj == kk) continue;
- k = d_neighbors(i,kk);
+ int k = d_neighbors_short(i,kk);
k &= NEIGHMASK;
- ktype = type(k);
+ const int ktype = type(k);
- delx2 = xtmp - x(k,0);
- dely2 = ytmp - x(k,1);
- delz2 = ztmp - x(k,2);
- rsq2 = delx2*delx2 + dely2*dely2 + delz2*delz2;
- cutsq2 = paramskk(itype,jtype,ktype).cutsq;
+ const F_FLOAT delx2 = xtmp - x(k,0);
+ const F_FLOAT dely2 = ytmp - x(k,1);
+ const F_FLOAT delz2 = ztmp - x(k,2);
+ const F_FLOAT rsq2 = delx2*delx2 + dely2*dely2 + delz2*delz2;
+ const F_FLOAT cutsq2 = paramskk(itype,jtype,ktype).cutsq;
if (rsq2 > cutsq2) continue;
- rik = sqrt(rsq2);
+ const F_FLOAT rik = sqrt(rsq2);
ters_dthb(itype,jtype,ktype,prefactor,rij,delx1,dely1,delz1,
- rik,delx2,dely2,delz2,fi,fj,fk);
-
- a_f(i,0) += fi[0];
- a_f(i,1) += fi[1];
- a_f(i,2) += fi[2];
- a_f(j,0) += fj[0];
- a_f(j,1) += fj[1];
- a_f(j,2) += fj[2];
+ rik,delx2,dely2,delz2,fi,fj,fk);
+
+ f_x += fi[0];
+ f_y += fi[1];
+ f_z += fi[2];
+ fj_x += fj[0];
+ fj_y += fj[1];
+ fj_z += fj[2];
a_f(k,0) += fk[0];
a_f(k,1) += fk[1];
a_f(k,2) += fk[2];
if (vflag_atom) {
- F_FLOAT delrij[3], delrik[3];
- delrij[0] = -delx1; delrij[1] = -dely1; delrij[2] = -delz1;
- delrik[0] = -delx2; delrik[1] = -dely2; delrik[2] = -delz2;
- if (vflag_either) this->template v_tally3<NEIGHFLAG>(ev,i,j,k,fj,fk,delrij,delrik);
+ F_FLOAT delrij[3], delrik[3];
+ delrij[0] = -delx1; delrij[1] = -dely1; delrij[2] = -delz1;
+ delrik[0] = -delx2; delrik[1] = -dely2; delrik[2] = -delz2;
+ if (vflag_either) this->template v_tally3<NEIGHFLAG>(ev,i,j,k,fj,fk,delrij,delrik);
}
}
+ a_f(j,0) += fj_x;
+ a_f(j,1) += fj_y;
+ a_f(j,2) += fj_z;
}
+ a_f(i,0) += f_x;
+ a_f(i,1) += f_y;
+ a_f(i,2) += f_z;
}
template<class DeviceType>
@@ -498,8 +546,11 @@ void PairTersoffZBLKokkos<DeviceType>::operator()(TagPairTersoffZBLComputeFullA<
// repulsive
+ F_FLOAT f_x = 0.0;
+ F_FLOAT f_y = 0.0;
+ F_FLOAT f_z = 0.0;
for (jj = 0; jj < jnum; jj++) {
- j = d_neighbors(i,jj);
+ j = d_neighbors_short(i,jj);
j &= NEIGHMASK;
const int jtype = type(j);
@@ -549,9 +600,9 @@ void PairTersoffZBLKokkos<DeviceType>::operator()(TagPairTersoffZBLComputeFullA<
eng = (1.0 - fermi_k(itype,jtype,jtype,r)) * eng_z +
fermi_k(itype,jtype,jtype,r) * eng_t;
- f(i,0) += delx*frep;
- f(i,1) += dely*frep;
- f(i,2) += delz*frep;
+ f_x += delx*frep;
+ f_y += dely*frep;
+ f_z += delz*frep;
if (EVFLAG) {
if (eflag)
@@ -564,7 +615,7 @@ void PairTersoffZBLKokkos<DeviceType>::operator()(TagPairTersoffZBLComputeFullA<
// attractive: bond order
for (jj = 0; jj < jnum; jj++) {
- j = d_neighbors(i,jj);
+ j = d_neighbors_short(i,jj);
j &= NEIGHMASK;
jtype = type(j);
@@ -580,7 +631,7 @@ void PairTersoffZBLKokkos<DeviceType>::operator()(TagPairTersoffZBLComputeFullA<
for (kk = 0; kk < jnum; kk++) {
if (jj == kk) continue;
- k = d_neighbors(i,kk);
+ k = d_neighbors_short(i,kk);
k &= NEIGHMASK;
ktype = type(k);
@@ -604,9 +655,9 @@ void PairTersoffZBLKokkos<DeviceType>::operator()(TagPairTersoffZBLComputeFullA<
const F_FLOAT prefactor = 0.5*fa * ters_dbij(itype,jtype,jtype,bo_ij);
const F_FLOAT eng = 0.5*bij * fa;
- f(i,0) += delx1*fatt;
- f(i,1) += dely1*fatt;
- f(i,2) += delz1*fatt;
+ f_x += delx1*fatt;
+ f_y += dely1*fatt;
+ f_z += delz1*fatt;
if (EVFLAG) {
if (eflag) ev.evdwl += 0.5*eng;
@@ -618,7 +669,7 @@ void PairTersoffZBLKokkos<DeviceType>::operator()(TagPairTersoffZBLComputeFullA<
for (kk = 0; kk < jnum; kk++) {
if (jj == kk) continue;
- k = d_neighbors(i,kk);
+ k = d_neighbors_short(i,kk);
k &= NEIGHMASK;
ktype = type(k);
@@ -633,9 +684,9 @@ void PairTersoffZBLKokkos<DeviceType>::operator()(TagPairTersoffZBLComputeFullA<
ters_dthb(itype,jtype,ktype,prefactor,rij,delx1,dely1,delz1,
rik,delx2,dely2,delz2,fi,fj,fk);
- f(i,0) += fi[0];
- f(i,1) += fi[1];
- f(i,2) += fi[2];
+ f_x += fi[0];
+ f_y += fi[1];
+ f_z += fi[2];
if (vflag_atom) {
F_FLOAT delrij[3], delrik[3];
@@ -645,6 +696,9 @@ void PairTersoffZBLKokkos<DeviceType>::operator()(TagPairTersoffZBLComputeFullA<
}
}
}
+ f(i,0) += f_x;
+ f(i,1) += f_y;
+ f(i,2) += f_z;
}
template<class DeviceType>
@@ -673,12 +727,16 @@ void PairTersoffZBLKokkos<DeviceType>::operator()(TagPairTersoffZBLComputeFullB<
F_FLOAT fj[3], fk[3];
X_FLOAT delx1, dely1, delz1, delx2, dely2, delz2;
- const int jnum = d_numneigh[i];
+ const int jnum = d_numneigh_short[i];
+
+ F_FLOAT f_x = 0.0;
+ F_FLOAT f_y = 0.0;
+ F_FLOAT f_z = 0.0;
// attractive: bond order
for (jj = 0; jj < jnum; jj++) {
- j = d_neighbors(i,jj);
+ j = d_neighbors_short(i,jj);
j &= NEIGHMASK;
if (j >= nlocal) continue;
jtype = type(j);
@@ -693,10 +751,10 @@ void PairTersoffZBLKokkos<DeviceType>::operator()(TagPairTersoffZBLComputeFullB<
if (rsq1 > cutsq1) continue;
rij = sqrt(rsq1);
- j_jnum = d_numneigh[j];
+ j_jnum = d_numneigh_short[j];
for (kk = 0; kk < j_jnum; kk++) {
- k = d_neighbors(j,kk);
+ k = d_neighbors_short(j,kk);
if (k == i) continue;
k &= NEIGHMASK;
ktype = type(k);
@@ -722,9 +780,9 @@ void PairTersoffZBLKokkos<DeviceType>::operator()(TagPairTersoffZBLComputeFullB<
const F_FLOAT prefactor = 0.5*fa * ters_dbij(jtype,itype,itype,bo_ij);
const F_FLOAT eng = 0.5*bij * fa;
- f(i,0) -= delx1*fatt;
- f(i,1) -= dely1*fatt;
- f(i,2) -= delz1*fatt;
+ f_x -= delx1*fatt;
+ f_y -= dely1*fatt;
+ f_z -= delz1*fatt;
if (EVFLAG) {
if (eflag)
@@ -736,7 +794,7 @@ void PairTersoffZBLKokkos<DeviceType>::operator()(TagPairTersoffZBLComputeFullB<
// attractive: three-body force
for (kk = 0; kk < j_jnum; kk++) {
- k = d_neighbors(j,kk);
+ k = d_neighbors_short(j,kk);
if (k == i) continue;
k &= NEIGHMASK;
ktype = type(k);
@@ -751,9 +809,9 @@ void PairTersoffZBLKokkos<DeviceType>::operator()(TagPairTersoffZBLComputeFullB<
rik = sqrt(rsq2);
ters_dthbj(jtype,itype,ktype,prefactor,rij,delx1,dely1,delz1,
rik,delx2,dely2,delz2,fj,fk);
- f(i,0) += fj[0];
- f(i,1) += fj[1];
- f(i,2) += fj[2];
+ f_x += fj[0];
+ f_y += fj[1];
+ f_z += fj[2];
if (vflag_atom) {
F_FLOAT delrji[3], delrjk[3];
@@ -766,11 +824,14 @@ void PairTersoffZBLKokkos<DeviceType>::operator()(TagPairTersoffZBLComputeFullB<
const F_FLOAT prefactor_jk = 0.5*fa_jk * ters_dbij(jtype,ktype,itype,bo_ij);
ters_dthbk(jtype,ktype,itype,prefactor_jk,rik,delx2,dely2,delz2,
rij,delx1,dely1,delz1,fk);
- f(i,0) += fk[0];
- f(i,1) += fk[1];
- f(i,2) += fk[2];
+ f_x += fk[0];
+ f_y += fk[1];
+ f_z += fk[2];
}
}
+ f(i,0) += f_x;
+ f(i,1) += f_y;
+ f(i,2) += f_z;
}
template<class DeviceType>
diff --git a/src/KOKKOS/pair_tersoff_zbl_kokkos.h b/src/KOKKOS/pair_tersoff_zbl_kokkos.h
index d580634197..136366d3fa 100644
--- a/src/KOKKOS/pair_tersoff_zbl_kokkos.h
+++ b/src/KOKKOS/pair_tersoff_zbl_kokkos.h
@@ -39,6 +39,8 @@ struct TagPairTersoffZBLComputeFullA{};
template<int NEIGHFLAG, int EVFLAG>
struct TagPairTersoffZBLComputeFullB{};
+struct TagPairTersoffZBLComputeShortNeigh{};
+
template<class DeviceType>
class PairTersoffZBLKokkos : public PairTersoffZBL {
public:
@@ -77,6 +79,8 @@ class PairTersoffZBLKokkos : public PairTersoffZBL {
KOKKOS_INLINE_FUNCTION
void operator()(TagPairTersoffZBLComputeFullB<NEIGHFLAG,EVFLAG>, const int&) const;
+ KOKKOS_INLINE_FUNCTION
+ void operator()(TagPairTersoffZBLComputeShortNeigh, const int&) const;
KOKKOS_INLINE_FUNCTION
double ters_fc_k(const int &i, const int &j, const int &k, const F_FLOAT &r) const;
@@ -190,6 +194,7 @@ class PairTersoffZBLKokkos : public PairTersoffZBL {
// hardwired to space for 15 atom types
//params_ters m_params[MAX_TYPES_STACKPARAMS+1][MAX_TYPES_STACKPARAMS+1][MAX_TYPES_STACKPARAMS+1];
+ int inum;
typename AT::t_x_array_randomread x;
typename AT::t_f_array f;
typename AT::t_int_1d_randomread type;
@@ -209,10 +214,12 @@ class PairTersoffZBLKokkos : public PairTersoffZBL {
typename ArrayTypes<DeviceType>::t_int_1d_randomread d_numneigh;
//NeighListKokkos<DeviceType> k_list;
- class AtomKokkos *atomKK;
int neighflag,newton_pair;
int nlocal,nall,eflag,vflag;
+ Kokkos::View<int**,DeviceType> d_neighbors_short;
+ Kokkos::View<int*,DeviceType> d_numneigh_short;
+
// ZBL
F_FLOAT global_a_0; // Bohr radius for Coulomb repulsion
F_FLOAT global_epsilon_0; // permittivity of vacuum for Coulomb repulsion
diff --git a/src/pair_morse.cpp b/src/pair_morse.cpp
index 4e2d47be5c..9170eb7a1d 100644
--- a/src/pair_morse.cpp
+++ b/src/pair_morse.cpp
@@ -176,7 +176,8 @@ void PairMorse::settings(int narg, char **arg)
void PairMorse::coeff(int narg, char **arg)
{
- if (narg < 5 || narg > 6) error->all(FLERR,"Incorrect args for pair coefficients");
+ if (narg < 5 || narg > 6)
+ error->all(FLERR,"Incorrect args for pair coefficients");
if (!allocated) allocate();
int ilo,ihi,jlo,jhi;
--
GitLab