Merge branch 'master' into doc-updates

This commit is contained in:
Axel Kohlmeyer 2016-11-17 20:44:07 -05:00
commit 81f68e06fd
67 changed files with 5686 additions and 118 deletions

BIN
doc/src/Eqs/pair_agni.jpg Normal file

Binary file not shown.

After

Width:  |  Height:  |  Size: 15 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 4.1 KiB

View File

@ -0,0 +1,10 @@
\documentclass[12pt]{article}
\pagestyle{empty}
\begin{document}
\begin{eqnarray*}
V_{ij} & = & f_C(r_{ij}) \left[ f_R(r_{ij}) + b_{ij} f_A(r_{ij}) + c_0 \right]
\end{eqnarray*}
\end{document}

Binary file not shown.

Before

Width:  |  Height:  |  Size: 4.0 KiB

After

Width:  |  Height:  |  Size: 4.2 KiB

View File

@ -3,7 +3,7 @@
\begin{document}
$$
P = \frac{N k_B T}{V} + \frac{\sum_{i}^{N} r_i \bullet f_i}{dV}
P = \frac{N k_B T}{V} + \frac{\sum_{i}^{N'} r_i \bullet f_i}{dV}
$$
\end{document}

Binary file not shown.

Before

Width:  |  Height:  |  Size: 4.9 KiB

After

Width:  |  Height:  |  Size: 5.3 KiB

View File

@ -4,7 +4,7 @@
$$
P_{IJ} = \frac{\sum_{k}^{N} m_k v_{k_I} v_{k_J}}{V} +
\frac{\sum_{k}^{N} r_{k_I} f_{k_J}}{V}
\frac{\sum_{k}^{N'} r_{k_I} f_{k_J}}{V}
$$
\end{document}

View File

@ -1,7 +1,7 @@
<!-- HTML_ONLY -->
<HEAD>
<TITLE>LAMMPS Users Manual</TITLE>
<META NAME="docnumber" CONTENT="9 Nov 2016 version">
<META NAME="docnumber" CONTENT="17 Nov 2016 version">
<META NAME="author" CONTENT="http://lammps.sandia.gov - Sandia National Laboratories">
<META NAME="copyright" CONTENT="Copyright (2003) Sandia Corporation. This software and manual is distributed under the GNU General Public License.">
</HEAD>
@ -21,7 +21,7 @@
<H1></H1>
LAMMPS Documentation :c,h3
9 Nov 2016 version :c,h4
17 Nov 2016 version :c,h4
Version info: :h4

View File

@ -981,11 +981,12 @@ KOKKOS, o = USER-OMP, t = OPT.
"table (gko)"_pair_table.html,
"tersoff (gkio)"_pair_tersoff.html,
"tersoff/mod (gko)"_pair_tersoff_mod.html,
"tersoff/mod/c (o)"_pair_tersoff_mod.html,
"tersoff/zbl (gko)"_pair_tersoff_zbl.html,
"tip4p/cut (o)"_pair_coul.html,
"tip4p/long (o)"_pair_coul.html,
"tri/lj"_pair_tri_lj.html,
"vashishta (o)"_pair_vashishta.html,
"vashishta (ko)"_pair_vashishta.html,
"vashishta/table (o)"_pair_vashishta.html,
"yukawa (go)"_pair_yukawa.html,
"yukawa/colloid (go)"_pair_yukawa_colloid.html,
@ -995,6 +996,7 @@ These are additional pair styles in USER packages, which can be used
if "LAMMPS is built with the appropriate
package"_Section_start.html#start_3.
"agni (o)"_pair_agni.html,
"awpmd/cut"_pair_awpmd.html,
"buck/mdf"_pair_mdf.html,
"coul/cut/soft (o)"_pair_lj_soft.html,

View File

@ -37,12 +37,18 @@ The pressure is computed by the formula
where N is the number of atoms in the system (see discussion of DOF
below), Kb is the Boltzmann constant, T is the temperature, d is the
dimensionality of the system (2 or 3 for 2d/3d), V is the system
volume (or area in 2d), and the second term is the virial, computed
within LAMMPS for all pairwise as well as 2-body, 3-body, and 4-body,
and long-range interactions. "Fixes"_fix.html that impose constraints
(e.g. the "fix shake"_fix_shake.html command) also contribute to the
virial term.
dimensionality of the system (2 or 3 for 2d/3d), and V is the system
volume (or area in 2d). The second term is the virial, equal to
-dU/dV, computed for all pairwise as well as 2-body, 3-body, 4-body,
manybody, and long-range interactions, where r_i and f_i are the
position and force vector of atom i, and the black dot indicates a dot
product. When periodic boundary conditions are used, N' necessarily
includes periodic image (ghost) atoms outside the central box, and the
position and force vectors of ghost atoms are thus included in the
summation. When periodic boundary conditions are not used, N' = N =
the number of atoms in the system. "Fixes"_fix.html that impose
constraints (e.g. the "fix shake"_fix_shake.html command) also
contribute to the virial term.
A symmetric pressure tensor, stored as a 6-element vector, is also
calculated by this compute. The 6 components of the vector are
@ -62,8 +68,9 @@ compute temperature or ke and/or the virial. The {virial} keyword
means include all terms except the kinetic energy {ke}.
Details of how LAMMPS computes the virial efficiently for the entire
system, including the effects of periodic boundary conditions is
discussed in "(Thompson)"_#Thompson.
system, including for manybody potentials and accounting for the
effects of periodic boundary conditions are discussed in
"(Thompson)"_#Thompson.
The temperature and kinetic energy tensor is not calculated by this
compute, but rather by the temperature compute specified with the

View File

@ -59,6 +59,7 @@ dump_h5md.html
dump_image.html
dump_modify.html
dump_molfile.html
dump_nc.html
echo.html
fix.html
fix_modify.html
@ -152,6 +153,7 @@ fix_colvars.html
fix_controller.html
fix_deform.html
fix_deposit.html
fix_dpd_energy.html
fix_drag.html
fix_drude.html
fix_drude_transform.html
@ -272,6 +274,7 @@ fix_viscosity.html
fix_viscous.html
fix_wall.html
fix_wall_gran.html
fix_wall_gran_region.html
fix_wall_piston.html
fix_wall_reflect.html
fix_wall_region.html
@ -390,6 +393,7 @@ compute_voronoi_atom.html
compute_xrd.html
pair_adp.html
pair_agni.html
pair_airebo.html
pair_awpmd.html
pair_beck.html
@ -622,3 +626,4 @@ USER/atc/man_unfix_flux.html
USER/atc/man_unfix_nodes.html
USER/atc/man_write_atom_weights.html
USER/atc/man_write_restart.html

128
doc/src/pair_agni.txt Normal file
View File

@ -0,0 +1,128 @@
"LAMMPS WWW Site"_lws - "LAMMPS Documentation"_ld - "LAMMPS Commands"_lc :c
:link(lws,http://lammps.sandia.gov)
:link(ld,Manual.html)
:link(lc,Section_commands.html#comm)
:line
pair_style agni command :h3
pair_style agni/omp command :h3
[Syntax:]
pair_style agni :pre
[Examples:]
pair_style agni
pair_coeff * * Al.agni Al
[Description:]
Style {agni} style computes the manybody vectorial force components for
an atom as
:c,image(Eqs/pair_agni.jpg)
{u} labels the individual components, i.e. x, y or z, and {V} is the
corresponding atomic fingerprint. {d} is the Euclidean distance between
any two atomic fingerprints. A total of N_t reference atomic
environments are considered to construct the force field file. {alpha_t}
and {l} are the weight coefficients and length scale parameter of the
non-linear regression model.
The method implements the recently proposed machine learning access to
atomic forces as discussed extensively in the following publications -
"(Botu1)"_#Botu2015adaptive and "(Botu2)"_#Botu2015learning. The premise
of the method is to map the atomic enviornment numerically into a
fingerprint, and use machine learning methods to create a mapping to the
vectorial atomic forces.
Only a single pair_coeff command is used with the {agni} style which
specifies an AGNI potential file containing the parameters of the
force field for the needed elements. These are mapped to LAMMPS atom
types by specifying N additional arguments after the filename in the
pair_coeff command, where N is the number of LAMMPS atom types:
filename
N element names = mapping of AGNI elements to atom types :ul
See the "pair_coeff"_pair_coeff.html doc page for alternate ways
to specify the path for the force field file.
An AGNI force field is fully specified by the filename which contains the
parameters of the force field, i.e., the reference training environments
used to construct the machine learning force field. Example force field
and input files are provided in the examples/USER/misc/agni directory.
:line
Styles with {omp} suffix is 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 5"_Section_accelerate.html
of the manual. The accelerated style takes the same arguments and
should produce the same results, except for round-off and precision
issues.
The accelerated style is part of the USER-OMP. 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 style 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 5"_Section_accelerate.html of the manual for
more instructions on how to use the accelerated styles effectively.
:line
[Mixing, shift, table, tail correction, restart, rRESPA info]:
This pair style does not support the "pair_modify"_pair_modify.html
shift, table, and tail options.
This pair style does not write its information to "binary restart
files"_restart.html, since it is stored in potential files. Thus, you
need to re-specify the pair_style and pair_coeff commands in an input
script that reads a restart file.
This pair style can only be used via the {pair} keyword of the
"run_style respa"_run_style.html command. It does not support the
{inner}, {middle}, {outer} keywords.
:line
[Restrictions:]
Currently, only elemental systems are implemented. Also, the method only
provides access to the forces and not energies or stresses. However, one
can access the energy via thermodynamic integration of the forces as
discussed in "(Botu3)"_#Botu2016construct. This pair style is part
of the USER-MISC package. It is only enabled if LAMMPS was built with
that package. See the "Making LAMMPS"_Section_start.html#start_3 section
for more info.
The AGNI force field files provided with LAMMPS (see the
potentials directory) are parameterized for metal "units"_units.html.
You can use the AGNI potential with any LAMMPS units, but you would need
to create your own AGNI potential file with coefficients listed in the
appropriate units if your simulation doesn't use "metal" units.
[Related commands:]
"pair_coeff"_pair_coeff.html
[Default:] none
:line
:link(Botu2015adaptive)
[(Botu1)] V. Botu and R. Ramprasad, Int. J. Quant. Chem., 115(16), 1074 (2015).
:link(Botu2015learning)
[(Botu2)] V. Botu and R. Ramprasad, Phys. Rev. B, 92(9), 094306 (2015).
:link(Botu2016construct)
[(Botu3)] V. Botu, R. Batra, J. Chapman and R. Ramprasad, https://arxiv.org/abs/1610.02098 (2016).

View File

@ -7,32 +7,43 @@
:line
pair_style tersoff/mod command :h3
pair_style tersoff/mod/c command :h3
pair_style tersoff/mod/gpu command :h3
pair_style tersoff/mod/kk command :h3
pair_style tersoff/mod/omp command :h3
pair_style tersoff/mod/c/omp command :h3
[Syntax:]
pair_style tersoff/mod :pre
pair_style tersoff/mod/c :pre
[Examples:]
pair_style tersoff/mod
pair_coeff * * Si.tersoff.mod Si Si :pre
pair_style tersoff/mod/c
pair_coeff * * Si.tersoff.modc Si Si :pre
[Description:]
The {tersoff/mod} style computes a bond-order type interatomic
potential "(Kumagai)"_#Kumagai based on a 3-body Tersoff potential
"(Tersoff_1)"_#Tersoff_1, "(Tersoff_2)"_#Tersoff_2 with modified
cutoff function and angular-dependent term, giving the energy E of a
system of atoms as
The {tersoff/mod} and {tersoff/mod/c} styles computes a bond-order type
interatomic potential "(Kumagai)"_#Kumagai based on a 3-body Tersoff
potential "(Tersoff_1)"_#Tersoff_1, "(Tersoff_2)"_#Tersoff_2 with
modified cutoff function and angular-dependent term, giving the energy
E of a system of atoms as
:c,image(Eqs/pair_tersoff_mod.jpg)
where f_R is a two-body term and f_A includes three-body interactions.
The summations in the formula are over all neighbors J and K of atom I
within a cutoff distance = R + D.
The {tersoff/mod/c} style differs from {tersoff/mod} only in the
formulation of the V_ij term, where it contains an additional c0 term.
:c,image(Eqs/pair_tersoff_mod_c.jpg)
The modified cutoff function f_C proposed by "(Murty)"_#Murty and
having a continuous second-order differential is employed. The
@ -69,10 +80,11 @@ are placeholders for atom types that will be used with other
potentials.
Tersoff/MOD file in the {potentials} directory of the LAMMPS
distribution have a ".tersoff.mod" suffix. Lines that are not blank
or comments (starting with #) define parameters for a triplet of
elements. The parameters in a single entry correspond to coefficients
in the formula above:
distribution have a ".tersoff.mod" suffix. Potential files for the
{tersoff/mod/c} style have the suffix ".tersoff.modc". Lines that are
not blank or comments (starting with #) define parameters for a triplet
of elements. The parameters in a single entry correspond to
coefficients in the formulae above:
element 1 (the center atom in a 3-body interaction)
element 2 (the atom bonded to the center atom)
@ -93,13 +105,15 @@ c1
c2
c3
c4
c5 :ul
c5
c0 (energy units, tersoff/mod/c only):ul
The n, eta, lambda2, B, lambda1, and A parameters are only used for
two-body interactions. The beta, alpha, c1, c2, c3, c4, c5, h
parameters are only used for three-body interactions. The R and D
parameters are used for both two-body and three-body interactions. The
non-annotated parameters are unitless.
parameters are used for both two-body and three-body interactions.
The c0 term applies to {tersoff/mod/c} only. The non-annotated
parameters are unitless.
The Tersoff/MOD potential file must contain entries for all the elements
listed in the pair_coeff command. It can also contain entries for

View File

@ -8,6 +8,7 @@
pair_style vashishta command :h3
pair_style vashishta/omp command :h3
pair_style vashishta/kk command :h3
pair_style vashishta/table command :h3
pair_style vashishta/table/omp command :h3

View File

@ -6,6 +6,7 @@ Pair Styles :h1
:maxdepth: 1
pair_adp
pair_agni
pair_airebo
pair_awpmd
pair_beck

File diff suppressed because it is too large Load Diff

View File

@ -0,0 +1,197 @@
Data File created from VASP POSCAR
181 atoms
1 atom types
0 17.121440767 xlo xhi
0 14.8276026536 ylo yhi
0 39.3197318979 zlo zhi
Masses
1 26.982
Atoms
1 1 0.0 1.64751140595 15.0
2 1 2.85357346116 1.64751140595 15.0
3 1 5.70714692232 1.64751140595 15.0
4 1 8.56072038348 1.64751140595 15.0
5 1 11.4142938446 1.64751140595 15.0
6 1 14.2678673058 1.64751140595 15.0
7 1 1.42678673058 4.11877851488 15.0
8 1 4.28036019174 4.11877851488 15.0
9 1 7.1339336529 4.11877851488 15.0
10 1 9.98750711406 4.11877851488 15.0
11 1 12.8410805752 4.11877851488 15.0
12 1 15.6946540364 4.11877851488 15.0
13 1 0.0 6.59004562381 15.0
14 1 2.85357346116 6.59004562381 15.0
15 1 5.70714692232 6.59004562381 15.0
16 1 8.56072038348 6.59004562381 15.0
17 1 11.4142938446 6.59004562381 15.0
18 1 14.2678673058 6.59004562381 15.0
19 1 1.42678673058 9.06131273274 15.0
20 1 4.28036019174 9.06131273274 15.0
21 1 7.1339336529 9.06131273274 15.0
22 1 9.98750711406 9.06131273274 15.0
23 1 12.8410805752 9.06131273274 15.0
24 1 15.6946540364 9.06131273274 15.0
25 1 0.0 11.5325798417 15.0
26 1 2.85357346116 11.5325798417 15.0
27 1 5.70714692232 11.5325798417 15.0
28 1 8.56072038348 11.5325798417 15.0
29 1 11.4142938446 11.5325798417 15.0
30 1 14.2678673058 11.5325798417 15.0
31 1 1.42678673058 14.0038469506 15.0
32 1 4.28036019174 14.0038469506 15.0
33 1 7.1339336529 14.0038469506 15.0
34 1 9.98750711406 14.0038469506 15.0
35 1 12.8410805752 14.0038469506 15.0
36 1 15.6946540364 14.0038469506 15.0
37 1 0.0 0.0 17.3299329745
38 1 2.85357346116 0.0 17.3299329745
39 1 5.70714692232 0.0 17.3299329745
40 1 8.56072038348 0.0 17.3299329745
41 1 11.4142938446 0.0 17.3299329745
42 1 14.2678673058 0.0 17.3299329745
43 1 1.42678673058 2.47126710893 17.3299329745
44 1 4.28036019174 2.47126710893 17.3299329745
45 1 7.1339336529 2.47126710893 17.3299329745
46 1 9.98750711406 2.47126710893 17.3299329745
47 1 12.8410805752 2.47126710893 17.3299329745
48 1 15.6946540364 2.47126710893 17.3299329745
49 1 0.0 4.94253421786 17.3299329745
50 1 2.85357346116 4.94253421786 17.3299329745
51 1 5.70714692232 4.94253421786 17.3299329745
52 1 8.56072038348 4.94253421786 17.3299329745
53 1 11.4142938446 4.94253421786 17.3299329745
54 1 14.2678673058 4.94253421786 17.3299329745
55 1 1.42678673058 7.41380132679 17.3299329745
56 1 4.28036019174 7.41380132679 17.3299329745
57 1 7.1339336529 7.41380132679 17.3299329745
58 1 9.98750711406 7.41380132679 17.3299329745
59 1 12.8410805752 7.41380132679 17.3299329745
60 1 15.6946540364 7.41380132679 17.3299329745
61 1 0.0 9.88506843572 17.3299329745
62 1 2.85357346116 9.88506843572 17.3299329745
63 1 5.70714692232 9.88506843572 17.3299329745
64 1 8.56072038348 9.88506843572 17.3299329745
65 1 11.4142938446 9.88506843572 17.3299329745
66 1 14.2678673058 9.88506843572 17.3299329745
67 1 1.42678673058 12.3563355446 17.3299329745
68 1 4.28036019174 12.3563355446 17.3299329745
69 1 7.1339336529 12.3563355446 17.3299329745
70 1 9.98750711406 12.3563355446 17.3299329745
71 1 12.8410805752 12.3563355446 17.3299329745
72 1 15.6946540364 12.3563355446 17.3299329745
73 1 1.42678673058 0.823755702976 19.6598659489
74 1 4.28036019174 0.823755702976 19.6598659489
75 1 7.1339336529 0.823755702976 19.6598659489
76 1 9.98750711406 0.823755702976 19.6598659489
77 1 12.8410805752 0.823755702976 19.6598659489
78 1 15.6946540364 0.823755702976 19.6598659489
79 1 0.0 3.29502281191 19.6598659489
80 1 2.85357346116 3.29502281191 19.6598659489
81 1 5.70714692232 3.29502281191 19.6598659489
82 1 8.56072038348 3.29502281191 19.6598659489
83 1 11.4142938446 3.29502281191 19.6598659489
84 1 14.2678673058 3.29502281191 19.6598659489
85 1 1.42678673058 5.76628992084 19.6598659489
86 1 4.28036019174 5.76628992084 19.6598659489
87 1 7.1339336529 5.76628992084 19.6598659489
88 1 9.98750711406 5.76628992084 19.6598659489
89 1 12.8410805752 5.76628992084 19.6598659489
90 1 15.6946540364 5.76628992084 19.6598659489
91 1 0.0 8.23755702976 19.6598659489
92 1 2.85357346116 8.23755702976 19.6598659489
93 1 5.70714692232 8.23755702976 19.6598659489
94 1 8.56072038348 8.23755702976 19.6598659489
95 1 11.4142938446 8.23755702976 19.6598659489
96 1 14.2678673058 8.23755702976 19.6598659489
97 1 1.42678673058 10.7088241387 19.6598659489
98 1 4.28036019174 10.7088241387 19.6598659489
99 1 7.1339336529 10.7088241387 19.6598659489
100 1 9.98750711406 10.7088241387 19.6598659489
101 1 12.8410805752 10.7088241387 19.6598659489
102 1 15.6946540364 10.7088241387 19.6598659489
103 1 0.0 13.1800912476 19.6598659489
104 1 2.85357346116 13.1800912476 19.6598659489
105 1 5.70714692232 13.1800912476 19.6598659489
106 1 8.56072038348 13.1800912476 19.6598659489
107 1 11.4142938446 13.1800912476 19.6598659489
108 1 14.2678673058 13.1800912476 19.6598659489
109 1 0.0 1.64751140595 21.9897989234
110 1 2.85357346116 1.64751140595 21.9897989234
111 1 5.70714692232 1.64751140595 21.9897989234
112 1 8.56072038348 1.64751140595 21.9897989234
113 1 11.4142938446 1.64751140595 21.9897989234
114 1 14.2678673058 1.64751140595 21.9897989234
115 1 1.42678673058 4.11877851488 21.9897989234
116 1 4.28036019174 4.11877851488 21.9897989234
117 1 7.1339336529 4.11877851488 21.9897989234
118 1 9.98750711406 4.11877851488 21.9897989234
119 1 12.8410805752 4.11877851488 21.9897989234
120 1 15.6946540364 4.11877851488 21.9897989234
121 1 0.0 6.59004562381 21.9897989234
122 1 2.85357346116 6.59004562381 21.9897989234
123 1 5.70714692232 6.59004562381 21.9897989234
124 1 8.56072038348 6.59004562381 21.9897989234
125 1 11.4142938446 6.59004562381 21.9897989234
126 1 14.2678673058 6.59004562381 21.9897989234
127 1 1.42678673058 9.06131273274 21.9897989234
128 1 4.28036019174 9.06131273274 21.9897989234
129 1 7.1339336529 9.06131273274 21.9897989234
130 1 9.98750711406 9.06131273274 21.9897989234
131 1 12.8410805752 9.06131273274 21.9897989234
132 1 15.6946540364 9.06131273274 21.9897989234
133 1 0.0 11.5325798417 21.9897989234
134 1 2.85357346116 11.5325798417 21.9897989234
135 1 5.70714692232 11.5325798417 21.9897989234
136 1 8.56072038348 11.5325798417 21.9897989234
137 1 11.4142938446 11.5325798417 21.9897989234
138 1 14.2678673058 11.5325798417 21.9897989234
139 1 1.42678673058 14.0038469506 21.9897989234
140 1 4.28036019174 14.0038469506 21.9897989234
141 1 7.1339336529 14.0038469506 21.9897989234
142 1 9.98750711406 14.0038469506 21.9897989234
143 1 12.8410805752 14.0038469506 21.9897989234
144 1 15.6946540364 14.0038469506 21.9897989234
145 1 0.0 0.0 24.3197318979
146 1 2.85357346116 0.0 24.3197318979
147 1 5.70714692232 0.0 24.3197318979
148 1 8.56072038348 0.0 24.3197318979
149 1 11.4142938446 0.0 24.3197318979
150 1 14.2678673058 0.0 24.3197318979
151 1 1.42678673058 2.47126710893 24.3197318979
152 1 4.28036019174 2.47126710893 24.3197318979
153 1 7.1339336529 2.47126710893 24.3197318979
154 1 9.98750711406 2.47126710893 24.3197318979
155 1 12.8410805752 2.47126710893 24.3197318979
156 1 15.6946540364 2.47126710893 24.3197318979
157 1 0.0 4.94253421786 24.3197318979
158 1 2.85357346116 4.94253421786 24.3197318979
159 1 5.70714692232 4.94253421786 24.3197318979
160 1 8.56072038348 4.94253421786 24.3197318979
161 1 11.4142938446 4.94253421786 24.3197318979
162 1 14.2678673058 4.94253421786 24.3197318979
163 1 1.42678673058 7.41380132679 24.3197318979
164 1 4.28036019174 7.41380132679 24.3197318979
165 1 7.1339336529 7.41380132679 24.3197318979
166 1 9.98750711406 7.41380132679 24.3197318979
167 1 12.8410805752 7.41380132679 24.3197318979
168 1 15.6946540364 7.41380132679 24.3197318979
169 1 0.0 9.88506843572 24.3197318979
170 1 2.85357346116 9.88506843572 24.3197318979
171 1 5.70714692232 9.88506843572 24.3197318979
172 1 8.56072038348 9.88506843572 24.3197318979
173 1 11.4142938446 9.88506843572 24.3197318979
174 1 14.2678673058 9.88506843572 24.3197318979
175 1 1.42678673058 12.3563355446 24.3197318979
176 1 4.28036019174 12.3563355446 24.3197318979
177 1 7.1339336529 12.3563355446 24.3197318979
178 1 9.98750711406 12.3563355446 24.3197318979
179 1 12.8410805752 12.3563355446 24.3197318979
180 1 15.6946540364 12.3563355446 24.3197318979
181 1 7.1339336529 4.11877851488 26.7197318979

View File

@ -0,0 +1,23 @@
processors * * 1
units metal
boundary p p f
read_data adatom.data
pair_style agni
pair_coeff * * Al_prb.agni Al
neighbor 0.3 bin
neigh_modify delay 2 check yes
timestep 0.0005
velocity all create 500 12345
fix 1 all nvt temp 250 250 0.2
fix 5 all momentum 1 linear 1 1 1
thermo 100
thermo_style custom step ke temp
# dump MyDump all custom 250 dump.atoms id type x y z vx vy vz fx fy fz
run 1000

View File

@ -0,0 +1,23 @@
units metal
boundary p p p
read_data vacancy.data
pair_style agni
pair_coeff * * Al_prb.agni Al
neighbor 0.3 bin
neigh_modify delay 2 check yes
timestep 0.0005
velocity all create 1000 12345
fix 1 all nvt temp 900 900 200
fix 5 all momentum 1 linear 1 1 1
thermo 100
thermo_style custom step ke etotal temp
# dump MyDump all custom 250 dump.atoms id type x y z vx vy vz fx fy fz
run 1000

View File

@ -0,0 +1,82 @@
LAMMPS (9 Nov 2016)
using 1 OpenMP thread(s) per MPI task
processors * * 1
units metal
boundary p p f
read_data adatom.data
orthogonal box = (0 0 0) to (17.1214 14.8276 39.3197)
1 by 1 by 1 MPI processor grid
reading atoms ...
181 atoms
pair_style agni
pair_coeff * * Al_prb.agni Al
Reading potential file Al_prb.agni with DATE: 2016-11-11
neighbor 0.3 bin
neigh_modify delay 2 check yes
timestep 0.0005
velocity all create 500 12345
fix 1 all nvt temp 250 250 0.2
fix 5 all momentum 1 linear 1 1 1
thermo 100
thermo_style custom step ke temp
# dump MyDump all custom 250 dump.atoms id type x y z vx vy vz fx fy fz
run 1000
Neighbor list info ...
1 neighbor list requests
update every 1 steps, delay 2 steps, check yes
max neighbors/atom: 2000, page size: 100000
master list distance cutoff = 8.3
ghost atom cutoff = 8.3
binsize = 4.15 -> bins = 5 4 10
Memory usage per processor = 2.37049 Mbytes
Step KinEng Temp
0 11.633413 500
100 4.049974 174.06646
200 5.8983472 253.50889
300 5.3667309 230.66021
400 4.9343935 212.0785
500 5.4054496 232.32432
600 6.1779127 265.52452
700 6.3749266 273.9921
800 6.0701481 260.89283
900 6.4582394 277.57286
1000 6.4047444 275.27366
Loop time of 20.8273 on 1 procs for 1000 steps with 181 atoms
Performance: 2.074 ns/day, 11.571 hours/ns, 48.014 timesteps/s
100.0% CPU use with 1 MPI tasks x 1 OpenMP threads
MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total
---------------------------------------------------------------
Pair | 20.79 | 20.79 | 20.79 | 0.0 | 99.82
Neigh | 0.022742 | 0.022742 | 0.022742 | 0.0 | 0.11
Comm | 0.0040836 | 0.0040836 | 0.0040836 | 0.0 | 0.02
Output | 0.00011086 | 0.00011086 | 0.00011086 | 0.0 | 0.00
Modify | 0.0089345 | 0.0089345 | 0.0089345 | 0.0 | 0.04
Other | | 0.001819 | | | 0.01
Nlocal: 181 ave 181 max 181 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Nghost: 562 ave 562 max 562 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Neighs: 0 ave 0 max 0 min
Histogram: 1 0 0 0 0 0 0 0 0 0
FullNghs: 18810 ave 18810 max 18810 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Total # of neighbors = 18810
Ave neighs/atom = 103.923
Neighbor list builds = 33
Dangerous builds = 0
Please see the log.cite file for references relevant to this simulation
Total wall time: 0:00:20

View File

@ -0,0 +1,82 @@
LAMMPS (9 Nov 2016)
using 1 OpenMP thread(s) per MPI task
processors * * 1
units metal
boundary p p f
read_data adatom.data
orthogonal box = (0 0 0) to (17.1214 14.8276 39.3197)
2 by 2 by 1 MPI processor grid
reading atoms ...
181 atoms
pair_style agni
pair_coeff * * Al_prb.agni Al
Reading potential file Al_prb.agni with DATE: 2016-11-11
neighbor 0.3 bin
neigh_modify delay 2 check yes
timestep 0.0005
velocity all create 500 12345
fix 1 all nvt temp 250 250 0.2
fix 5 all momentum 1 linear 1 1 1
thermo 100
thermo_style custom step ke temp
# dump MyDump all custom 250 dump.atoms id type x y z vx vy vz fx fy fz
run 1000
Neighbor list info ...
1 neighbor list requests
update every 1 steps, delay 2 steps, check yes
max neighbors/atom: 2000, page size: 100000
master list distance cutoff = 8.3
ghost atom cutoff = 8.3
binsize = 4.15 -> bins = 5 4 10
Memory usage per processor = 2.48695 Mbytes
Step KinEng Temp
0 11.633413 500
100 4.049974 174.06646
200 5.8983472 253.50889
300 5.3667309 230.66021
400 4.9343935 212.0785
500 5.4054496 232.32432
600 6.1779127 265.52451
700 6.3749266 273.9921
800 6.0701481 260.89283
900 6.4582394 277.57286
1000 6.4047444 275.27366
Loop time of 5.96868 on 4 procs for 1000 steps with 181 atoms
Performance: 7.238 ns/day, 3.316 hours/ns, 167.541 timesteps/s
99.9% CPU use with 4 MPI tasks x 1 OpenMP threads
MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total
---------------------------------------------------------------
Pair | 5.4176 | 5.4602 | 5.5505 | 2.3 | 91.48
Neigh | 0.0056074 | 0.0060464 | 0.0062635 | 0.3 | 0.10
Comm | 0.39544 | 0.48696 | 0.53111 | 7.9 | 8.16
Output | 0.0001545 | 0.00015736 | 0.0001595 | 0.0 | 0.00
Modify | 0.010492 | 0.011565 | 0.012588 | 0.9 | 0.19
Other | | 0.003794 | | | 0.06
Nlocal: 45.25 ave 47 max 42 min
Histogram: 1 0 0 0 0 0 0 0 2 1
Nghost: 374.75 ave 380 max 373 min
Histogram: 3 0 0 0 0 0 0 0 0 1
Neighs: 0 ave 0 max 0 min
Histogram: 4 0 0 0 0 0 0 0 0 0
FullNghs: 4702.5 ave 4868 max 4389 min
Histogram: 1 0 0 0 0 0 0 1 1 1
Total # of neighbors = 18810
Ave neighs/atom = 103.923
Neighbor list builds = 33
Dangerous builds = 0
Please see the log.cite file for references relevant to this simulation
Total wall time: 0:00:05

View File

@ -0,0 +1,82 @@
LAMMPS (9 Nov 2016)
using 1 OpenMP thread(s) per MPI task
units metal
boundary p p p
read_data vac.data
orthogonal box = (0 0 0) to (8.07113 8.07113 8.07113)
1 by 1 by 1 MPI processor grid
reading atoms ...
31 atoms
pair_style agni
pair_coeff * * Al_prb.agni Al
Reading potential file Al_prb.agni with DATE: 2016-11-11
neighbor 0.3 bin
neigh_modify delay 2 check yes
timestep 0.0005
velocity all create 1000 12345
fix 1 all nvt temp 900 900 200
fix 5 all momentum 1 linear 1 1 1
thermo 100
thermo_style custom step ke etotal temp
# dump MyDump all custom 250 dump.atoms id type x y z vx vy vz fx fy fz
run 1000
Neighbor list info ...
1 neighbor list requests
update every 1 steps, delay 2 steps, check yes
max neighbors/atom: 2000, page size: 100000
master list distance cutoff = 8.3
ghost atom cutoff = 8.3
binsize = 4.15 -> bins = 2 2 2
Memory usage per processor = 2.40842 Mbytes
Step KinEng TotEng Temp
0 3.8778043 3.8778043 1000
100 2.8126642 2.8126642 725.32391
200 3.7110413 3.7110413 956.9955
300 3.2361084 3.2361084 834.52081
400 3.4625769 3.4625769 892.92201
500 3.4563307 3.4563307 891.31126
600 2.8486344 2.8486344 734.59982
700 3.1183057 3.1183057 804.14208
800 2.9164818 2.9164818 752.09618
900 3.464416 3.464416 893.39629
1000 3.5954546 3.5954546 927.18825
Loop time of 3.86777 on 1 procs for 1000 steps with 31 atoms
Performance: 11.169 ns/day, 2.149 hours/ns, 258.547 timesteps/s
99.9% CPU use with 1 MPI tasks x 1 OpenMP threads
MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total
---------------------------------------------------------------
Pair | 3.8463 | 3.8463 | 3.8463 | 0.0 | 99.44
Neigh | 0.011294 | 0.011294 | 0.011294 | 0.0 | 0.29
Comm | 0.0057271 | 0.0057271 | 0.0057271 | 0.0 | 0.15
Output | 0.00014257 | 0.00014257 | 0.00014257 | 0.0 | 0.00
Modify | 0.0029459 | 0.0029459 | 0.0029459 | 0.0 | 0.08
Other | | 0.001395 | | | 0.04
Nlocal: 31 ave 31 max 31 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Nghost: 878 ave 878 max 878 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Neighs: 0 ave 0 max 0 min
Histogram: 1 0 0 0 0 0 0 0 0 0
FullNghs: 4334 ave 4334 max 4334 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Total # of neighbors = 4334
Ave neighs/atom = 139.806
Neighbor list builds = 51
Dangerous builds = 0
Please see the log.cite file for references relevant to this simulation
Total wall time: 0:00:03

View File

@ -0,0 +1,82 @@
LAMMPS (9 Nov 2016)
using 1 OpenMP thread(s) per MPI task
units metal
boundary p p p
read_data vac.data
orthogonal box = (0 0 0) to (8.07113 8.07113 8.07113)
1 by 2 by 2 MPI processor grid
reading atoms ...
31 atoms
pair_style agni
pair_coeff * * Al_prb.agni Al
Reading potential file Al_prb.agni with DATE: 2016-11-11
neighbor 0.3 bin
neigh_modify delay 2 check yes
timestep 0.0005
velocity all create 1000 12345
fix 1 all nvt temp 900 900 200
fix 5 all momentum 1 linear 1 1 1
thermo 100
thermo_style custom step ke etotal temp
# dump MyDump all custom 250 dump.atoms id type x y z vx vy vz fx fy fz
run 1000
Neighbor list info ...
1 neighbor list requests
update every 1 steps, delay 2 steps, check yes
max neighbors/atom: 2000, page size: 100000
master list distance cutoff = 8.3
ghost atom cutoff = 8.3
binsize = 4.15 -> bins = 2 2 2
Memory usage per processor = 2.3974 Mbytes
Step KinEng TotEng Temp
0 3.8778044 3.8778044 1000
100 2.8126642 2.8126642 725.32391
200 3.7110413 3.7110413 956.99549
300 3.2361084 3.2361084 834.52081
400 3.4625769 3.4625769 892.92201
500 3.4563307 3.4563307 891.31126
600 2.8486343 2.8486343 734.59981
700 3.1183057 3.1183057 804.14208
800 2.9164819 2.9164819 752.09618
900 3.4644161 3.4644161 893.39631
1000 3.5954546 3.5954546 927.18824
Loop time of 1.11007 on 4 procs for 1000 steps with 31 atoms
Performance: 38.916 ns/day, 0.617 hours/ns, 900.843 timesteps/s
99.8% CPU use with 4 MPI tasks x 1 OpenMP threads
MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total
---------------------------------------------------------------
Pair | 0.906 | 0.99469 | 1.0291 | 5.1 | 89.61
Neigh | 0.0026186 | 0.0027665 | 0.0028622 | 0.2 | 0.25
Comm | 0.066125 | 0.10079 | 0.1896 | 16.2 | 9.08
Output | 0.00012875 | 0.00014615 | 0.00018787 | 0.2 | 0.01
Modify | 0.0080338 | 0.0083079 | 0.00861 | 0.2 | 0.75
Other | | 0.003372 | | | 0.30
Nlocal: 7.75 ave 8 max 7 min
Histogram: 1 0 0 0 0 0 0 0 0 3
Nghost: 623.5 ave 630 max 616 min
Histogram: 1 0 0 0 0 2 0 0 0 1
Neighs: 0 ave 0 max 0 min
Histogram: 4 0 0 0 0 0 0 0 0 0
FullNghs: 1083.5 ave 1131 max 988 min
Histogram: 1 0 0 0 0 0 0 1 1 1
Total # of neighbors = 4334
Ave neighs/atom = 139.806
Neighbor list builds = 51
Dangerous builds = 0
Please see the log.cite file for references relevant to this simulation
Total wall time: 0:00:01

View File

@ -0,0 +1,46 @@
Data File created from VASP POSCAR
31 atoms
1 atom types
0 8.071125 xlo xhi
0 8.071125 ylo yhi
0 8.071125 zlo zhi
Masses
1 26.9815
Atoms
1 1 8.05711986217 3.20498589607 7.09652861184
2 1 8.05262255028 3.62006786258 3.16719841667
3 1 2.08891866821 1.38430927213 3.14852514324
4 1 4.25446836692 3.27689661974 3.35678388118
5 1 7.92524269451 7.20500664579 3.03232792051
6 1 6.04056771113 7.24499020906 1.11223380379
7 1 2.32585852889 5.29910389395 7.31500292009
8 1 2.09613190567 1.27658214906 7.44277603054
9 1 3.96852985867 7.2805082905 3.37568009522
10 1 0.0773420461671 1.29964047903 5.27451616984
11 1 7.96501442334 1.24471347504 1.17853896176
12 1 2.13035246804 5.36148411996 3.3817805118
13 1 2.06211525033 7.25482811482 1.52039033766
14 1 3.99735704234 7.4099829467 7.05753768668
15 1 3.84113228596 5.1855444403 1.41642147402
16 1 0.231862769544 5.38528175164 5.51171817022
17 1 0.12718452785 5.35814065671 1.11669573581
18 1 8.05303937039 7.38861123542 7.41398359808
19 1 1.88506066609 3.17578974033 1.20929473631
20 1 4.33739926831 1.37976783613 5.28141762358
21 1 2.23200994743 3.12419127088 5.36881641316
22 1 6.22871004896 1.34968648416 7.24032447626
23 1 6.08380394159 1.16222146146 3.30535465675
24 1 6.16629028099 5.22806528503 3.7675179617
25 1 4.30194966153 1.14526017671 1.45054175732
26 1 6.24221620153 5.05377575942 7.17573714759
27 1 3.92820642281 2.9627641757 7.71515743722
28 1 4.33686872315 4.73096617728 5.57649231331
29 1 6.05033104136 3.51389714904 1.34127903322
30 1 6.27311587476 7.19257797516 5.46814369382
31 1 1.81274009101 7.47392095028 5.35484578074

View File

@ -114,3 +114,35 @@ neighbor 1.0 bin
neigh_modify every 1 delay 10 check yes
run 100
# Test Tersoff/Mod model for Si
clear
read_restart restart.equil
pair_style tersoff/mod
pair_coeff * * Si.tersoff.mod 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
# Test Tersoff/Mod/C model for Si
clear
read_restart restart.equil
pair_style tersoff/mod/c
pair_coeff * * Si.tersoff.modc 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

1019
potentials/Al_prb.agni Normal file

File diff suppressed because it is too large Load Diff

View File

@ -0,0 +1,10 @@
# DATE: 2016-11-09 CONTRIBUTOR: Ganga P Purja Pun (George Mason University, Fairfax) CITATION: Unknown
#
# Format:
# element1 element2 element3
# beta alpha h eta
# beta_ters lam2 B R D lam1 A
# n c1 c2 c3 c4 c5 C
Si Si Si 3.00000000 1.80536502 -0.38136087 2.16152496
1 1.39343356 117.78072440 2.87478837 0.33090566 3.18011795 3198.51383127
1.98633876 0.20123243 614230.04310619 996439.09714140 3.33560562 25.20963770 -0.00592042

2
src/.gitignore vendored
View File

@ -539,6 +539,8 @@
/pair_adp.cpp
/pair_adp.h
/pair_agni.cpp
/pair_agni.h
/pair_airebo.cpp
/pair_airebo.h
/pair_airebo_morse.cpp

View File

@ -53,6 +53,7 @@ fi
if (test $1 = "CLASS2") then
depend GPU
depend KOKKOS
depend USER-OMP
fi

View File

@ -169,6 +169,8 @@ action pair_reax_c_kokkos.cpp pair_reax_c.cpp
action pair_reax_c_kokkos.h pair_reax_c.h
action pair_sw_kokkos.cpp pair_sw.cpp
action pair_sw_kokkos.h pair_sw.h
action pair_vashishta_kokkos.cpp pair_vashishta.cpp
action pair_vashishta_kokkos.h pair_vashishta.h
action pair_table_kokkos.cpp
action pair_table_kokkos.h
action pair_tersoff_kokkos.cpp pair_tersoff.cpp

View File

@ -83,8 +83,13 @@ class AtomVecKokkos : public AtomVec {
std::is_same<typename ViewType::execution_space,LMPDeviceType>::value,
Kokkos::CudaHostPinnedSpace,typename ViewType::memory_space>::type,
Kokkos::MemoryTraits<Kokkos::Unmanaged> > mirror_type;
if(buffer_size < src.capacity())
if (buffer_size == 0) {
buffer = Kokkos::kokkos_malloc<Kokkos::CudaHostPinnedSpace>(src.capacity());
buffer_size = src.capacity();
} else if (buffer_size < src.capacity()) {
buffer = Kokkos::kokkos_realloc<Kokkos::CudaHostPinnedSpace>(buffer,src.capacity());
buffer_size = src.capacity();
}
return mirror_type( buffer ,
src.dimension_0() ,
src.dimension_1() ,
@ -104,8 +109,13 @@ class AtomVecKokkos : public AtomVec {
std::is_same<typename ViewType::execution_space,LMPDeviceType>::value,
Kokkos::CudaHostPinnedSpace,typename ViewType::memory_space>::type,
Kokkos::MemoryTraits<Kokkos::Unmanaged> > mirror_type;
if(buffer_size < src.capacity())
if (buffer_size == 0) {
buffer = Kokkos::kokkos_malloc<Kokkos::CudaHostPinnedSpace>(src.capacity()*sizeof(typename ViewType::value_type));
buffer_size = src.capacity();
} else if (buffer_size < src.capacity()) {
buffer = Kokkos::kokkos_realloc<Kokkos::CudaHostPinnedSpace>(buffer,src.capacity()*sizeof(typename ViewType::value_type));
buffer_size = src.capacity();
}
mirror_type tmp_view( (typename ViewType::value_type*)buffer ,
src.dimension_0() ,
src.dimension_1() ,

View File

@ -21,6 +21,11 @@
#include "atom_masks.h"
#include "error.h"
#include "kokkos.h"
#include "force.h"
#include "bond.h"
#include "angle.h"
#include "dihedral.h"
#include "improper.h"
using namespace LAMMPS_NS;
@ -601,7 +606,19 @@ void NeighborKokkos::build_topology_kokkos() {
k_anglelist.modify<LMPDeviceType>();
k_dihedrallist.modify<LMPDeviceType>();
k_improperlist.modify<LMPDeviceType>();
} else {
// Transfer topology neighbor lists to Host for non-Kokkos styles
if (force->bond && force->bond->execution_space == Host)
k_bondlist.sync<LMPHostType>();
if (force->angle && force->angle->execution_space == Host)
k_anglelist.sync<LMPHostType>();
if (force->dihedral && force->dihedral->execution_space == Host)
k_dihedrallist.sync<LMPHostType>();
if (force->improper && force->improper->execution_space == Host)
k_improperlist.sync<LMPHostType>();
} else {
neighbond_host.build_topology_kk();
k_bondlist = neighbond_host.k_bondlist;

View File

@ -126,26 +126,26 @@ void PairTersoffKokkos<DeviceType>::setup_params()
for (i = 1; i <= n; i++)
for (j = 1; j <= n; j++)
for (k = 1; k <= n; k++) {
m = elem2param[i-1][j-1][k-1];
k_params.h_view(i,j,k).powerm = params[m].powerm;
k_params.h_view(i,j,k).gamma = params[m].gamma;
k_params.h_view(i,j,k).lam3 = params[m].lam3;
k_params.h_view(i,j,k).c = params[m].c;
k_params.h_view(i,j,k).d = params[m].d;
k_params.h_view(i,j,k).h = params[m].h;
k_params.h_view(i,j,k).powern = params[m].powern;
k_params.h_view(i,j,k).beta = params[m].beta;
k_params.h_view(i,j,k).lam2 = params[m].lam2;
k_params.h_view(i,j,k).bigb = params[m].bigb;
k_params.h_view(i,j,k).bigr = params[m].bigr;
k_params.h_view(i,j,k).bigd = params[m].bigd;
k_params.h_view(i,j,k).lam1 = params[m].lam1;
k_params.h_view(i,j,k).biga = params[m].biga;
k_params.h_view(i,j,k).cutsq = params[m].cutsq;
k_params.h_view(i,j,k).c1 = params[m].c1;
k_params.h_view(i,j,k).c2 = params[m].c2;
k_params.h_view(i,j,k).c3 = params[m].c3;
k_params.h_view(i,j,k).c4 = params[m].c4;
m = elem2param[map[i]][map[j]][map[k]];
k_params.h_view(i,j,k).powerm = params[m].powerm;
k_params.h_view(i,j,k).gamma = params[m].gamma;
k_params.h_view(i,j,k).lam3 = params[m].lam3;
k_params.h_view(i,j,k).c = params[m].c;
k_params.h_view(i,j,k).d = params[m].d;
k_params.h_view(i,j,k).h = params[m].h;
k_params.h_view(i,j,k).powern = params[m].powern;
k_params.h_view(i,j,k).beta = params[m].beta;
k_params.h_view(i,j,k).lam2 = params[m].lam2;
k_params.h_view(i,j,k).bigb = params[m].bigb;
k_params.h_view(i,j,k).bigr = params[m].bigr;
k_params.h_view(i,j,k).bigd = params[m].bigd;
k_params.h_view(i,j,k).lam1 = params[m].lam1;
k_params.h_view(i,j,k).biga = params[m].biga;
k_params.h_view(i,j,k).cutsq = params[m].cutsq;
k_params.h_view(i,j,k).c1 = params[m].c1;
k_params.h_view(i,j,k).c2 = params[m].c2;
k_params.h_view(i,j,k).c3 = params[m].c3;
k_params.h_view(i,j,k).c4 = params[m].c4;
}
k_params.template modify<LMPHostType>();

View File

@ -125,27 +125,27 @@ void PairTersoffMODKokkos<DeviceType>::setup_params()
for (i = 1; i <= n; i++)
for (j = 1; j <= n; j++)
for (k = 1; k <= n; k++) {
m = elem2param[i-1][j-1][k-1];
k_params.h_view(i,j,k).powerm = params[m].powerm;
k_params.h_view(i,j,k).lam3 = params[m].lam3;
k_params.h_view(i,j,k).h = params[m].h;
k_params.h_view(i,j,k).powern = params[m].powern;
k_params.h_view(i,j,k).beta = params[m].beta;
k_params.h_view(i,j,k).lam2 = params[m].lam2;
k_params.h_view(i,j,k).bigb = params[m].bigb;
k_params.h_view(i,j,k).bigr = params[m].bigr;
k_params.h_view(i,j,k).bigd = params[m].bigd;
k_params.h_view(i,j,k).lam1 = params[m].lam1;
k_params.h_view(i,j,k).biga = params[m].biga;
k_params.h_view(i,j,k).cutsq = params[m].cutsq;
k_params.h_view(i,j,k).c1 = params[m].c1;
k_params.h_view(i,j,k).c2 = params[m].c2;
k_params.h_view(i,j,k).c3 = params[m].c3;
k_params.h_view(i,j,k).c4 = params[m].c4;
k_params.h_view(i,j,k).c5 = params[m].c5;
k_params.h_view(i,j,k).ca1 = params[m].ca1;
k_params.h_view(i,j,k).ca4 = params[m].ca4;
k_params.h_view(i,j,k).powern_del = params[m].powern_del;
m = elem2param[map[i]][map[j]][map[k]];
k_params.h_view(i,j,k).powerm = params[m].powerm;
k_params.h_view(i,j,k).lam3 = params[m].lam3;
k_params.h_view(i,j,k).h = params[m].h;
k_params.h_view(i,j,k).powern = params[m].powern;
k_params.h_view(i,j,k).beta = params[m].beta;
k_params.h_view(i,j,k).lam2 = params[m].lam2;
k_params.h_view(i,j,k).bigb = params[m].bigb;
k_params.h_view(i,j,k).bigr = params[m].bigr;
k_params.h_view(i,j,k).bigd = params[m].bigd;
k_params.h_view(i,j,k).lam1 = params[m].lam1;
k_params.h_view(i,j,k).biga = params[m].biga;
k_params.h_view(i,j,k).cutsq = params[m].cutsq;
k_params.h_view(i,j,k).c1 = params[m].c1;
k_params.h_view(i,j,k).c2 = params[m].c2;
k_params.h_view(i,j,k).c3 = params[m].c3;
k_params.h_view(i,j,k).c4 = params[m].c4;
k_params.h_view(i,j,k).c5 = params[m].c5;
k_params.h_view(i,j,k).ca1 = params[m].ca1;
k_params.h_view(i,j,k).ca4 = params[m].ca4;
k_params.h_view(i,j,k).powern_del = params[m].powern_del;
}
k_params.template modify<LMPHostType>();

View File

@ -136,30 +136,30 @@ void PairTersoffZBLKokkos<DeviceType>::setup_params()
for (i = 1; i <= n; i++)
for (j = 1; j <= n; j++)
for (k = 1; k <= n; k++) {
m = elem2param[i-1][j-1][k-1];
k_params.h_view(i,j,k).powerm = params[m].powerm;
k_params.h_view(i,j,k).gamma = params[m].gamma;
k_params.h_view(i,j,k).lam3 = params[m].lam3;
k_params.h_view(i,j,k).c = params[m].c;
k_params.h_view(i,j,k).d = params[m].d;
k_params.h_view(i,j,k).h = params[m].h;
k_params.h_view(i,j,k).powern = params[m].powern;
k_params.h_view(i,j,k).beta = params[m].beta;
k_params.h_view(i,j,k).lam2 = params[m].lam2;
k_params.h_view(i,j,k).bigb = params[m].bigb;
k_params.h_view(i,j,k).bigr = params[m].bigr;
k_params.h_view(i,j,k).bigd = params[m].bigd;
k_params.h_view(i,j,k).lam1 = params[m].lam1;
k_params.h_view(i,j,k).biga = params[m].biga;
k_params.h_view(i,j,k).cutsq = params[m].cutsq;
k_params.h_view(i,j,k).c1 = params[m].c1;
k_params.h_view(i,j,k).c2 = params[m].c2;
k_params.h_view(i,j,k).c3 = params[m].c3;
k_params.h_view(i,j,k).c4 = params[m].c4;
k_params.h_view(i,j,k).Z_i = params[m].Z_i;
k_params.h_view(i,j,k).Z_j = params[m].Z_j;
k_params.h_view(i,j,k).ZBLcut = params[m].ZBLcut;
k_params.h_view(i,j,k).ZBLexpscale = params[m].ZBLexpscale;
m = elem2param[map[i]][map[j]][map[k]];
k_params.h_view(i,j,k).powerm = params[m].powerm;
k_params.h_view(i,j,k).gamma = params[m].gamma;
k_params.h_view(i,j,k).lam3 = params[m].lam3;
k_params.h_view(i,j,k).c = params[m].c;
k_params.h_view(i,j,k).d = params[m].d;
k_params.h_view(i,j,k).h = params[m].h;
k_params.h_view(i,j,k).powern = params[m].powern;
k_params.h_view(i,j,k).beta = params[m].beta;
k_params.h_view(i,j,k).lam2 = params[m].lam2;
k_params.h_view(i,j,k).bigb = params[m].bigb;
k_params.h_view(i,j,k).bigr = params[m].bigr;
k_params.h_view(i,j,k).bigd = params[m].bigd;
k_params.h_view(i,j,k).lam1 = params[m].lam1;
k_params.h_view(i,j,k).biga = params[m].biga;
k_params.h_view(i,j,k).cutsq = params[m].cutsq;
k_params.h_view(i,j,k).c1 = params[m].c1;
k_params.h_view(i,j,k).c2 = params[m].c2;
k_params.h_view(i,j,k).c3 = params[m].c3;
k_params.h_view(i,j,k).c4 = params[m].c4;
k_params.h_view(i,j,k).Z_i = params[m].Z_i;
k_params.h_view(i,j,k).Z_j = params[m].Z_j;
k_params.h_view(i,j,k).ZBLcut = params[m].ZBLcut;
k_params.h_view(i,j,k).ZBLexpscale = params[m].ZBLexpscale;
}
k_params.template modify<LMPHostType>();

View File

@ -0,0 +1,961 @@
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
Contributing author: Stan Moore (SNL)
Anders Hafreager (UiO), andershaf@gmail.com
------------------------------------------------------------------------- */
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include "pair_vashishta_kokkos.h"
#include "kokkos.h"
#include "pair_kokkos.h"
#include "atom_kokkos.h"
#include "neighbor.h"
#include "neigh_request.h"
#include "force.h"
#include "comm.h"
#include "memory.h"
#include "neighbor.h"
#include "neigh_list_kokkos.h"
#include "memory.h"
#include "error.h"
#include "atom_masks.h"
#include "math_const.h"
using namespace LAMMPS_NS;
using namespace MathConst;
#define MAXLINE 1024
#define DELTA 4
/* ---------------------------------------------------------------------- */
template<class DeviceType>
PairVashishtaKokkos<DeviceType>::PairVashishtaKokkos(LAMMPS *lmp) : PairVashishta(lmp)
{
respa_enable = 0;
atomKK = (AtomKokkos *) atom;
execution_space = ExecutionSpaceFromDevice<DeviceType>::space;
datamask_read = X_MASK | F_MASK | TAG_MASK | TYPE_MASK | ENERGY_MASK | VIRIAL_MASK;
datamask_modify = F_MASK | ENERGY_MASK | VIRIAL_MASK;
}
/* ----------------------------------------------------------------------
check if allocated, since class can be destructed when incomplete
------------------------------------------------------------------------- */
template<class DeviceType>
PairVashishtaKokkos<DeviceType>::~PairVashishtaKokkos()
{
if (!copymode) {
memory->destroy_kokkos(k_eatom,eatom);
memory->destroy_kokkos(k_vatom,vatom);
eatom = NULL;
vatom = NULL;
}
}
/* ---------------------------------------------------------------------- */
template<class DeviceType>
void PairVashishtaKokkos<DeviceType>::compute(int eflag_in, int vflag_in)
{
eflag = eflag_in;
vflag = vflag_in;
if (neighflag == FULL) no_virial_fdotr_compute = 1;
if (eflag || vflag) ev_setup(eflag,vflag);
else evflag = vflag_fdotr = 0;
// reallocate per-atom arrays if necessary
if (eflag_atom) {
memory->destroy_kokkos(k_eatom,eatom);
memory->create_kokkos(k_eatom,eatom,maxeatom,"pair:eatom");
d_eatom = k_eatom.d_view;
}
if (vflag_atom) {
memory->destroy_kokkos(k_vatom,vatom);
memory->create_kokkos(k_vatom,vatom,maxvatom,6,"pair:vatom");
d_vatom = k_vatom.d_view;
}
atomKK->sync(execution_space,datamask_read);
if (eflag || vflag) atomKK->modified(execution_space,datamask_modify);
else atomKK->modified(execution_space,F_MASK);
x = atomKK->k_x.view<DeviceType>();
f = atomKK->k_f.view<DeviceType>();
tag = atomKK->k_tag.view<DeviceType>();
type = atomKK->k_type.view<DeviceType>();
nlocal = atom->nlocal;
newton_pair = force->newton_pair;
nall = atom->nlocal + atom->nghost;
inum = list->inum;
const int ignum = inum + list->gnum;
NeighListKokkos<DeviceType>* k_list = static_cast<NeighListKokkos<DeviceType>*>(list);
d_ilist = k_list->d_ilist;
d_numneigh = k_list->d_numneigh;
d_neighbors = k_list->d_neighbors;
k_list->clean_copy();
copymode = 1;
EV_FLOAT ev;
EV_FLOAT ev_all;
int max_neighs = d_neighbors.dimension_1();
if ((d_neighbors_short_2body.dimension_1() != max_neighs) ||
(d_neighbors_short_2body.dimension_0() != ignum)) {
d_neighbors_short_2body = Kokkos::View<int**,DeviceType>("Vashishta::neighbors_short_2body",ignum,max_neighs);
}
if (d_numneigh_short_2body.dimension_0()!=ignum) {
d_numneigh_short_2body = Kokkos::View<int*,DeviceType>("Vashishta::numneighs_short_2body",ignum);
}
if ((d_neighbors_short_3body.dimension_1() != max_neighs) ||
(d_neighbors_short_3body.dimension_0() != ignum)) {
d_neighbors_short_3body = Kokkos::View<int**,DeviceType>("Vashishta::neighbors_short_3body",ignum,max_neighs);
}
if (d_numneigh_short_3body.dimension_0()!=ignum) {
d_numneigh_short_3body = Kokkos::View<int*,DeviceType>("Vashishta::numneighs_short_3body",ignum);
}
Kokkos::parallel_for(Kokkos::RangePolicy<DeviceType,TagPairVashishtaComputeShortNeigh>(0,neighflag==FULL?ignum:inum), *this);
// loop over neighbor list of my atoms
if (neighflag == HALF) {
if (evflag)
Kokkos::parallel_reduce(Kokkos::RangePolicy<DeviceType, TagPairVashishtaComputeHalf<HALF,1> >(0,inum),*this,ev);
else
Kokkos::parallel_for(Kokkos::RangePolicy<DeviceType, TagPairVashishtaComputeHalf<HALF,0> >(0,inum),*this);
ev_all += ev;
} else if (neighflag == HALFTHREAD) {
if (evflag)
Kokkos::parallel_reduce(Kokkos::RangePolicy<DeviceType, TagPairVashishtaComputeHalf<HALFTHREAD,1> >(0,inum),*this,ev);
else
Kokkos::parallel_for(Kokkos::RangePolicy<DeviceType, TagPairVashishtaComputeHalf<HALFTHREAD,0> >(0,inum),*this);
ev_all += ev;
} else if (neighflag == FULL) {
if (evflag)
Kokkos::parallel_reduce(Kokkos::RangePolicy<DeviceType, TagPairVashishtaComputeFullA<FULL,1> >(0,inum),*this,ev);
else
Kokkos::parallel_for(Kokkos::RangePolicy<DeviceType, TagPairVashishtaComputeFullA<FULL,0> >(0,inum),*this);
ev_all += ev;
if (evflag)
Kokkos::parallel_reduce(Kokkos::RangePolicy<DeviceType, TagPairVashishtaComputeFullB<FULL,1> >(0,ignum),*this,ev);
else
Kokkos::parallel_for(Kokkos::RangePolicy<DeviceType, TagPairVashishtaComputeFullB<FULL,0> >(0,ignum),*this);
ev_all += ev;
}
if (eflag_global) eng_vdwl += ev_all.evdwl;
if (vflag_global) {
virial[0] += ev_all.v[0];
virial[1] += ev_all.v[1];
virial[2] += ev_all.v[2];
virial[3] += ev_all.v[3];
virial[4] += ev_all.v[4];
virial[5] += ev_all.v[5];
}
if (vflag_fdotr) pair_virial_fdotr_compute(this);
if (eflag_atom) {
k_eatom.template modify<DeviceType>();
k_eatom.template sync<LMPHostType>();
}
if (vflag_atom) {
k_vatom.template modify<DeviceType>();
k_vatom.template sync<LMPHostType>();
}
copymode = 0;
}
template<class DeviceType>
KOKKOS_INLINE_FUNCTION
void PairVashishtaKokkos<DeviceType>::operator()(TagPairVashishtaComputeShortNeigh, const int& ii) const {
const int i = d_ilist[ii];
const int itype = d_map[type[i]];
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_2body = 0;
int inside_3body = 0;
for (int jj = 0; jj < jnum; jj++) {
int j = d_neighbors(i,jj);
j &= NEIGHMASK;
const int jtype = d_map[type[j]];
const int ijparam = d_elem2param(itype,jtype,jtype);
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 < d_params[ijparam].cutsq) {
d_neighbors_short_2body(i,inside_2body) = j;
inside_2body++;
}
if (rsq < d_params[ijparam].cutsq2) {
d_neighbors_short_3body(i,inside_3body) = j;
inside_3body++;
}
}
d_numneigh_short_2body(i) = inside_2body;
d_numneigh_short_3body(i) = inside_3body;
}
template<class DeviceType>
template<int NEIGHFLAG, int EVFLAG>
KOKKOS_INLINE_FUNCTION
void PairVashishtaKokkos<DeviceType>::operator()(TagPairVashishtaComputeHalf<NEIGHFLAG,EVFLAG>, const int &ii, EV_FLOAT& ev) const {
// The f array is atomic
Kokkos::View<F_FLOAT*[3], typename DAT::t_f_array::array_layout,DeviceType,Kokkos::MemoryTraits<AtomicF<NEIGHFLAG>::value> > a_f = f;
F_FLOAT delr1[3],delr2[3],fj[3],fk[3];
F_FLOAT evdwl = 0.0;
F_FLOAT fpair = 0.0;
const int i = d_ilist[ii];
const tagint itag = tag[i];
const int itype = d_map[type[i]];
const X_FLOAT xtmp = x(i,0);
const X_FLOAT ytmp = x(i,1);
const X_FLOAT ztmp = x(i,2);
// two-body interactions, skip half of them
const int jnum = d_numneigh_short_2body[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_short_2body(i,jj);
j &= NEIGHMASK;
const tagint jtag = tag[j];
if (itag > jtag) {
if ((itag+jtag) % 2 == 0) continue;
} else if (itag < jtag) {
if ((itag+jtag) % 2 == 1) continue;
} else {
if (x(j,2) < ztmp) continue;
if (x(j,2) == ztmp && x(j,1) < ytmp) continue;
if (x(j,2) == ztmp && x(j,1) == ytmp && x(j,0) < xtmp) continue;
}
const int jtype = d_map[type[j]];
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;
const int ijparam = d_elem2param(itype,jtype,jtype);
twobody(d_params[ijparam],rsq,fpair,eflag,evdwl);
fxtmpi += delx*fpair;
fytmpi += dely*fpair;
fztmpi += delz*fpair;
a_f(j,0) -= delx*fpair;
a_f(j,1) -= dely*fpair;
a_f(j,2) -= delz*fpair;
if (EVFLAG) {
if (eflag) ev.evdwl += evdwl;
if (vflag_either || eflag_atom) this->template ev_tally<NEIGHFLAG>(ev,i,j,evdwl,fpair,delx,dely,delz);
}
}
const int jnumm1 = d_numneigh_short_3body[i];
for (int jj = 0; jj < jnumm1-1; jj++) {
int j = d_neighbors_short_3body(i,jj);
j &= NEIGHMASK;
const int jtype = d_map[type[j]];
const int ijparam = d_elem2param(itype,jtype,jtype);
delr1[0] = x(j,0) - xtmp;
delr1[1] = x(j,1) - ytmp;
delr1[2] = x(j,2) - ztmp;
const F_FLOAT rsq1 = delr1[0]*delr1[0] + delr1[1]*delr1[1] + delr1[2]*delr1[2];
F_FLOAT fxtmpj = 0.0;
F_FLOAT fytmpj = 0.0;
F_FLOAT fztmpj = 0.0;
for (int kk = jj+1; kk < jnumm1; kk++) {
int k = d_neighbors_short_3body(i,kk);
k &= NEIGHMASK;
const int ktype = d_map[type[k]];
const int ikparam = d_elem2param(itype,ktype,ktype);
const int ijkparam = d_elem2param(itype,jtype,ktype);
delr2[0] = x(k,0) - xtmp;
delr2[1] = x(k,1) - ytmp;
delr2[2] = x(k,2) - ztmp;
const F_FLOAT rsq2 = delr2[0]*delr2[0] + delr2[1]*delr2[1] + delr2[2]*delr2[2];
threebody(d_params[ijparam],d_params[ikparam],d_params[ijkparam],
rsq1,rsq2,delr1,delr2,fj,fk,eflag,evdwl);
fxtmpi -= fj[0] + fk[0];
fytmpi -= fj[1] + fk[1];
fztmpi -= fj[2] + fk[2];
fxtmpj += fj[0];
fytmpj += fj[1];
fztmpj += fj[2];
a_f(k,0) += fk[0];
a_f(k,1) += fk[1];
a_f(k,2) += fk[2];
if (EVFLAG) {
if (eflag) ev.evdwl += evdwl;
if (vflag_either || eflag_atom) this->template ev_tally3<NEIGHFLAG>(ev,i,j,k,evdwl,0.0,fj,fk,delr1,delr2);
}
}
a_f(j,0) += fxtmpj;
a_f(j,1) += fytmpj;
a_f(j,2) += fztmpj;
}
a_f(i,0) += fxtmpi;
a_f(i,1) += fytmpi;
a_f(i,2) += fztmpi;
}
template<class DeviceType>
template<int NEIGHFLAG, int EVFLAG>
KOKKOS_INLINE_FUNCTION
void PairVashishtaKokkos<DeviceType>::operator()(TagPairVashishtaComputeHalf<NEIGHFLAG,EVFLAG>, const int &ii) const {
EV_FLOAT ev;
this->template operator()<NEIGHFLAG,EVFLAG>(TagPairVashishtaComputeHalf<NEIGHFLAG,EVFLAG>(), ii, ev);
}
/* ---------------------------------------------------------------------- */
template<class DeviceType>
template<int NEIGHFLAG, int EVFLAG>
KOKKOS_INLINE_FUNCTION
void PairVashishtaKokkos<DeviceType>::operator()(TagPairVashishtaComputeFullA<NEIGHFLAG,EVFLAG>, const int &ii, EV_FLOAT& ev) const {
F_FLOAT delr1[3],delr2[3],fj[3],fk[3];
F_FLOAT evdwl = 0.0;
F_FLOAT fpair = 0.0;
const int i = d_ilist[ii];
const tagint itag = tag[i];
const int itype = d_map[type[i]];
const X_FLOAT xtmp = x(i,0);
const X_FLOAT ytmp = x(i,1);
const X_FLOAT ztmp = x(i,2);
// two-body interactions
const int jnum = d_numneigh_short_2body[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_short_2body(i,jj);
j &= NEIGHMASK;
const tagint jtag = tag[j];
const int jtype = d_map[type[j]];
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;
const int ijparam = d_elem2param(itype,jtype,jtype);
twobody(d_params[ijparam],rsq,fpair,eflag,evdwl);
fxtmpi += delx*fpair;
fytmpi += dely*fpair;
fztmpi += delz*fpair;
if (EVFLAG) {
if (eflag) ev.evdwl += 0.5*evdwl;
if (vflag_either || eflag_atom) this->template ev_tally<NEIGHFLAG>(ev,i,j,evdwl,fpair,delx,dely,delz);
}
}
const int jnumm1 = d_numneigh_short_3body[i];
for (int jj = 0; jj < jnumm1-1; jj++) {
int j = d_neighbors_short_3body(i,jj);
j &= NEIGHMASK;
const int jtype = d_map[type[j]];
const int ijparam = d_elem2param(itype,jtype,jtype);
delr1[0] = x(j,0) - xtmp;
delr1[1] = x(j,1) - ytmp;
delr1[2] = x(j,2) - ztmp;
const F_FLOAT rsq1 = delr1[0]*delr1[0] + delr1[1]*delr1[1] + delr1[2]*delr1[2];
for (int kk = jj+1; kk < jnumm1; kk++) {
int k = d_neighbors_short_3body(i,kk);
k &= NEIGHMASK;
const int ktype = d_map[type[k]];
const int ikparam = d_elem2param(itype,ktype,ktype);
const int ijkparam = d_elem2param(itype,jtype,ktype);
delr2[0] = x(k,0) - xtmp;
delr2[1] = x(k,1) - ytmp;
delr2[2] = x(k,2) - ztmp;
const F_FLOAT rsq2 = delr2[0]*delr2[0] + delr2[1]*delr2[1] + delr2[2]*delr2[2];
threebody(d_params[ijparam],d_params[ikparam],d_params[ijkparam],
rsq1,rsq2,delr1,delr2,fj,fk,eflag,evdwl);
fxtmpi -= fj[0] + fk[0];
fytmpi -= fj[1] + fk[1];
fztmpi -= fj[2] + fk[2];
if (EVFLAG) {
if (eflag) ev.evdwl += evdwl;
if (vflag_either || eflag_atom) this->template ev_tally3<NEIGHFLAG>(ev,i,j,k,evdwl,0.0,fj,fk,delr1,delr2);
}
}
}
f(i,0) += fxtmpi;
f(i,1) += fytmpi;
f(i,2) += fztmpi;
}
template<class DeviceType>
template<int NEIGHFLAG, int EVFLAG>
KOKKOS_INLINE_FUNCTION
void PairVashishtaKokkos<DeviceType>::operator()(TagPairVashishtaComputeFullA<NEIGHFLAG,EVFLAG>, const int &ii) const {
EV_FLOAT ev;
this->template operator()<NEIGHFLAG,EVFLAG>(TagPairVashishtaComputeFullA<NEIGHFLAG,EVFLAG>(), ii, ev);
}
/* ---------------------------------------------------------------------- */
template<class DeviceType>
template<int NEIGHFLAG, int EVFLAG>
KOKKOS_INLINE_FUNCTION
void PairVashishtaKokkos<DeviceType>::operator()(TagPairVashishtaComputeFullB<NEIGHFLAG,EVFLAG>, const int &ii, EV_FLOAT& ev) const {
F_FLOAT delr1[3],delr2[3],fj[3],fk[3];
F_FLOAT evdwl = 0.0;
F_FLOAT fpair = 0.0;
const int i = d_ilist[ii];
const int itype = d_map[type[i]];
const X_FLOAT xtmpi = x(i,0);
const X_FLOAT ytmpi = x(i,1);
const X_FLOAT ztmpi = x(i,2);
const int jnum = d_numneigh_short_3body[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_short_3body(i,jj);
j &= NEIGHMASK;
if (j >= nlocal) continue;
const int jtype = d_map[type[j]];
const int jiparam = d_elem2param(jtype,itype,itype);
const X_FLOAT xtmpj = x(j,0);
const X_FLOAT ytmpj = x(j,1);
const X_FLOAT ztmpj = x(j,2);
delr1[0] = xtmpi - xtmpj;
delr1[1] = ytmpi - ytmpj;
delr1[2] = ztmpi - ztmpj;
const F_FLOAT rsq1 = delr1[0]*delr1[0] + delr1[1]*delr1[1] + delr1[2]*delr1[2];
const int j_jnum = d_numneigh_short_3body[j];
for (int kk = 0; kk < j_jnum; kk++) {
int k = d_neighbors_short_3body(j,kk);
k &= NEIGHMASK;
if (k == i) continue;
const int ktype = d_map[type[k]];
const int jkparam = d_elem2param(jtype,ktype,ktype);
const int jikparam = d_elem2param(jtype,itype,ktype);
delr2[0] = x(k,0) - xtmpj;
delr2[1] = x(k,1) - ytmpj;
delr2[2] = x(k,2) - ztmpj;
const F_FLOAT rsq2 = delr2[0]*delr2[0] + delr2[1]*delr2[1] + delr2[2]*delr2[2];
if (vflag_atom)
threebody(d_params[jiparam],d_params[jkparam],d_params[jikparam],
rsq1,rsq2,delr1,delr2,fj,fk,eflag,evdwl);
else
threebodyj(d_params[jiparam],d_params[jkparam],d_params[jikparam],
rsq1,rsq2,delr1,delr2,fj);
fxtmpi += fj[0];
fytmpi += fj[1];
fztmpi += fj[2];
if (EVFLAG)
if (vflag_atom || eflag_atom) ev_tally3_atom(ev,i,evdwl,0.0,fj,fk,delr1,delr2);
}
}
f(i,0) += fxtmpi;
f(i,1) += fytmpi;
f(i,2) += fztmpi;
}
template<class DeviceType>
template<int NEIGHFLAG, int EVFLAG>
KOKKOS_INLINE_FUNCTION
void PairVashishtaKokkos<DeviceType>::operator()(TagPairVashishtaComputeFullB<NEIGHFLAG,EVFLAG>, const int &ii) const {
EV_FLOAT ev;
this->template operator()<NEIGHFLAG,EVFLAG>(TagPairVashishtaComputeFullB<NEIGHFLAG,EVFLAG>(), ii, ev);
}
/* ----------------------------------------------------------------------
set coeffs for one or more type pairs
------------------------------------------------------------------------- */
template<class DeviceType>
void PairVashishtaKokkos<DeviceType>::coeff(int narg, char **arg)
{
PairVashishta::coeff(narg,arg);
// sync map
int n = atom->ntypes;
DAT::tdual_int_1d k_map = DAT::tdual_int_1d("pair:map",n+1);
HAT::t_int_1d h_map = k_map.h_view;
for (int i = 1; i <= n; i++)
h_map[i] = map[i];
k_map.template modify<LMPHostType>();
k_map.template sync<DeviceType>();
d_map = k_map.d_view;
}
/* ----------------------------------------------------------------------
init specific to this pair style
------------------------------------------------------------------------- */
template<class DeviceType>
void PairVashishtaKokkos<DeviceType>::init_style()
{
PairVashishta::init_style();
// irequest = neigh request made by parent class
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;
// always request a full neighbor list
if (neighflag == FULL || neighflag == HALF || neighflag == HALFTHREAD) {
neighbor->requests[irequest]->full = 1;
neighbor->requests[irequest]->half = 0;
neighbor->requests[irequest]->full_cluster = 0;
if (neighflag == FULL)
neighbor->requests[irequest]->ghost = 1;
else
neighbor->requests[irequest]->ghost = 0;
} else {
error->all(FLERR,"Cannot use chosen neighbor list style with pair vashishta/kk");
}
}
/* ---------------------------------------------------------------------- */
template<class DeviceType>
void PairVashishtaKokkos<DeviceType>::setup_params()
{
PairVashishta::setup_params();
// sync elem2param and params
tdual_int_3d k_elem2param = tdual_int_3d("pair:elem2param",nelements,nelements,nelements);
t_host_int_3d h_elem2param = k_elem2param.h_view;
tdual_param_1d k_params = tdual_param_1d("pair:params",nparams);
t_host_param_1d h_params = k_params.h_view;
for (int i = 0; i < nelements; i++)
for (int j = 0; j < nelements; j++)
for (int k = 0; k < nelements; k++)
h_elem2param(i,j,k) = elem2param[i][j][k];
for (int m = 0; m < nparams; m++)
h_params[m] = params[m];
k_elem2param.template modify<LMPHostType>();
k_elem2param.template sync<DeviceType>();
k_params.template modify<LMPHostType>();
k_params.template sync<DeviceType>();
d_elem2param = k_elem2param.d_view;
d_params = k_params.d_view;
}
/* ---------------------------------------------------------------------- */
template<class DeviceType>
KOKKOS_INLINE_FUNCTION
void PairVashishtaKokkos<DeviceType>::twobody(const Param& param, const F_FLOAT& rsq, F_FLOAT& fforce,
const int& eflag, F_FLOAT& eng) const
{
F_FLOAT r,rinvsq,r4inv,r6inv,reta,lam1r,lam4r,vc2,vc3;
r = sqrt(rsq);
rinvsq = 1.0/rsq;
r4inv = rinvsq*rinvsq;
r6inv = rinvsq*r4inv;
reta = pow(r,-param.eta);
lam1r = r*param.lam1inv;
lam4r = r*param.lam4inv;
vc2 = param.zizj * exp(-lam1r)/r;
vc3 = param.mbigd * r4inv*exp(-lam4r);
fforce = (param.dvrc*r
- (4.0*vc3 + lam4r*vc3+param.big6w*r6inv
- param.heta*reta - vc2 - lam1r*vc2)
) * rinvsq;
if (eflag) eng = param.bigh*reta + vc2 - vc3 - param.bigw*r6inv - r*param.dvrc + param.c0;
}
/* ---------------------------------------------------------------------- */
template<class DeviceType>
KOKKOS_INLINE_FUNCTION
void PairVashishtaKokkos<DeviceType>::threebody(const Param& paramij, const Param& paramik, const Param& paramijk,
const F_FLOAT& rsq1, const F_FLOAT& rsq2,
F_FLOAT *delr1, F_FLOAT *delr2,
F_FLOAT *fj, F_FLOAT *fk, const int& eflag, F_FLOAT& eng) const
{
F_FLOAT r1,rinvsq1,rainv1,gsrainv1,gsrainvsq1,expgsrainv1;
F_FLOAT r2,rinvsq2,rainv2,gsrainv2,gsrainvsq2,expgsrainv2;
F_FLOAT rinv12,cs,delcs,delcssq,facexp,facrad,frad1,frad2,pcsinv,pcsinvsq,pcs;
F_FLOAT facang,facang12,csfacang,csfac1,csfac2;
r1 = sqrt(rsq1);
rinvsq1 = 1.0/rsq1;
rainv1 = 1.0/(r1 - paramij.r0);
gsrainv1 = paramij.gamma * rainv1;
gsrainvsq1 = gsrainv1*rainv1/r1;
expgsrainv1 = exp(gsrainv1);
r2 = sqrt(rsq2);
rinvsq2 = 1.0/rsq2;
rainv2 = 1.0/(r2 - paramik.r0);
gsrainv2 = paramik.gamma * rainv2;
gsrainvsq2 = gsrainv2*rainv2/r2;
expgsrainv2 = exp(gsrainv2);
rinv12 = 1.0/(r1*r2);
cs = (delr1[0]*delr2[0] + delr1[1]*delr2[1] + delr1[2]*delr2[2]) * rinv12;
delcs = cs - paramijk.costheta;
delcssq = delcs*delcs;
pcsinv = paramijk.bigc*delcssq + 1.0;
pcsinvsq = pcsinv*pcsinv;
pcs = delcssq/pcsinv;
facexp = expgsrainv1*expgsrainv2;
facrad = paramijk.bigb * facexp * pcs;
frad1 = facrad*gsrainvsq1;
frad2 = facrad*gsrainvsq2;
facang = paramijk.big2b * facexp * delcs/pcsinvsq;
facang12 = rinv12*facang;
csfacang = cs*facang;
csfac1 = rinvsq1*csfacang;
fj[0] = delr1[0]*(frad1+csfac1)-delr2[0]*facang12;
fj[1] = delr1[1]*(frad1+csfac1)-delr2[1]*facang12;
fj[2] = delr1[2]*(frad1+csfac1)-delr2[2]*facang12;
csfac2 = rinvsq2*csfacang;
fk[0] = delr2[0]*(frad2+csfac2)-delr1[0]*facang12;
fk[1] = delr2[1]*(frad2+csfac2)-delr1[1]*facang12;
fk[2] = delr2[2]*(frad2+csfac2)-delr1[2]*facang12;
if (eflag) eng = facrad;
}
/* ---------------------------------------------------------------------- */
template<class DeviceType>
KOKKOS_INLINE_FUNCTION
void PairVashishtaKokkos<DeviceType>::threebodyj(const Param& paramij, const Param& paramik, const Param& paramijk,
const F_FLOAT& rsq1, const F_FLOAT& rsq2, F_FLOAT *delr1, F_FLOAT *delr2, F_FLOAT *fj) const
{
F_FLOAT r1,rinvsq1,rainv1,gsrainv1,gsrainvsq1,expgsrainv1;
F_FLOAT r2,rinvsq2,rainv2,gsrainv2,gsrainvsq2,expgsrainv2;
F_FLOAT rinv12,cs,delcs,delcssq,facexp,facrad,frad1,frad2,pcsinv,pcsinvsq,pcs;
F_FLOAT facang,facang12,csfacang,csfac1,csfac2;
r1 = sqrt(rsq1);
rinvsq1 = 1.0/rsq1;
rainv1 = 1.0/(r1 - paramij.r0);
gsrainv1 = paramij.gamma * rainv1;
gsrainvsq1 = gsrainv1*rainv1/r1;
expgsrainv1 = exp(gsrainv1);
r2 = sqrt(rsq2);
rinvsq2 = 1.0/rsq2;
rainv2 = 1.0/(r2 - paramik.r0);
gsrainv2 = paramik.gamma * rainv2;
gsrainvsq2 = gsrainv2*rainv2/r2;
expgsrainv2 = exp(gsrainv2);
rinv12 = 1.0/(r1*r2);
cs = (delr1[0]*delr2[0] + delr1[1]*delr2[1] + delr1[2]*delr2[2]) * rinv12;
delcs = cs - paramijk.costheta;
delcssq = delcs*delcs;
pcsinv = paramijk.bigc*delcssq + 1.0;
pcsinvsq = pcsinv*pcsinv;
pcs = delcssq/pcsinv;
facexp = expgsrainv1*expgsrainv2;
facrad = paramijk.bigb * facexp * pcs;
frad1 = facrad*gsrainvsq1;
frad2 = facrad*gsrainvsq2;
facang = paramijk.big2b * facexp * delcs/pcsinvsq;
facang12 = rinv12*facang;
csfacang = cs*facang;
csfac1 = rinvsq1*csfacang;
fj[0] = delr1[0]*(frad1+csfac1)-delr2[0]*facang12;
fj[1] = delr1[1]*(frad1+csfac1)-delr2[1]*facang12;
fj[2] = delr1[2]*(frad1+csfac1)-delr2[2]*facang12;
}
/* ---------------------------------------------------------------------- */
template<class DeviceType>
template<int NEIGHFLAG>
KOKKOS_INLINE_FUNCTION
void PairVashishtaKokkos<DeviceType>::ev_tally(EV_FLOAT &ev, const int &i, const int &j,
const F_FLOAT &epair, const F_FLOAT &fpair, const F_FLOAT &delx,
const F_FLOAT &dely, const F_FLOAT &delz) const
{
const int VFLAG = vflag_either;
// The eatom and vatom arrays are atomic for half/thread neighbor list
Kokkos::View<E_FLOAT*, typename DAT::t_efloat_1d::array_layout,DeviceType,Kokkos::MemoryTraits<AtomicF<NEIGHFLAG>::value> > v_eatom = k_eatom.view<DeviceType>();
Kokkos::View<F_FLOAT*[6], typename DAT::t_virial_array::array_layout,DeviceType,Kokkos::MemoryTraits<AtomicF<NEIGHFLAG>::value> > v_vatom = k_vatom.view<DeviceType>();
if (eflag_atom) {
const E_FLOAT epairhalf = 0.5 * epair;
v_eatom[i] += epairhalf;
if (NEIGHFLAG != FULL)
v_eatom[j] += epairhalf;
}
if (VFLAG) {
const E_FLOAT v0 = delx*delx*fpair;
const E_FLOAT v1 = dely*dely*fpair;
const E_FLOAT v2 = delz*delz*fpair;
const E_FLOAT v3 = delx*dely*fpair;
const E_FLOAT v4 = delx*delz*fpair;
const E_FLOAT v5 = dely*delz*fpair;
if (vflag_global) {
if (NEIGHFLAG != FULL) {
ev.v[0] += v0;
ev.v[1] += v1;
ev.v[2] += v2;
ev.v[3] += v3;
ev.v[4] += v4;
ev.v[5] += v5;
} else {
ev.v[0] += 0.5*v0;
ev.v[1] += 0.5*v1;
ev.v[2] += 0.5*v2;
ev.v[3] += 0.5*v3;
ev.v[4] += 0.5*v4;
ev.v[5] += 0.5*v5;
}
}
if (vflag_atom) {
v_vatom(i,0) += 0.5*v0;
v_vatom(i,1) += 0.5*v1;
v_vatom(i,2) += 0.5*v2;
v_vatom(i,3) += 0.5*v3;
v_vatom(i,4) += 0.5*v4;
v_vatom(i,5) += 0.5*v5;
if (NEIGHFLAG != FULL) {
v_vatom(j,0) += 0.5*v0;
v_vatom(j,1) += 0.5*v1;
v_vatom(j,2) += 0.5*v2;
v_vatom(j,3) += 0.5*v3;
v_vatom(j,4) += 0.5*v4;
v_vatom(j,5) += 0.5*v5;
}
}
}
}
/* ----------------------------------------------------------------------
tally eng_vdwl and virial into global and per-atom accumulators
called by SW and hbond potentials, newton_pair is always on
virial = riFi + rjFj + rkFk = (rj-ri) Fj + (rk-ri) Fk = drji*fj + drki*fk
------------------------------------------------------------------------- */
template<class DeviceType>
template<int NEIGHFLAG>
KOKKOS_INLINE_FUNCTION
void PairVashishtaKokkos<DeviceType>::ev_tally3(EV_FLOAT &ev, const int &i, const int &j, int &k,
const F_FLOAT &evdwl, const F_FLOAT &ecoul,
F_FLOAT *fj, F_FLOAT *fk, F_FLOAT *drji, F_FLOAT *drki) const
{
F_FLOAT epairthird,v[6];
const int VFLAG = vflag_either;
// The eatom and vatom arrays are atomic for half/thread neighbor list
Kokkos::View<E_FLOAT*, typename DAT::t_efloat_1d::array_layout,DeviceType,Kokkos::MemoryTraits<AtomicF<NEIGHFLAG>::value> > v_eatom = k_eatom.view<DeviceType>();
Kokkos::View<F_FLOAT*[6], typename DAT::t_virial_array::array_layout,DeviceType,Kokkos::MemoryTraits<AtomicF<NEIGHFLAG>::value> > v_vatom = k_vatom.view<DeviceType>();
if (eflag_atom) {
epairthird = THIRD * (evdwl + ecoul);
v_eatom[i] += epairthird;
if (NEIGHFLAG != FULL) {
v_eatom[j] += epairthird;
v_eatom[k] += epairthird;
}
}
if (VFLAG) {
v[0] = drji[0]*fj[0] + drki[0]*fk[0];
v[1] = drji[1]*fj[1] + drki[1]*fk[1];
v[2] = drji[2]*fj[2] + drki[2]*fk[2];
v[3] = drji[0]*fj[1] + drki[0]*fk[1];
v[4] = drji[0]*fj[2] + drki[0]*fk[2];
v[5] = drji[1]*fj[2] + drki[1]*fk[2];
if (vflag_global) {
ev.v[0] += v[0];
ev.v[1] += v[1];
ev.v[2] += v[2];
ev.v[3] += v[3];
ev.v[4] += v[4];
ev.v[5] += v[5];
}
if (vflag_atom) {
v_vatom(i,0) += THIRD*v[0]; v_vatom(i,1) += THIRD*v[1];
v_vatom(i,2) += THIRD*v[2]; v_vatom(i,3) += THIRD*v[3];
v_vatom(i,4) += THIRD*v[4]; v_vatom(i,5) += THIRD*v[5];
if (NEIGHFLAG != FULL) {
v_vatom(j,0) += THIRD*v[0]; v_vatom(j,1) += THIRD*v[1];
v_vatom(j,2) += THIRD*v[2]; v_vatom(j,3) += THIRD*v[3];
v_vatom(j,4) += THIRD*v[4]; v_vatom(j,5) += THIRD*v[5];
v_vatom(k,0) += THIRD*v[0]; v_vatom(k,1) += THIRD*v[1];
v_vatom(k,2) += THIRD*v[2]; v_vatom(k,3) += THIRD*v[3];
v_vatom(k,4) += THIRD*v[4]; v_vatom(k,5) += THIRD*v[5];
}
}
}
}
/* ----------------------------------------------------------------------
tally eng_vdwl and virial into global and per-atom accumulators
called by SW and hbond potentials, newton_pair is always on
virial = riFi + rjFj + rkFk = (rj-ri) Fj + (rk-ri) Fk = drji*fj + drki*fk
------------------------------------------------------------------------- */
template<class DeviceType>
KOKKOS_INLINE_FUNCTION
void PairVashishtaKokkos<DeviceType>::ev_tally3_atom(EV_FLOAT &ev, const int &i,
const F_FLOAT &evdwl, const F_FLOAT &ecoul,
F_FLOAT *fj, F_FLOAT *fk, F_FLOAT *drji, F_FLOAT *drki) const
{
F_FLOAT epairthird,v[6];
const int VFLAG = vflag_either;
if (eflag_atom) {
epairthird = THIRD * (evdwl + ecoul);
d_eatom[i] += epairthird;
}
if (VFLAG) {
v[0] = drji[0]*fj[0] + drki[0]*fk[0];
v[1] = drji[1]*fj[1] + drki[1]*fk[1];
v[2] = drji[2]*fj[2] + drki[2]*fk[2];
v[3] = drji[0]*fj[1] + drki[0]*fk[1];
v[4] = drji[0]*fj[2] + drki[0]*fk[2];
v[5] = drji[1]*fj[2] + drki[1]*fk[2];
if (vflag_atom) {
d_vatom(i,0) += THIRD*v[0]; d_vatom(i,1) += THIRD*v[1];
d_vatom(i,2) += THIRD*v[2]; d_vatom(i,3) += THIRD*v[3];
d_vatom(i,4) += THIRD*v[4]; d_vatom(i,5) += THIRD*v[5];
}
}
}
namespace LAMMPS_NS {
template class PairVashishtaKokkos<LMPDeviceType>;
#ifdef KOKKOS_HAVE_CUDA
template class PairVashishtaKokkos<LMPHostType>;
#endif
}

View File

@ -0,0 +1,162 @@
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#ifdef PAIR_CLASS
PairStyle(vashishta/kk,PairVashishtaKokkos<LMPDeviceType>)
PairStyle(vashishta/kk/device,PairVashishtaKokkos<LMPDeviceType>)
PairStyle(vashishta/kk/host,PairVashishtaKokkos<LMPHostType>)
#else
#ifndef LMP_PAIR_VASHISHTA_KOKKOS_H
#define LMP_PAIR_VASHISHTA_KOKKOS_H
#include "pair_vashishta.h"
#include "pair_kokkos.h"
template<int NEIGHFLAG, int EVFLAG>
struct TagPairVashishtaComputeHalf{};
template<int NEIGHFLAG, int EVFLAG>
struct TagPairVashishtaComputeFullA{};
template<int NEIGHFLAG, int EVFLAG>
struct TagPairVashishtaComputeFullB{};
struct TagPairVashishtaComputeShortNeigh{};
namespace LAMMPS_NS {
template<class DeviceType>
class PairVashishtaKokkos : public PairVashishta {
public:
enum {EnabledNeighFlags=FULL};
enum {COUL_FLAG=0};
typedef DeviceType device_type;
typedef ArrayTypes<DeviceType> AT;
typedef EV_FLOAT value_type;
PairVashishtaKokkos(class LAMMPS *);
virtual ~PairVashishtaKokkos();
virtual void compute(int, int);
virtual void coeff(int, char **);
virtual void init_style();
template<int NEIGHFLAG, int EVFLAG>
KOKKOS_INLINE_FUNCTION
void operator()(TagPairVashishtaComputeHalf<NEIGHFLAG,EVFLAG>, const int&, EV_FLOAT&) const;
template<int NEIGHFLAG, int EVFLAG>
KOKKOS_INLINE_FUNCTION
void operator()(TagPairVashishtaComputeHalf<NEIGHFLAG,EVFLAG>, const int&) const;
template<int NEIGHFLAG, int EVFLAG>
KOKKOS_INLINE_FUNCTION
void operator()(TagPairVashishtaComputeFullA<NEIGHFLAG,EVFLAG>, const int&, EV_FLOAT&) const;
template<int NEIGHFLAG, int EVFLAG>
KOKKOS_INLINE_FUNCTION
void operator()(TagPairVashishtaComputeFullA<NEIGHFLAG,EVFLAG>, const int&) const;
template<int NEIGHFLAG, int EVFLAG>
KOKKOS_INLINE_FUNCTION
void operator()(TagPairVashishtaComputeFullB<NEIGHFLAG,EVFLAG>, const int&, EV_FLOAT&) const;
template<int NEIGHFLAG, int EVFLAG>
KOKKOS_INLINE_FUNCTION
void operator()(TagPairVashishtaComputeFullB<NEIGHFLAG,EVFLAG>, const int&) const;
KOKKOS_INLINE_FUNCTION
void operator()(TagPairVashishtaComputeShortNeigh, const int&) const;
template<int NEIGHFLAG>
KOKKOS_INLINE_FUNCTION
void ev_tally(EV_FLOAT &ev, const int &i, const int &j,
const F_FLOAT &epair, const F_FLOAT &fpair, const F_FLOAT &delx,
const F_FLOAT &dely, const F_FLOAT &delz) const;
template<int NEIGHFLAG>
KOKKOS_INLINE_FUNCTION
void ev_tally3(EV_FLOAT &ev, const int &i, const int &j, int &k,
const F_FLOAT &evdwl, const F_FLOAT &ecoul,
F_FLOAT *fj, F_FLOAT *fk, F_FLOAT *drji, F_FLOAT *drki) const;
KOKKOS_INLINE_FUNCTION
void ev_tally3_atom(EV_FLOAT &ev, const int &i,
const F_FLOAT &evdwl, const F_FLOAT &ecoul,
F_FLOAT *fj, F_FLOAT *fk, F_FLOAT *drji, F_FLOAT *drki) const;
protected:
typedef Kokkos::DualView<int***,DeviceType> tdual_int_3d;
typedef typename tdual_int_3d::t_dev_const_randomread t_int_3d_randomread;
typedef typename tdual_int_3d::t_host t_host_int_3d;
t_int_3d_randomread d_elem2param;
DAT::t_int_1d_randomread d_map;
typedef Kokkos::DualView<Param*,DeviceType> tdual_param_1d;
typedef typename tdual_param_1d::t_dev t_param_1d;
typedef typename tdual_param_1d::t_host t_host_param_1d;
t_param_1d d_params;
virtual void setup_params();
void twobody(const Param&, const F_FLOAT&, F_FLOAT&, const int&, F_FLOAT&) const;
void threebody(const Param&, const Param&, const Param&, const F_FLOAT&, const F_FLOAT&, F_FLOAT *, F_FLOAT *,
F_FLOAT *, F_FLOAT *, const int&, F_FLOAT&) const;
void threebodyj(const Param&, const Param&, const Param&, const F_FLOAT&, const F_FLOAT&, F_FLOAT *, F_FLOAT *,
F_FLOAT *) const;
typename ArrayTypes<DeviceType>::t_x_array_randomread x;
typename ArrayTypes<DeviceType>::t_f_array f;
typename ArrayTypes<DeviceType>::t_tagint_1d tag;
typename ArrayTypes<DeviceType>::t_int_1d_randomread type;
DAT::tdual_efloat_1d k_eatom;
DAT::tdual_virial_array k_vatom;
DAT::t_efloat_1d d_eatom;
DAT::t_virial_array d_vatom;
DAT::t_int_1d_randomread d_type2frho;
DAT::t_int_2d_randomread d_type2rhor;
DAT::t_int_2d_randomread d_type2z2r;
typename ArrayTypes<DeviceType>::t_neighbors_2d d_neighbors;
typename ArrayTypes<DeviceType>::t_int_1d_randomread d_ilist;
typename ArrayTypes<DeviceType>::t_int_1d_randomread d_numneigh;
//NeighListKokkos<DeviceType> k_list;
int neighflag,newton_pair;
int nlocal,nall,eflag,vflag;
int inum;
Kokkos::View<int**,DeviceType> d_neighbors_short_2body;
Kokkos::View<int*,DeviceType> d_numneigh_short_2body;
Kokkos::View<int**,DeviceType> d_neighbors_short_3body;
Kokkos::View<int*,DeviceType> d_numneigh_short_3body;
friend void pair_virial_fdotr_compute<PairVashishtaKokkos>(PairVashishtaKokkos*);
};
}
#endif
#endif
/* ERROR/WARNING messages:
E: Cannot use chosen neighbor list style with pair vashishta/kk
Self-explanatory.
*/

View File

@ -33,10 +33,10 @@ template<class DeviceType>
class RegBlockKokkos : public RegBlock {
friend class FixPour;
public:
typedef DeviceType device_type;
typedef ArrayTypes<DeviceType> AT;
public:
RegBlockKokkos(class LAMMPS *, int, char **);
~RegBlockKokkos();
void match_all_kokkos(int, DAT::t_int_1d);

View File

@ -7,7 +7,7 @@ SHELL = /bin/sh
# specify flags and libraries needed for your compiler
CC = mpicxx -cxx=icc
CCFLAGS = -g -O3
CCFLAGS = -g -O3 -restrict
SHFLAGS = -fPIC
DEPFLAGS = -M

View File

@ -7,7 +7,7 @@ SHELL = /bin/sh
# specify flags and libraries needed for your compiler
CC = icc
CCFLAGS = -g -O3
CCFLAGS = -g -O3 -restrict
SHFLAGS = -fPIC
DEPFLAGS = -M

View File

@ -8,7 +8,7 @@ SHELL = /bin/sh
export OMPI_CXX = icc
CC = mpicxx
CCFLAGS = -g -O3
CCFLAGS = -g -O3 -restrict
SHFLAGS = -fPIC
DEPFLAGS = -M

View File

@ -7,7 +7,7 @@ SHELL = /bin/sh
# specify flags and libraries needed for your compiler
CC = icc
CCFLAGS = -g -O3
CCFLAGS = -g -O3 -restrict
SHFLAGS = -fPIC
DEPFLAGS = -M

View File

@ -7,7 +7,7 @@ SHELL = /bin/sh
# specify flags and libraries needed for your compiler
CC = icc
CCFLAGS = -g -O3
CCFLAGS = -g -O3 -restrict
SHFLAGS = -fPIC
DEPFLAGS = -M

View File

@ -49,6 +49,7 @@ class PairTersoff : public Pair {
double ZBLcut,ZBLexpscale;
double c5,ca1,ca4; // added for TersoffMOD
double powern_del;
double c0; // added for TersoffMODC
};
Param *params; // parameter set for an I-J-K interaction

View File

@ -30,7 +30,7 @@ class PairTersoffMOD : public PairTersoff {
~PairTersoffMOD() {}
protected:
void read_file(char *);
virtual void read_file(char *);
virtual void setup_params();
double zeta(Param *, double, double, double *, double *);

View File

@ -0,0 +1,196 @@
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
Contributing author: Ganga P Purja Pun (George Mason University, Fairfax)
------------------------------------------------------------------------- */
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include "pair_tersoff_mod_c.h"
#include "atom.h"
#include "neighbor.h"
#include "neigh_list.h"
#include "neigh_request.h"
#include "force.h"
#include "comm.h"
#include "memory.h"
#include "error.h"
#include "math_const.h"
using namespace LAMMPS_NS;
using namespace MathConst;
#define MAXLINE 1024
#define DELTA 4
/* ---------------------------------------------------------------------- */
void PairTersoffMODC::read_file(char *file)
{
int params_per_line = 21;
char **words = new char*[params_per_line+1];
memory->sfree(params);
params = NULL;
nparams = maxparam = 0;
// open file on proc 0
FILE *fp;
if (comm->me == 0) {
fp = force->open_potential(file);
if (fp == NULL) {
char str[128];
sprintf(str,"Cannot open Tersoff potential file %s",file);
error->one(FLERR,str);
}
}
// read each line out of file, skipping blank lines or leading '#'
// store line of params if all 3 element tags are in element list
int n,nwords,ielement,jelement,kelement;
char line[MAXLINE],*ptr;
int eof = 0;
while (1) {
if (comm->me == 0) {
ptr = fgets(line,MAXLINE,fp);
if (ptr == NULL) {
eof = 1;
fclose(fp);
} else n = strlen(line) + 1;
}
MPI_Bcast(&eof,1,MPI_INT,0,world);
if (eof) break;
MPI_Bcast(&n,1,MPI_INT,0,world);
MPI_Bcast(line,n,MPI_CHAR,0,world);
// strip comment, skip line if blank
if ((ptr = strchr(line,'#'))) *ptr = '\0';
nwords = atom->count_words(line);
if (nwords == 0) continue;
// concatenate additional lines until have params_per_line words
while (nwords < params_per_line) {
n = strlen(line);
if (comm->me == 0) {
ptr = fgets(&line[n],MAXLINE-n,fp);
if (ptr == NULL) {
eof = 1;
fclose(fp);
} else n = strlen(line) + 1;
}
MPI_Bcast(&eof,1,MPI_INT,0,world);
if (eof) break;
MPI_Bcast(&n,1,MPI_INT,0,world);
MPI_Bcast(line,n,MPI_CHAR,0,world);
if ((ptr = strchr(line,'#'))) *ptr = '\0';
nwords = atom->count_words(line);
}
if (nwords != params_per_line)
error->all(FLERR,"Incorrect format in Tersoff potential file");
// words = ptrs to all words in line
nwords = 0;
words[nwords++] = strtok(line," \t\n\r\f");
while ((words[nwords++] = strtok(NULL," \t\n\r\f"))) continue;
// ielement,jelement,kelement = 1st args
// if all 3 args are in element list, then parse this line
// else skip to next line
for (ielement = 0; ielement < nelements; ielement++)
if (strcmp(words[0],elements[ielement]) == 0) break;
if (ielement == nelements) continue;
for (jelement = 0; jelement < nelements; jelement++)
if (strcmp(words[1],elements[jelement]) == 0) break;
if (jelement == nelements) continue;
for (kelement = 0; kelement < nelements; kelement++)
if (strcmp(words[2],elements[kelement]) == 0) break;
if (kelement == nelements) continue;
// load up parameter settings and error check their values
if (nparams == maxparam) {
maxparam += DELTA;
params = (Param *) memory->srealloc(params,maxparam*sizeof(Param),
"pair:params");
}
params[nparams].ielement = ielement;
params[nparams].jelement = jelement;
params[nparams].kelement = kelement;
params[nparams].powerm = atof(words[3]);
params[nparams].lam3 = atof(words[4]);
params[nparams].h = atof(words[5]);
params[nparams].powern = atof(words[6]);
params[nparams].beta = atof(words[7]);
params[nparams].lam2 = atof(words[8]);
params[nparams].bigb = atof(words[9]);
params[nparams].bigr = atof(words[10]);
params[nparams].bigd = atof(words[11]);
params[nparams].lam1 = atof(words[12]);
params[nparams].biga = atof(words[13]);
params[nparams].powern_del = atof(words[14]);
params[nparams].c1 = atof(words[15]);
params[nparams].c2 = atof(words[16]);
params[nparams].c3 = atof(words[17]);
params[nparams].c4 = atof(words[18]);
params[nparams].c5 = atof(words[19]);
params[nparams].c0 = atof(words[20]);
// currently only allow m exponent of 1
params[nparams].powermint = int(params[nparams].powerm);
if (
params[nparams].lam3 < 0.0 || params[nparams].powern < 0.0 ||
params[nparams].beta < 0.0 || params[nparams].lam2 < 0.0 ||
params[nparams].bigb < 0.0 || params[nparams].bigr < 0.0 ||
params[nparams].bigd < 0.0 ||
params[nparams].bigd > params[nparams].bigr ||
params[nparams].lam3 < 0.0 || params[nparams].biga < 0.0 ||
params[nparams].powerm - params[nparams].powermint != 0.0 ||
(params[nparams].powermint != 3 && params[nparams].powermint != 1))
error->all(FLERR,"Illegal Tersoff parameter");
nparams++;
}
delete [] words;
}
/* ---------------------------------------------------------------------- */
void PairTersoffMODC::repulsive(Param *param, double rsq, double &fforce,
int eflag, double &eng)
{
double r,tmp_fc,tmp_fc_d,tmp_exp;
r = sqrt(rsq);
tmp_fc = ters_fc(r,param);
tmp_fc_d = ters_fc_d(r,param);
tmp_exp = exp(-param->lam1 * r);
fforce = -param->biga * tmp_exp * (tmp_fc_d - tmp_fc * param->lam1) / r - param->c0 * tmp_fc_d / r;
if (eflag) eng = tmp_fc * param->biga * tmp_exp + param->c0 * tmp_fc;
}

View File

@ -0,0 +1,66 @@
/* -*- c++ -*- ----------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#ifdef PAIR_CLASS
PairStyle(tersoff/mod/c,PairTersoffMODC)
#else
#ifndef LMP_PAIR_TERSOFF_MOD_C_H
#define LMP_PAIR_TERSOFF_MOD_C_H
#include "pair_tersoff_mod.h"
namespace LAMMPS_NS {
class PairTersoffMODC : public PairTersoffMOD {
public:
PairTersoffMODC(class LAMMPS *lmp) : PairTersoffMOD(lmp) {};
~PairTersoffMODC() {}
protected:
void read_file(char *);
void repulsive(Param *, double, double &, int, double &);
};
}
#endif
#endif
/* ERROR/WARNING messages:
E: Cannot open Tersoff potential file %s
The specified potential file cannot be opened. Check that the path
and name are correct.
E: Incorrect format in Tersoff potential file
Incorrect number of words per line in the potential file.
E: Illegal Tersoff parameter
One or more of the coefficients defined in the potential file is
invalid.
E: Potential file has duplicate entry
The potential file has more than one entry for the same element.
E: Potential file is missing an entry
The potential file does not have a needed entry.
*/

View File

@ -64,6 +64,8 @@ PairVashishta::PairVashishta(LAMMPS *lmp) : Pair(lmp)
PairVashishta::~PairVashishta()
{
if (copymode) return;
if (elements)
for (int i = 0; i < nelements; i++) delete [] elements[i];
delete [] elements;

View File

@ -34,7 +34,6 @@ class PairVashishta : public Pair {
double init_one(int, int);
void init_style();
protected:
struct Param {
double bigb,gamma,r0,bigc,costheta;
double bigh,eta,zi,zj;
@ -45,7 +44,7 @@ class PairVashishta : public Pair {
double lam1rc,lam4rc,vrcc2,vrcc3,vrc,dvrc,c0;
int ielement,jelement,kelement;
};
protected:
double cutmax; // max cutoff for all elements
int nelements; // # of unique elements
char **elements; // names of unique elements

View File

@ -55,7 +55,8 @@ static const char cite_fix_orient_fcc[] =
FixOrientFCC::FixOrientFCC(LAMMPS *lmp, int narg, char **arg) :
Fix(lmp, narg, arg),
xifilename(NULL), chifilename(NULL), order(NULL), nbr(NULL), sort(NULL), list(NULL)
xifilename(NULL), chifilename(NULL), order(NULL), nbr(NULL),
sort(NULL), list(NULL)
{
if (lmp->citeme) lmp->citeme->add(cite_fix_orient_fcc);

View File

@ -62,7 +62,12 @@ static const char cite_fix_poems[] =
------------------------------------------------------------------------- */
FixPOEMS::FixPOEMS(LAMMPS *lmp, int narg, char **arg) :
Fix(lmp, narg, arg)
Fix(lmp, narg, arg), step_respa(NULL), natom2body(NULL),
atom2body(NULL), displace(NULL), nrigid(NULL), masstotal(NULL),
xcm(NULL), vcm(NULL), fcm(NULL), inertia(NULL), ex_space(NULL),
ey_space(NULL), ez_space(NULL), angmom(NULL), omega(NULL),
torque(NULL), sum(NULL), all(NULL), jointbody(NULL),
xjoint(NULL), freelist(NULL), poems(NULL)
{
if (lmp->citeme) lmp->citeme->add(cite_fix_poems);

View File

@ -28,12 +28,10 @@
#ifndef LMP_INTEL_PREPROCESS_H
#define LMP_INTEL_PREPROCESS_H
#ifndef LAMMPS_MEMALIGN
#error Please set -DLAMMPS_MEMALIGN=64 in CCFLAGS for your LAMMPS makefile.
#else
#if (LAMMPS_MEMALIGN != 64)
#error Please set -DLAMMPS_MEMALIGN=64 in CCFLAGS for your LAMMPS makefile.
#endif
// LAMMPS_MEMALIGN is set to 64 by default for -DLMP_USER_INTEL
// so we only need to error out in case of a different alignment
#if LAMMPS_MEMALIGN && (LAMMPS_MEMALIGN != 64)
#error Please set -DLAMMPS_MEMALIGN=64 in CCFLAGS of your LAMMPS makefile for USER-INTEL package
#endif
#if defined(_OPENMP)

View File

@ -49,6 +49,7 @@ improper_style cossq, Georgios Vogiatzis, gvog at chemeng.ntua.gr, 25 May 12
improper_style fourier, Loukas Peristeras, loukas.peristeras at scienomics.com, 27 Oct 12
improper_style ring, Georgios Vogiatzis, gvog at chemeng.ntua.gr, 25 May 12
improper_style distance, Paolo Raiteri, p.raiteri at curtin.edu.au, 2 Dec 15
pair_style agni, Axel Kohlmeyer, akohlmey at gmail.com, 9 Nov 16
pair_style buck/mdf, Paolo Raiteri, p.raiteri at curtin.edu.au, 2 Dec 15
pair_style coul/diel, Axel Kohlmeyer, akohlmey at gmail.com, 1 Dec 11
pair_style dipole/sf, Mario Orsi, orsimario at gmail.com, 8 Aug 11

508
src/USER-MISC/pair_agni.cpp Normal file
View File

@ -0,0 +1,508 @@
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
Contributing authors: Axel Kohlmeyer (Temple U), Venkatesh Botu
------------------------------------------------------------------------- */
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include "pair_agni.h"
#include "atom.h"
#include "neighbor.h"
#include "neigh_request.h"
#include "force.h"
#include "comm.h"
#include "memory.h"
#include "neighbor.h"
#include "neigh_list.h"
#include "memory.h"
#include "error.h"
#include "citeme.h"
#include "math_special.h"
#include "math_const.h"
using namespace LAMMPS_NS;
using namespace MathSpecial;
static const char cite_pair_agni[] =
"pair agni command:\n\n"
"@article{botu2015adaptive,\n"
" author = {Botu, Venkatesh and Ramprasad, Rampi},\n"
" title = {Adaptive machine learning framework to"
" accelerate ab initio molecular dynamics},\n"
" journal = {International Journal of Quantum Chemistry},\n"
" volume = {115},\n"
" number = {16},\n"
" pages = {1074--1083},\n"
" year = {2015},\n"
" publisher = {Wiley Online Library}\n"
"}\n\n"
"@article{botu2015learning,\n"
" author = {Botu, Venkatesh and Ramprasad, Rampi},\n"
" title = {Learning scheme to predict atomic forces"
" and accelerate materials simulations},\n"
" journal = {Physical Review B},\n"
" volume = {92},\n"
" number = {9},\n"
" pages = {094306},\n"
" year = {2015},\n"
" publisher = {APS}\n"
"}\n\n";
#define AGNI_VERSION 1
#define MAXLINE 10240
#define MAXWORD 40
/* ---------------------------------------------------------------------- */
PairAGNI::PairAGNI(LAMMPS *lmp) : Pair(lmp)
{
if (lmp->citeme) lmp->citeme->add(cite_pair_agni);
single_enable = 0;
restartinfo = 0;
one_coeff = 1;
manybody_flag = 1;
no_virial_fdotr_compute = 1;
nelements = 0;
elements = NULL;
elem2param = NULL;
nparams = 0;
params = NULL;
map = NULL;
}
/* ----------------------------------------------------------------------
check if allocated, since class can be destructed when incomplete
------------------------------------------------------------------------- */
PairAGNI::~PairAGNI()
{
if (elements)
for (int i = 0; i < nelements; i++) delete [] elements[i];
delete [] elements;
if (params) {
for (int i = 0; i < nparams; ++i) {
int n = params[i].numeta;
for (int j = 0; j < n; ++j) {
delete [] params[i].xU[j];
}
delete [] params[i].eta;
delete [] params[i].alpha;
delete [] params[i].xU;
delete [] params[i].yU;
}
memory->destroy(params);
params = NULL;
}
if (allocated) {
memory->destroy(setflag);
memory->destroy(cutsq);
delete [] map;
}
}
/* ---------------------------------------------------------------------- */
void PairAGNI::compute(int eflag, int vflag)
{
int i,j,k,ii,jj,inum,jnum,itype;
double xtmp,ytmp,ztmp,delx,dely,delz;
double rsq;
int *ilist,*jlist,*numneigh,**firstneigh;
if (eflag || vflag) ev_setup(eflag,vflag);
else evflag = vflag_fdotr = 0;
double **x = atom->x;
double **f = atom->f;
int *type = atom->type;
inum = list->inum;
ilist = list->ilist;
numneigh = list->numneigh;
firstneigh = list->firstneigh;
double fxtmp,fytmp,fztmp;
double *Vx, *Vy, *Vz;
// loop over full neighbor list of my atoms
for (ii = 0; ii < inum; ii++) {
i = ilist[ii];
itype = map[type[i]];
xtmp = x[i][0];
ytmp = x[i][1];
ztmp = x[i][2];
fxtmp = fytmp = fztmp = 0.0;
const Param &iparam = params[elem2param[itype]];
Vx = new double[iparam.numeta];
Vy = new double[iparam.numeta];
Vz = new double[iparam.numeta];
memset(Vx,0,iparam.numeta*sizeof(double));
memset(Vy,0,iparam.numeta*sizeof(double));
memset(Vz,0,iparam.numeta*sizeof(double));
jlist = firstneigh[i];
jnum = numneigh[i];
for (jj = 0; jj < jnum; jj++) {
j = jlist[jj];
j &= NEIGHMASK;
delx = xtmp - x[j][0];
dely = ytmp - x[j][1];
delz = ztmp - x[j][2];
rsq = delx*delx + dely*dely + delz*delz;
if ((rsq > 0.0) && (rsq < iparam.cutsq)) {
const double r = sqrt(rsq);
const double cF = 0.5*(cos((MathConst::MY_PI*r)/iparam.cut)+1.0);
const double wX = cF*delx/r;
const double wY = cF*dely/r;
const double wZ = cF*delz/r;
for (k = 0; k < iparam.numeta; ++k) {
const double e = exp(-(iparam.eta[k]*rsq));
Vx[k] += wX*e;
Vy[k] += wY*e;
Vz[k] += wZ*e;
}
}
}
for (j = 0; j < iparam.numtrain; ++j) {
double kx = 0.0;
double ky = 0.0;
double kz = 0.0;
for(int k = 0; k < iparam.numeta; ++k) {
const double xu = iparam.xU[k][j];
kx += square(Vx[k] - xu);
ky += square(Vy[k] - xu);
kz += square(Vz[k] - xu);
}
const double e = -0.5/(square(iparam.sigma));
fxtmp += iparam.alpha[j]*exp(kx*e);
fytmp += iparam.alpha[j]*exp(ky*e);
fztmp += iparam.alpha[j]*exp(kz*e);
}
fxtmp += iparam.b;
fytmp += iparam.b;
fztmp += iparam.b;
f[i][0] += fxtmp;
f[i][1] += fytmp;
f[i][2] += fztmp;
if (evflag) ev_tally_xyz_full(i,0.0,0.0,fxtmp,fytmp,fztmp,delx,dely,delz);
delete [] Vx;
delete [] Vy;
delete [] Vz;
}
if (vflag_fdotr) virial_fdotr_compute();
}
/* ---------------------------------------------------------------------- */
void PairAGNI::allocate()
{
allocated = 1;
int n = atom->ntypes;
memory->create(setflag,n+1,n+1,"pair:setflag");
memory->create(cutsq,n+1,n+1,"pair:cutsq");
map = new int[n+1];
}
/* ----------------------------------------------------------------------
global settings
------------------------------------------------------------------------- */
void PairAGNI::settings(int narg, char **arg)
{
if (narg != 0) error->all(FLERR,"Illegal pair_style command");
}
/* ----------------------------------------------------------------------
set coeffs for one or more type pairs
------------------------------------------------------------------------- */
void PairAGNI::coeff(int narg, char **arg)
{
int i,j,n;
if (!allocated) allocate();
if (narg != 3 + atom->ntypes)
error->all(FLERR,"Incorrect args for pair coefficients");
// insure I,J args are * *
if (strcmp(arg[0],"*") != 0 || strcmp(arg[1],"*") != 0)
error->all(FLERR,"Incorrect args for pair coefficients");
// read args that map atom types to elements in potential file
// map[i] = which element the Ith atom type is, -1 if NULL
// nelements = # of unique elements
// elements = list of element names
if (elements) {
for (i = 0; i < nelements; i++) delete [] elements[i];
delete [] elements;
}
elements = new char*[atom->ntypes];
for (i = 0; i < atom->ntypes; i++) elements[i] = NULL;
nelements = 0;
for (i = 3; i < narg; i++) {
if (strcmp(arg[i],"NULL") == 0) {
map[i-2] = -1;
continue;
}
for (j = 0; j < nelements; j++)
if (strcmp(arg[i],elements[j]) == 0) break;
map[i-2] = j;
if (j == nelements) {
n = strlen(arg[i]) + 1;
elements[j] = new char[n];
strcpy(elements[j],arg[i]);
nelements++;
}
}
// read potential file and initialize potential parameters
read_file(arg[2]);
setup_params();
// clear setflag since coeff() called once with I,J = * *
n = atom->ntypes;
for (int i = 1; i <= n; i++)
for (int j = i; j <= n; j++)
setflag[i][j] = 0;
// set setflag i,j for type pairs where both are mapped to elements
int count = 0;
for (int i = 1; i <= n; i++)
for (int j = i; j <= n; j++)
if (map[i] >= 0 && map[j] >= 0) {
setflag[i][j] = 1;
count++;
}
if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients");
}
/* ----------------------------------------------------------------------
init specific to this pair style
------------------------------------------------------------------------- */
void PairAGNI::init_style()
{
// need a full neighbor list
int irequest = neighbor->request(this,instance_me);
neighbor->requests[irequest]->half = 0;
neighbor->requests[irequest]->full = 1;
}
/* ----------------------------------------------------------------------
init for one type pair i,j and corresponding j,i
------------------------------------------------------------------------- */
double PairAGNI::init_one(int i, int j)
{
if (setflag[i][j] == 0) error->all(FLERR,"All pair coeffs are not set");
return cutmax;
}
/* ---------------------------------------------------------------------- */
void PairAGNI::read_file(char *file)
{
memory->sfree(params);
params = NULL;
nparams = 0;
// open file on proc 0 only
// then read line by line and broadcast the line to all MPI ranks
FILE *fp;
if (comm->me == 0) {
fp = force->open_potential(file);
if (fp == NULL) {
char str[128];
sprintf(str,"Cannot open AGNI potential file %s",file);
error->one(FLERR,str);
}
}
int i,j,n,nwords,curparam,wantdata;
char line[MAXLINE],*ptr;
int eof = 0;
char **words = new char*[MAXWORD+1];
while (1) {
n = 0;
if (comm->me == 0) {
ptr = fgets(line,MAXLINE,fp);
if (ptr == NULL) {
eof = 1;
fclose(fp);
} else n = strlen(line) + 1;
}
MPI_Bcast(&eof,1,MPI_INT,0,world);
if (eof) break;
MPI_Bcast(&n,1,MPI_INT,0,world);
MPI_Bcast(line,n,MPI_CHAR,0,world);
// strip comment, skip line if blank
if ((ptr = strchr(line,'#'))) *ptr = '\0';
nwords = atom->count_words(line);
if (nwords == 0) continue;
if (nwords > MAXWORD)
error->all(FLERR,"Increase MAXWORD and recompile");
// words = ptrs to all words in line
nwords = 0;
words[nwords++] = strtok(line," \t\n\r\f");
while ((words[nwords++] = strtok(NULL," \t\n\r\f"))) continue;
--nwords;
if ((nwords == 2) && (strcmp(words[0],"generation") == 0)) {
int ver = atoi(words[1]);
if (ver != AGNI_VERSION)
error->all(FLERR,"Incompatible AGNI potential file version");
if ((ver == 1) && (nelements != 1))
error->all(FLERR,"Cannot handle multi-element systems with this potential");
} else if ((nwords == 2) && (strcmp(words[0],"n_elements") == 0)) {
nparams = atoi(words[1]);
if ((nparams < 1) || params) // sanity check
error->all(FLERR,"Invalid AGNI potential file");
params = memory->create(params,nparams,"pair:params");
memset(params,0,nparams*sizeof(Param));
curparam = -1;
} else if (params && (nwords == nparams+1) && (strcmp(words[0],"element") == 0)) {
wantdata = -1;
for (i = 0; i < nparams; ++i) {
for (j = 0; j < nelements; ++j)
if (strcmp(words[i+1],elements[j]) == 0) break;
if (j == nelements)
error->all(FLERR,"No suitable parameters for requested element found");
else params[i].ielement = j;
}
} else if (params && (nwords == 2) && (strcmp(words[0],"interaction") == 0)) {
for (i = 0; i < nparams; ++i)
if (strcmp(words[1],elements[params[i].ielement]) == 0) curparam = i;
} else if ((curparam >=0) && (nwords == 1) && (strcmp(words[0],"endVar") == 0)) {
int numtrain = params[curparam].numtrain;
int numeta = params[curparam].numeta;
params[curparam].alpha = new double[numtrain];
params[curparam].yU = new double[numtrain];
params[curparam].xU = new double*[numeta];
for (i = 0; i < numeta; ++i)
params[curparam].xU[i] = new double[numtrain];
wantdata = curparam;
curparam = -1;
} else if ((curparam >=0) && (nwords == 2) && (strcmp(words[0],"Rc") == 0)) {
params[curparam].cut = atof(words[1]);
} else if ((curparam >=0) && (nwords == 2) && (strcmp(words[0],"Rs") == 0)) {
; // ignored
} else if ((curparam >=0) && (nwords == 2) && (strcmp(words[0],"neighbors") == 0)) {
; // ignored
} else if ((curparam >=0) && (nwords == 2) && (strcmp(words[0],"sigma") == 0)) {
params[curparam].sigma = atof(words[1]);
} else if ((curparam >=0) && (nwords == 2) && (strcmp(words[0],"lambda") == 0)) {
params[curparam].lambda = atof(words[1]);
} else if ((curparam >=0) && (nwords == 2) && (strcmp(words[0],"b") == 0)) {
params[curparam].b = atof(words[1]);
} else if ((curparam >=0) && (nwords == 2) && (strcmp(words[0],"n_train") == 0)) {
params[curparam].numtrain = atoi(words[1]);
} else if ((curparam >=0) && (nwords > 1) && (strcmp(words[0],"eta") == 0)) {
params[curparam].numeta = nwords-1;
params[curparam].eta = new double[nwords-1];
for (i = 0, j = 1 ; j < nwords; ++i, ++j)
params[curparam].eta[i] = atof(words[j]);
} else if (params && (wantdata >=0) && (nwords == params[wantdata].numeta+3)) {
n = (int) atof(words[0]);
for (i = 0; i < params[wantdata].numeta; ++i) {
params[wantdata].xU[i][n] = atof(words[i+1]);
}
params[wantdata].yU[n] = atof(words[params[wantdata].numeta+1]);
params[wantdata].alpha[n] = atof(words[params[wantdata].numeta+2]);
} else {
if (comm->me == 0)
error->warning(FLERR,"Ignoring unknown content in AGNI potential file.");
}
}
delete [] words;
}
/* ---------------------------------------------------------------------- */
void PairAGNI::setup_params()
{
int i,m,n;
double rtmp;
// set elem2param for all elements
memory->destroy(elem2param);
memory->create(elem2param,nelements,"pair:elem2param");
for (i = 0; i < nelements; i++) {
n = -1;
for (m = 0; m < nparams; m++) {
if (i == params[m].ielement) {
if (n >= 0) error->all(FLERR,"Potential file has duplicate entry");
n = m;
}
}
if (n < 0) error->all(FLERR,"Potential file is missing an entry");
elem2param[i] = n;
}
// compute parameter values derived from inputs
// set cutsq using shortcut to reduce neighbor list for accelerated
// calculations. cut must remain unchanged as it is a potential parameter
// (cut = a*sigma)
cutmax = 0.0;
for (m = 0; m < nparams; m++) {
rtmp = params[m].cut;
params[m].cutsq = rtmp * rtmp;
if (rtmp > cutmax) cutmax = rtmp;
}
}

85
src/USER-MISC/pair_agni.h Normal file
View File

@ -0,0 +1,85 @@
/* -*- c++ -*- ----------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#ifdef PAIR_CLASS
PairStyle(agni,PairAGNI)
#else
#ifndef LMP_PAIR_AGNI_H
#define LMP_PAIR_AGNI_H
#include "pair.h"
namespace LAMMPS_NS {
class PairAGNI : public Pair {
public:
PairAGNI(class LAMMPS *);
virtual ~PairAGNI();
virtual void compute(int, int);
void settings(int, char **);
virtual void coeff(int, char **);
virtual double init_one(int, int);
virtual void init_style();
struct Param {
double cut,cutsq;
double *eta,**xU,*yU,*alpha;
double sigma,lambda,b;
int numeta,numtrain,ielement;
};
protected:
double cutmax; // max cutoff for all elements
int nelements; // # of unique atom type labels
char **elements; // names of unique elements
int *elem2param; // mapping from element pairs to parameters
int *map; // mapping from atom types to elements
int nparams; // # of stored parameter sets
Param *params; // parameter set for an I-J interaction
virtual void allocate();
void read_file(char *);
virtual void setup_params();
};
}
#endif
#endif
/* ERROR/WARNING messages:
E: Illegal ... command
Self-explanatory. Check the input script syntax and compare to the
documentation for the command. You can use -echo screen as a
command-line option when running LAMMPS to see the offending line.
E: Incorrect args for pair coefficients
Self-explanatory. Check the input script or data file.
E: Cannot open AGNI potential file %s
The specified AGNI potential file cannot be opened. Check that the path
and name are correct.
E: Incorrect format in AGNI potential file
The potential file is not compatible with the AGNI pair style
implementation in this LAMMPS version.
*/

View File

@ -0,0 +1,299 @@
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
This software is distributed under the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
Contributing author: Axel Kohlmeyer (Temple U)
------------------------------------------------------------------------- */
#include <math.h>
#include <string.h>
#include <stdint.h>
#include "pair_agni_omp.h"
#include "atom.h"
#include "comm.h"
#include "force.h"
#include "memory.h"
#include "neighbor.h"
#include "neigh_list.h"
#include "math_special.h"
#include "math_const.h"
#include "suffix.h"
using namespace LAMMPS_NS;
using namespace MathSpecial;
/*
Copyright (c) 2012,2013 Axel Kohlmeyer <akohlmey@gmail.com>
All rights reserved.
Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions
are met:
* Redistributions of source code must retain the above copyright
notice, this list of conditions and the following disclaimer.
* Redistributions in binary form must reproduce the above copyright
notice, this list of conditions and the following disclaimer in the
documentation and/or other materials provided with the distribution.
* Neither the name of the <organization> nor the
names of its contributors may be used to endorse or promote products
derived from this software without specific prior written permission.
THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
ARE DISCLAIMED. IN NO EVENT SHALL <COPYRIGHT HOLDER> BE LIABLE FOR ANY
DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
(INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
*/
/* faster versions of 2**x, e**x, and 10**x in single and double precision.
*
* Based on the Cephes math library 2.8
*/
/* internal definitions for the fastermath library */
/* IEEE 754 double precision floating point data manipulation */
typedef union
{
double f;
uint64_t u;
struct {int32_t i0,i1;};
} udi_t;
#define FM_DOUBLE_BIAS 1023
#define FM_DOUBLE_EMASK 2146435072
#define FM_DOUBLE_MBITS 20
#define FM_DOUBLE_MMASK 1048575
#define FM_DOUBLE_EZERO 1072693248
/* generate 2**num in floating point by bitshifting */
#define FM_DOUBLE_INIT_EXP(var,num) \
var.i0 = 0; \
var.i1 = (((int) num) + FM_DOUBLE_BIAS) << 20
/* double precision constants */
#define FM_DOUBLE_LOG2OFE 1.4426950408889634074
#define FM_DOUBLE_LOGEOF2 6.9314718055994530942e-1
#define FM_DOUBLE_LOG2OF10 3.32192809488736234789
#define FM_DOUBLE_LOG10OF2 3.0102999566398119521e-1
#define FM_DOUBLE_LOG10OFE 4.3429448190325182765e-1
#define FM_DOUBLE_SQRT2 1.41421356237309504880
#define FM_DOUBLE_SQRTH 0.70710678118654752440
/* optimizer friendly implementation of exp2(x).
*
* strategy:
*
* split argument into an integer part and a fraction:
* ipart = floor(x+0.5);
* fpart = x - ipart;
*
* compute exp2(ipart) from setting the ieee754 exponent
* compute exp2(fpart) using a pade' approximation for x in [-0.5;0.5[
*
* the result becomes: exp2(x) = exp2(ipart) * exp2(fpart)
*/
static const double fm_exp2_q[] = {
/* 1.00000000000000000000e0, */
2.33184211722314911771e2,
4.36821166879210612817e3
};
static const double fm_exp2_p[] = {
2.30933477057345225087e-2,
2.02020656693165307700e1,
1.51390680115615096133e3
};
static double fm_exp2(double x)
{
double ipart, fpart, px, qx;
udi_t epart;
ipart = floor(x+0.5);
fpart = x - ipart;
FM_DOUBLE_INIT_EXP(epart,ipart);
x = fpart*fpart;
px = fm_exp2_p[0];
px = px*x + fm_exp2_p[1];
qx = x + fm_exp2_q[0];
px = px*x + fm_exp2_p[2];
qx = qx*x + fm_exp2_q[1];
px = px * fpart;
x = 1.0 + 2.0*(px/(qx-px));
return epart.f*x;
}
static double fm_exp(double x)
{
#if defined(__BYTE_ORDER__)
#if __BYTE_ORDER__ == __ORDER_LITTLE_ENDIAN__
return fm_exp2(FM_DOUBLE_LOG2OFE * (x));
#endif
#endif
return exp(x);
}
/* ---------------------------------------------------------------------- */
PairAGNIOMP::PairAGNIOMP(LAMMPS *lmp) :
PairAGNI(lmp), ThrOMP(lmp, THR_PAIR)
{
suffix_flag |= Suffix::OMP;
respa_enable = 0;
}
/* ---------------------------------------------------------------------- */
void PairAGNIOMP::compute(int eflag, int vflag)
{
if (eflag || vflag) {
ev_setup(eflag,vflag);
} else evflag = vflag_fdotr = 0;
const int nall = atom->nlocal + atom->nghost;
const int nthreads = comm->nthreads;
const int inum = list->inum;
#if defined(_OPENMP)
#pragma omp parallel default(none) shared(eflag,vflag)
#endif
{
int ifrom, ito, tid;
loop_setup_thr(ifrom, ito, tid, inum, nthreads);
ThrData *thr = fix->get_thr(tid);
thr->timer(Timer::START);
ev_setup_thr(eflag, vflag, nall, eatom, vatom, thr);
if (evflag) eval<1>(ifrom, ito, thr);
else eval<0>(ifrom, ito, thr);
thr->timer(Timer::PAIR);
reduce_thr(this, eflag, vflag, thr);
} // end of omp parallel region
}
template <int EVFLAG>
void PairAGNIOMP::eval(int iifrom, int iito, ThrData * const thr)
{
int i,j,k,ii,jj,itype,jnum;
double xtmp,ytmp,ztmp,delx,dely,delz;
double rsq;
int *ilist,*jlist,*numneigh,**firstneigh;
const dbl3_t * _noalias const x = (dbl3_t *) atom->x[0];
dbl3_t * _noalias const f = (dbl3_t *) thr->get_f()[0];
const int * _noalias const type = atom->type;
ilist = list->ilist;
numneigh = list->numneigh;
firstneigh = list->firstneigh;
double fxtmp,fytmp,fztmp;
double *Vx, *Vy, *Vz;
// loop over full neighbor list of my atoms
for (ii = iifrom; ii < iito; ++ii) {
i = ilist[ii];
itype = map[type[i]];
xtmp = x[i].x;
ytmp = x[i].y;
ztmp = x[i].z;
fxtmp = fytmp = fztmp = 0.0;
const Param &iparam = params[elem2param[itype]];
Vx = new double[iparam.numeta];
Vy = new double[iparam.numeta];
Vz = new double[iparam.numeta];
memset(Vx,0,iparam.numeta*sizeof(double));
memset(Vy,0,iparam.numeta*sizeof(double));
memset(Vz,0,iparam.numeta*sizeof(double));
jlist = firstneigh[i];
jnum = numneigh[i];
for (jj = 0; jj < jnum; jj++) {
j = jlist[jj];
j &= NEIGHMASK;
delx = xtmp - x[j].x;
dely = ytmp - x[j].y;
delz = ztmp - x[j].z;
rsq = delx*delx + dely*dely + delz*delz;
if ((rsq > 0.0) && (rsq < iparam.cutsq)) {
const double r = sqrt(rsq);
const double cF = 0.5*(cos((MathConst::MY_PI*r)/iparam.cut)+1.0);
const double wX = cF*delx/r;
const double wY = cF*dely/r;
const double wZ = cF*delz/r;
for (k = 0; k < iparam.numeta; ++k) {
const double e = fm_exp(-(iparam.eta[k]*rsq));
Vx[k] += wX*e;
Vy[k] += wY*e;
Vz[k] += wZ*e;
}
}
}
for (j = 0; j < iparam.numtrain; ++j) {
double kx = 0.0;
double ky = 0.0;
double kz = 0.0;
for(int k = 0; k < iparam.numeta; ++k) {
const double xu = iparam.xU[k][j];
kx += square(Vx[k] - xu);
ky += square(Vy[k] - xu);
kz += square(Vz[k] - xu);
}
const double e = -0.5/(square(iparam.sigma));
fxtmp += iparam.alpha[j]*fm_exp(kx*e);
fytmp += iparam.alpha[j]*fm_exp(ky*e);
fztmp += iparam.alpha[j]*fm_exp(kz*e);
}
fxtmp += iparam.b;
fytmp += iparam.b;
fztmp += iparam.b;
f[i].x += fxtmp;
f[i].y += fytmp;
f[i].z += fztmp;
if (EVFLAG) ev_tally_xyz_full_thr(this,i,0.0,0.0,
fxtmp,fytmp,fztmp,
delx,dely,delz,thr);
delete [] Vx;
delete [] Vy;
delete [] Vz;
}
}
/* ---------------------------------------------------------------------- */
double PairAGNIOMP::memory_usage()
{
double bytes = memory_usage_thr();
bytes += PairAGNI::memory_usage();
return bytes;
}

View File

@ -0,0 +1,48 @@
/* -*- c++ -*- ----------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
Contributing author: Axel Kohlmeyer (Temple U)
------------------------------------------------------------------------- */
#ifdef PAIR_CLASS
PairStyle(agni/omp,PairAGNIOMP)
#else
#ifndef LMP_PAIR_AGNI_OMP_H
#define LMP_PAIR_AGNI_OMP_H
#include "pair_agni.h"
#include "thr_omp.h"
namespace LAMMPS_NS {
class PairAGNIOMP : public PairAGNI, public ThrOMP {
public:
PairAGNIOMP(class LAMMPS *);
virtual void compute(int, int);
virtual double memory_usage();
private:
template <int EVFLAG>
void eval(int ifrom, int ito, ThrData * const thr);
};
}
#endif
#endif

View File

@ -0,0 +1,253 @@
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
This software is distributed under the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
Contributing author: Axel Kohlmeyer (Temple U)
------------------------------------------------------------------------- */
#include <math.h>
#include "pair_tersoff_mod_c_omp.h"
#include "atom.h"
#include "comm.h"
#include "force.h"
#include "neighbor.h"
#include "neigh_list.h"
#include "suffix.h"
using namespace LAMMPS_NS;
/* ---------------------------------------------------------------------- */
PairTersoffMODCOMP::PairTersoffMODCOMP(LAMMPS *lmp) :
PairTersoffMODC(lmp), ThrOMP(lmp, THR_PAIR)
{
suffix_flag |= Suffix::OMP;
respa_enable = 0;
}
/* ---------------------------------------------------------------------- */
void PairTersoffMODCOMP::compute(int eflag, int vflag)
{
if (eflag || vflag) {
ev_setup(eflag,vflag);
} else evflag = vflag_fdotr = vflag_atom = 0;
const int nall = atom->nlocal + atom->nghost;
const int nthreads = comm->nthreads;
const int inum = list->inum;
#if defined(_OPENMP)
#pragma omp parallel default(none) shared(eflag,vflag)
#endif
{
int ifrom, ito, tid;
loop_setup_thr(ifrom, ito, tid, inum, nthreads);
ThrData *thr = fix->get_thr(tid);
thr->timer(Timer::START);
ev_setup_thr(eflag, vflag, nall, eatom, vatom, thr);
if (evflag) {
if (eflag) {
if (vflag_atom) eval<1,1,1>(ifrom, ito, thr);
else eval<1,1,0>(ifrom, ito, thr);
} else {
if (vflag_atom) eval<1,0,1>(ifrom, ito, thr);
else eval<1,0,0>(ifrom, ito, thr);
}
} else eval<0,0,0>(ifrom, ito, thr);
thr->timer(Timer::PAIR);
reduce_thr(this, eflag, vflag, thr);
} // end of omp parallel region
}
template <int EVFLAG, int EFLAG, int VFLAG_ATOM>
void PairTersoffMODCOMP::eval(int iifrom, int iito, ThrData * const thr)
{
int i,j,k,ii,jj,kk,jnum;
tagint itag,jtag;
int itype,jtype,ktype,iparam_ij,iparam_ijk;
double xtmp,ytmp,ztmp,delx,dely,delz,evdwl,fpair;
double rsq,rsq1,rsq2;
double delr1[3],delr2[3],fi[3],fj[3],fk[3];
double zeta_ij,prefactor;
int *ilist,*jlist,*numneigh,**firstneigh;
evdwl = 0.0;
const dbl3_t * _noalias const x = (dbl3_t *) atom->x[0];
dbl3_t * _noalias const f = (dbl3_t *) thr->get_f()[0];
const tagint * _noalias const tag = atom->tag;
const int * _noalias const type = atom->type;
const int nlocal = atom->nlocal;
ilist = list->ilist;
numneigh = list->numneigh;
firstneigh = list->firstneigh;
double fxtmp,fytmp,fztmp;
// loop over full neighbor list of my atoms
for (ii = iifrom; ii < iito; ++ii) {
i = ilist[ii];
itag = tag[i];
itype = map[type[i]];
xtmp = x[i].x;
ytmp = x[i].y;
ztmp = x[i].z;
fxtmp = fytmp = fztmp = 0.0;
// two-body interactions, skip half of them
jlist = firstneigh[i];
jnum = numneigh[i];
for (jj = 0; jj < jnum; jj++) {
j = jlist[jj];
j &= NEIGHMASK;
jtag = tag[j];
if (itag > jtag) {
if ((itag+jtag) % 2 == 0) continue;
} else if (itag < jtag) {
if ((itag+jtag) % 2 == 1) continue;
} else {
if (x[j].z < ztmp) continue;
if (x[j].z == ztmp && x[j].y < ytmp) continue;
if (x[j].z == ztmp && x[j].y == ytmp && x[j].x < xtmp) continue;
}
jtype = map[type[j]];
delx = xtmp - x[j].x;
dely = ytmp - x[j].y;
delz = ztmp - x[j].z;
rsq = delx*delx + dely*dely + delz*delz;
iparam_ij = elem2param[itype][jtype][jtype];
if (rsq > params[iparam_ij].cutsq) continue;
repulsive(&params[iparam_ij],rsq,fpair,EFLAG,evdwl);
fxtmp += delx*fpair;
fytmp += dely*fpair;
fztmp += delz*fpair;
f[j].x -= delx*fpair;
f[j].y -= dely*fpair;
f[j].z -= delz*fpair;
if (EVFLAG) ev_tally_thr(this,i,j,nlocal,/* newton_pair */ 1,
evdwl,0.0,fpair,delx,dely,delz,thr);
}
// three-body interactions
// skip immediately if I-J is not within cutoff
double fjxtmp,fjytmp,fjztmp;
for (jj = 0; jj < jnum; jj++) {
j = jlist[jj];
j &= NEIGHMASK;
jtype = map[type[j]];
iparam_ij = elem2param[itype][jtype][jtype];
delr1[0] = x[j].x - xtmp;
delr1[1] = x[j].y - ytmp;
delr1[2] = x[j].z - ztmp;
rsq1 = delr1[0]*delr1[0] + delr1[1]*delr1[1] + delr1[2]*delr1[2];
if (rsq1 > params[iparam_ij].cutsq) continue;
// accumulate bondorder zeta for each i-j interaction via loop over k
fjxtmp = fjytmp = fjztmp = 0.0;
zeta_ij = 0.0;
for (kk = 0; kk < jnum; kk++) {
if (jj == kk) continue;
k = jlist[kk];
k &= NEIGHMASK;
ktype = map[type[k]];
iparam_ijk = elem2param[itype][jtype][ktype];
delr2[0] = x[k].x - xtmp;
delr2[1] = x[k].y - ytmp;
delr2[2] = x[k].z - ztmp;
rsq2 = delr2[0]*delr2[0] + delr2[1]*delr2[1] + delr2[2]*delr2[2];
if (rsq2 > params[iparam_ijk].cutsq) continue;
zeta_ij += zeta(&params[iparam_ijk],rsq1,rsq2,delr1,delr2);
}
// pairwise force due to zeta
force_zeta(&params[iparam_ij],rsq1,zeta_ij,fpair,prefactor,EFLAG,evdwl);
fxtmp += delr1[0]*fpair;
fytmp += delr1[1]*fpair;
fztmp += delr1[2]*fpair;
fjxtmp -= delr1[0]*fpair;
fjytmp -= delr1[1]*fpair;
fjztmp -= delr1[2]*fpair;
if (EVFLAG) ev_tally_thr(this,i,j,nlocal,/* newton_pair */ 1,evdwl,0.0,
-fpair,-delr1[0],-delr1[1],-delr1[2],thr);
// attractive term via loop over k
for (kk = 0; kk < jnum; kk++) {
if (jj == kk) continue;
k = jlist[kk];
k &= NEIGHMASK;
ktype = map[type[k]];
iparam_ijk = elem2param[itype][jtype][ktype];
delr2[0] = x[k].x - xtmp;
delr2[1] = x[k].y - ytmp;
delr2[2] = x[k].z - ztmp;
rsq2 = delr2[0]*delr2[0] + delr2[1]*delr2[1] + delr2[2]*delr2[2];
if (rsq2 > params[iparam_ijk].cutsq) continue;
attractive(&params[iparam_ijk],prefactor,
rsq1,rsq2,delr1,delr2,fi,fj,fk);
fxtmp += fi[0];
fytmp += fi[1];
fztmp += fi[2];
fjxtmp += fj[0];
fjytmp += fj[1];
fjztmp += fj[2];
f[k].x += fk[0];
f[k].y += fk[1];
f[k].z += fk[2];
if (VFLAG_ATOM) v_tally3_thr(i,j,k,fj,fk,delr1,delr2,thr);
}
f[j].x += fjxtmp;
f[j].y += fjytmp;
f[j].z += fjztmp;
}
f[i].x += fxtmp;
f[i].y += fytmp;
f[i].z += fztmp;
}
}
/* ---------------------------------------------------------------------- */
double PairTersoffMODCOMP::memory_usage()
{
double bytes = memory_usage_thr();
bytes += PairTersoffMOD::memory_usage();
return bytes;
}

View File

@ -0,0 +1,43 @@
/* -*- c++ -*- ----------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
Contributing author: Axel Kohlmeyer (Temple U)
------------------------------------------------------------------------- */
#ifdef PAIR_CLASS
PairStyle(tersoff/mod/c/omp,PairTersoffMODCOMP)
#else
#ifndef LMP_PAIR_TERSOFF_MOD_C_OMP_H
#define LMP_PAIR_TERSOFF_MOD_C_OMP_H
#include "pair_tersoff_mod_c.h"
#include "thr_omp.h"
namespace LAMMPS_NS {
class PairTersoffMODCOMP : public PairTersoffMODC, public ThrOMP {
public:
PairTersoffMODCOMP(class LAMMPS *);
virtual void compute(int, int);
virtual double memory_usage();
private:
template <int EVFLAG, int EFLAG, int VFLAG_ATOM>
void eval(int ifrom, int ito, ThrData * const thr);
};
}
#endif
#endif

View File

@ -557,6 +557,37 @@ void ThrOMP::ev_tally_xyz_thr(Pair * const pair, const int i, const int j,
}
}
/* ----------------------------------------------------------------------
tally eng_vdwl and virial into global and per-atom accumulators
for virial, have delx,dely,delz and fx,fy,fz
called when using full neighbor lists
------------------------------------------------------------------------- */
void ThrOMP::ev_tally_xyz_full_thr(Pair * const pair, const int i,
const double evdwl, const double ecoul,
const double fx, const double fy,
const double fz, const double delx,
const double dely, const double delz,
ThrData * const thr)
{
if (pair->eflag_either)
e_tally_thr(pair,i,i,i+1,0,0.5*evdwl,ecoul,thr);
if (pair->vflag_either) {
double v[6];
v[0] = 0.5*delx*fx;
v[1] = 0.5*dely*fy;
v[2] = 0.5*delz*fz;
v[3] = 0.5*delx*fy;
v[4] = 0.5*delx*fz;
v[5] = 0.5*dely*fz;
v_tally_thr(pair,i,i,i+1,0,v,thr);
}
}
/* ----------------------------------------------------------------------
tally eng_vdwl and virial into global and per-atom accumulators
called by SW and hbond potentials, newton_pair is always on

View File

@ -118,6 +118,9 @@ class ThrOMP {
const int, const double, const double, const double,
const double, const double, const double,
const double, const double, ThrData * const);
void ev_tally_xyz_full_thr(Pair * const, const int, const double, const double,
const double, const double, const double,
const double, const double, const double, ThrData * const);
void ev_tally3_thr(Pair * const, const int, const int, const int, const double,
const double, const double * const, const double * const,
const double * const, const double * const, ThrData * const);

View File

@ -1423,7 +1423,7 @@ void PairSMTBQ::tabqeq()
{
gam = dgam = dza = dzb = d2zaa = d2zab =
d2zbb = d2zra = d2zrb = d2gamr2 = 0.0 ;
aCoeff = bCoeff = dij = 0.0 ;
dij = 0.0 ;
s = static_cast<double>(k)*ds ; r = sqrt(s) ;
if (k==0) r=10e-30;
@ -1438,7 +1438,6 @@ void PairSMTBQ::tabqeq()
// Cutting Fonction
if (dij < 0.01 && ii==0)
{
ii=2;
@ -1452,7 +1451,6 @@ void PairSMTBQ::tabqeq()
ddij = aCoeff*(r- rc-nang) *(2+bCoeff*(r-rc-nang))*exp(bCoeff*r);
}
if (r > (rc+nang)) {dij = 0.0 ; ddij = 0.0;}
fafb[k][m] = potqn[k] - dij ;
@ -1471,7 +1469,6 @@ void PairSMTBQ::tabqeq()
rb = ROxSurf;
zb = (2.0*params[j].ne + 1.0)/(4.0*rb); }
ii = 0 ; nang =cang= 5.0 ;
// --------------------------
for (k = 0; k < kmax+5; k++)

View File

@ -42,7 +42,10 @@ using namespace voro;
/* ---------------------------------------------------------------------- */
ComputeVoronoi::ComputeVoronoi(LAMMPS *lmp, int narg, char **arg) :
Compute(lmp, narg, arg)
Compute(lmp, narg, arg), con_mono(NULL), con_poly(NULL),
radstr(NULL), voro(NULL), edge(NULL), sendvector(NULL),
rfield(NULL), tags(NULL), occvec(NULL), sendocc(NULL),
lroot(NULL), lnext(NULL), faces(NULL)
{
int sgroup;

View File

@ -44,8 +44,9 @@ enum{PERATOM,LOCAL};
ComputeReduce::ComputeReduce(LAMMPS *lmp, int narg, char **arg) :
Compute(lmp, narg, arg),
nvalues(0), which(NULL), argindex(NULL), flavor(NULL), value2index(NULL), ids(NULL),
onevec(NULL), replace(NULL), indices(NULL), owner(NULL), idregion(NULL), varatom(NULL)
nvalues(0), which(NULL), argindex(NULL), flavor(NULL),
value2index(NULL), ids(NULL), onevec(NULL), replace(NULL), indices(NULL),
owner(NULL), idregion(NULL), varatom(NULL)
{
int iarg = 0;
if (strcmp(style,"reduce") == 0) {

View File

@ -22,6 +22,10 @@
#include "tbb/scalable_allocator.h"
#endif
#if defined(LMP_USER_INTEL) && !defined(LAMMPS_MEMALIGN)
#define LAMMPS_MEMALIGN 64
#endif
using namespace LAMMPS_NS;
/* ---------------------------------------------------------------------- */

View File

@ -48,6 +48,10 @@ methods:
#ifndef LAMMPS_MY_PAGE_H
#define LAMMPS_MY_PAGE_H
#if defined(LMP_USER_INTEL) && !defined(LAMMPS_MEMALIGN)
#define LAMMPS_MEMALIGN 64
#endif
#include <stdlib.h>
namespace LAMMPS_NS {

View File

@ -1 +1 @@
#define LAMMPS_VERSION "9 Nov 2016"
#define LAMMPS_VERSION "17 Nov 2016"