Merge pull request #1151 from akohlmey/next-patch-release

update version number for next patch release

.github/CODEOWNERS

@@ -30,6 +30,7 @@ src/USER-MOFFF/*      @hheenen
 src/USER-MOLFILE/*    @akohlmey
 src/USER-NETCDF/*     @pastewka
 src/USER-PHONON/*     @lingtikong
+src/USER-PTM/*        @pmla
 src/USER-OMP/*        @akohlmey
 src/USER-QMMM/*       @akohlmey
 src/USER-REAXC/*      @hasanmetin

cmake/CMakeLists.txt

@@ -164,12 +164,13 @@ set(LAMMPS_DEPS)
 set(LAMMPS_API_DEFINES)
 
 set(DEFAULT_PACKAGES ASPHERE BODY CLASS2 COLLOID COMPRESS DIPOLE GRANULAR
-  KSPACE MANYBODY MC MEAM MESSAGE MISC MOLECULE PERI REAX REPLICA RIGID SHOCK SPIN SNAP
-  SRD KIM PYTHON MSCG MPIIO VORONOI POEMS LATTE USER-ATC USER-AWPMD USER-BOCS
-  USER-CGDNA USER-MESO USER-CGSDK USER-COLVARS USER-DIFFRACTION USER-DPD USER-DRUDE
-  USER-EFF USER-FEP USER-H5MD USER-LB USER-MANIFOLD USER-MEAMC USER-MGPT USER-MISC
-  USER-MOFFF USER-MOLFILE USER-NETCDF USER-PHONON USER-QTB USER-REAXC USER-SCAFACOS 
-  USER-SMD USER-SMTBQ USER-SPH USER-TALLY USER-UEF USER-VTK USER-QUIP USER-QMMM)
+  KSPACE MANYBODY MC MEAM MESSAGE MISC MOLECULE PERI REAX REPLICA RIGID SHOCK
+  SPIN SNAP SRD KIM PYTHON MSCG MPIIO VORONOI POEMS LATTE USER-ATC USER-AWPMD
+  USER-BOCS USER-CGDNA USER-MESO USER-CGSDK USER-COLVARS USER-DIFFRACTION
+  USER-DPD USER-DRUDE USER-EFF USER-FEP USER-H5MD USER-LB USER-MANIFOLD
+  USER-MEAMC USER-MGPT USER-MISC USER-MOFFF USER-MOLFILE USER-NETCDF
+  USER-PHONON USER-PTM USER-QTB USER-REAXC USER-SCAFACOS USER-SMD USER-SMTBQ
+  USER-SPH USER-TALLY USER-UEF USER-VTK USER-QUIP USER-QMMM)
 set(ACCEL_PACKAGES USER-OMP KOKKOS OPT USER-INTEL GPU)
 set(OTHER_PACKAGES CORESHELL QEQ)
 foreach(PKG ${DEFAULT_PACKAGES})

doc/Makefile

@@ -38,7 +38,7 @@ OBJECTS=$(SOURCES:src/%.txt=$(RSTDIR)/%.rst)
 help:
 	@echo "Please use \`make <target>' where <target> is one of"
 	@echo "  html       create HTML doc pages in html dir"
-	@echo "  pdf        create Manual.pdf and Developer.pdf in this dir"
+	@echo "  pdf        create Developer.pdf and Manual.pdf in this dir"
 	@echo "  old        create old-style HTML doc pages in old dir"
 	@echo "  fetch      fetch HTML and PDF files from LAMMPS web site"
 	@echo "  epub       create ePUB format manual for e-book readers"
@@ -116,17 +116,17 @@ mobi: epub
 pdf: utils/txt2html/txt2html.exe
 	@(\
 		set -e; \
-		cd src; \
+		cd src/Developer; \
+		pdflatex developer; \
+		pdflatex developer; \
+		mv developer.pdf ../../Developer.pdf; \
+		cd ..; \
 		../utils/txt2html/txt2html.exe -b *.txt; \
 		htmldoc --batch lammps.book;          \
 		for s in `echo *.txt | sed -e 's,\.txt,\.html,g'` ; \
 			do grep -q $$s lammps.book || \
 			echo doc file $$s missing in src/lammps.book; done; \
 		rm *.html; \
-		cd Developer; \
-		pdflatex developer; \
-		pdflatex developer; \
-		mv developer.pdf ../../Developer.pdf; \
 	)
 
 old: utils/txt2html/txt2html.exe

doc/src/Commands_bond.txt

@@ -70,7 +70,7 @@ OPT.
 "fourier/simple (o)"_angle_fourier_simple.html,
 "harmonic (iko)"_angle_harmonic.html,
 "quartic (o)"_angle_quartic.html,
-"sdk"_angle_sdk.html,
+"sdk (o)"_angle_sdk.html,
 "table (o)"_angle_table.html :tb(c=4,ea=c)
 
 :line

doc/src/Commands_compute.txt

@@ -96,6 +96,7 @@ KOKKOS, o = USER-OMP, t = OPT.
 "property/atom"_compute_property_atom.html,
 "property/chunk"_compute_property_chunk.html,
 "property/local"_compute_property_local.html,
+"ptm/atom"_compute_ptm_atom.html,
 "rdf"_compute_rdf.html,
 "reduce"_compute_reduce.html,
 "reduce/chunk"_compute_reduce_chunk.html,
@@ -117,7 +118,7 @@ KOKKOS, o = USER-OMP, t = OPT.
 "smd/tlsph/strain"_compute_smd_tlsph_strain.html,
 "smd/tlsph/strain/rate"_compute_smd_tlsph_strain_rate.html,
 "smd/tlsph/stress"_compute_smd_tlsph_stress.html,
-"smd/triangle/mesh/vertices"_compute_smd_triangle_mesh_vertices.html,
+"smd/triangle/mesh/vertices"_compute_smd_triangle_vertices.html,
 "smd/ulsph/num/neighs"_compute_smd_ulsph_num_neighs.html,
 "smd/ulsph/strain"_compute_smd_ulsph_strain.html,
 "smd/ulsph/strain/rate"_compute_smd_ulsph_strain_rate.html,

doc/src/Commands_fix.txt

@@ -65,6 +65,7 @@ OPT.
 "eos/table/rx (k)"_fix_eos_table_rx.html,
 "evaporate"_fix_evaporate.html,
 "external"_fix_external.html,
+"ffl"_fix_ffl.html,
 "filter/corotate"_fix_filter_corotate.html,
 "flow/gauss"_fix_flow_gauss.html,
 "freeze"_fix_freeze.html,
@@ -216,7 +217,7 @@ OPT.
 "wall/body/polyhedron"_fix_wall_body_polyhedron.html,
 "wall/colloid"_fix_wall.html,
 "wall/ees"_fix_wall_ees.html,
-"wall/gran"_fix_wall_gran.html,
+"wall/gran (o)"_fix_wall_gran.html,
 "wall/gran/region"_fix_wall_gran_region.html,
 "wall/harmonic"_fix_wall.html,
 "wall/lj1043"_fix_wall.html,

doc/src/Commands_pair.txt

@@ -26,7 +26,7 @@ OPT.
 
 "none"_pair_none.html,
 "zero"_pair_zero.html,
-"hybrid"_pair_hybrid.html,
+"hybrid (k)"_pair_hybrid.html,
 "hybrid/overlay (k)"_pair_hybrid.html :tb(c=4,ea=c)
 
 "adp (o)"_pair_adp.html,
@@ -81,6 +81,7 @@ OPT.
 "eam (gikot)"_pair_eam.html,
 "eam/alloy (gikot)"_pair_eam.html,
 "eam/cd (o)"_pair_eam.html,
+"eam/cd/old (o)"_pair_eam.html,
 "eam/fs (gikot)"_pair_eam.html,
 "edip (o)"_pair_edip.html,
 "edip/multi"_pair_edip.html,
@@ -167,7 +168,7 @@ OPT.
 "morse/soft"_pair_morse.html,
 "multi/lucy"_pair_multi_lucy.html,
 "multi/lucy/rx (k)"_pair_multi_lucy_rx.html,
-"nb3b/harmonic (o)"_pair_nb3b_harmonic.html,
+"nb3b/harmonic"_pair_nb3b_harmonic.html,
 "nm/cut (o)"_pair_nm.html,
 "nm/cut/coul/cut (o)"_pair_nm.html,
 "nm/cut/coul/long (o)"_pair_nm.html,

doc/src/Eqs/ptm_rmsd.tex

@@ -0,0 +1,21 @@
+\documentclass[12pt,article]{article}
+
+\usepackage{indentfirst}
+\usepackage{amsmath}
+
+\newcommand{\set}[1]{\ensuremath{\mathbf{#1}}}
+\newcommand{\mean}[1]{\ensuremath{\overline{#1}}}
+\newcommand{\norm}[1]{\ensuremath{\left|\left|{#1}\right|\right|}}
+
+\begin{document}
+
+\begin{equation*}
+\text{RMSD}(\set{u}, \set{v}) = \min_{s, \set{Q}} \sqrt{\frac{1}{N} \sum\limits_{i=1}^{N}
+\norm{
+s[\vec{u_i} - \mean{\set{u}}]
+-
+\set{Q} \vec{v_i}
+}^2}
+\end{equation*}
+
+\end{document}

doc/src/Errors_messages.txt

@@ -1092,11 +1092,6 @@ correct. :dd
 The specified file cannot be opened.  Check that the path and name are
 correct. :dd
 
-{Cannot open fix ave/spatial file %s} :dt
-
-The specified file cannot be opened.  Check that the path and name are
-correct. :dd
-
 {Cannot open fix ave/time file %s} :dt
 
 The specified file cannot be opened.  Check that the path and name are
@@ -1677,10 +1672,6 @@ provided by an atom map.  An atom map does not exist (by default) for
 non-molecular problems.  Using the atom_modify map command will force
 an atom map to be created. :dd
 
-{Cannot use fix ave/spatial z for 2 dimensional model} :dt
-
-Self-explanatory. :dd
-
 {Cannot use fix bond/break with non-molecular systems} :dt
 
 Only systems with bonds that can be changed can be used.  Atom_style
@@ -2425,10 +2416,6 @@ Self-explanatory. :dd
 
 Self-explanatory. :dd
 
-{Compute ID for fix ave/spatial does not exist} :dt
-
-Self-explanatory. :dd
-
 {Compute ID for fix ave/time does not exist} :dt
 
 Self-explanatory. :dd
@@ -4074,10 +4061,6 @@ Self-explanatory. :dd
 
 Self-explanatory. :dd
 
-{Fix ID for fix ave/spatial does not exist} :dt
-
-Self-explanatory. :dd
-
 {Fix ID for fix ave/time does not exist} :dt
 
 Self-explanatory. :dd
@@ -4379,51 +4362,6 @@ same style. :dd
 
 Self-explanatory. :dd
 
-{Fix ave/spatial compute does not calculate a per-atom array} :dt
-
-Self-explanatory. :dd
-
-{Fix ave/spatial compute does not calculate a per-atom vector} :dt
-
-A compute used by fix ave/spatial must generate per-atom values. :dd
-
-{Fix ave/spatial compute does not calculate per-atom values} :dt
-
-A compute used by fix ave/spatial must generate per-atom values. :dd
-
-{Fix ave/spatial compute vector is accessed out-of-range} :dt
-
-The index for the vector is out of bounds. :dd
-
-{Fix ave/spatial fix does not calculate a per-atom array} :dt
-
-Self-explanatory. :dd
-
-{Fix ave/spatial fix does not calculate a per-atom vector} :dt
-
-A fix used by fix ave/spatial must generate per-atom values. :dd
-
-{Fix ave/spatial fix does not calculate per-atom values} :dt
-
-A fix used by fix ave/spatial must generate per-atom values. :dd
-
-{Fix ave/spatial fix vector is accessed out-of-range} :dt
-
-The index for the vector is out of bounds. :dd
-
-{Fix ave/spatial for triclinic boxes requires units reduced} :dt
-
-Self-explanatory. :dd
-
-{Fix ave/spatial settings invalid with changing box size} :dt
-
-If the box size changes, only the units reduced option can be
-used. :dd
-
-{Fix ave/spatial variable is not atom-style variable} :dt
-
-A variable used by fix ave/spatial must generate per-atom values. :dd
-
 {Fix ave/time cannot set output array intensive/extensive from these inputs} :dt
 
 One of more of the vector inputs has individual elements which are

doc/src/Errors_warnings.txt

@@ -291,24 +291,6 @@ This may cause accuracy problems. :dd
 
 This may cause accuracy problems. :dd
 
-{Fix thermal/conductivity comes before fix ave/spatial} :dt
-
-The order of these 2 fixes in your input script is such that fix
-thermal/conductivity comes first.  If you are using fix ave/spatial to
-measure the temperature profile induced by fix viscosity, then this
-may cause a glitch in the profile since you are averaging immediately
-after swaps have occurred.  Flipping the order of the 2 fixes
-typically helps. :dd
-
-{Fix viscosity comes before fix ave/spatial} :dt
-
-The order of these 2 fixes in your input script is such that
-fix viscosity comes first.  If you are using fix ave/spatial
-to measure the velocity profile induced by fix viscosity, then
-this may cause a glitch in the profile since you are averaging
-immediately after swaps have occurred.  Flipping the order
-of the 2 fixes typically helps. :dd
-
 {Fixes cannot send data in Kokkos communication, switching to classic communication} :dt
 
 This is current restriction with Kokkos. :dd

doc/src/Install_linux.txt

@@ -9,39 +9,16 @@ Documentation"_ld - "LAMMPS Commands"_lc :c
 
 Download an executable for Linux :h3
 
-Binaries are available for many different versions of Linux:
+Binaries are available for different versions of Linux:
 
-"Pre-built binary RPMs for Fedora/RedHat/CentOS/openSUSE"_#rpm
 "Pre-built Ubuntu Linux executables"_#ubuntu
+"Pre-built Fedora Linux executables"_#fedora
+"Pre-built EPEL Linux executables (RHEL, CentOS)"_#epel
+"Pre-built OpenSuse Linux executables"_#opensuse
 "Pre-built Gentoo Linux executable"_#gentoo :all(b)
 
 :line
 
-Pre-built binary RPMs for Fedora/RedHat/CentOS/openSUSE :h4,link(rpm)
-
-Pre-built LAMMPS executables for various Linux distributions
-can be downloaded as binary RPM files from this site:
-
-"http://rpm.lammps.org"_http://rpm.lammps.org
-
-There are multiple package variants supporting serial, parallel and
-Python wrapper versions.  The LAMMPS binaries contain all optional
-packages included in the source distribution except: GPU, KIM, REAX,
-and USER-INTEL.
-
-Installation instructions for the various versions are here:
-
-"http://rpm.lammps.org/install.html"_http://rpm.lammps.org/install.html
-
-The instructions show how to enable the repository in the respective
-system's package management system.  Installing and updating are then
-straightforward and automatic.  
-
-Thanks to Axel Kohlmeyer (Temple U, akohlmey at gmail.com) for setting
-up this RPM capability.
-
-:line
-
 Pre-built Ubuntu Linux executables :h4,link(ubuntu)
 
 A pre-built LAMMPS executable suitable for running on the latest
@@ -60,10 +37,10 @@ To install LAMMPS do the following once:
 
 sudo apt-get install lammps-daily :pre
 
-This downloads an executable named "lammps-daily" to your box, which
+This downloads an executable named "lmp_daily" to your box, which
 can then be used in the usual way to run input scripts:
 
-lammps-daily < in.lj :pre
+lmp_daily -in in.lj :pre
 
 To update LAMMPS to the most current version, do the following:
 
@@ -99,6 +76,80 @@ Ubuntu package capability.
 
 :line
 
+Pre-built Fedora Linux executables :h4,link(fedora)
+
+Pre-built LAMMPS packages for stable releases are available
+in the Fedora Linux distribution as of version 28. The packages
+can be installed via the dnf package manager. There are 3 basic
+varieties (lammps = no MPI, lammps-mpich = MPICH MPI library,
+lammps-openmpi = OpenMPI MPI library) and for each support for
+linking to the C library interface (lammps-devel, lammps-mpich-devel,
+lammps-openmpi-devel), the header for compiling programs using
+the C library interface (lammps-headers), and the LAMMPS python
+module for Python 3. All packages can be installed at the same
+time and the name of the LAMMPS executable is {lmp} in all 3 cases.
+By default, {lmp} will refer to the serial executable, unless
+one of the MPI environment modules is loaded 
+("module load mpi/mpich-x86_64" or "module load mpi/openmpi-x86_64").
+Then the corresponding parallel LAMMPS executable is used.
+The same mechanism applies when loading the LAMMPS python module.
+
+To install LAMMPS with OpenMPI and run an input in.lj with 2 CPUs do:
+
+dnf install lammps-openmpi
+module load mpi/openmpi-x86_64
+mpirun -np 2 lmp -in in.lj :pre
+
+The "dnf install" command is needed only once. In case of a new LAMMPS
+stable release, "dnf update" will automatically update to the newer
+version as soon at the RPM files are built and uploaded to the download
+mirrors. The "module load" command is needed once per (shell) session
+or shell terminal instance, unless it is automatically loaded from the
+shell profile.
+
+Please use "lmp -help" to see which compilation options, packages,
+and styles are included in the binary.
+
+Thanks to Christoph Junghans (LANL) for making LAMMPS available in Fedora.
+
+:line
+
+Pre-built EPEL Linux executable :h4,link(epel)
+
+Pre-built LAMMPS packages for stable releases are available
+in the "Extra Packages for Enterprise Linux (EPEL) repository"_https://fedoraproject.org/wiki/EPEL
+for use with Red Hat Enterprise Linux (RHEL) or CentOS version 7.x
+and compatible Linux distributions. Names of packages, executable,
+and content are the same as described above for Fedora Linux.
+But RHEL/CentOS 7.x uses the "yum" package manager instead of "dnf"
+in Fedora 28.
+
+Please use "lmp -help" to see which compilation options, packages,
+and styles are included in the binary.
+
+Thanks to Christoph Junghans (LANL) for making LAMMPS available in EPEL.
+
+:line
+
+Pre-built OpenSuse Linux executable :h4,link(opensuse)
+
+A pre-built LAMMPS package for stable releases is available
+in OpenSuse as of Leap 15.0. You can install the package with:
+
+zypper install lammps :pre
+
+This includes support for OpenMPI. The name of the LAMMPS executable
+is {lmp}. Thus to run an input in parallel on 2 CPUs you would do:
+
+mpirun -np 2 lmp -in in.lj :pre
+
+Please use "lmp -help" to see which compilation options, packages,
+and styles are included in the binary.
+
+Thanks to Christoph Junghans (LANL) for making LAMMPS available in OpenSuse.
+
+:line
+
 Pre-built Gentoo Linux executable :h4,link(gentoo)
 
 LAMMPS is part of Gentoo's main package tree and can be installed by

doc/src/Manual.txt

@@ -1,7 +1,7 @@
 <!-- HTML_ONLY -->
 <HEAD>
 <TITLE>LAMMPS Users Manual</TITLE>
-<META NAME="docnumber" CONTENT="18 Sep 2018 version">
+<META NAME="docnumber" CONTENT="10 Oct 2018 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 @@
 :line
 
 LAMMPS Documentation :c,h1
-18 Sep 2018 version :c,h2
+10 Oct 2018 version :c,h2
 
 "What is a LAMMPS version?"_Manual_version.html
 

doc/src/Packages_details.txt

@@ -89,6 +89,7 @@ as contained in the file name.
 "USER-NETCDF"_#PKG-USER-NETCDF,
 "USER-OMP"_#PKG-USER-OMP,
 "USER-PHONON"_#PKG-USER-PHONON,
+"USER-PTM"_#PKG-USER-PTM,
 "USER-QMMM"_#PKG-USER-QMMM,
 "USER-QTB"_#PKG-USER-QTB,
 "USER-QUIP"_#PKG-USER-QUIP,
@@ -1744,6 +1745,24 @@ examples/USER/phonon :ul
 
 :line
 
+USER-PTM package :link(PKG-USER-PTM),h4
+
+[Contents:]
+
+A "compute ptm/atom"_compute_ptm_atom.html command that calculates
+local structure characterization using the Polyhedral Template
+Matching methodology.
+
+[Author:] Peter Mahler Larsen (MIT).
+
+[Supporting info:]
+
+src/USER-PTM: filename starting with ptm_ -> supporting code, other filenames -> commands
+src/USER-PTM/LICENSE
+"compute ptm/atom"_compute_ptm_atom.html :ul
+
+:line
+
 USER-QMMM package :link(PKG-USER-QMMM),h4
 
 [Contents:]

doc/src/Packages_user.txt

@@ -62,6 +62,7 @@ Package, Description, Doc page, Example, Library
 "USER-NETCDF"_Packages_details.html#PKG-USER-NETCDF, dump output via NetCDF,"dump netcdf"_dump_netcdf.html, n/a, ext
 "USER-OMP"_Packages_details.html#PKG-USER-OMP, OpenMP-enabled styles,"Speed omp"_Speed_omp.html, "Benchmarks"_http://lammps.sandia.gov/bench.html, no
 "USER-PHONON"_Packages_details.html#PKG-USER-PHONON, phonon dynamical matrix,"fix phonon"_fix_phonon.html, USER/phonon, no
+"USER-PTM"_Packages_details.html#PKG-USER-PTM, Polyhedral Template Matching,"compute ptm/atom"_compute_ptm_atom.html, n/a, no
 "USER-QMMM"_Packages_details.html#PKG-USER-QMMM, QM/MM coupling,"fix qmmm"_fix_qmmm.html, USER/qmmm, ext
 "USER-QTB"_Packages_details.html#PKG-USER-QTB, quantum nuclear effects,"fix qtb"_fix_qtb.html "fix qbmsst"_fix_qbmsst.html, qtb, no
 "USER-QUIP"_Packages_details.html#PKG-USER-QUIP, QUIP/libatoms interface,"pair_style quip"_pair_quip.html, USER/quip, ext

doc/src/angle_sdk.txt

@@ -7,6 +7,7 @@
 :line
 
 angle_style sdk command :h3
+angle_style sdk/omp command :h3
 
 [Syntax:]
 
@@ -43,6 +44,30 @@ internally; hence the units of K are in energy/radian^2.
 The also required {lj/sdk} parameters will be extracted automatically
 from the pair_style.
 
+:line
+
+Styles with a {gpu}, {intel}, {kk}, {omp}, or {opt} suffix are
+functionally the same as the corresponding style without the suffix.
+They have been optimized to run faster, depending on your available
+hardware, as discussed on the "Speed packages"_Speed_packages.html doc
+page.  The accelerated styles take the same arguments and should
+produce the same results, except for round-off and precision issues.
+
+These accelerated styles are part of the GPU, USER-INTEL, KOKKOS,
+USER-OMP and OPT packages, respectively.  They are only enabled if
+LAMMPS was built with those packages.  See the "Build
+package"_Build_package.html doc page for more info.
+
+You can specify the accelerated styles explicitly in your input script
+by including their suffix, or you can use the "-suffix command-line
+switch"_Run_options.html when you invoke LAMMPS, or you can use the
+"suffix"_suffix.html command in your input script.
+
+See the "Speed packages"_Speed_packages.html doc page for more
+instructions on how to use the accelerated styles effectively.
+
+:line
+
 [Restrictions:]
 
 This angle style can only be used if LAMMPS was built with the

doc/src/compute_msd_chunk.txt

@@ -90,12 +90,12 @@ This is so that the fix this compute creates to store per-chunk
 quantities will also have the same ID, and thus be initialized
 correctly with chunk reference positions from the restart file.
 
-The simplest way to output the results of the compute com/msd
+The simplest way to output the results of the compute msd/chunk
 calculation to a file is to use the "fix ave/time"_fix_ave_time.html
 command, for example:
 
 compute cc1 all chunk/atom molecule
-compute myChunk all com/msd cc1
+compute myChunk all msd/chunk cc1
 fix 1 all ave/time 100 1 100 c_myChunk\[*\] file tmp.out mode vector :pre
 
 [Output info:]

doc/src/compute_pair.txt

@@ -10,17 +10,20 @@ compute pair command :h3
 
 [Syntax:]
 
-compute ID group-ID pair pstyle evalue :pre
+compute ID group-ID pair pstyle \[nstyle\] \[evalue\]  :pre
 
-ID, group-ID are documented in "compute"_compute.html command
-pair = style name of this compute command
-pstyle = style name of a pair style that calculates additional values
-evalue = {epair} or {evdwl} or {ecoul} or blank (optional setting) :ul
+ID, group-ID are documented in "compute"_compute.html command :ulb,l
+pair = style name of this compute command :l
+pstyle = style name of a pair style that calculates additional values :l
+nsub = {n}-instance of a substyle, if a pair style is used multiple times in a hybrid style :l
+{evalue} = {epair} or {evdwl} or {ecoul} or blank (optional) :l
+:ule
 
 [Examples:]
 
 compute 1 all pair gauss
 compute 1 all pair lj/cut/coul/cut ecoul
+compute 1 all pair tersoff 2 epair
 compute 1 all pair reax :pre
 
 [Description:]
@@ -33,15 +36,19 @@ NOTE: The group specified for this command is [ignored].
 
 The specified {pstyle} must be a pair style used in your simulation
 either by itself or as a sub-style in a "pair_style hybrid or
-hybrid/overlay"_pair_hybrid.html command.
+hybrid/overlay"_pair_hybrid.html command. If the sub-style is
+used more than once, an additional number {nsub} has to be specified
+in order to choose which instance of the sub-style will be used by
+the compute. Not specifying the number in this case will cause the
+compute to fail.
 
-The {evalue} setting is optional; it may be left off the command.  All
+The {evalue} setting is optional.  All
 pair styles tally a potential energy {epair} which may be broken into
 two parts: {evdwl} and {ecoul} such that {epair} = {evdwl} + {ecoul}.
 If the pair style calculates Coulombic interactions, their energy will
 be tallied in {ecoul}.  Everything else (whether it is a Lennard-Jones
 style van der Waals interaction or not) is tallied in {evdwl}.  If
-{evalue} is specified as {epair} or left out, then {epair} is stored
+{evalue} is blank or specified as {epair}, then {epair} is stored
 as a global scalar by this compute.  This is useful when using
 "pair_style hybrid"_pair_hybrid.html if you want to know the portion
 of the total energy contributed by one sub-style.  If {evalue} is
@@ -82,4 +89,4 @@ the doc page for the pair style for details.
 
 [Default:]
 
-The default for {evalue} is {epair}.
+The keyword defaults are {evalue} = {epair}, nsub = 0.

doc/src/compute_pressure_cylinder.txt

@@ -0,0 +1,81 @@
+"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
+
+compute pressure/cylinder command :h3
+
+[Syntax:]
+
+compute ID group-ID pressure/cylinder zlo zhi Rmax bin_width :pre
+
+ID, group-ID are documented in "compute"_compute.html command
+pressure/cylinder = style name of this compute command
+zlo = minimum z-boundary for cylinder
+zhi = maximum z-boundary for cylinder
+Rmax = maximum radius to perform calculation to
+bin_width = width of radial bins to use for calculation :ul
+
+[Examples:]
+
+compute 1 all pressure/cylinder -10.0 10.0 15.0 0.25 :pre
+
+[Description:]
+
+Define a computation that calculates the pressure tensor of a system in
+cylindrical coordinates, as discussed in "(Addington)"_#Addington1.
+This is useful for systems with a single axis of rotational symmetry,
+such as cylindrical micelles or carbon nanotubes. The compute splits the
+system into radial, cylindrical-shell-type bins of width bin_width,
+centered at x=0,y=0, and calculates the radial (P_rhorho), azimuthal
+(P_phiphi), and axial (P_zz) components of the configurational pressure
+tensor. The local density is also calculated for each bin, so that the
+true pressure can be recovered as P_kin+P_conf=density*k*T+P_conf.  The
+output is a global array with 5 columns; one each for bin radius, local
+number density, P_rhorho, P_phiphi, and P_zz. The number of rows is
+governed by the values of Rmax and bin_width. Pressure tensor values are
+output in pressure units.
+
+[Output info:]
+
+This compute calculates a global array with 5 columns and Rmax/bin_width
+rows. The output columns are: R (distance units), number density (inverse 
+volume units), configurational radial pressure (pressure units), 
+configurational azimuthal pressure (pressure units), and configurational
+axial pressure (pressure units).
+
+The values calculated by this compute are
+"intensive".  The pressure values will be in pressure
+"units"_units.html. The number density values will be in 
+inverse volume "units"_units.html.
+
+[Restrictions:]
+
+This compute currently calculates the pressure tensor contributions
+for pair styles only (i.e. no bond, angle, dihedral, etc. contributions
+and in the presence of bonded interactions, the result will be incorrect
+due to exclusions for special bonds)  and requires pair-wise force
+calculations not available for most manybody pair styles. K-space
+calculations are also excluded. Note that this pressure compute outputs
+the configurational terms only; the kinetic contribution is not included
+and may be calculated from the number density output by P_kin=density*k*T.
+
+This compute is part of the USER-MISC package.  It is only enabled
+if LAMMPS was built with that package.  See the "Build
+package"_Build_package.html doc page for more info.
+
+[Related commands:]
+
+"compute temp"_compute_temp.html, "compute
+stress/atom"_compute_stress_atom.html,
+"thermo_style"_thermo_style.html,
+
+[Default:] none
+
+:line
+
+:link(Addington1)
+[(Addington)] Addington, Long, Gubbins, J Chem Phys, 149, 084109 (2018).

doc/src/compute_ptm_atom.txt

@@ -0,0 +1,121 @@
+"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
+
+compute ptm/atom command :h3
+
+[Syntax:]
+
+compute ID group-ID ptm/atom structures threshold :pre
+
+ID, group-ID are documented in "compute"_compute.html command
+ptm/atom = style name of this compute command
+structures = structure types to search for
+threshold = lattice distortion threshold (RMSD) :ul
+
+[Examples:]
+
+compute 1 all ptm/atom default 0.1
+compute 1 all ptm/atom fcc-hcp-dcub-dhex 0.15
+compute 1 all ptm/atom all 0 :pre
+
+[Description:]
+
+Define a computation that determines the local lattice structure
+around an atom using the PTM (Polyhedral Template Matching) method.
+The PTM method is described in "(Larsen)"_#Larsen.
+
+Currently, there are seven lattice structures PTM recognizes:
+
+fcc = 1
+hcp = 2
+bcc = 3
+ico (icosahedral) = 4
+sc (simple cubic) = 5
+dcub (diamond cubic) = 6
+dhex (diamond hexagonal) = 7
+other = 8 :ul
+
+The value of the PTM structure will be 0 for atoms not in the specified
+compute group.  The choice of structures to search for can be specified using the "structures"
+argument, which is a hyphen-separated list of structure keywords.
+Two convenient pre-set options are provided:
+
+default: fcc-hcp-bcc-ico
+all: fcc-hcp-bcc-ico-sc-dcub-dhex :ul
+
+The 'default' setting detects the same structures as the Common Neighbor Analysis method.
+The 'all' setting searches for all structure types.  A small performance penalty is
+incurred for the diamond structures, so it is not recommended to use this option if
+it is known that the simulation does not contain diamond structures.
+
+
+PTM identifies structures using two steps.  First, a graph isomorphism test is used
+to identify potential structure matches.  Next, the deviation is computed between the
+local structure (in the simulation) and a template of the ideal lattice structure.
+The deviation is calculated as:
+
+:c,image(Eqs/ptm_rmsd.jpg)
+
+Here, u and v contain the coordinates of the local and ideal structures respectively,
+s is a scale factor, and Q is a rotation.  The best match is identified by the
+lowest RMSD value, using the optimal scaling, rotation, and correspondence between the
+points.
+
+The 'threshold' keyword sets an upper limit on the maximum permitted deviation before
+a local structure is identified as disordered.  Typical values are in the range 0.1-0.15,
+but larger values may be desirable at higher temperatures.
+A value of 0 is equivalent to infinity and can be used if no threshold is desired.
+
+
+The neighbor list needed to compute this quantity is constructed each
+time the calculation is performed (e.g. each time a snapshot of atoms
+is dumped).  Thus it can be inefficient to compute/dump this quantity
+too frequently or to have multiple compute/dump commands, each with a
+{ptm/atom} style.
+
+[Output info:]
+
+This compute calculates a per-atom arry, which can be accessed by
+any command that uses per-atom values from a compute as input.  See
+the "Howto output"_Howto_output.html doc page for an overview of
+LAMMPS output options.
+
+Results are stored in the per-atom array in the following order:
+
+type
+rmsd
+interatomic distance
+qw
+qx
+qy
+qw :ul
+
+The type is a number from 0 to 8.  The rmsd is a positive real number.
+The interatomic distance is computed from the scale factor in the RMSD equation.
+The (qw,qx,qy,qz) parameters represent the orientation of the local structure
+in quaternion form.  The reference coordinates for each template (from which the
+orientation is determined) can be found in the {ptm_constants.h} file in the PTM source directory.
+
+[Restrictions:]
+
+This fix is part of the USER-PTM package.  It is only enabled if
+LAMMPS was built with that package.  See the "Build
+package"_Build_package.html doc page for more info.
+
+[Related commands:]
+
+"compute centro/atom"_compute_centro_atom.html
+"compute cna/atom"_compute_cna_atom.html
+
+[Default:] none
+
+:line
+
+:link(Larsen)
+[(Larsen)] Larsen, Schmidt, Schiøtz, Modelling Simul Mater Sci Eng, 24, 055007 (2016).
+

doc/src/compute_smd_triangle_vertices.txt

@@ -6,14 +6,14 @@
 
 :line
 
-compute smd/triangle/mesh/vertices :h3
+compute smd/triangle/vertices command :h3
 
 [Syntax:]
 
-compute ID group-ID smd/triangle/mesh/vertices :pre
+compute ID group-ID smd/triangle/vertices :pre
 
 ID, group-ID are documented in "compute"_compute.html command
-smd/triangle/mesh/vertices = style name of this compute command :ul
+smd/triangle/vertices = style name of this compute command :ul
 
 [Examples:]
 

doc/src/compute_spin.txt

@@ -10,14 +10,14 @@ compute spin command :h3
 
 [Syntax:]
 
-compute ID group-ID compute/spin :pre
+compute ID group-ID spin :pre
 
 ID, group-ID are documented in "compute"_compute.html command
-compute/spin = style name of this compute command :ul
+spin = style name of this compute command :ul
 
 [Examples:]
 
-compute out_mag all compute/spin :pre
+compute out_mag all spin :pre
 
 [Description:]
 
@@ -26,7 +26,8 @@ of atoms having spins.
 
 This compute calculates 6 magnetic quantities.
 
-The three first quantities are the x,y and z coordinates of the total magnetization.
+The three first quantities are the x,y and z coordinates of the total
+magnetization.
 
 The fourth quantity is the norm of the total magnetization.
 
@@ -39,7 +40,7 @@ The simplest way to output the results of the compute spin calculation
 is to define some of the quantities as variables, and to use the thermo and
 thermo_style commands, for example:
 
-compute out_mag		all compute/spin :pre
+compute out_mag		all spin :pre
 
 variable mag_z      	equal c_out_mag\[3\]
 variable mag_norm	equal c_out_mag\[4\]
@@ -53,7 +54,6 @@ the total magnetization, and the magnetic temperature. Three variables are
 assigned to those quantities. The thermo and thermo_style commands print them 
 every 10 timesteps.
 
-
 [Output info:]
 
 The array values are "intensive".  The array values will be in
@@ -68,7 +68,6 @@ has to be "spin" for this compute to be valid.
 
 [Related commands:] none
 
-
 [Default:] none
 
 :line

doc/src/computes.txt

@@ -67,10 +67,12 @@ Computes :h1
    compute_pe_atom
    compute_plasticity_atom
    compute_pressure
+   compute_pressure_cylinder
    compute_pressure_uef
    compute_property_atom
    compute_property_chunk
    compute_property_local
+   compute_ptm_atom
    compute_rdf
    compute_reduce
    compute_reduce_chunk
@@ -91,7 +93,7 @@ Computes :h1
    compute_smd_tlsph_strain
    compute_smd_tlsph_strain_rate
    compute_smd_tlsph_stress
-   compute_smd_triangle_mesh_vertices
+   compute_smd_triangle_vertices
    compute_smd_ulsph_num_neighs
    compute_smd_ulsph_strain
    compute_smd_ulsph_strain_rate

doc/src/dihedral_nharmonic.txt

@@ -16,7 +16,7 @@ dihedral_style nharmonic :pre
 [Examples:]
 
 dihedral_style nharmonic
-dihedral_coeff 3 10.0 20.0 30.0 :pre
+dihedral_coeff * 3 10.0 20.0 30.0 :pre
 
 [Description:]
 

doc/src/dump_image.txt

@@ -384,12 +384,7 @@ change this via the "dump_modify"_dump_modify.html command.
 :line
 
 The {fix} keyword can be used with a "fix"_fix.html that produces
-objects to be drawn.  An example is the "fix
-surface/global"_fix_surface_global.html command which can draw lines
-or triangles for 2d/3d simulations.
-
-NOTE: Aug 2016 - The fix surface/global command is not yet added to
-LAMMPS.
+objects to be drawn.
 
 The {fflag1} and {fflag2} settings are numerical values which are
 passed to the fix to affect how the drawing of its objects is done.

doc/src/fix_ffl.txt

@@ -0,0 +1,124 @@
+<script type="text/javascript"
+  src="https://cdn.mathjax.org/mathjax/latest/MathJax.js?config=TeX-AMS-MML_HTMLorMML">
+</script>
+<script type="text/x-mathjax-config">
+  MathJax.Hub.Config({ TeX: { equationNumbers: {autoNumber: "AMS"} } });
+</script>
+
+"LAMMPS WWW Site"_lws - "LAMMPS Documentation"_ld - "LAMMPS Commands"_lc :c
+
+:link(lws,http://lammps.sandia.gov)
+:link(ld,Manual.html)
+:link(lc,Commands_all.html)
+
+:line
+
+fix ffl command :h3
+
+[Syntax:]
+
+fix ID id-group ffl tau Tstart Tstop seed \[flip-type\]  :pre
+
+ID, group-ID are documented in "fix"_fix.html command :ulb,l
+ffl = style name of this fix command :l
+tau = thermostat parameter (positive real) :l
+Tstart, Tstop = temperature ramp during the run :l
+seed = random number seed to use for generating noise (positive integer) :l
+one more value may be appended :l
+    flip-type  = determines the flipping type, can be chosen between rescale - no_flip - hard - soft, if no flip type is given, rescale will be chosen by default :pre
+:ule
+
+[Examples:]
+
+fix 3 boundary ffl 10 300 300 31415
+fix 1 all ffl 100 500 500 9265 soft :pre
+
+[Description:]
+
+Apply a Fast-Forward Langevin Equation (FFL) thermostat as described
+in "(Hijazi)"_#Hijazi. Contrary to
+"fix langevin"_fix_langevin.html, this fix performs both
+thermostatting and evolution of the Hamiltonian equations of motion, so it
+should not be used together with "fix nve"_fix_nve.html -- at least not
+on the same atom groups.
+
+The time-evolution of a single particle undergoing Langevin dynamics is described
+by the equations
+
+\begin\{equation\} \frac \{dq\}\{dt\} = \frac\{p\}\{m\}, \end\{equation\}
+
+\begin\{equation\} \frac \{dp\}\{dt\} = -\gamma p + W + F, \end\{equation\}
+
+where \(F\) is the physical force, \(\gamma\) is the friction coefficient, and \(W\) is a
+Gaussian random force.
+
+The friction coefficient is the inverse of the thermostat parameter : \(\gamma = 1/\tau\), with \(\tau\) the thermostat parameter {tau}.
+The thermostat parameter is given in the time units, \(\gamma\) is in inverse time units.
+
+Equilibrium sampling a temperature T is obtained by specifying the
+target value as the {Tstart} and {Tstop} arguments, so that the internal
+constants depending on the temperature are computed automatically.
+
+The random number {seed} must be a positive integer.  A Marsaglia random
+number generator is used.  Each processor uses the input seed to
+generate its own unique seed and its own stream of random numbers.
+Thus the dynamics of the system will not be identical on two runs on
+different numbers of processors.
+
+The flipping type {flip-type} can be chosen between 4 types described in
+"(Hijazi)"_#Hijazi. The flipping operation occurs during the thermostatting
+step and it flips the momenta of the atoms. If no_flip is chosen, no flip
+will be executed and the integration will be the same as a standard
+Langevin thermostat "(Bussi)"_#Bussi3. The other flipping types are : rescale - hard - soft.
+
+[Restart, fix_modify, output, run start/stop, minimize info:]
+
+The instantaneous values of the extended variables are written to
+"binary restart files"_restart.html.  Because the state of the random
+number generator is not saved in restart files, this means you cannot
+do "exact" restarts with this fix, where the simulation continues on
+the same as if no restart had taken place. However, in a statistical
+sense, a restarted simulation should produce the same behavior.
+Note however that you should use a different seed each time you
+restart, otherwise the same sequence of random numbers will be used
+each time, which might lead to stochastic synchronization and
+subtle artefacts in the sampling.
+
+This fix can ramp its target temperature over multiple runs, using the
+{start} and {stop} keywords of the "run"_run.html command.  See the
+"run"_run.html command for details of how to do this.
+
+The "fix_modify"_fix_modify.html {energy} option is supported by this
+fix to add the energy change induced by Langevin thermostatting to the
+system's potential energy as part of "thermodynamic
+output"_thermo_style.html.
+
+This fix computes a global scalar which can be accessed by various
+"output commands"_Howto_output.html.  The scalar is the cumulative
+energy change due to this fix.  The scalar value calculated by this
+fix is "extensive".
+
+[Restrictions:]
+
+In order to perform constant-pressure simulations please use
+"fix press/berendsen"_fix_press_berendsen.html, rather than
+"fix npt"_fix_nh.html, to avoid duplicate integration of the
+equations of motion.
+
+This fix is part of the USER-MISC package.  It is only enabled if
+LAMMPS was built with that package.  See the "Build
+package"_Build_package.html doc page for more info.
+
+[Related commands:]
+
+"fix nvt"_fix_nh.html, "fix temp/rescale"_fix_temp_rescale.html, "fix
+viscous"_fix_viscous.html, "fix nvt"_fix_nh.html, "pair_style
+dpd/tstat"_pair_dpd.html, "fix gld"_fix_gld.html, "fix gle"_fix_gle.html
+
+:line
+
+:link(Hijazi)
+[(Hijazi)] M. Hijazi, D. M. Wilkins, M. Ceriotti, J. Chem. Phys. 148, 184109 (2018)
+:link(Bussi3)
+[(Bussi)] G. Bussi, M. Parrinello, Phs. Rev. E 75, 056707 (2007)
+

doc/src/fix_msst.txt

@@ -6,7 +6,7 @@
 
 :line
 
- fix msst command :h3
+fix msst command :h3
 
 [Syntax:]
 

doc/src/fix_nve_awpmd.txt

@@ -0,0 +1,56 @@
+"LAMMPS WWW Site"_lws - "LAMMPS Documentation"_ld - "LAMMPS Commands"_lc :c
+
+:link(lws,http://lammps.sandia.gov)
+:link(ld,Manual.html)
+:link(lc,Commands_all.html)
+
+:line
+
+fix nve/awpmd command :h3
+
+[Syntax:]
+
+fix ID group-ID nve/awpmd :pre
+
+ID, group-ID are documented in "fix"_fix.html command
+nve/awpmd = style name of this fix command :ul
+
+[Examples:]
+
+fix 1 all nve/awpmd :pre
+
+[Description:]
+
+Perform constant NVE integration to update position and velocity for
+nuclei and electrons in the group for the "Antisymmetrized Wave Packet
+Molecular Dynamics"_pair_awpmd.html model.  V is volume; E is energy.
+This creates a system trajectory consistent with the microcanonical
+ensemble.
+
+The operation of this fix is exactly like that described by the "fix
+nve"_fix_nve.html command, except that the width and width-velocity of
+the electron wavefunctions are also updated.
+
+:line
+
+[Restart, fix_modify, output, run start/stop, minimize info:]
+
+No information about this fix is written to "binary restart
+files"_restart.html.  None of the "fix_modify"_fix_modify.html options
+are relevant to this fix.  No global or per-atom quantities are stored
+by this fix for access by various "output commands"_Howto_output.html.
+No parameter of this fix can be used with the {start/stop} keywords of
+the "run"_run.html command.  This fix is not invoked during "energy
+minimization"_minimize.html.
+
+[Restrictions:]
+
+This fix is part of the USER-AWPMD package.  It is only enabled if
+LAMMPS was built with that package.  See the "Build
+package"_Build_package.html doc page for more info.
+
+[Related commands:]
+
+"fix nve"_fix_nve.html
+
+[Default:] none

doc/src/fix_poems.txt

@@ -6,7 +6,7 @@
 
 :line
 
-fix poems :h3
+fix poems command :h3
 
 Syntax:
 

doc/src/fix_property_atom.txt

@@ -7,6 +7,7 @@
 :line
 
 fix property/atom command :h3
+fix property/atom/kk command :h3
 
 [Syntax:]
 
@@ -201,6 +202,7 @@ added classes.
 :line
 
 :link(isotopes)
+
 Example for using per-atom masses with TIP4P water to
 study isotope effects. When setting up simulations with the "TIP4P
 pair styles"_Howto_tip4p.html for water, you have to provide exactly
@@ -238,6 +240,28 @@ set group hwat mass 2.0141018 :pre
 
 :line
 
+Styles with a {gpu}, {intel}, {kk}, {omp}, or {opt} suffix are
+functionally the same as the corresponding style without the suffix.
+They have been optimized to run faster, depending on your available
+hardware, as discussed on the "Speed packages"_Speed_packages.html doc
+page.  The accelerated styles take the same arguments and should
+produce the same results, except for round-off and precision issues.
+
+These accelerated styles are part of the GPU, USER-INTEL, KOKKOS,
+USER-OMP and OPT packages, respectively.  They are only enabled if
+LAMMPS was built with those packages.  See the "Build
+package"_Build_package.html doc page for more info.
+
+You can specify the accelerated styles explicitly in your input script
+by including their suffix, or you can use the "-suffix command-line
+switch"_Run_options.html when you invoke LAMMPS, or you can use the
+"suffix"_suffix.html command in your input script.
+
+See the "Speed packages"_Speed_packages.html doc page for more
+instructions on how to use the accelerated styles effectively.
+
+:line
+
 [Restart, fix_modify, output, run start/stop, minimize info:]
 
 This fix writes the per-atom values it stores to "binary restart

doc/src/fix_smd_move_triangulated_surface.txt

@@ -73,7 +73,7 @@ package"_Build_package.html doc page for more info.
 
 [Related commands:]
 
-"smd/triangle_mesh_vertices"_compute_smd_triangle_mesh_vertices.html,
+"smd/triangle_mesh_vertices"_compute_smd_triangle_vertices.html,
 "smd/wall_surface"_fix_smd_wall_surface.html
 
 [Default:] none

doc/src/fix_smd_wall_surface.txt

@@ -64,7 +64,7 @@ multiple objects in one file.
 
 [Related commands:]
 
-"smd/triangle_mesh_vertices"_compute_smd_triangle_mesh_vertices.html,
+"smd/triangle_mesh_vertices"_compute_smd_triangle_vertices.html,
 "smd/move_tri_surf"_fix_smd_move_triangulated_surface.html,
 "smd/tri_surface"_pair_smd_triangulated_surface.html
 

doc/src/fix_surface_global.txt

@@ -1,19 +0,0 @@
-"LAMMPS WWW Site"_lws - "LAMMPS Documentation"_ld - "LAMMPS Commands"_lc :c
-
-:link(lws,http://lammps.sandia.gov)
-:link(ld,Manual.html)
-:link(lc,Commands_all.html)
-
-:line
-
-fix wall/surface/global command :h3
-
-[Description:]
-
-This feature is not yet implemented.
-
-[Related commands:]
-
-"dump image"_dump_image.html
-
-

doc/src/fix_wall_gran.txt

@@ -7,6 +7,7 @@
 :line
 
 fix wall/gran command :h3
+fix wall/gran/omp command :h3
 
 [Syntax:]
 
@@ -136,6 +137,28 @@ the clockwise direction for {vshear} > 0 or counter-clockwise for
 {vshear} < 0.  In this case, {vshear} is the tangential velocity of
 the wall at whatever {radius} has been defined.
 
+:line
+
+Styles with a {gpu}, {intel}, {kk}, {omp}, or {opt} suffix are
+functionally the same as the corresponding style without the suffix.
+They have been optimized to run faster, depending on your available
+hardware, as discussed on the "Speed packages"_Speed_packages.html doc
+page.  The accelerated styles take the same arguments and should
+produce the same results, except for round-off and precision issues.
+
+These accelerated styles are part of the GPU, USER-INTEL, KOKKOS,
+USER-OMP and OPT packages, respectively.  They are only enabled if
+LAMMPS was built with those packages.  See the "Build
+package"_Build_package.html doc page for more info.
+
+You can specify the accelerated styles explicitly in your input script
+by including their suffix, or you can use the "-suffix command-line
+switch"_Run_options.html when you invoke LAMMPS, or you can use the
+"suffix"_suffix.html command in your input script.
+
+See the "Speed packages"_Speed_packages.html doc page for more
+instructions on how to use the accelerated styles effectively.
+
 [Restart, fix_modify, output, run start/stop, minimize info:]
 
 This fix writes the shear friction state of atoms interacting with the

doc/src/fixes.txt

@@ -46,6 +46,7 @@ Fixes :h1
    fix_eos_table_rx
    fix_evaporate
    fix_external
+   fix_ffl
    fix_filter_corotate
    fix_flow_gauss
    fix_freeze
@@ -92,6 +93,7 @@ Fixes :h1
    fix_nve
    fix_nve_asphere
    fix_nve_asphere_noforce
+   fix_nve_awpmd
    fix_nve_body
    fix_nve_dot
    fix_nve_dotc_langevin
@@ -154,7 +156,6 @@ Fixes :h1
    fix_srd
    fix_store_force
    fix_store_state
-   fix_surface_global
    fix_temp_berendsen
    fix_temp_csvr
    fix_temp_rescale

doc/src/lammps.book

@@ -266,6 +266,7 @@ fix_eos_table.html
 fix_eos_table_rx.html
 fix_evaporate.html
 fix_external.html
+fix_ffl.html
 fix_filter_corotate.html
 fix_flow_gauss.html
 fix_freeze.html
@@ -311,6 +312,7 @@ fix_npt_sphere.html
 fix_nve.html
 fix_nve_asphere.html
 fix_nve_asphere_noforce.html
+fix_nve_awpmd.html
 fix_nve_body.html
 fix_nve_dot.html
 fix_nve_dotc_langevin.html
@@ -374,7 +376,6 @@ fix_spring_self.html
 fix_srd.html
 fix_store_force.html
 fix_store_state.html
-fix_surface_global.html
 fix_temp_berendsen.html
 fix_temp_csvr.html
 fix_temp_rescale.html
@@ -464,10 +465,12 @@ compute_pe.html
 compute_pe_atom.html
 compute_plasticity_atom.html
 compute_pressure.html
+compute_pressure_cylinder.html
 compute_pressure_uef.html
 compute_property_atom.html
 compute_property_chunk.html
 compute_property_local.html
+compute_ptm_atom.html
 compute_rdf.html
 compute_reduce.html
 compute_reduce_chunk.html
@@ -488,7 +491,7 @@ compute_smd_tlsph_shape.html
 compute_smd_tlsph_strain.html
 compute_smd_tlsph_strain_rate.html
 compute_smd_tlsph_stress.html
-compute_smd_triangle_mesh_vertices.html
+compute_smd_triangle_vertices.html
 compute_smd_ulsph_num_neighs.html
 compute_smd_ulsph_strain.html
 compute_smd_ulsph_strain_rate.html

doc/src/pair_body_nparticle.txt

@@ -6,7 +6,7 @@
 
 :line
 
-pair_style body command :h3
+pair_style body/nparticle command :h3
 
 [Syntax:]
 

doc/src/pair_born.txt

@@ -11,19 +11,14 @@ pair_style born command :h3
 pair_style born/omp command :h3
 pair_style born/gpu command :h3
 pair_style born/coul/long command :h3
-pair_style born/coul/long/cs command :h3
-pair_style born/coul/long/cs/gpu command :h3
 pair_style born/coul/long/gpu command :h3
 pair_style born/coul/long/omp command :h3
 pair_style born/coul/msm command :h3
 pair_style born/coul/msm/omp command :h3
 pair_style born/coul/wolf command :h3
-pair_style born/coul/wolf/cs command :h3
-pair_style born/coul/wolf/cs/gpu command :h3
 pair_style born/coul/wolf/gpu command :h3
 pair_style born/coul/wolf/omp command :h3
 pair_style born/coul/dsf command :h3
-pair_style born/coul/dsf/cs command :h3
 
 [Syntax:]
 
@@ -55,9 +50,7 @@ pair_coeff * * 6.08 0.317 2.340 24.18 11.51
 pair_coeff 1 1 6.08 0.317 2.340 24.18 11.51 :pre
 
 pair_style born/coul/long 10.0
-pair_style born/coul/long/cs 10.0
-pair_style born/coul/long 10.0 8.0
-pair_style born/coul/long/cs 10.0 8.0
+pair_style born/coul/long 10.0 8.
 pair_coeff * * 6.08 0.317 2.340 24.18 11.51
 pair_coeff 1 1 6.08 0.317 2.340 24.18 11.51 :pre
 
@@ -68,7 +61,6 @@ pair_coeff 1 1 6.08 0.317 2.340 24.18 11.51 :pre
 
 pair_style born/coul/wolf 0.25 10.0
 pair_style born/coul/wolf 0.25 10.0 9.0
-pair_style born/coul/wolf/cs 0.25 10.0 9.0
 pair_coeff * * 6.08 0.317 2.340 24.18 11.51
 pair_coeff 1 1 6.08 0.317 2.340 24.18 11.51 :pre
 
@@ -107,13 +99,6 @@ Wolf potential in the "coul/wolf"_pair_coul.html pair style.
 The {born/coul/dsf} style computes the Coulomb contribution with the
 damped shifted force model as in the "coul/dsf"_pair_coul.html style.
 
-Style {born/coul/long/cs} is identical to {born/coul/long} except that
-a term is added for the "core/shell model"_Howto_coreshell.html to
-allow charges on core and shell particles to be separated by r = 0.0.
-The same correction is introduced for the {born/coul/dsf/cs} style
-which is identical to {born/coul/dsf}.  And likewise for
-{born/coul/wolf/cs} style which is identical to {born/coul/wolf}.
-
 Note that these potentials are related to the "Buckingham
 potential"_pair_buck.html.
 
@@ -174,7 +159,7 @@ for the energy of the exp(), 1/r^6, and 1/r^8 portion of the pair
 interaction.
 
 The {born/coul/long} pair style supports the
-"pair_modify"_pair_modify.html table option ti tabulate the
+"pair_modify"_pair_modify.html table option to tabulate the
 short-range portion of the long-range Coulombic interaction.
 
 These styles support the pair_modify tail option for adding long-range

doc/src/pair_buck.txt

@@ -17,7 +17,6 @@ pair_style buck/coul/cut/intel command :h3
 pair_style buck/coul/cut/kk command :h3
 pair_style buck/coul/cut/omp command :h3
 pair_style buck/coul/long command :h3
-pair_style buck/coul/long/cs command :h3
 pair_style buck/coul/long/gpu command :h3
 pair_style buck/coul/long/intel command :h3
 pair_style buck/coul/long/kk command :h3
@@ -29,14 +28,14 @@ pair_style buck/coul/msm/omp command :h3
 
 pair_style style args :pre
 
-style = {buck} or {buck/coul/cut} or {buck/coul/long} or {buck/coul/long/cs} or {buck/coul/msm}
+style = {buck} or {buck/coul/cut} or {buck/coul/long} or {buck/coul/msm}
 args = list of arguments for a particular style :ul
   {buck} args = cutoff
     cutoff = global cutoff for Buckingham interactions (distance units)
   {buck/coul/cut} args = cutoff (cutoff2)
     cutoff = global cutoff for Buckingham (and Coulombic if only 1 arg) (distance units)
     cutoff2 = global cutoff for Coulombic (optional) (distance units)
-  {buck/coul/long} or {buck/coul/long/cs} args = cutoff (cutoff2)
+  {buck/coul/long} args = cutoff (cutoff2)
     cutoff = global cutoff for Buckingham (and Coulombic if only 1 arg) (distance units)
     cutoff2 = global cutoff for Coulombic (optional) (distance units)
   {buck/coul/msm} args = cutoff (cutoff2)
@@ -56,9 +55,7 @@ pair_coeff 1 1 100.0 1.5 200.0 9.0
 pair_coeff 1 1 100.0 1.5 200.0 9.0 8.0 :pre
 
 pair_style buck/coul/long 10.0
-pair_style buck/coul/long/cs 10.0
 pair_style buck/coul/long 10.0 8.0
-pair_style buck/coul/long/cs 10.0 8.0
 pair_coeff * * 100.0 1.5 200.0
 pair_coeff 1 1 100.0 1.5 200.0 9.0 :pre
 
@@ -92,10 +89,6 @@ A,C and Coulombic terms.  If two cutoffs are specified, the first is
 used as the cutoff for the A,C terms, and the second is the cutoff for
 the Coulombic term.
 
-Style {buck/coul/long/cs} is identical to {buck/coul/long} except that
-a term is added for the "core/shell model"_Howto_coreshell.html to
-allow charges on core and shell particles to be separated by r = 0.0.
-
 Note that these potentials are related to the "Born-Mayer-Huggins
 potential"_pair_born.html.
 
@@ -184,8 +177,7 @@ respa"_run_style.html command.  They do not support the {inner},
 
 [Restrictions:]
 
-The {buck/coul/long} style is part of the KSPACE package.  The
-{buck/coul/long/cs} style is part of the CORESHELL package.  They are
+The {buck/coul/long} style is part of the KSPACE package.  They are
 only enabled if LAMMPS was built with that package.  See the "Build
 package"_Build_package.html doc page for more info.
 

doc/src/pair_buck6d_coul_gauss.txt

@@ -6,8 +6,8 @@
 
 :line
 
-pair_style buck6d/coul/gauss/dsf :h3
-pair_style buck6d/coul/gauss/long :h3
+pair_style buck6d/coul/gauss/dsf command :h3
+pair_style buck6d/coul/gauss/long command :h3
 
 [Syntax:]
 

doc/src/pair_charmm.txt

@@ -8,12 +8,15 @@
 
 pair_style lj/charmm/coul/charmm command :h3
 pair_style lj/charmm/coul/charmm/intel command :h3
+pair_style lj/charmm/coul/charmm/kk command :h3
 pair_style lj/charmm/coul/charmm/omp command :h3
 pair_style lj/charmm/coul/charmm/implicit command :h3
+pair_style lj/charmm/coul/charmm/implicit/kk command :h3
 pair_style lj/charmm/coul/charmm/implicit/omp command :h3
 pair_style lj/charmm/coul/long command :h3
 pair_style lj/charmm/coul/long/gpu command :h3
 pair_style lj/charmm/coul/long/intel command :h3
+pair_style lj/charmm/coul/long/kk command :h3
 pair_style lj/charmm/coul/long/opt command :h3
 pair_style lj/charmm/coul/long/omp command :h3
 pair_style lj/charmm/coul/msm command :h3

doc/src/pair_coul_shield.txt

@@ -38,7 +38,8 @@ charge and molecule ID information is included.
 
 Where Tap(r_ij) is the taper function which provides a continuous cutoff
 (up to third derivative) for inter-atomic separations larger than r_c
-"(Maaravi)"_#Maaravi1. Here {lambda} is the shielding parameter that
+"(Leven1)"_#Leven3, "(Leven2)"_#Leven4 and "(Maaravi)"_#Maaravi1.
+Here {lambda} is the shielding parameter that
 eliminates the short-range singularity of the classical mono-polar
 electrostatic interaction expression "(Maaravi)"_#Maaravi1.
 
@@ -82,5 +83,11 @@ package"_Build_package.html doc page for more info.
 
 :line
 
+:link(Leven3)
+[(Leven1)] I. Leven, I. Azuri, L. Kronik and O. Hod, J. Chem. Phys. 140, 104106 (2014).
+
+:link(Leven4)
+[(Leven2)] I. Leven et al, J. Chem.Theory Comput. 12, 2896-905 (2016).
+
 :link(Maaravi1)
 [(Maaravi)] T. Maaravi et al, J. Phys. Chem. C 121, 22826-22835 (2017).

doc/src/pair_cs.txt

@@ -7,9 +7,11 @@
 :line
 
 pair_style born/coul/long/cs command :h3
+pair_style born/coul/long/cs/gpu command :h3
 pair_style buck/coul/long/cs command :h3
 pair_style born/coul/dsf/cs command :h3
 pair_style born/coul/wolf/cs command :h3
+pair_style born/coul/wolf/cs/gpu command :h3
 
 [Syntax:]
 
@@ -97,6 +99,38 @@ a long-range solver, thus the only correction is the addition of a
 minimal distance to avoid the possible r = 0.0 case for a
 core/shell pair.
 
+:line
+
+Styles with a {gpu}, {intel}, {kk}, {omp}, or {opt} suffix are
+functionally the same as the corresponding style without the suffix.
+They have been optimized to run faster, depending on your available
+hardware, as discussed on the "Speed packages"_Speed_packages.html doc
+page.  The accelerated styles take the same arguments and should
+produce the same results, except for round-off and precision issues.
+
+These accelerated styles are part of the GPU, USER-INTEL, KOKKOS,
+USER-OMP and OPT packages, respectively.  They are only enabled if
+LAMMPS was built with those packages.  See the "Build
+package"_Build_package.html doc page for more info.
+
+You can specify the accelerated styles explicitly in your input script
+by including their suffix, or you can use the "-suffix command-line
+switch"_Run_options.html when you invoke LAMMPS, or you can use the
+"suffix"_suffix.html command in your input script.
+
+See the "Speed packages"_Speed_packages.html doc page for more
+instructions on how to use the accelerated styles effectively.
+
+:line
+
+[Mixing, shift, table, tail correction, restart, rRESPA info]:
+
+See the corresponding doc pages for pair styles without the "cs"
+suffix to see how mixing, shifting, tabulation, tail correction,
+restarting, and rRESPA are handled by theses pair styles.
+
+:line
+
 [Restrictions:]
 
 These pair styles are part of the CORESHELL package.  They are only

doc/src/pair_dipole.txt

@@ -13,6 +13,7 @@ pair_style lj/sf/dipole/sf command :h3
 pair_style lj/sf/dipole/sf/gpu command :h3
 pair_style lj/sf/dipole/sf/omp command :h3
 pair_style lj/cut/dipole/long command :h3
+pair_style lj/cut/dipole/long/gpu command :h3
 pair_style lj/long/dipole/long command :h3
 
 [Syntax:]

doc/src/pair_eam.txt

@@ -20,6 +20,8 @@ pair_style eam/alloy/omp command :h3
 pair_style eam/alloy/opt command :h3
 pair_style eam/cd command :h3
 pair_style eam/cd/omp command :h3
+pair_style eam/cd/old command :h3
+pair_style eam/cd/old/omp command :h3
 pair_style eam/fs command :h3
 pair_style eam/fs/gpu command :h3
 pair_style eam/fs/intel command :h3
@@ -31,7 +33,7 @@ pair_style eam/fs/opt command :h3
 
 pair_style style :pre
 
-style = {eam} or {eam/alloy} or {eam/cd} or {eam/fs} :ul
+style = {eam} or {eam/alloy} or {eam/cd} or {eam/cd/old} or {eam/fs} :ul
 
 [Examples:]
 
@@ -268,7 +270,8 @@ Style {eam/cd} is similar to the {eam/alloy} style, except that it
 computes alloy pairwise interactions using the concentration-dependent
 embedded-atom method (CD-EAM).  This model can reproduce the enthalpy
 of mixing of alloys over the full composition range, as described in
-"(Stukowski)"_#Stukowski.
+"(Stukowski)"_#Stukowski. Style {eam/cd/old} is an older, slightly
+different and slower two-site formulation of the model "(Caro)"_#Caro.
 
 The pair_coeff command is specified the same as for the {eam/alloy}
 style.  However the DYNAMO {setfl} file must has two
@@ -442,3 +445,6 @@ Daw, Baskes, Phys Rev B, 29, 6443 (1984).
 :link(Stukowski)
 [(Stukowski)] Stukowski, Sadigh, Erhart, Caro; Modeling Simulation
 Materials Science & Engineering, 7, 075005 (2009).
+
+:link(Caro)
+[(Caro)] A Caro, DA Crowson, M Caro; Phys Rev Lett, 95, 075702 (2005)

doc/src/pair_edip.txt

@@ -7,6 +7,7 @@
 :line
 
 pair_style edip command :h3
+pair_style edip/omp command :h3
 pair_style edip/multi command :h3
 
 [Syntax:]

doc/src/pair_gran.txt

@@ -7,7 +7,7 @@
 :line
 
 pair_style gran/hooke command :h3
-pair_style gran/omp command :h3
+pair_style gran/hooke/omp command :h3
 pair_style gran/hooke/history command :h3
 pair_style gran/hooke/history/omp command :h3
 pair_style gran/hertz/history command :h3

doc/src/pair_gromacs.txt

@@ -8,8 +8,10 @@
 
 pair_style lj/gromacs command :h3
 pair_style lj/gromacs/gpu command :h3
+pair_style lj/gromacs/kk command :h3
 pair_style lj/gromacs/omp command :h3
 pair_style lj/gromacs/coul/gromacs command :h3
+pair_style lj/gromacs/coul/gromacs/kk command :h3
 pair_style lj/gromacs/coul/gromacs/omp command :h3
 
 [Syntax:]

doc/src/pair_hybrid.txt

@@ -7,9 +7,8 @@
 :line
 
 pair_style hybrid command :h3
-pair_style hybrid/omp command :h3
+pair_style hybrid/kk command :h3
 pair_style hybrid/overlay command :h3
-pair_style hybrid/overlay/omp command :h3
 pair_style hybrid/overlay/kk command :h3
 
 [Syntax:]

doc/src/pair_ilp_graphene_hbn.txt

@@ -25,22 +25,24 @@ pair_coeff  * * rebo                 CH.airebo   NULL NULL C
 pair_coeff  * * tersoff              BNC.tersoff B    N    NULL
 pair_coeff  * * ilp/graphene/hbn     BNCH.ILP    B    N    C
 pair_coeff  1 1 coul/shield 0.70
-pair_coeff  1 2 coul/shield 0.69498201415576216335
+pair_coeff  1 2 coul/shield 0.695
 pair_coeff  2 2 coul/shield 0.69 :pre
 
 [Description:]
 
 The {ilp/graphene/hbn} style computes the registry-dependent interlayer
-potential (ILP) potential as described in "(Leven)"_#Leven and
-"(Maaravi)"_#Maaravi2. The normals are calculated in the way as described
+potential (ILP) potential as described in "(Leven1)"_#Leven1,
+"(Leven2)"_#Leven2 and "(Maaravi)"_#Maaravi2.
+The normals are calculated in the way as described
 in "(Kolmogorov)"_#Kolmogorov2.
 
 :c,image(Eqs/pair_ilp_graphene_hbn.jpg)
 
 Where Tap(r_ij) is the taper function which provides a continuous
 cutoff (up to third derivative) for interatomic separations larger than
-r_c "(Maaravi)"_#Maaravi2. The definitions of each parameter in the above
-equation can be found in "(Leven)"_#Leven and "(Maaravi)"_#Maaravi2.
+r_c "(Maaravi)"_#Maaravi2. The definitons of each parameter in the above
+equation can be found in "(Leven1)"_#Leven1 and "(Maaravi)"_#Maaravi2.
+
 
 It is important to include all the pairs to build the neighbor list for
 calculating the normals.
@@ -61,13 +63,15 @@ NOTE: The parameters presented in the parameter file (e.g. BNCH.ILP),
 are fitted with taper function by setting the cutoff equal to 16.0
 Angstrom.  Using different cutoff or taper function should be careful.
 
-NOTE: Two new sets of parameters of ILP for two-dimensional hexagonal Materials
-are presented in "(Ouyang)"_#Ouyang1. These parameters provide a good description
-in both short- and long-range interaction regime, while the old ILP parameters
-published in "(Leven)"_#Leven and "(Maaravi)"_#Maaravi2 are only suitable for
-long-range interaction regime. This feature is essential for simulations in
-high-pressure regime (i.e., the interlayer distance smaller than the equilibrium distance).
-The benchmark tests and comparison of these parameters can be found in "(Ouyang)"_#Ouyang1.
+NOTE: Two new sets of parameters of ILP for two-dimensional hexagonal
+Materials are presented in "(Ouyang)"_#Ouyang.  These parameters provide
+a good description in both short- and long-range interaction regimes.
+While the old ILP parameters published in "(Leven2)"_#Leven2 and
+"(Maaravi)"_#Maaravi2 are only suitable for long-range interaction
+regime. This feature is essential for simulations in high pressure
+regime (i.e., the interlayer distance is smaller than the equilibrium
+distance).  The benchmark tests and comparison of these parameters can
+be found in "(Ouyang)"_#Ouyang.
 
 This potential must be used in combination with hybrid/overlay.
 Other interactions can be set to zero using pair_style {none}.
@@ -112,14 +116,17 @@ units, if your simulation does not use {metal} units.
 
 :line
 
-:link(Leven)
-[(Leven)] I. Leven et al, J. Chem.Theory Comput. 12, 2896-905 (2016)
+:link(Leven1)
+[(Leven1)] I. Leven, I. Azuri, L. Kronik and O. Hod, J. Chem. Phys. 140, 104106 (2014).
+
+:link(Leven2)
+[(Leven2)] I. Leven et al, J. Chem.Theory Comput. 12, 2896-905 (2016).
 
 :link(Maaravi2)
 [(Maaravi)] T. Maaravi et al, J. Phys. Chem. C 121, 22826-22835 (2017).
 
 :link(Kolmogorov2)
-[(Kolmogorov)] A. N. Kolmogorov, V. H. Crespi, Phys. Rev. B 71, 235415 (2005)
+[(Kolmogorov)] A. N. Kolmogorov, V. H. Crespi, Phys. Rev. B 71, 235415 (2005).
 
-:link(Ouyang1)
-[(Ouyang)] W. Ouyang, D. Mandelli,  M. Urbakh, O. Hod, arXiv:1806.09555 (2018).
+:link(Ouyang)
+[(Ouyang)] W. Ouyang, D. Mandelli, M. Urbakh and O. Hod, Nano Lett. 18, 6009-6016 (2018).

doc/src/pair_kolmogorov_crespi_full.txt

@@ -19,11 +19,11 @@ tap_flag = 0/1 to turn off/on the taper function
 
 pair_style hybrid/overlay kolmogorov/crespi/full 20.0 0
 pair_coeff * * none
-pair_coeff * * kolmogorov/crespi/full  CC.KC   C C :pre
+pair_coeff * * kolmogorov/crespi/full  CH.KC   C C :pre
 
-pair_style hybrid/overlay rebo kolmogorov/crespi/full 16.0
-pair_coeff * * rebo                    CH.airebo  C C
-pair_coeff * * kolmogorov/crespi/full  CC.KC      C C :pre
+pair_style hybrid/overlay rebo kolmogorov/crespi/full 16.0 1
+pair_coeff * * rebo                    CH.airebo    C H
+pair_coeff * * kolmogorov/crespi/full  CH_taper.KC  C H :pre
 
 [Description:]
 
@@ -38,27 +38,32 @@ forces and to include all the pairs to build the neighbor list for
 calculating the normals.  Energies are shifted so that they go
 continuously to zero at the cutoff assuming that the exponential part of
 {Vij} (first term) decays sufficiently fast.  This shift is achieved by
-the last term in the equation for {Vij} above.
+the last term in the equation for {Vij} above. This is essential only
+when the tapper function is turned off. The formula of taper function
+can be found in pair style "ilp/graphene/hbn"_pair_ilp_graphene_hbn.html.
 
 NOTE: This potential is intended for interactions between two different
 graphene layers. Therefore, to avoid interaction within the same layers,
 each layer should have a separate molecule id and is recommended to use
 "full" atom style in the data file.
 
-The parameter file (e.g. CC.KC), is intended for use with {metal}
+The parameter file (e.g. CH.KC), is intended for use with {metal}
 "units"_units.html, with energies in meV. Two additional parameters, {S},
 and {rcut} are included in the parameter file. {S} is designed to
 facilitate scaling of energies. {rcut} is designed to build the neighbor
 list for calculating the normals for each atom pair.
 
-NOTE: A new set of parameters of KC potential for hydrocarbons (CH.KC)
-is presented in "(Ouyang)"_#Ouyang2. The parameters in CH.KC provides
-a good description in both short- and long-range interaction regime,
-while the original parameters (CC.KC) published in "(Kolmogorov)"_#Kolmogorov1
-are only suitable for long-range interaction regime. 
-This feature is essential for simulations in high-pressure regime
-(i.e., the interlayer distance smaller than the equilibrium distance).
-The benchmark tests and comparison of these parameters can be found in "(Ouyang)"_#Ouyang2.
+NOTE: Two new sets of parameters of KC potential for hydrocarbons, CH.KC
+(without the taper function) and CH_taper.KC (with the taper function)
+are presented in "(Ouyang)"_#Ouyang1.  The energy for the KC potential
+with the taper function goes continuously to zero at the cutoff.  The
+parameters in both CH.KC and CH_taper.KC provide a good description in
+both short- and long-range interaction regimes. While the original
+parameters (CC.KC) published in "(Kolmogorov)"_#Kolmogorov1 are only
+suitable for long-range interaction regime.  This feature is essential
+for simulations in high pressure regime (i.e., the interlayer distance
+is smaller than the equilibrium distance).  The benchmark tests and
+comparison of these parameters can be found in "(Ouyang)"_#Ouyang1.
 
 This potential must be used in combination with hybrid/overlay.
 Other interactions can be set to zero using pair_style {none}.
@@ -84,7 +89,7 @@ package"_Build_package.html doc page for more info.
 This pair potential requires the newton setting to be {on} for pair
 interactions.
 
-The CC.KC potential file provided with LAMMPS (see the potentials
+The CH.KC potential file provided with LAMMPS (see the potentials
 folder) are parameterized for metal units.  You can use this potential
 with any LAMMPS units, but you would need to create your own custom
 CC.KC potential file with all coefficients converted to the appropriate
@@ -105,5 +110,5 @@ units.
 :link(Kolmogorov1)
 [(Kolmogorov)] A. N. Kolmogorov, V. H. Crespi, Phys. Rev. B 71, 235415 (2005)
 
-:link(Ouyang2)
-[(Ouyang)] W. Ouyang, D. Mandelli,  M. Urbakh, O. Hod, arXiv:1806.09555 (2018).
+:link(Ouyang1)
+[(Ouyang)] W. Ouyang, D. Mandelli, M. Urbakh and O. Hod, Nano Lett. 18, 6009-6016 (2018).

doc/src/pair_lj.txt

@@ -14,6 +14,7 @@ pair_style lj/cut/opt command :h3
 pair_style lj/cut/omp command :h3
 pair_style lj/cut/coul/cut command :h3
 pair_style lj/cut/coul/cut/gpu command :h3
+pair_style lj/cut/coul/cut/kk command :h3
 pair_style lj/cut/coul/cut/omp command :h3
 pair_style lj/cut/coul/debye command :h3
 pair_style lj/cut/coul/debye/gpu command :h3
@@ -26,6 +27,7 @@ pair_style lj/cut/coul/dsf/omp command :h3
 pair_style lj/cut/coul/long command :h3
 pair_style lj/cut/coul/long/cs command :h3
 pair_style lj/cut/coul/long/gpu command :h3
+pair_style lj/cut/coul/long/kk command :h3
 pair_style lj/cut/coul/long/intel command :h3
 pair_style lj/cut/coul/long/opt command :h3
 pair_style lj/cut/coul/long/omp command :h3

doc/src/pair_lj_expand.txt

@@ -8,7 +8,10 @@
 
 pair_style lj/expand command :h3
 pair_style lj/expand/gpu command :h3
+pair_style lj/expand/kk command :h3
 pair_style lj/expand/omp command :h3
+pair_style lj/expand/coul/long command :h3
+pair_style lj/expand/coul/long/gpu command :h3
 
 [Syntax:]
 
@@ -22,6 +25,11 @@ pair_style lj/expand 2.5
 pair_coeff * * 1.0 1.0 0.5
 pair_coeff 1 1 1.0 1.0 -0.2 2.0 :pre
 
+pair_style lj/expand/coul/long 2.5
+pair_style lj/expand/coul/long 2.5 4.0
+pair_coeff * * 1.0 1.0 0.5
+pair_coeff 1 1 1.0 1.0 -0.2 3.0 :pre
+
 [Description:]
 
 Style {lj/expand} computes a LJ interaction with a distance shifted by
@@ -34,11 +42,12 @@ formula:
 Rc is the cutoff which does not include the delta distance.  I.e. the
 actual force cutoff is the sum of cutoff + delta.
 
-The following coefficients must be defined for each pair of atoms
-types via the "pair_coeff"_pair_coeff.html command as in the examples
-above, or in the data file or restart files read by the
-"read_data"_read_data.html or "read_restart"_read_restart.html
-commands, or by mixing as described below:
+For all of the {lj/expand} pair styles, the following coefficients must
+be defined for each pair of atoms types via the
+"pair_coeff"_pair_coeff.html command as in the examples above, or in
+the data file or restart files read by the "read_data"_read_data.html
+or "read_restart"_read_restart.html commands, or by mixing as
+described below:
 
 epsilon (energy units)
 sigma (distance units)
@@ -48,6 +57,11 @@ cutoff (distance units) :ul
 The delta values can be positive or negative.  The last coefficient is
 optional.  If not specified, the global LJ cutoff is used.
 
+For {lj/expand/coul/long} only the LJ cutoff can be specified since a
+Coulombic cutoff cannot be specified for an individual I,J type pair.
+All type pairs use the same global Coulombic cutoff specified in the
+pair_style command.
+
 :line
 
 Styles with a {gpu}, {intel}, {kk}, {omp}, or {opt} suffix are

doc/src/pair_lj_long.txt

@@ -11,6 +11,7 @@ pair_style lj/long/coul/long/intel command :h3
 pair_style lj/long/coul/long/omp command :h3
 pair_style lj/long/coul/long/opt command :h3
 pair_style lj/long/tip4p/long command :h3
+pair_style lj/long/tip4p/long/omp command :h3
 
 [Syntax:]
 

doc/src/pair_meam_spline.txt

@@ -6,8 +6,8 @@
 
 :line
 
-pair_style meam/spline :h3
-pair_style meam/spline/omp :h3
+pair_style meam/spline command :h3
+pair_style meam/spline/omp command :h3
 
 [Syntax:]
 

doc/src/pair_meam_sw_spline.txt

@@ -6,8 +6,7 @@
 
 :line
 
-pair_style meam/sw/spline :h3
-pair_style meam/sw/spline/omp :h3
+pair_style meam/sw/spline command :h3
 
 [Syntax:]
 

doc/src/pair_nb3b_harmonic.txt

@@ -7,7 +7,6 @@
 :line
 
 pair_style nb3b/harmonic command :h3
-pair_style nb3b/harmonic/omp command :h3
 
 [Syntax:]
 
@@ -89,28 +88,6 @@ a particular simulation; LAMMPS ignores those entries.
 
 :line
 
-Styles with a {gpu}, {intel}, {kk}, {omp}, or {opt} suffix are
-functionally the same as the corresponding style without the suffix.
-They have been optimized to run faster, depending on your available
-hardware, as discussed on the "Speed packages"_Speed_packages.html doc
-page.  The accelerated styles take the same arguments and should
-produce the same results, except for round-off and precision issues.
-
-These accelerated styles are part of the GPU, USER-INTEL, KOKKOS,
-USER-OMP and OPT packages, respectively.  They are only enabled if
-LAMMPS was built with those packages.  See the "Build
-package"_Build_package.html doc page for more info.
-
-You can specify the accelerated styles explicitly in your input script
-by including their suffix, or you can use the "-suffix command-line
-switch"_Run_options.html when you invoke LAMMPS, or you can use the
-"suffix"_suffix.html command in your input script.
-
-See the "Speed packages"_Speed_packages.html doc page for more
-instructions on how to use the accelerated styles effectively.
-
-:line
-
 [Restrictions:]
 
 This pair style can only be used if LAMMPS was built with the MANYBODY

doc/src/pair_sdk.txt

@@ -13,6 +13,8 @@ pair_style lj/sdk/omp command :h3
 pair_style lj/sdk/coul/long command :h3
 pair_style lj/sdk/coul/long/gpu command :h3
 pair_style lj/sdk/coul/long/omp command :h3
+pair_style lj/sdk/coul/msm command :h3
+pair_style lj/sdk/coul/msm/omp command :h3
 
 [Syntax:]
 
@@ -35,6 +37,10 @@ pair_style lj/sdk/coul/long 10.0
 pair_style lj/sdk/coul/long 10.0 12.0
 pair_coeff 1 1 lj9_6 100.0 3.5 12.0 :pre
 
+pair_style lj/sdk/coul/msm 10.0
+pair_style lj/sdk/coul/msm 10.0 12.0
+pair_coeff 1 1 lj9_6 100.0 3.5 12.0 :pre
+
 [Description:]
 
 The {lj/sdk} styles compute a 9/6, 12/4, or 12/6 Lennard-Jones potential,
@@ -75,10 +81,10 @@ and Coulombic interactions for this type pair.  If both coefficients
 are specified, they are used as the LJ and Coulombic cutoffs for this
 type pair.
 
-For {lj/sdk/coul/long} only the LJ cutoff can be specified since a
-Coulombic cutoff cannot be specified for an individual I,J type pair.
-All type pairs use the same global Coulombic cutoff specified in the
-pair_style command.
+For {lj/sdk/coul/long} and {lj/sdk/coul/msm} only the LJ cutoff can be
+specified since a Coulombic cutoff cannot be specified for an
+individual I,J type pair.  All type pairs use the same global
+Coulombic cutoff specified in the pair_style command.
 
 :line
 

doc/src/pair_spin_magelec.txt

@@ -6,11 +6,11 @@
 
 :line
 
-pair_style spin/me command :h3
+pair_style spin/magelec command :h3
 
 [Syntax:]
 
-pair_style spin/me cutoff :pre
+pair_style spin/magelec cutoff :pre
 
 cutoff = global cutoff pair (distance in metal units) :ulb,l
 
@@ -18,8 +18,8 @@ cutoff = global cutoff pair (distance in metal units) :ulb,l
 
 [Examples:]
 
-pair_style spin/me 4.5
-pair_coeff * * me 4.5 0.00109 1.0 1.0 1.0 :pre
+pair_style spin/magelec 4.5
+pair_coeff * * magelec 4.5 0.00109 1.0 1.0 1.0 :pre
 
 [Description:]
 

doc/src/pair_spin_neel.txt

@@ -29,34 +29,36 @@ between pairs of magnetic spins:
 
 :c,image(Eqs/pair_spin_neel_interaction.jpg)
 
-where si and sj are two neighboring magnetic spins of two particles, 
+where si and sj are two neighboring magnetic spins of two particles,
 rij = ri - rj is the inter-atomic distance between the two particles,
-eij = (ri - rj)/|ri-rj| is their normalized separation vector 
-and g1, q1 and q2 are three functions defining the intensity of the
-dipolar and quadrupolar contributions, with:
+eij = (ri - rj)/|ri-rj| is their normalized separation vector and g1,
+q1 and q2 are three functions defining the intensity of the dipolar
+and quadrupolar contributions, with:
 
 :c,image(Eqs/pair_spin_neel_functions.jpg)
 
-With the functions g(rij) and q(rij) defined and fitted according to the same
-Bethe-Slater function used to fit the exchange interaction: 
+With the functions g(rij) and q(rij) defined and fitted according to
+the same Bethe-Slater function used to fit the exchange interaction:
 
 :c,image(Eqs/pair_spin_exchange_function.jpg)
 
-where a, b and d are the three constant coefficients defined in the associated 
-"pair_coeff" command.   
+where a, b and d are the three constant coefficients defined in the
+associated "pair_coeff" command.
 
-The coefficients a, b, and d need to be fitted so that the function above matches with 
-the values of the magneto-elastic constant of the materials at stake. 
+The coefficients a, b, and d need to be fitted so that the function
+above matches with the values of the magneto-elastic constant of the
+materials at stake.
 
-Examples and more explanations about this function and its parametrization are reported 
-in "(Tranchida)"_#Tranchida6. More examples of parametrization will be provided in 
-future work.
+Examples and more explanations about this function and its
+parametrization are reported in "(Tranchida)"_#Tranchida6. More
+examples of parametrization will be provided in future work.
 
-From this DM interaction, each spin i will be submitted to a magnetic torque
-omega and its associated atom to a force F (for spin-lattice calculations only). 
+From this DM interaction, each spin i will be submitted to a magnetic
+torque omega and its associated atom to a force F (for spin-lattice
+calculations only).
 
-More details about the derivation of these torques/forces are reported in
-"(Tranchida)"_#Tranchida6.
+More details about the derivation of these torques/forces are reported
+in "(Tranchida)"_#Tranchida6.
 
 :line
 

doc/src/pair_tersoff.txt

@@ -8,10 +8,10 @@
 
 pair_style tersoff command :h3
 pair_style tersoff/table command :h3
-pair_style tersoff/gpu :h3
-pair_style tersoff/intel :h3
-pair_style tersoff/kk :h3
-pair_style tersoff/omp :h3
+pair_style tersoff/gpu command :h3
+pair_style tersoff/intel command :h3
+pair_style tersoff/kk command :h3
+pair_style tersoff/omp command :h3
 pair_style tersoff/table/omp command :h3
 
 [Syntax:]

doc/src/special_bonds.txt

@@ -13,7 +13,7 @@ special_bonds command :h3
 special_bonds keyword values ... :pre
 
 one or more keyword/value pairs may be appended :ulb,l
-keyword = {amber} or {charmm} or {dreiding} or {fene} or {lj/coul} or {lj} or {coul} or {angle} or {dihedral} or {extra} :l
+keyword = {amber} or {charmm} or {dreiding} or {fene} or {lj/coul} or {lj} or {coul} or {angle} or {dihedral} :l
   {amber} values = none
   {charmm} values = none
   {dreiding} values = none
@@ -219,7 +219,7 @@ each time the command is issued.
 [Default:]
 
 All 3 Lennard-Jones and 3 Coulombic weighting coefficients = 0.0,
-angle = no, dihedral = no, and extra = 0.
+angle = no, dihedral = no.
 
 :line
 

examples/SPIN/bfo/in.spin.bfo

@@ -36,7 +36,7 @@ fix 		3 all nve/spin lattice no
 
 timestep	0.0002
 
-compute 	out_mag	 all compute/spin
+compute 	out_mag	 all spin
 compute 	out_pe	 all pe
 compute 	out_ke	 all ke
 compute 	out_temp all temp

examples/SPIN/cobalt_fcc/in.spin.cobalt_fcc

@@ -40,7 +40,7 @@ timestep	0.0001
 
 # compute and output options
 
-compute 	out_mag    all compute/spin
+compute 	out_mag    all spin
 compute 	out_pe     all pe
 compute 	out_ke     all ke
 compute 	out_temp   all temp

examples/SPIN/cobalt_hcp/in.spin.cobalt_hcp

@@ -40,7 +40,7 @@ fix 		3 all nve/spin lattice yes
 timestep	0.0001
 
 
-compute 	out_mag    all compute/spin
+compute 	out_mag    all spin
 compute 	out_pe     all pe
 compute 	out_ke     all ke
 compute 	out_temp   all temp

examples/SPIN/iron/in.spin.iron

@@ -37,7 +37,7 @@ timestep	0.0001
 
 # compute and output options
 
-compute 	out_mag    all compute/spin
+compute 	out_mag    all spin
 compute 	out_pe     all pe
 compute 	out_ke     all ke
 compute 	out_temp   all temp

examples/SPIN/nickel/in.spin.nickel

@@ -38,7 +38,7 @@ timestep	0.0001
 
 # compute and output options
 
-compute 	out_mag    all compute/spin
+compute 	out_mag    all spin
 compute 	out_pe     all pe
 compute 	out_ke     all ke
 compute 	out_temp   all temp

examples/SPIN/read_restart/in.spin.read_data

@@ -25,7 +25,7 @@ timestep	0.0001
 
 # define outputs and computes
 
-compute 	out_mag    all compute/spin
+compute 	out_mag    all spin
 compute 	out_pe     all pe
 compute 	out_ke     all ke
 compute 	out_temp   all temp

examples/SPIN/read_restart/in.spin.restart

@@ -29,7 +29,7 @@ timestep	0.0001
 
 # define outputs
 
-compute 	out_mag    all compute/spin
+compute 	out_mag    all spin
 compute 	out_pe     all pe
 compute 	out_ke     all ke
 compute 	out_temp   all temp

examples/SPIN/read_restart/in.spin.write_restart

@@ -34,7 +34,7 @@ timestep	0.0001
 
 # compute and output options
 
-compute 	out_mag    all compute/spin
+compute 	out_mag    all spin
 compute 	out_pe     all pe
 compute 	out_ke     all ke
 compute 	out_temp   all temp

examples/USER/smd/fluid_structure_interaction/fluid_structure_interaction.lmp

@@ -129,7 +129,7 @@ compute         contact_radius all smd/contact/radius
 compute         S solids smd/tlsph/stress
 compute         nn water smd/ulsph/num/neighs
 compute         epl solids smd/plastic/strain
-compute         vol all smd/volume
+compute         vol all smd/vol
 compute         rho all smd/rho
 
 dump            dump_id all custom 100 dump.LAMMPS id type x y &

potentials/BNCH.ILP

@@ -1,6 +1,6 @@
 # Interlayer Potential (ILP) for graphene/graphene, graphene/hBN and hBN/hBN junctions
 #
-# Cite as Wengen Ouyang,Davide Mandelli, Michael Urbakh, Oded Hod, arXiv:1806.09555 (2018).
+# Cite as W. Ouyang, D. Mandelli, M. Urbakh and O. Hod, Nano Letters 18, 6009-6016 (2018).
 #
 # ----------------- Repulsion Potential ------------------++++++++++++++ Vdw Potential ++++++++++++++++************  
 #    beta(A)    alpha    delta(A)  epsilon(meV)  C(meV)        d         sR      reff(A)   C6(meV*A^6)    S    rcut

potentials/CH.KC

@@ -1,6 +1,6 @@
-# Refined parameters for Kolmogorov-Crespi Potential
+# Refined parameters for Kolmogorov-Crespi Potential without taper function
 #
-# Cite as Wengen Ouyang,Davide Mandelli, Michael Urbakh, Oded Hod, arXiv:1806.09555 (2018).
+# Cite as W. Ouyang, D. Mandelli, M. Urbakh and O. Hod, Nano Letters 18, 6009-6016 (2018).
 # 
 #        z0         C0             C2          C4         C          delta     lambda       A          S    rcut
 C  C  3.328819  21.847167     12.060173     4.711099  6.678908e-4  0.7718101  3.143921  12.660270     1.0   2.0

potentials/CH_taper.KC

@@ -0,0 +1,10 @@
+# Refined parameters for Kolmogorov-Crespi Potential with taper function
+#
+# Cite as W. Ouyang, D. Mandelli, M. Urbakh and O. Hod, Nano Letters 18, 6009-6016 (2018).
+# 
+#        z0         C0          C2             C4            C            delta        lambda       A           S    rcut
+C C  3.416084	20.021583    10.9055107	    4.2756354	1.0010836E-2	0.8447122    2.9360584	 14.3132588    1.0    2.0
+C H  2.849054	72.557245    1.0164169E-2   65.923312	8.7962504E-5	0.3349237    3.0402632	 14.7533201    1.0    1.5
+H H  2.187478	3.915802E-5  5.0896431E-5   3.6657827	1.5373722446	0.9633581    0.4249989	 1.570737E-4   1.0    1.2
+H C  2.849054	72.557245    1.0164169E-2   65.923312	8.7962504E-5	0.3349237    3.0402632	 14.7533201    1.0    2.2
+

src/.gitignore

@@ -41,6 +41,11 @@
 /pair_meamc.cpp
 /pair_meamc.h
 
+/ptm_*.cpp
+/ptm_*.h
+/compute_ptm.cpp
+/compute_ptm.h
+
 /fix_qeq*.cpp
 /fix_qeq*.h
 
@@ -296,6 +301,8 @@
 /compute_plasticity_atom.h
 /compute_pressure_bocs.cpp
 /compute_pressure_bocs.h
+/compute_pressure_cylinder.cpp
+/compute_pressure_cylinder.h
 /compute_pressure_grem.cpp
 /compute_pressure_grem.h
 /compute_rigid_local.cpp
@@ -424,6 +431,8 @@
 /fix_eos_table.h
 /fix_evaporate.cpp
 /fix_evaporate.h
+/fix_ffl.cpp
+/fix_ffl.h
 /fix_filter_corotate.cpp
 /fix_filter_corotate.h
 /fix_viscosity.cpp

src/KOKKOS/pair_eam_alloy_kokkos.cpp

@@ -165,6 +165,9 @@ void PairEAMAlloyKokkos<DeviceType>::compute(int eflag_in, int vflag_in)
       }
     }
 
+    if (need_dup)
+      Kokkos::Experimental::contribute(d_rho, dup_rho);
+
     // communicate and sum densities (on the host)
 
     if (newton_pair) {

src/KOKKOS/pair_eam_fs_kokkos.cpp

@@ -165,6 +165,9 @@ void PairEAMFSKokkos<DeviceType>::compute(int eflag_in, int vflag_in)
       }
     }
 
+    if (need_dup)
+      Kokkos::Experimental::contribute(d_rho, dup_rho);
+
     // communicate and sum densities (on the host)
 
     if (newton_pair) {

src/KOKKOS/pair_eam_kokkos.cpp

@@ -159,6 +159,9 @@ void PairEAMKokkos<DeviceType>::compute(int eflag_in, int vflag_in)
       }
     }
 
+    if (need_dup)
+      Kokkos::Experimental::contribute(d_rho, dup_rho);
+
     // communicate and sum densities (on the host)
 
     if (newton_pair) {

src/MC/fix_atom_swap.cpp

@@ -377,7 +377,7 @@ int FixAtomSwap::attempt_semi_grand()
   int success = 0;
   if (i >= 0)
     if (random_unequal->uniform() <
-      exp(-beta*(energy_after - energy_before
+      exp(beta*(energy_before - energy_after
             + mu[jtype] - mu[itype]))) success = 1;
 
   int success_all = 0;

src/MEAM/pair_meam.cpp

@@ -48,8 +48,8 @@ static const char *keywords[] = {
 PairMEAM::PairMEAM(LAMMPS *lmp) : Pair(lmp)
 {
   if (comm->me == 0)
-    error->warning(FLERR,"The pair_style meam command is unsupported. "
-                   "Please use pair_style meam/c instead");
+    error->warning(FLERR,"THE pair_style meam COMMAND IS OBSOLETE AND "
+                   "WILL BE REMOVED VERY SOON. PLEASE USE meam/c");
 
   single_enable = 0;
   restartinfo = 0;

src/Makefile

@@ -61,7 +61,7 @@ PACKUSER = user-atc user-awpmd user-bocs user-cgdna user-cgsdk user-colvars \
 	   user-diffraction user-dpd user-drude user-eff user-fep user-h5md \
 	   user-intel user-lb user-manifold user-meamc user-meso \
 	   user-mgpt user-misc user-mofff user-molfile \
-	   user-netcdf user-omp user-phonon user-qmmm user-qtb \
+	   user-netcdf user-omp user-phonon user-ptm user-qmmm user-qtb \
 	   user-quip user-reaxc user-scafacos user-smd user-smtbq \
 	   user-sph user-tally user-uef user-vtk
 

src/PYTHON/fix_python_invoke.h

@@ -13,8 +13,8 @@
 
 #ifdef FIX_CLASS
 
-FixStyle(python,FixPythonInvoke)
 FixStyle(python/invoke,FixPythonInvoke)
+FixStyle(python,FixPythonInvoke)
 
 #else
 

src/Purge.list

@@ -24,6 +24,9 @@ style_nstencil.h
 style_ntopo.h
 # other auto-generated files
 lmpinstalledpkgs.h
+# renamed on 25 September 2018
+compute_smd_triangle_mesh_vertices.h
+compute_smd_triangle_mesh_vertices.cpp
 # renamed on 6 September 2018
 pair_cdeam_omp.h
 pair_cdeam_omp.cpp

src/REAX/pair_reax.cpp

@@ -46,8 +46,8 @@ using namespace LAMMPS_NS;
 PairREAX::PairREAX(LAMMPS *lmp) : Pair(lmp)
 {
   if (comm->me == 0)
-    error->warning(FLERR,"The pair_style reax command is unsupported. "
-                   "Please switch to pair_style reax/c instead");
+    error->warning(FLERR,"THE pair_style reax COMMAND IS OBSOLETE AND "
+                   "WILL BE REMOVED VERY SOON. PLEASE USE reax/c");
 
   single_enable = 0;
   restartinfo = 0;

src/SPIN/compute_spin.h

@@ -13,7 +13,7 @@
 
 #ifdef COMPUTE_CLASS
 
-ComputeStyle(compute/spin,ComputeSpin)
+ComputeStyle(spin,ComputeSpin)
 
 #else
 

src/SPIN/pair_spin.h

@@ -11,12 +11,6 @@
    See the README file in the top-level LAMMPS directory.
 ------------------------------------------------------------------------- */
 
-#ifdef PAIR_CLASS
-
-PairStyle(pair/spin,PairSpin)
-
-#else
-
 #ifndef LMP_PAIR_SPIN_H
 #define LMP_PAIR_SPIN_H
 
@@ -46,7 +40,6 @@ friend class FixNVESpin;
 
 }
 
-#endif
 #endif
 
 /* ERROR/WARNING messages:

src/USER-INTEL/fix_intel.cpp

@@ -532,8 +532,13 @@ void FixIntel::check_neighbor_intel()
       _offload_noghost = 0;
     }
     #endif
+
+    // avoid flagging a neighbor list as both USER-INTEL and USER-OMP
+    if (neighbor->requests[i]->intel)
+      neighbor->requests[i]->omp = 0;
+
     if (neighbor->requests[i]->skip)
-      error->all(FLERR, "Cannot yet use hybrid styles with Intel package.");
+      error->all(FLERR, "Hybrid styles with Intel package are unsupported.");
   }
 }
 

src/USER-MEAMC/meam_dens_init.cpp

@@ -101,7 +101,7 @@ MEAM::getscreen(int i, double* scrfcn, double* dscrfcn, double* fcpair, double**
   double xjtmp, yjtmp, zjtmp, delxik, delyik, delzik, rik2 /*,rik*/;
   double xktmp, yktmp, zktmp, delxjk, delyjk, delzjk, rjk2 /*,rjk*/;
   double xik, xjk, sij, fcij, sfcij, dfcij, sikj, dfikj, cikj;
-  double Cmin, Cmax, delc, /*ebound,*/ rbound, a, coef1, coef2;
+  double Cmin, Cmax, delc, /*ebound,*/ a, coef1, coef2;
   double dCikj;
   double rnorm, fc, dfc, drinv;
 
@@ -129,6 +129,7 @@ MEAM::getscreen(int i, double* scrfcn, double* dscrfcn, double* fcpair, double**
     rij2 = delxij * delxij + delyij * delyij + delzij * delzij;
     rij = sqrt(rij2);
 
+    const double rbound = this->ebound_meam[elti][eltj] * rij2;
     if (rij > this->rc_meam) {
       fcij = 0.0;
       dfcij = 0.0;
@@ -138,7 +139,6 @@ MEAM::getscreen(int i, double* scrfcn, double* dscrfcn, double* fcpair, double**
       sij = 1.0;
 
       //     if rjk2 > ebound*rijsq, atom k is definitely outside the ellipse
-      const double rbound = this->ebound_meam[elti][eltj] * rij2;
       for (kn = 0; kn < numneigh_full; kn++) {
         k = firstneigh_full[kn];
         eltk = fmap[type[k]];
@@ -193,7 +193,6 @@ MEAM::getscreen(int i, double* scrfcn, double* dscrfcn, double* fcpair, double**
     if (iszero(sfcij) || iszero(sfcij - 1.0))
       goto LABEL_100;
 
-    rbound = this->ebound_meam[elti][eltj] * rij2;
     for (kn = 0; kn < numneigh_full; kn++) {
       k = firstneigh_full[kn];
       if (k == j) continue;

src/USER-MEAMC/meam_funcs.cpp

@@ -87,7 +87,7 @@ MEAM::dG_gam(const double gamma, const int ibar, double& dG) const
         //         e.g. gsmooth_factor is 99, {:
         //         gsmooth_switchpoint = -0.99
         //         G = 0.01*(-0.99/gamma)**99
-        double G = 1 / (gsmooth_factor + 1) * pow((gsmooth_switchpoint / gamma), gsmooth_factor);
+        G = 1 / (gsmooth_factor + 1) * pow((gsmooth_switchpoint / gamma), gsmooth_factor);
         G = sqrt(G);
         dG = -gsmooth_factor * G / (2.0 * gamma);
         return G;

src/USER-MEAMC/meam_impl.cpp

@@ -34,6 +34,19 @@ MEAM::MEAM(Memory* mem)
 
   maxneigh = 0;
   scrfcn = dscrfcn = fcpair = NULL;
+  
+  neltypes = 0;
+  for (int i = 0; i < maxelt; i++) {
+    Omega_meam[i] = Z_meam[i] = A_meam[i] = rho0_meam[i] = beta0_meam[i] = 
+      beta1_meam[i]= beta2_meam[i] = beta3_meam[i] = 
+      t0_meam[i] = t1_meam[i] = t2_meam[i] = t3_meam[i] = 
+      rho_ref_meam[i] = ibar_meam[i] = ielt_meam[i] = 0.0;
+    for (int j = 0; j < maxelt; j++) {
+      lattce_meam[i][j] = FCC;
+      Ec_meam[i][j] = re_meam[i][j] = alpha_meam[i][j] = delta_meam[i][j] = Ec_meam[i][j] = ebound_meam[i][j] = attrac_meam[i][j] = repuls_meam[i][j] = 0.0;
+      nn2_meam[i][j] = zbl_meam[i][j] = eltind[i][j] = 0;
+    }
+  }
 }
 
 MEAM::~MEAM()

src/USER-MEAMC/pair_meamc.cpp

@@ -206,7 +206,7 @@ void PairMEAMC::settings(int narg, char **/*arg*/)
 
 void PairMEAMC::coeff(int narg, char **arg)
 {
-  int i,j,m,n;
+  int m,n;
 
   if (!allocated) allocate();
 
@@ -222,7 +222,7 @@ void PairMEAMC::coeff(int narg, char **arg)
   // elements = list of unique element names
 
   if (nelements) {
-    for (i = 0; i < nelements; i++) delete [] elements[i];
+    for (int i = 0; i < nelements; i++) delete [] elements[i];
     delete [] elements;
     delete [] mass;
   }
@@ -231,7 +231,7 @@ void PairMEAMC::coeff(int narg, char **arg)
   elements = new char*[nelements];
   mass = new double[nelements];
 
-  for (i = 0; i < nelements; i++) {
+  for (int i = 0; i < nelements; i++) {
     n = strlen(arg[i+3]) + 1;
     elements[i] = new char[n];
     strcpy(elements[i],arg[i+3]);
@@ -247,8 +247,9 @@ void PairMEAMC::coeff(int narg, char **arg)
   // read args that map atom types to MEAM elements
   // map[i] = which element the Ith atom type is, -1 if not mapped
 
-  for (i = 4 + nelements; i < narg; i++) {
+  for (int i = 4 + nelements; i < narg; i++) {
     m = i - (4+nelements) + 1;
+    int j;
     for (j = 0; j < nelements; j++)
       if (strcmp(arg[i],elements[j]) == 0) break;
     if (j < nelements) map[m] = j;

src/USER-MISC/README

@@ -30,6 +30,7 @@ compute ackland/atom, Gerolf Ziegenhain, gerolf at ziegenhain.com, 4 Oct 2007
 compute basal/atom, Christopher Barrett, cdb333 at cavs.msstate.edu, 3 Mar 2013
 compute cnp/atom, Paulo Branicio (USC), branicio at usc.edu, 31 May 2017
 compute entropy/atom, Pablo Piaggi (EPFL), pablo.piaggi at epfl.ch, 15 June 2018
+compute pressure/cylinder, Cody K. Addington (NCSU), , 2 Oct 2018
 compute stress/mop, Romain Vermorel (U Pau) & Laurent Joly (U Lyon), romain.vermorel at univ-pau.fr & ljoly.ulyon at gmail.com, 5 Sep 18
 compute stress/mop/profile, Romain Vermorel (U Pau) & Laurent Joly (U Lyon), romain.vermorel at univ-pau.fr & ljoly.ulyon at gmail.com, 5 Sep 18
 compute temp/rotate, Laurent Joly (U Lyon), ljoly.ulyon at gmail.com, 8 Aug 11
@@ -44,6 +45,7 @@ dihedral_style table/cut, Mike Salerno, ksalerno@pha.jhu.edu, 11 May 18
 fix addtorque, Laurent Joly (U Lyon), ljoly.ulyon at gmail.com, 8 Aug 11
 fix ave/correlate/long, Jorge Ramirez (UPM Madrid), jorge.ramirez at upm.es, 21 Oct 2015
 fix bond/react, Jacob Gissinger (CU Boulder), info at disarmmd.org, 24 Feb 2018
+fix ffl, David Wilkins (EPFL Lausanne), david.wilkins @ epfl.ch, 28 Sep 2018
 fix filter/corotate, Lukas Fath (KIT), lukas.fath at kit.edu, 15 Mar 2017
 fix flow/gauss, Joel Eaves (CU Boulder), Joel.Eaves@Colorado.edu, 23 Aug 2016
 fix gle, Michele Ceriotti (EPFL Lausanne), michele.ceriotti at gmail.com, 24 Nov 2014

src/USER-MISC/compute_pressure_cylinder.cpp

@@ -0,0 +1,484 @@
+/* ----------------------------------------------------------------------
+   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.
+------------------------------------------------------------------------- */
+
+#include <cmath>
+#include <cstring>
+#include <cstdlib>
+#include "compute_pressure_cylinder.h"
+#include "atom.h"
+#include "update.h"
+#include "force.h"
+#include "pair.h"
+#include "neighbor.h"
+#include "neigh_request.h"
+#include "neigh_list.h"
+#include "group.h"
+#include "memory.h"
+#include "error.h"
+#include "citeme.h"
+#include "domain.h"
+#include "math_const.h"
+
+using namespace LAMMPS_NS;
+using namespace MathConst;
+
+static const char cite_compute_pressure_cylinder[] =
+  "compute pressure/cylinder:\n\n"
+  "@Article{Addington,\n"
+  " author = {C. K. Addington, Y. Long, K. E. Gubbins},\n"
+  " title = {The pressure in interfaces having cylindrical geometry},\n"
+  " journal = {J.~Chem.~Phys.},\n"
+  " year =    2018,\n"
+  " volume =  149,\n"
+  " pages =   {084109}\n"
+  "}\n\n";
+
+/*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+  Calculate the configurational components of the pressure tensor in
+  cylindrical geometry, according to the formulation of Addington et al. (2018)
+  +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
+
+ComputePressureCyl::ComputePressureCyl(LAMMPS *lmp, int narg, char **arg) :
+  Compute(lmp, narg, arg),
+  Pr_temp(NULL), Pr_all(NULL), Pz_temp(NULL), Pz_all(NULL), Pphi_temp(NULL),
+  Pphi_all(NULL), R(NULL), Rinv(NULL), R2(NULL), PrAinv(NULL), PzAinv(NULL),
+  R2kin(NULL), density_temp(NULL), invVbin(NULL), density_all(NULL),
+  tangent(NULL), ephi_x(NULL), ephi_y(NULL), binz(NULL)
+{
+  if (lmp->citeme) lmp->citeme->add(cite_compute_pressure_cylinder);
+  if (narg != 7) error->all(FLERR,"Illegal compute pressure/cylinder command");
+
+  zlo=force->numeric(FLERR,arg[3]);
+  zhi=force->numeric(FLERR,arg[4]);
+  Rmax=force->numeric(FLERR,arg[5]);
+  bin_width=force->numeric(FLERR,arg[6]);
+
+  if ((bin_width <= 0.0) || (bin_width < Rmax))
+    error->all(FLERR,"Illegal compute pressure/cylinder command");
+  if ((zhi < zlo) || ((zhi-zlo) < bin_width))
+    error->all(FLERR,"Illegal compute pressure/cylinder command");
+  if ((zhi > domain->boxhi[2]) || (zlo < domain->boxlo[2]))
+    error->all(FLERR,"Illegal compute pressure/cylinder command");
+
+  nbins=(int)(Rmax/bin_width);
+  nzbins=(int)((zhi-zlo)/bin_width);
+
+  // NOTE: at 2^22 = 4.2M bins, we will be close to exhausting allocatable
+  // memory on a 32-bit environment. so we use this as an upper limit.
+
+  if ((nbins < 1) || (nzbins < 1) || (nbins > 2>>22) || (nbins > 2>>22))
+    error->all(FLERR,"Illegal compute pressure/cylinder command");
+
+  array_flag=1;
+  vector_flag=0;
+  extarray=0;
+  size_array_cols = 5;  // r, number density, Pr, Pphi, Pz
+  size_array_rows = nbins;
+
+  Pr_temp = new double[nbins];
+  Pr_all = new double[nbins];
+  Pz_temp = new double[nbins];
+  Pz_all = new double[nbins];
+  Pphi_temp = new double[nbins];
+  Pphi_all = new double[nbins];
+  R  = new double[nbins];
+  R2 = new double[nbins];
+  PrAinv = new double[nbins];
+  PzAinv = new double[nbins];
+  Rinv = new double[nbins];
+  binz = new double[nzbins];
+
+  R2kin = new double[nbins];
+  density_temp = new double[nbins];
+  invVbin = new double[nbins];
+  density_all = new double[nbins];
+
+  memory->create(array,nbins,5,"PN:array");
+
+  nphi=360;
+  tangent = new double[nphi];
+  ephi_x = new double[nphi];
+  ephi_y = new double[nphi];
+
+  nktv2p = force->nktv2p;
+
+}
+
+/* ---------------------------------------------------------------------- */
+
+ComputePressureCyl::~ComputePressureCyl()
+{
+  // count all of these for memory usage
+  memory->destroy(array);
+  delete [] R;
+  delete [] Rinv;
+  delete [] R2;
+  delete [] R2kin;
+  delete [] invVbin;
+  delete [] density_temp;
+  delete [] density_all;
+  delete [] tangent;
+  delete [] ephi_x;
+  delete [] ephi_y;
+  delete [] Pr_temp;
+  delete [] Pr_all;
+  delete [] Pz_temp;
+  delete [] Pz_all;
+  delete [] Pphi_temp;
+  delete [] Pphi_all;
+  delete [] PrAinv;
+  delete [] PzAinv;
+  delete [] binz;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void ComputePressureCyl::init()
+{
+  if (force->pair == NULL)
+    error->all(FLERR,"No pair style is defined for compute pressure/cylinder");
+  if (force->pair->single_enable == 0)
+    error->all(FLERR,"Pair style does not support compute pressure/cylinder");
+
+  double phi;
+
+  for (int iphi = 0; iphi < nphi; iphi++) {
+    phi=((double)iphi)*MY_PI/180.0;
+    tangent[iphi]=tan(phi);
+    ephi_x[iphi]=-sin(phi);
+    ephi_y[iphi]=cos(phi);
+  }
+
+  for (int iq = 0; iq < nbins; iq++) {
+    R[iq]=((double)iq+0.5)*bin_width;
+    Rinv[iq]=1.0/R[iq];
+    R2[iq]=R[iq]*R[iq];
+    R2kin[iq]=(((double)iq)+1.0)*bin_width;
+    R2kin[iq]*=R2kin[iq];
+    PrAinv[iq]=1.0/(2.0*MY_PI*(zhi-zlo)*R[iq]);
+  }
+  PphiAinv=1.0/((zhi-zlo)*bin_width*2.0*(double)nphi);
+
+  invVbin[0]=1.0/((zhi-zlo)*MY_PI*R2kin[0]);
+  PzAinv[0]=1.0/(MY_PI*R2kin[0]*((double)nzbins));
+
+  for (int jq = 1; jq < nbins; jq++) {
+    invVbin[jq]=1.0/((zhi-zlo)*MY_PI*(R2kin[jq]-R2kin[jq-1]));
+    PzAinv[jq]=1.0/(MY_PI*(R2kin[jq]-R2kin[jq-1])*((double)nzbins));
+  }
+
+  // need an occasional half neighbor list
+  int irequest = neighbor->request(this,instance_me);
+  neighbor->requests[irequest]->pair = 0;
+  neighbor->requests[irequest]->compute = 1;
+  neighbor->requests[irequest]->occasional = 1;
+
+  for (int zzz = 0; zzz < nzbins; zzz++) binz[zzz]=(((double)zzz)+0.5)*bin_width+zlo;
+
+}
+
+/* ---------------------------------------------------------------------- */
+
+void ComputePressureCyl::init_list(int /* id */, NeighList *ptr)
+{
+  list = ptr;
+}
+
+/* ---------------------------------------------------------------------- */
+
+
+/* ----------------------------------------------------------------------
+   count pairs and compute pair info on this proc
+   only count pair once if newton_pair is off
+   both atom I,J must be in group
+   if flag is set, compute requested info about pair
+------------------------------------------------------------------------- */
+
+void ComputePressureCyl::compute_array()
+{
+  invoked_array = update->ntimestep;
+
+  int ibin;
+
+  // clear pressures
+  for (ibin = 0; ibin < nbins; ibin++) {
+    density_temp[ibin]=0.0;
+    density_all[ibin]=0.0;
+    Pr_temp[ibin]=0.0;
+    Pr_all[ibin]=0.0;
+    Pphi_temp[ibin]=0.0;
+    Pphi_all[ibin]=0.0;
+    Pz_temp[ibin]=0.0;
+    Pz_all[ibin]=0.0;
+  }
+
+  // what processor am I?
+  int me;
+  MPI_Comm_rank(world,&me);
+
+  int i,j,ii,jj,inum,jnum,itype,jtype;
+  tagint itag,jtag;
+  double xtmp,ytmp,ztmp,delx,dely,delz;
+  double rsq,fpair,factor_coul,factor_lj;
+  int *ilist,*jlist,*numneigh,**firstneigh;
+
+  double **x = atom->x;
+  tagint *tag = atom->tag;
+  int *type = atom->type;
+  int *mask = atom->mask;
+  int nlocal = atom->nlocal;
+  double *special_coul = force->special_coul;
+  double *special_lj = force->special_lj;
+  int newton_pair = force->newton_pair;
+
+  // invoke half neighbor list (will copy or build if necessary)
+  neighbor->build_one(list);
+
+  inum = list->inum;
+  ilist = list->ilist;
+  numneigh = list->numneigh;
+  firstneigh = list->firstneigh;
+
+  // calculate number density (by radius)
+  double temp_R2;
+  for (i = 0; i < nlocal; i++) if ((x[i][2] < zhi) && (x[i][2] > zlo)) {
+    temp_R2=x[i][0]*x[i][0]+x[i][1]*x[i][1];
+    if (temp_R2 > R2kin[nbins-1]) continue; // outside of Rmax
+
+    for (j = 0; j < nbins; j++) if (temp_R2 < R2kin[j]) break;
+
+    density_temp[j]+=invVbin[j];
+  }
+  MPI_Allreduce(density_temp,density_all,nbins,MPI_DOUBLE,MPI_SUM,world);
+  for (i = 0; i < nbins; i++) array[i][1]=density_all[i]; // NEW
+
+  // loop over neighbors of my atoms
+  // skip if I or J are not in group
+  // for newton = 0 and J = ghost atom,
+  //   need to insure I,J pair is only output by one proc
+  //   use same itag,jtag logic as in Neighbor::neigh_half_nsq()
+  // for flag = 0, just count pair interactions within force cutoff
+  // for flag = 1, calculate requested output fields
+
+  Pair *pair = force->pair;
+  double **cutsq = force->pair->cutsq;
+
+  double r1=0.0;
+  double r2=0.0;
+  double risq,rjsq;
+  double A,B,C,D;
+  double alpha1,alpha2;
+  double xi,yi,zi,dx,dy,dz;
+  double xR,yR,zR,fn;
+  double alpha,xL,yL,zL,L2,ftphi,ftz;
+  double sqrtD;
+
+  for (ii = 0; ii < inum; ii++) {
+    i = ilist[ii];
+    if (!(mask[i] & groupbit)) continue;
+
+    xtmp = x[i][0];
+    ytmp = x[i][1];
+    ztmp = x[i][2];
+    itag = tag[i];
+    itype = type[i];
+    jlist = firstneigh[i];
+    jnum = numneigh[i];
+
+    r1=x[i][0]*x[i][0]+x[i][1]*x[i][1];
+
+    for (jj = 0; jj < jnum; jj++) {
+      j = jlist[jj];
+      factor_lj = special_lj[sbmask(j)];
+      factor_coul = special_coul[sbmask(j)];
+      j &= NEIGHMASK;
+
+      if (!(mask[j] & groupbit)) continue;
+
+      // itag = jtag is possible for long cutoffs that include images of self
+      // do calculation only on appropriate processor
+      if (newton_pair == 0 && j >= nlocal) {
+        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) {
+            if (x[j][1] < ytmp) continue;
+            if (x[j][1] == ytmp && x[j][0] < xtmp) continue;
+          }
+        }
+      }
+
+      delx = xtmp - x[j][0];
+      dely = ytmp - x[j][1];
+      delz = ztmp - x[j][2];
+
+      r2=x[j][0]*x[j][0]+x[j][1]*x[j][1];
+
+      // ri is smaller of r1 and r2
+      if (r2 < r1) {
+        risq=r2;
+        rjsq=r1;
+        xi=x[j][0];
+        yi=x[j][1];
+        zi=x[j][2];
+        dx=x[i][0]-x[j][0];
+        dy=x[i][1]-x[j][1];
+        dz=x[i][2]-x[j][2];
+      } else {
+        risq=r1;
+        rjsq=r2;
+        xi=x[i][0];
+        yi=x[i][1];
+        zi=x[i][2];
+        dx=x[j][0]-x[i][0];
+        dy=x[j][1]-x[i][1];
+        dz=x[j][2]-x[i][2];
+      }
+
+      rsq = delx*delx + dely*dely + delz*delz;
+      jtype = type[j];
+      if (rsq >= cutsq[itype][jtype]) continue;
+
+      pair->single(i,j,itype,jtype,rsq,factor_coul,factor_lj,fpair);
+
+      A=dx*dx+dy*dy;
+      B=2.0*(xi*dx+yi*dy);
+
+      // normal pressure contribution P_rhorho
+      for (ibin = 0; ibin < nbins; ibin++) {
+        // completely inside of R
+        if (rjsq < R2[ibin]) continue;
+
+        C=risq-R2[ibin];
+        D=B*B-4.0*A*C;
+
+        // completely outside of R
+        if (D < 0.0) continue;
+
+        sqrtD=sqrt(D);
+        alpha1=0.5*(-B+sqrtD)/A;
+        alpha2=0.5*(-B-sqrtD)/A;
+
+        if ((alpha1 > 0.0) && (alpha1 < 1.0)) {
+          zR=zi+alpha1*dz;
+          if ((zR < zhi) && (zR > zlo))
+          {
+            xR=xi+alpha1*dx;
+            yR=yi+alpha1*dy;
+            fn=fpair*fabs(xR*dx+yR*dy);
+
+            Pr_temp[ibin]+=fn;
+          }
+        }
+        if ((alpha2 > 0.0) && (alpha2 < 1.0)) {
+          zR=zi+alpha2*dz;
+          if ((zR < zhi) && (zR > zlo)) {
+            xR=xi+alpha2*dx;
+            yR=yi+alpha2*dy;
+            fn=fpair*fabs(xR*dx+yR*dy);
+
+            Pr_temp[ibin]+=fn;
+          }
+        }
+      }
+
+      // azimuthal pressure contribution (P_phiphi)
+      for (int iphi = 0; iphi < nphi; iphi++) {
+        alpha=(yi-xi*tangent[iphi])/(dx*tangent[iphi]-dy);
+
+        // no intersection with phi surface
+        if ((alpha >= 1.0) || (alpha <= 0.0)) continue;
+
+        // no contribution (outside of averaging region)
+        zL=zi+alpha*dz;
+        if ((zL > zhi) || (zL < zlo)) continue;
+
+        xL=xi+alpha*dx;
+        yL=yi+alpha*dy;
+
+        L2=xL*xL+yL*yL;
+
+        // no intersection (outside of Rmax)
+        if (L2 > R2kin[nbins-1]) continue;
+
+        ftphi=fabs(dx*ephi_x[iphi]+dy*ephi_y[iphi])*fpair;
+
+        // add to appropriate bin
+        for (ibin = 0; ibin < nbins; ibin++) if (L2 < R2kin[ibin]) {
+          Pphi_temp[ibin]+=ftphi;
+          break;
+        }
+      }
+
+      // z pressure contribution (P_zz)
+      for (int zbin = 0; zbin < nzbins; zbin++) {
+        // check if interaction contributes
+        if ((x[i][2] > binz[zbin]) && (x[j][2] > binz[zbin])) continue;
+        if ((x[i][2] < binz[zbin]) && (x[j][2] < binz[zbin])) continue;
+
+        alpha=(binz[zbin]-zi)/dz;
+
+        xL=xi+alpha*dx;
+        yL=yi+alpha*dy;
+
+        L2=xL*xL+yL*yL;
+
+        if (L2 > R2kin[nbins-1]) continue;
+
+        ftz=fabs(dz)*fpair;
+
+        // add to appropriate bin
+        for (ibin = 0; ibin < nbins; ibin++) if (L2 < R2kin[ibin]) {
+          Pz_temp[ibin]+=ftz;
+          break;
+        }
+      }
+    }
+  }
+
+  // calculate pressure (force over area)
+  for (ibin = 0; ibin < nbins; ibin++) {
+    Pr_temp[ibin]*=PrAinv[ibin]*Rinv[ibin];
+    Pphi_temp[ibin]*=PphiAinv;
+    Pz_temp[ibin]*=PzAinv[ibin];
+  }
+
+  // communicate these values across processors
+  MPI_Allreduce(Pr_temp,Pr_all,nbins,MPI_DOUBLE,MPI_SUM,world);
+  MPI_Allreduce(Pphi_temp,Pphi_all,nbins,MPI_DOUBLE,MPI_SUM,world);
+  MPI_Allreduce(Pz_temp,Pz_all,nbins,MPI_DOUBLE,MPI_SUM,world);
+
+  // populate array
+  for (ibin = 0; ibin < nbins; ibin++) {
+    array[ibin][0]=R[ibin];
+    array[ibin][2]=Pr_all[ibin]*nktv2p;
+    array[ibin][3]=Pphi_all[ibin]*nktv2p;
+    array[ibin][4]=Pz_all[ibin]*nktv2p;
+  }
+
+}
+
+/* ----------------------------------------------------------------------
+memory usage of data
+------------------------------------------------------------------------- */
+
+double ComputePressureCyl::memory_usage()
+{
+  double bytes =
+  (3.0*(double)nphi + 16.0*(double)nbins+5.0*(double)nbins) * sizeof(double);
+  return bytes;
+}

src/USER-MISC/compute_pressure_cylinder.h

@@ -0,0 +1,74 @@
+/* -*- 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 COMPUTE_CLASS
+
+ComputeStyle(pressure/cylinder,ComputePressureCyl)
+
+#else
+
+#ifndef LMP_COMPUTE_PRESSURE_CYLINDER
+#define LMP_COMPUTE_PRESSURE_CYLINDER
+
+#include "compute.h"
+
+namespace LAMMPS_NS {
+
+class ComputePressureCyl : public Compute {
+ public:
+  ComputePressureCyl(class LAMMPS *, int, char **);
+  ~ComputePressureCyl();
+  void init();
+  void init_list(int, class NeighList *);
+  void compute_array();
+  double memory_usage();
+
+ private:
+  int nbins,nphi,nzbins;
+  double *Pr_temp,*Pr_all,*Pz_temp,*Pz_all,*Pphi_temp,*Pphi_all;
+  double *R,*Rinv,*R2,*PrAinv,*PzAinv,PphiAinv;
+  double Rmax,bin_width,nktv2p;
+  double *R2kin,*density_temp,*invVbin,*density_all;
+  double *tangent,*ephi_x,*ephi_y;
+  double *binz;
+
+  double zlo,zhi;
+
+  class NeighList *list;
+
+};
+
+}
+
+#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: No pair style is defined for compute pressure/cylinder
+
+Self-explanatory.
+
+E: Pair style does not support compute pressure/cylinder
+
+The pair style does not have a single() function, so it can
+not be invoked by compute pressure/cylinder.
+
+*/
+

src/USER-MISC/fix_ffl.cpp

@@ -0,0 +1,483 @@
+/* ----------------------------------------------------------------------
+   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.
+------------------------------------------------------------------------- */
+
+/* ----------------------------------------------------------------------
+   Fast-forward Langevin thermostat, see
+   M. Hijazi, D. M. Wilkins, M. Ceriotti, J. Chem. Phys. 148, 184109 (2018)
+------------------------------------------------------------------------- */
+
+/* ----------------------------------------------------------------------
+   Contributing authors: Lionel Constantin (EPFL), David M. Wilkins (EPFL),
+			 Michele Ceriotti (EPFL)
+------------------------------------------------------------------------- */
+
+#include <mpi.h>
+#include <cmath>
+#include <cstring>
+#include <cstdlib>
+#include "fix_ffl.h"
+#include "math_extra.h"
+#include "atom.h"
+#include "atom_vec_ellipsoid.h"
+#include "force.h"
+#include "update.h"
+#include "modify.h"
+#include "compute.h"
+#include "domain.h"
+#include "region.h"
+#include "respa.h"
+#include "comm.h"
+#include "input.h"
+#include "variable.h"
+#include "random_mars.h"
+#include "memory.h"
+#include "error.h"
+#include "group.h"
+
+using namespace LAMMPS_NS;
+using namespace FixConst;
+
+enum {NOBIAS,BIAS};
+enum {CONSTANT,EQUAL,ATOM};
+enum {NO_FLIP, FLIP_RESCALE, FLIP_HARD, FLIP_SOFT};
+//#define FFL_DEBUG 1
+
+#define MAXLINE 1024
+
+/* syntax for fix_ffl:
+ * fix nfix id-group ffl tau Tstart Tstop seed [flip_type]
+ *
+ *                                                                        */
+
+/* ---------------------------------------------------------------------- */
+
+
+FixFFL::FixFFL(LAMMPS *lmp, int narg, char **arg) :
+  Fix(lmp, narg, arg) {
+
+
+  if (narg < 7)
+    error->all(FLERR,"Illegal fix ffl command. Expecting: fix <fix-ID>"
+               " <group-ID> ffl <tau> <Tstart> <Tstop> <seed>  ");
+
+  restart_peratom = 1;
+  time_integrate = 1;
+  scalar_flag = 1;
+
+  //gamma = 1/ time constant(tau)
+  if (force->numeric(FLERR,arg[3]) <= 0)
+    error->all(FLERR,"Illegal fix ffl tau value, should be greater than 0");
+  gamma = 1.0/force->numeric(FLERR,arg[3]);
+  ffl_every=1;
+  ffl_step=0;
+
+  // start temperature (t ramp)
+  t_start = force->numeric(FLERR,arg[4]);
+
+  // final temperature (t ramp)
+  t_stop = force->numeric(FLERR,arg[5]);
+
+  // PRNG seed
+  int seed = force->inumeric(FLERR,arg[6]);
+
+  // Flip type used, uses rescale if no flip is given
+  if (narg == 8) {
+    if (strcmp(arg[7],"no_flip") == 0) {
+      flip_int = NO_FLIP;
+    } else if (strcmp(arg[7],"rescale") == 0) {
+      flip_int = FLIP_RESCALE;
+    } else if (strcmp(arg[7],"hard") == 0) {
+      flip_int = FLIP_HARD;
+    } else if (strcmp(arg[7],"soft") == 0) {
+      flip_int = FLIP_SOFT;
+    } else {
+      error->all(FLERR,"Illegal fix ffl flip type, only accepts : rescale - hard - soft - no_flip");
+    }
+  } else {
+    flip_int = FLIP_RESCALE;
+  }
+
+  t_target=t_start;
+  const double kT = t_target * force->boltz / force->mvv2e;
+
+
+  // initialize Marsaglia RNG with processor-unique seed
+  // NB: this means runs will not be the same with different numbers of processors
+  if (seed <= 0) error->all(FLERR,"Illegal fix ffl command");
+  random = new RanMars(lmp,seed + comm->me);
+
+  // allocate per-type arrays for mass-scaling
+  sqrt_m=NULL;
+  memory->grow(sqrt_m, atom->ntypes+1,"ffl:sqrt_m");
+
+  // allocates space for temporaries
+  ffl_tmp1=ffl_tmp2=NULL;
+
+  grow_arrays(atom->nmax);
+
+  // add callbacks to enable restarts
+  atom->add_callback(0);
+  atom->add_callback(1);
+
+  energy = 0.0;
+}
+
+
+/* --- Frees up memory used by temporaries and buffers ------------------ */
+
+FixFFL::~FixFFL() {
+  delete random;
+
+  memory->destroy(sqrt_m);
+  memory->destroy(ffl_tmp1);
+  memory->destroy(ffl_tmp2);
+}
+
+/* ---------------------------------------------------------------------- */
+
+int FixFFL::setmask() {
+  int mask = 0;
+
+  mask |= INITIAL_INTEGRATE;
+  mask |= FINAL_INTEGRATE;
+  mask |= INITIAL_INTEGRATE_RESPA;
+  mask |= FINAL_INTEGRATE_RESPA;
+  mask |= THERMO_ENERGY;
+
+
+  return mask;
+}
+
+/* ------- Initializes one-time quantities for FFL ---------------------- */
+
+void FixFFL::init() {
+  doffl = 1;
+  dtv = update->dt;
+  dtf = 0.5 * update->dt * force->ftm2v;
+
+  // set force prefactors
+  if (!atom->rmass) {
+    for (int i = 1; i <= atom->ntypes; i++) {
+      sqrt_m[i] = sqrt(atom->mass[i]);
+    }
+  }
+
+  if (strstr(update->integrate_style,"respa")) {
+    nlevels_respa = ((Respa *) update->integrate)->nlevels;
+    step_respa = ((Respa *) update->integrate)->step;
+  }
+
+  init_ffl();
+}
+
+/* ------- Initializes constants for FFL (change with T and dt) ------- */
+
+void FixFFL::init_ffl() {
+  const double kT = t_target * force->boltz / force->mvv2e;
+
+  // compute constants for FFL
+
+  c1 = exp ( - gamma * 0.5 * dtv );
+  c2 = sqrt( (1.0 - c1*c1)* kT ); //without the mass term
+
+
+}
+
+
+
+/* ---------------------------------------------------------------------- */
+
+void FixFFL::setup(int vflag) {
+  if (strstr(update->integrate_style,"verlet"))
+    post_force(vflag);
+  else {
+    ((Respa *) update->integrate)->copy_flevel_f(nlevels_respa-1);
+    post_force_respa(vflag,nlevels_respa-1,0);
+    ((Respa *) update->integrate)->copy_f_flevel(nlevels_respa-1);
+  }
+}
+
+void FixFFL::ffl_integrate() {
+  double **v = atom->v;
+  double *rmass = atom->rmass, smi, ismi;
+  double factor;
+  int *type = atom->type;
+  int *mask = atom->mask;
+  int nlocal = atom->nlocal;
+  if (igroup == atom->firstgroup) nlocal = atom->nfirst;
+
+  // loads momentum data (mass-scaled) into the temporary vectors for the propagation
+  int nk=0;
+  double deltae=0.0;
+  for (int i = 0; i < nlocal; i++) {
+    if (mask[i] & groupbit) {
+      if (rmass) smi = sqrt(rmass[i]);
+      else smi = sqrt_m[type[i]];
+
+      for (int k = 0; k<3; k++) {
+        // first loads velocities and accumulates conserved quantity
+        ffl_tmp2[nk] = v[i][k] * smi;
+        deltae += ffl_tmp2[nk] * ffl_tmp2[nk];
+        nk++;
+      }
+    }
+  }
+
+  //fills up a vector of random numbers
+  for (int i = 0; i < nk; i++) ffl_tmp1[i] = random->gaussian();
+
+
+  // unloads momentum data (mass-scaled) from the temporary vectors
+  nk=0;
+  for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) {
+      if (rmass) ismi = 1.0 / sqrt(rmass[i]);
+      else ismi = 1.0/ sqrt_m[type[i]];
+
+      for (int k = 0; k<3; k++) {
+        // fetches new velocities and completes computation of the conserved quantity change
+        v[i][k]= c1*v[i][k] + c2*ffl_tmp1[nk]*ismi;
+
+        deltae-= v[i][k]*v[i][k] /ismi /ismi;
+
+        //flips the sign of the momentum (HARD FLIP)
+        if ( flip_int == FLIP_HARD) {
+          if (v[i][k]*ffl_tmp2[nk] < 0.0) v[i][k] = -v[i][k];
+        }
+
+        nk++;
+      }
+    }
+
+  //rescale operation (RESCALE FLIP)
+  if (flip_int == FLIP_RESCALE) {
+    nk=0;
+    for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) {
+      factor = sqrt ((v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2]) /
+                       (ffl_tmp2[nk]*ffl_tmp2[nk] + ffl_tmp2[nk+1]*ffl_tmp2[nk+1]
+                        + ffl_tmp2[nk+2]*ffl_tmp2[nk+2]));
+
+      for (int k = 0; k<3; k++) {
+        v[i][k]= factor * ffl_tmp2[nk];
+        nk++;
+      }
+    }
+  }
+
+
+  //soft flip operation (SOFT FLIP)
+  if (flip_int == FLIP_SOFT) {
+    nk=0;
+    for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) {
+      factor = v[i][0]*ffl_tmp2[nk] + v[i][1]*ffl_tmp2[nk+1] + v[i][2]*ffl_tmp2[nk+2];
+      if (factor < 0) {
+        factor =  factor / (ffl_tmp2[nk]*ffl_tmp2[nk] + ffl_tmp2[nk+1]*ffl_tmp2[nk+1]
+                            + ffl_tmp2[nk+2]*ffl_tmp2[nk+2]);
+
+        for (int k = 0; k<3; k++) {
+          v[i][k] -= 2.0 * factor * ffl_tmp2[nk];
+          nk++;
+        }
+      } else {
+        nk += 3;
+      }
+    }
+
+  }
+
+  energy += deltae*0.5*force->mvv2e;
+
+}
+
+void FixFFL::initial_integrate(int vflag) {
+  double dtfm;
+
+  // update v and x of atoms in group
+  double **x = atom->x;
+  double **v = atom->v;
+  double **f = atom->f;
+  double *rmass = atom->rmass;
+  double *mass = atom->mass;
+  int *type = atom->type;
+  int *mask = atom->mask;
+  int nlocal = atom->nlocal;
+  if (igroup == atom->firstgroup) nlocal = atom->nfirst;
+
+  ffl_step--;
+  if (doffl && ffl_step<1) ffl_integrate();
+
+  if (rmass) {
+    for (int i = 0; i < nlocal; i++)
+      if (mask[i] & groupbit) {
+        dtfm = dtf / rmass[i];
+        v[i][0] += dtfm * f[i][0];
+        v[i][1] += dtfm * f[i][1];
+        v[i][2] += dtfm * f[i][2];
+        x[i][0] += dtv * v[i][0];
+        x[i][1] += dtv * v[i][1];
+        x[i][2] += dtv * v[i][2];
+      }
+
+  } else {
+    for (int i = 0; i < nlocal; i++)
+      if (mask[i] & groupbit) {
+        dtfm = dtf / mass[type[i]];
+        v[i][0] += dtfm * f[i][0];
+        v[i][1] += dtfm * f[i][1];
+        v[i][2] += dtfm * f[i][2];
+        x[i][0] += dtv * v[i][0];
+        x[i][1] += dtv * v[i][1];
+        x[i][2] += dtv * v[i][2];
+      }
+  }
+}
+
+void FixFFL::final_integrate() {
+  double dtfm;
+
+  // update v of atoms in group
+
+  double **v = atom->v;
+  double **f = atom->f;
+  double *rmass = atom->rmass;
+  double *mass = atom->mass;
+  int *type = atom->type;
+  int *mask = atom->mask;
+  int nlocal = atom->nlocal;
+  if (igroup == atom->firstgroup) nlocal = atom->nfirst;
+
+  if (rmass) {
+    for (int i = 0; i < nlocal; i++)
+      if (mask[i] & groupbit) {
+        dtfm = dtf / rmass[i];
+        v[i][0] += dtfm * f[i][0];
+        v[i][1] += dtfm * f[i][1];
+        v[i][2] += dtfm * f[i][2];
+      }
+
+  } else {
+    for (int i = 0; i < nlocal; i++)
+      if (mask[i] & groupbit) {
+        dtfm = dtf / mass[type[i]];
+        v[i][0] += dtfm * f[i][0];
+        v[i][1] += dtfm * f[i][1];
+        v[i][2] += dtfm * f[i][2];
+      }
+  }
+
+  if (doffl && ffl_step<1) {
+    ffl_integrate();
+    ffl_step = ffl_every;
+  }
+
+  // Change the temperature for the next step
+  double delta = update->ntimestep - update->beginstep;
+  delta /= update->endstep - update->beginstep;
+  t_target = t_start + delta * (t_stop - t_start);
+  if (t_stop != t_start) {
+    // only updates if it is really necessary
+    init_ffl();
+  }
+
+}
+/* ---------------------------------------------------------------------- */
+
+void FixFFL::initial_integrate_respa(int vflag, int ilevel, int iloop) {
+  dtv = step_respa[ilevel];
+  dtf = 0.5 * step_respa[ilevel] * force->ftm2v;
+
+  // innermost level - NVE update of v and x
+  // all other levels - NVE update of v
+
+  if (ilevel==nlevels_respa-1) ffl_integrate();
+  doffl=0;
+  if (ilevel == 0) initial_integrate(vflag);
+  else {
+    final_integrate();
+  }
+}
+
+void FixFFL::final_integrate_respa(int ilevel, int iloop) {
+
+  dtv = step_respa[ilevel];
+  dtf = 0.5 * step_respa[ilevel] * force->ftm2v;
+  doffl=0;
+  final_integrate();
+  if (ilevel==nlevels_respa-1) ffl_integrate();
+}
+
+
+double FixFFL::compute_scalar() {
+
+  double energy_me = energy;
+  double energy_all;
+  MPI_Allreduce(&energy_me,&energy_all,1,MPI_DOUBLE,MPI_SUM,world);
+
+  return energy_all;
+}
+
+/* ----------------------------------------------------------------------
+   extract thermostat properties
+------------------------------------------------------------------------- */
+
+void *FixFFL::extract(const char *str, int &dim) {
+  dim = 0;
+  if (strcmp(str,"t_target") == 0) {
+    return &t_target;
+  }
+  return NULL;
+}
+
+
+/* ----------------------------------------------------------------------
+   Called when a change to the target temperature is requested mid-run
+------------------------------------------------------------------------- */
+
+void FixFFL::reset_target(double t_new) {
+
+  t_target = t_start = t_stop = t_new;
+}
+
+/* ----------------------------------------------------------------------
+   Called when a change to the timestep is requested mid-run
+------------------------------------------------------------------------- */
+
+void FixFFL::reset_dt() {
+  // set the time integration constants
+  dtv = update->dt;
+  dtf = 0.5 * update->dt * (force->ftm2v);
+  init_ffl();
+}
+
+/* ----------------------------------------------------------------------
+   memory usage of local atom-based arrays
+------------------------------------------------------------------------- */
+
+double FixFFL::memory_usage() {
+  double bytes = atom->nmax*(3*2)*sizeof(double);
+  return bytes;
+}
+
+
+/* ----------------------------------------------------------------------
+   allocate local atom-based arrays
+------------------------------------------------------------------------- */
+
+void FixFFL::grow_arrays(int nmax) {
+  memory->grow(ffl_tmp1, nmax*3,"ffl:tmp1");
+  memory->grow(ffl_tmp2, nmax*3,"ffl:tmp2");
+  //zeroes out temporary buffers
+  for (int i=0; i< nmax*3; ++i) ffl_tmp1[i] = 0.0;
+  for (int i=0; i< nmax*3; ++i) ffl_tmp2[i] = 0.0;
+}
+
+

src/USER-MISC/fix_ffl.h

@@ -0,0 +1,67 @@
+/* -*- 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 FIX_CLASS
+
+FixStyle(ffl,FixFFL)
+
+#else
+
+#ifndef LMP_FIX_FFL_H
+#define LMP_FIX_FFL_H
+
+#include "fix.h"
+
+namespace LAMMPS_NS {
+
+class FixFFL : public Fix {
+ public:
+  FixFFL(class LAMMPS *, int, char **);
+  virtual ~FixFFL();
+  int setmask();
+  void init();
+  void setup(int);
+  void ffl_integrate();
+  void initial_integrate_respa(int vflag, int ilevel, int iloop);
+  void final_integrate_respa(int ilevel, int iloop);
+  void initial_integrate(int vflag);
+  void final_integrate();
+  double compute_scalar();
+  void reset_target(double);
+  virtual void reset_dt();
+  double memory_usage();
+  void grow_arrays(int);
+
+  virtual void *extract(const char *, int &);
+
+  void init_ffl();
+ protected:
+  double *ffl_tmp1, *ffl_tmp2;
+  double t_start, t_stop, t_target;
+  double dtv, dtf, c1, c2, gamma;
+  char flip_type[10];
+
+  int doffl, ffl_every, ffl_step, flip_int;
+  class RanMars *random;
+  double *sqrt_m;
+  double *step_respa;
+  double energy;
+  int nlevels_respa;
+
+  double **vaux;
+};
+
+}
+
+#endif
+#endif

src/USER-MISC/pair_ilp_graphene_hbn.cpp

@@ -249,11 +249,9 @@ void PairILPGrapheneHBN::compute(int eflag, int vflag)
         f[i][0] += fkcx - fprod1[0]*Tap;
         f[i][1] += fkcy - fprod1[1]*Tap;
         f[i][2] += fkcz - fprod1[2]*Tap;
-        if (newton_pair || j < nlocal) {
-          f[j][0] -= fkcx + fprod2[0]*Tap;
-          f[j][1] -= fkcy + fprod2[1]*Tap;
-          f[j][2] -= fkcz + fprod2[2]*Tap;
-        }
+        f[j][0] -= fkcx + fprod2[0]*Tap;
+        f[j][1] -= fkcy + fprod2[1]*Tap;
+        f[j][2] -= fkcz + fprod2[2]*Tap;
 
         // calculate the forces acted on the neighbors of atom i from atom j
         ILP_neighs_i = ILP_firstneigh[i];
@@ -274,15 +272,13 @@ void PairILPGrapheneHBN::compute(int eflag, int vflag)
         for (ll = 0; ll < ILP_numneigh[j]; ll++) {
           l = ILP_neighs_j[ll];
           if (l == j) continue;
-          if (newton_pair || l < nlocal) {
-            // derivatives of the product of rji and nj respect to rl, l=0,1,2, where atom l is the neighbors of atom j
-            dprodnorm2[0] = dnormal[0][0][ll][j]*delx + dnormal[1][0][ll][j]*dely + dnormal[2][0][ll][j]*delz;
-            dprodnorm2[1] = dnormal[0][1][ll][j]*delx + dnormal[1][1][ll][j]*dely + dnormal[2][1][ll][j]*delz;
-            dprodnorm2[2] = dnormal[0][2][ll][j]*delx + dnormal[1][2][ll][j]*dely + dnormal[2][2][ll][j]*delz;
-            f[l][0] += (-prodnorm2*dprodnorm2[0]*fpair2)*Tap;
-            f[l][1] += (-prodnorm2*dprodnorm2[1]*fpair2)*Tap;
-            f[l][2] += (-prodnorm2*dprodnorm2[2]*fpair2)*Tap;
-          }
+          // derivatives of the product of rji and nj respect to rl, l=0,1,2, where atom l is the neighbors of atom j
+          dprodnorm2[0] = dnormal[0][0][ll][j]*delx + dnormal[1][0][ll][j]*dely + dnormal[2][0][ll][j]*delz;
+          dprodnorm2[1] = dnormal[0][1][ll][j]*delx + dnormal[1][1][ll][j]*dely + dnormal[2][1][ll][j]*delz;
+          dprodnorm2[2] = dnormal[0][2][ll][j]*delx + dnormal[1][2][ll][j]*dely + dnormal[2][2][ll][j]*delz;
+          f[l][0] += (-prodnorm2*dprodnorm2[0]*fpair2)*Tap;
+          f[l][1] += (-prodnorm2*dprodnorm2[1]*fpair2)*Tap;
+          f[l][2] += (-prodnorm2*dprodnorm2[2]*fpair2)*Tap;
         }
 
         if (eflag) {
@@ -729,7 +725,8 @@ void PairILPGrapheneHBN::ILP_neigh()
 
     ILP_firstneigh[i] = neighptr;
     ILP_numneigh[i] = n;
-    if (n > 3) error->all(FLERR,"There are too many neighbors for some atoms, please reduce the cutoff for normals");
+    if (n == 0) error->all(FLERR,"Could not build neighbor list to calculate normals, please check your configuration");
+    if (n > 3) error->all(FLERR,"There are too many neighbors for some atoms, please check your configuration");
     ipage->vgot(n);
     if (ipage->status())
       error->one(FLERR,"Neighbor list overflow, boost neigh_modify one");

src/USER-MISC/pair_kolmogorov_crespi_full.cpp

@@ -249,11 +249,9 @@ void PairKolmogorovCrespiFull::compute(int eflag, int vflag)
         f[i][0] += fkcx - fprod1[0]*Tap;
         f[i][1] += fkcy - fprod1[1]*Tap;
         f[i][2] += fkcz - fprod1[2]*Tap;
-        if (newton_pair || j < nlocal) {
-          f[j][0] -= fkcx + fprod2[0]*Tap;
-          f[j][1] -= fkcy + fprod2[1]*Tap;
-          f[j][2] -= fkcz + fprod2[2]*Tap;
-        }
+        f[j][0] -= fkcx + fprod2[0]*Tap;
+        f[j][1] -= fkcy + fprod2[1]*Tap;
+        f[j][2] -= fkcz + fprod2[2]*Tap;
 
 	// calculate the forces acted on the neighbors of atom i from atom j
 	KC_neighs_i = KC_firstneigh[i];
@@ -274,15 +272,13 @@ void PairKolmogorovCrespiFull::compute(int eflag, int vflag)
   	for (ll = 0; ll < KC_numneigh[j]; ll++) {
 	  l = KC_neighs_j[ll];
           if (l == j) continue;
-          if (newton_pair || l < nlocal) {
-            // derivatives of the product of rji and nj respect to rl, l=0,1,2, where atom l is the neighbors of atom j
-            dprodnorm2[0] = dnormal[0][0][ll][j]*delx + dnormal[1][0][ll][j]*dely + dnormal[2][0][ll][j]*delz;
-            dprodnorm2[1] = dnormal[0][1][ll][j]*delx + dnormal[1][1][ll][j]*dely + dnormal[2][1][ll][j]*delz;
-            dprodnorm2[2] = dnormal[0][2][ll][j]*delx + dnormal[1][2][ll][j]*dely + dnormal[2][2][ll][j]*delz;
-	    f[l][0] += (-prodnorm2*dprodnorm2[0]*fpair2)*Tap;
-	    f[l][1] += (-prodnorm2*dprodnorm2[1]*fpair2)*Tap;
-	    f[l][2] += (-prodnorm2*dprodnorm2[2]*fpair2)*Tap;
-	  }
+          // derivatives of the product of rji and nj respect to rl, l=0,1,2, where atom l is the neighbors of atom j
+          dprodnorm2[0] = dnormal[0][0][ll][j]*delx + dnormal[1][0][ll][j]*dely + dnormal[2][0][ll][j]*delz;
+          dprodnorm2[1] = dnormal[0][1][ll][j]*delx + dnormal[1][1][ll][j]*dely + dnormal[2][1][ll][j]*delz;
+          dprodnorm2[2] = dnormal[0][2][ll][j]*delx + dnormal[1][2][ll][j]*dely + dnormal[2][2][ll][j]*delz;
+	  f[l][0] += (-prodnorm2*dprodnorm2[0]*fpair2)*Tap;
+	  f[l][1] += (-prodnorm2*dprodnorm2[1]*fpair2)*Tap;
+	  f[l][2] += (-prodnorm2*dprodnorm2[2]*fpair2)*Tap;
 	}
 
         if (eflag) {
@@ -734,7 +730,8 @@ void PairKolmogorovCrespiFull::KC_neigh()
 
     KC_firstneigh[i] = neighptr;
     KC_numneigh[i] = n;
-    if (n > 3) error->all(FLERR,"There are too many neighbors for some atoms, please reduce the cutoff for normals");
+    if (n == 0) error->all(FLERR,"Could not build neighbor list to calculate normals, please check your configuration");
+    if (n > 3) error->all(FLERR,"There are too many neighbors for some atoms, please check your configuration");
     ipage->vgot(n);
     if (ipage->status())
       error->one(FLERR,"Neighbor list overflow, boost neigh_modify one");

src/USER-OMP/fix_omp.cpp

@@ -318,8 +318,11 @@ void FixOMP::set_neighbor_omp()
   const int neigh_omp = _neighbor ? 1 : 0;
   const int nrequest = neighbor->nrequest;
 
+  // flag *all* neighbor list requests as USER-OMP threaded,
+  // but skip lists already flagged as USER-INTEL threaded
   for (int i = 0; i < nrequest; ++i)
-    neighbor->requests[i]->omp = neigh_omp;
+    if (! neighbor->requests[i]->intel)
+      neighbor->requests[i]->omp = neigh_omp;
 }
 
 /* ---------------------------------------------------------------------- */

src/USER-OMP/pair_nb3b_harmonic_omp.cpp

@@ -1,170 +0,0 @@
-/* ----------------------------------------------------------------------
-   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 <cmath>
-#include "pair_nb3b_harmonic_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;
-
-/* ---------------------------------------------------------------------- */
-
-PairNb3bHarmonicOMP::PairNb3bHarmonicOMP(LAMMPS *lmp) :
-  PairNb3bHarmonic(lmp), ThrOMP(lmp, THR_PAIR)
-{
-  suffix_flag |= Suffix::OMP;
-  respa_enable = 0;
-}
-
-/* ---------------------------------------------------------------------- */
-
-void PairNb3bHarmonicOMP::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) {
-      if (eflag) {
-        eval<1,1>(ifrom, ito, thr);
-      } else {
-        eval<1,0>(ifrom, ito, thr);
-      }
-    } else eval<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>
-void PairNb3bHarmonicOMP::eval(int iifrom, int iito, ThrData * const thr)
-{
-  int i,j,k,ii,jj,kk,jnum,jnumm1;
-  int itype,jtype,ktype,ijparam,ikparam,ijkparam;
-  double xtmp,ytmp,ztmp,evdwl;
-  double rsq1,rsq2;
-  double delr1[3],delr2[3],fj[3],fk[3];
-  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 int * _noalias const type = atom->type;
-
-  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];
-    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];
-    jnumm1 = jnum - 1;
-
-    for (jj = 0; jj < jnumm1; jj++) {
-      j = jlist[jj];
-      j &= NEIGHMASK;
-      jtype = map[type[j]];
-      ijparam = 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[ijparam].cutsq) continue;
-
-      double fjxtmp,fjytmp,fjztmp;
-      fjxtmp = fjytmp = fjztmp = 0.0;
-
-      for (kk = jj+1; kk < jnum; kk++) {
-        k = jlist[kk];
-        k &= NEIGHMASK;
-        ktype = map[type[k]];
-        ikparam = elem2param[itype][ktype][ktype];
-        ijkparam = 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[ikparam].cutsq) continue;
-
-        threebody(&params[ijparam],&params[ikparam],&params[ijkparam],
-                  rsq1,rsq2,delr1,delr2,fj,fk,EFLAG,evdwl);
-
-        fxtmp -= fj[0] + fk[0];
-        fytmp -= fj[1] + fk[1];
-        fztmp -= fj[2] + fk[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 (EVFLAG) ev_tally3_thr(this,i,j,k,evdwl,0.0,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 PairNb3bHarmonicOMP::memory_usage()
-{
-  double bytes = memory_usage_thr();
-  bytes += PairNb3bHarmonic::memory_usage();
-
-  return bytes;
-}