Merge branch 'master' into lammps_gjf

cmake/CMakeLists.txt

@@ -589,6 +589,7 @@ if(BUILD_TOOLS)
     enable_language(Fortran)
     add_executable(chain.x ${LAMMPS_TOOLS_DIR}/chain.f)
     target_link_libraries(chain.x ${CMAKE_Fortran_IMPLICIT_LINK_LIBRARIES})
+    install(TARGETS chain.x DESTINATION ${CMAKE_INSTALL_BINDIR})
   endif()
 
   enable_language(C)

cmake/Modules/FindKIM-API.cmake

@@ -37,6 +37,26 @@
 # KIM-API-CMAKE_Fortran_COMPILER
 #
 
+function(_KIMAPI_GET_VERSION _OUT_ver _version_hdr)
+  if(NOT EXISTS ${_version_hdr})
+    message(FATAL_ERROR "Header file ${_version_hdr} not found (check value of KIM-API_INCLUDE_DIR)")
+  endif()
+  foreach(_var KIM_VERSION_MAJOR KIM_VERSION_MINOR KIM_VERSION_PATCH)  
+    file(STRINGS ${_version_hdr} _contents REGEX "#define ${_var}[ \t]+")
+    if(_contents)
+      string(REGEX REPLACE ".*#define ${_var}[ \t]+([0-9]+).*" "\\1" _${_var} "${_contents}")
+      if(${_${_var}} STREQUAL "")
+        message(FATAL_ERROR "Version parsing failed for ${_var} in ${_version_hdr}, got empty return!")
+      elseif(NOT ${_${_var}} MATCHES "^[0-9]+$")
+        message(FATAL_ERROR "Version parsing failed for ${_var} in ${_version_hdr}, excepted a number but got ${_${_var}}!")
+      endif()
+    else()
+      message(FATAL_ERROR "No ${_var} line found in include file ${_version_hdr}")
+    endif()
+  endforeach()
+  set(${_OUT_ver} ${_KIM_VERSION_MAJOR}.${_KIM_VERSION_MINOR}.${_KIM_VERSION_PATCH} PARENT_SCOPE)
+endfunction()
+
 if(KIM-API_FIND_QUIETLY)
   set(REQ_OR_QUI "QUIET")
 else()
@@ -54,6 +74,12 @@ if(KIM-API_FOUND)
   pkg_get_variable(KIM-API_CMAKE_Fortran_COMPILER libkim-api CMAKE_Fortran_COMPILER)
 endif()
 
+if(KIM-API_INCLUDEDIR)
+  _KIMAPI_GET_VERSION(KIM-API_VERSION ${KIM-API_INCLUDEDIR}/KIM_Version.h)
+else()
+  set(KIM-API_VERSION 0)
+endif()
+
 # handle the QUIETLY and REQUIRED arguments and set KIM-API_FOUND to TRUE
 # if all listed variables are TRUE
-find_package_handle_standard_args(KIM-API REQUIRED_VARS KIM-API_LIBRARIES)
+find_package_handle_standard_args(KIM-API REQUIRED_VARS KIM-API_LIBRARIES VERSION_VAR KIM-API_VERSION)

cmake/Modules/Packages/KIM.cmake

@@ -1,4 +1,5 @@
 if(PKG_KIM)
+  set(KIM-API_MIN_VERSION 2.1)
   find_package(CURL)
   if(CURL_FOUND)
     include_directories(${CURL_INCLUDE_DIRS})
@@ -7,10 +8,17 @@ if(PKG_KIM)
   endif()
   find_package(KIM-API QUIET)
   if(KIM-API_FOUND)
-    set(DOWNLOAD_KIM_DEFAULT OFF)
+    if (KIM-API_VERSION VERSION_LESS ${KIM-API_MIN_VERSION})
+      if ("${DOWNLOAD_KIM}" STREQUAL "")
+        message(WARNING "Unsuitable KIM-API version \"${KIM-API_VERSION}\" found, but required is at least \"${KIM-API_MIN_VERSION}\".  Default behavior set to download and build our own.")
+      endif()
+      set(DOWNLOAD_KIM_DEFAULT ON)
+    else()
+      set(DOWNLOAD_KIM_DEFAULT OFF)
+    endif()
   else()
-    if (NOT DOWNLOAD_KIM)
+    if ("${DOWNLOAD_KIM}" STREQUAL "")
-      message(WARNING "KIM-API package not found.  We will download and build our own")
+      message(WARNING "KIM-API package not found.  Default behavior set to download and build our own")
     endif()
     set(DOWNLOAD_KIM_DEFAULT ON)
   endif()
@@ -28,8 +36,8 @@ if(PKG_KIM)
       message(FATAL_ERROR "Compiling the KIM-API library requires a Fortran compiler")
     endif()
     ExternalProject_Add(kim_build
-      URL https://s3.openkim.org/kim-api/kim-api-2.1.2.txz
+      URL https://s3.openkim.org/kim-api/kim-api-2.1.3.txz
-      URL_MD5 6ac52e14ef52967fc7858220b208cba5
+      URL_MD5 6ee829a1bbba5f8b9874c88c4c4ebff8
       BINARY_DIR build
       CMAKE_ARGS -DCMAKE_C_COMPILER=${CMAKE_C_COMPILER}
                  -DCMAKE_CXX_COMPILER=${CMAKE_CXX_COMPILER}
@@ -42,7 +50,7 @@ if(PKG_KIM)
     set(KIM-API_LDFLAGS ${INSTALL_DIR}/${CMAKE_INSTALL_LIBDIR}/libkim-api${CMAKE_SHARED_LIBRARY_SUFFIX})
     list(APPEND LAMMPS_DEPS kim_build)
   else()
-    find_package(KIM-API REQUIRED)
+    find_package(KIM-API ${KIM-API_MIN_VERSION} REQUIRED)
   endif()
   list(APPEND LAMMPS_LINK_LIBS "${KIM-API_LDFLAGS}")
   include_directories(${KIM-API_INCLUDE_DIRS})

cmake/Modules/Packages/KOKKOS.cmake

@@ -23,6 +23,7 @@ if(PKG_KOKKOS)
                          ${KOKKOS_PKG_SOURCES_DIR}/fix_nh_kokkos.cpp
                          ${KOKKOS_PKG_SOURCES_DIR}/nbin_kokkos.cpp
                          ${KOKKOS_PKG_SOURCES_DIR}/npair_kokkos.cpp
+                         ${KOKKOS_PKG_SOURCES_DIR}/npair_halffull_kokkos.cpp
                          ${KOKKOS_PKG_SOURCES_DIR}/domain_kokkos.cpp
                          ${KOKKOS_PKG_SOURCES_DIR}/modify_kokkos.cpp)
 
@@ -38,6 +39,7 @@ if(PKG_KOKKOS)
   # register kokkos-only styles
   RegisterNBinStyle(${KOKKOS_PKG_SOURCES_DIR}/nbin_kokkos.h)
   RegisterNPairStyle(${KOKKOS_PKG_SOURCES_DIR}/npair_kokkos.h)
+  RegisterNPairStyle(${KOKKOS_PKG_SOURCES_DIR}/npair_halffull_kokkos.h)
 
   if(PKG_USER-DPD)
     get_property(KOKKOS_PKG_SOURCES GLOBAL PROPERTY KOKKOS_PKG_SOURCES)

doc/src/Commands_fix.txt

@@ -232,6 +232,7 @@ OPT.
 "wall/lj1043"_fix_wall.html,
 "wall/lj126"_fix_wall.html,
 "wall/lj93 (k)"_fix_wall.html,
+"wall/morse"_fix_wall.html,
 "wall/piston"_fix_wall_piston.html,
 "wall/reflect (k)"_fix_wall_reflect.html,
 "wall/region"_fix_wall_region.html,

doc/src/Commands_pair.txt

@@ -65,6 +65,7 @@ OPT.
 "colloid (go)"_pair_colloid.html,
 "comb (o)"_pair_comb.html,
 "comb3"_pair_comb.html,
+"cosine/squared"_pair_cosine_squared.html,
 "coul/cut (gko)"_pair_coul.html,
 "coul/cut/soft (o)"_pair_fep_soft.html,
 "coul/debye (gko)"_pair_coul.html,

doc/src/Eqs/pair_cosine_squared.tex

@@ -0,0 +1,16 @@
+\documentclass[12pt]{article}
+
+\usepackage{amsmath}
+
+\begin{document}
+
+\begin{align*}
+E = 
+\begin{cases}
+-\epsilon& \quad r < \sigma \\
+-\epsilon\cos\left(\frac{\pi\left(r - \sigma\right)}{2\left(r_c - \sigma\right)}\right)&\quad \sigma \leq r < r_c \\
+0& \quad r \geq r_c
+\end{cases}
+\end{align*}
+
+\end{document}

doc/src/Eqs/pair_cosine_squared_wca.tex

@@ -0,0 +1,11 @@
+\documentstyle[12pt]{article}
+
+\begin{document}
+
+$$
+ E = \epsilon \left[ \left(\frac{\sigma}{r}\right)^{12} - 
+                       2\left(\frac{\sigma}{r}\right)^6 + 1\right]
+                       , \quad r < \sigma
+$$
+
+\end{document}

doc/src/Howto_chunk.txt

@@ -136,7 +136,9 @@ The "compute chunk/spread/atom"_compute_chunk_spread_atom.html command
 spreads per-chunk values to each atom in the chunk, producing per-atom
 values as its output.  This can be useful for outputting per-chunk
 values to a per-atom "dump file"_dump.html.  Or for using an atom's
-associated chunk value in an "atom-style variable"_variable.html.
+associated chunk value in an "atom-style variable"_variable.html.  Or
+as input to the "fix ave/chunk"_fix_ave_chunk.html command to
+spatially average per-chunk values calculated by a per-chunk compute.
 
 The "compute reduce/chunk"_compute_reduce_chunk.html command reduces a
 peratom value across the atoms in each chunk to produce a value per
@@ -184,12 +186,20 @@ compute cc1 all chunk/atom c_cluster compress yes
 compute size all property/chunk cc1 count
 fix 1 all ave/histo 100 1 100 0 20 20 c_size mode vector ave running beyond ignore file tmp.histo :pre
 
-(6) An example of using a per-chunk value to apply per-atom forces to
+(6) An example for using a per-chunk value to apply per-atom forces to
 compress individual polymer chains (molecules) in a mixture, is
 explained on the "compute
 chunk/spread/atom"_compute_chunk_spread_atom.html command doc page.
 
-(7) An example of using one set of per-chunk values for molecule
+(7) An example for using one set of per-chunk values for molecule
 chunks, to create a 2nd set of micelle-scale chunks (clustered
 molecules, due to hydrophobicity), is explained on the "compute
 chunk/reduce"_compute_reduce_chunk.html command doc page.
+
+(8) An example for using one set of per-chunk values (dipole moment
+vectors) for molecule chunks, spreading the values to each atom in
+each chunk, then defining a second set of chunks as spatial bins, and
+using the "fix ave/chunk"_fix_ave_chunk.html command to calculate an
+average dipole moment vector for each bin.  This example is explained
+on the "compute chunk/spread/atom"_compute_chunk_spread_atom.html
+command doc page.

doc/src/Howto_walls.txt

@@ -59,14 +59,15 @@ granular particles; all the other commands create smooth walls.
 "fix wall/lj126"_fix_wall.html - flat walls, with Lennard-Jones 12/6 potential
 "fix wall/colloid"_fix_wall.html - flat walls, with "pair_style colloid"_pair_colloid.html potential
 "fix wall/harmonic"_fix_wall.html - flat walls, with repulsive harmonic spring potential
+"fix wall/morse"_fix_wall.html - flat walls, with Morse potential
 "fix wall/region"_fix_wall_region.html - use region surface as wall
 "fix wall/gran"_fix_wall_gran.html - flat or curved walls with "pair_style granular"_pair_gran.html potential :ul
 
-The {lj93}, {lj126}, {colloid}, and {harmonic} styles all allow the
+The {lj93}, {lj126}, {colloid}, {harmonic}, and {morse} styles all
-flat walls to move with a constant velocity, or oscillate in time.
+allow the flat walls to move with a constant velocity, or oscillate in
-The "fix wall/region"_fix_wall_region.html command offers the most
+time.  The "fix wall/region"_fix_wall_region.html command offers the
-generality, since the region surface is treated as a wall, and the
+most generality, since the region surface is treated as a wall, and
-geometry of the region can be a simple primitive volume (e.g. a
+the geometry of the region can be a simple primitive volume (e.g. a
 sphere, or cube, or plane), or a complex volume made from the union
 and intersection of primitive volumes.  "Regions"_region.html can also
 specify a volume "interior" or "exterior" to the specified primitive

doc/src/compute.txt

@@ -217,6 +217,7 @@ compute"_Commands_compute.html doc page are followed by one or more of
 "heat/flux"_compute_heat_flux.html - heat flux through a group of atoms
 "heat/flux/tally"_compute_tally.html -
 "hexorder/atom"_compute_hexorder_atom.html - bond orientational order parameter q6
+"hma"_compute_hma.html - harmonically mapped averaging for atomic crystals
 "improper"_compute_improper.html - energy of each improper sub-style
 "improper/local"_compute_improper_local.html - angle of each improper
 "inertia/chunk"_compute_inertia_chunk.html - inertia tensor for each chunk

doc/src/compute_chunk_spread_atom.txt

@@ -30,11 +30,18 @@ compute 1 all chunk/spread/atom mychunk c_com[*] c_gyration :pre
 [Description:]
 
 Define a calculation that "spreads" one or more per-chunk values to
-each atom in the chunk.  This can be useful for creating a "dump
+each atom in the chunk.  This can be useful in several scenarios:
-file"_dump.html where each atom lists info about the chunk it is in,
+
-e.g. for post-processing purposes.  It can also be used in "atom-style
+For creating a "dump file"_dump.html where each atom lists info about
-variables"_variable.html that need info about the chunk each atom is
+the chunk it is in, e.g. for post-processing purposes. :ulb,l
-in.  Examples are given below.
+
+To access chunk value in "atom-style variables"_variable.html that
+need info about the chunk each atom is in. :l
+
+To use the "fix ave/chunk"_fix_ave_chunk.html command to spatially
+average per-chunk values calculated by a per-chunk compute. :l,ule
+
+Examples are given below.
 
 In LAMMPS, chunks are collections of atoms defined by a "compute
 chunk/atom"_compute_chunk_atom.html command, which assigns each atom
@@ -148,6 +155,28 @@ thermo_style    custom step etotal press v_ave :pre
 
 :line
 
+Here is an example for using one set of chunks, defined for molecules,
+to compute the dipole moment vector for each chunk.  E.g. for water
+molecules. Then spreading those values to each atom in each chunk.
+Then defining a second set of chunks based on spatial bins.  And
+finally, using the "fix ave/chunk"_fix_ave_chunk.html command to
+calculate an average dipole moment vector per spatial bin.
+
+compute       cmol all chunk/atom molecule
+compute       dipole all dipole/chunk cmol
+compute       spread all chunk/spread/atom cmol c_dipole\[1\] c_dipole\[2\] c_dipole\[3\]
+compute       cspatial all chunk/atom bin/1d z lower 0.1 units reduced
+fix           ave all ave/chunk 100 10 1000 cspatial c_spread\[*\] :pre
+
+Note that the "fix ave/chunk"_fix_ave_chunk.html command requires
+per-atom values as input.  That is why the compute chunk/spread/atom
+command is used to assign per-chunk values to each atom in the chunk.
+If a molecule straddles bin boundaries, each of its atoms contributes
+in a weighted manner to the average dipole moment of the spatial bin
+it is in.
+
+:line
+
 [Output info:]
 
 This compute calculates a per-atom vector or array, which can be

doc/src/compute_hma.txt

@@ -0,0 +1,184 @@
+"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 hma command :h3
+
+[Syntax:]
+
+compute ID group-ID hma temp-ID keyword ... :pre
+
+ID, group-ID are documented in "compute"_compute.html command :l
+hma = style name of this compute command :l
+temp-ID = ID of fix that specifies the set temperature during canonical simulation :l
+keyword = {anharmonic} {u} {p Pharm} {cv} :l
+  {anharmonic} = compute will return anharmonic property values
+  {u} = compute will return potential energy
+  {p} = compute will return pressure.  the following keyword must be the difference between the harmonic pressure and lattice pressure as described below
+  {cv} = compute will return the heat capacity :pre
+:ule
+
+[Examples:]
+
+compute 2 all hma 1 u
+compute 2 all hma 1 anharmonic u p 0.9
+compute 2 all hma 1 u cv :pre
+
+
+
+[Description:]
+
+Define a computation that calculates the properties of a solid (potential
+energy, pressure or heat capacity), using the harmonically-mapped averaging
+(HMA) method. 
+This command yields much higher precision than the equivalent compute commands
+("compute pe"_compute_pe.html, "compute pressure"_compute_pressure.html, etc.)
+commands during a canonical simulation of an atomic crystal. Specifically,
+near melting HMA can yield averages of a given precision an order of magnitude
+faster than conventional methods, and this only improves as the temperatures is
+lowered.  This is particularly important for evaluating the free energy by
+thermodynamic integration, where the low-temperature contributions are the
+greatest source of statistical uncertainty.  Moreover, HMA has other
+advantages, including smaller potential-truncation effects, finite-size
+effects, smaller timestep inaccuracy, faster equilibration and shorter
+decorrelation time.
+
+HMA should not be used if atoms are expected to diffuse.  It is also
+restricted to simulations in the NVT ensemble.  While this compute may be
+used with any potential in LAMMPS, it will provide inaccurate results
+for potentials that do not go to 0 at the truncation distance;
+"pair_lj_smooth_linear"_pair_lj_smooth_linear.html and Ewald summation should
+work fine, while "pair_lj"_pair_lj.html will perform poorly unless 
+the potential is shifted (via "pair_modify"_pair_modify.html shift) or the cutoff is large.  Furthermore, computation of the heat capacity with
+this compute is restricted to those that implement the single_hessian method
+in Pair.  Implementing single_hessian in additional pair styles is simple.
+Please contact Andrew Schultz (ajs42 at buffalo.edu) and David Kofke (kofke at
+buffalo.edu) if your desired pair style does not have this method.  This is
+the list of pair styles that currently implement pair_hessian:
+
+"lj_smooth_linear"_pair_lj_smooth_linear.html :l
+:ule
+
+In this method, the analytically known harmonic behavior of a crystal is removed from the traditional ensemble
+averages, which leads to an accurate and precise measurement of the anharmonic contributions without contamination 
+by noise produced by the already-known harmonic behavior. 
+A detailed description of this method can be found in ("Moustafa"_#hma-Moustafa). The potential energy is computed by the formula:
+
+\begin\{equation\}
+\left< U\right>_\{HMA\} = \frac\{d\}\{2\} (N-1) k_B T  + \left< U + \frac\{1\}\{2\} F\bullet\Delta r \right>
+\end\{equation\}
+
+where \(N\) is the number of atoms in the system, \(k_B\) is Boltzmann's
+constant, \(T\) is the temperature, \(d\) is the
+dimensionality of the system (2 or 3 for 2d/3d), \(F\bullet\Delta r\) is the sum of dot products of the 
+atomic force vectors and displacement (from lattice sites) vectors, and \(U\) is the sum of 
+pair, bond, angle, dihedral, improper, kspace (long-range), and fix energies. 
+
+The pressure is computed by the formula:
+
+\begin\{equation\}
+\left< P\right>_\{HMA\} = \Delta \hat P + \left< P_\{vir\} + \frac\{\beta \Delta \hat P - \rho\}\{d(N-1)\} F\bullet\Delta r \right>
+\end\{equation\}
+
+where \(\rho\) is the number density of the system, \(\Delta \hat P\) is the
+difference between the harmonic and lattice pressure, \(P_\{vir\}\) is
+the virial pressure computed as the sum of pair, bond, angle, dihedral,
+improper, kspace (long-range), and fix contributions to the force on each
+atom, and \(k_B=1/k_B T\).  Although the method will work for any value of \(\Delta \hat P\)
+specified (use pressure "units"_units.html), the precision of the resultant
+pressure is sensitive to \(\Delta \hat P\); the precision tends to be
+best when \(\Delta \hat P\) is the actual the difference between the lattice
+pressure and harmonic pressure.
+
+\begin\{equation\}
+\left<C_V \right>_\{HMA\} = \frac\{d\}\{2\} (N-1) k_B + \frac\{1\}\{k_B T^2\} \left( \left<
+U_\{HMA\}^2 \right> - \left<U_\{HMA\}\right>^2 \right) + \frac\{1\}\{4 T\}
+\left< F\bullet\Delta r + \Delta r \bullet \Phi \bullet \Delta r \right>
+\end\{equation\}
+
+where \(\Phi\) is the Hessian matrix. The compute hma command
+computes the full expression for \(C_V\) except for the
+\(\left<U_\{HMA\}^2\right>^2\) in the variance term, which can be obtained by
+passing the {u} keyword; you must add this extra contribution to the \(C_V\)
+value reported by this compute.  The variance term can cause significant
+round-off error when computing \(C_V\).  To address this, the {anharmonic}
+keyword can be passed and/or the output format can be specified with more
+digits.
+
+thermo_modify format float '%22.15e' :pre
+
+The {anharmonic} keyword will instruct the compute to return anharmonic
+properties rather than the full properties, which include lattice, harmonic
+and anharmonic contributions.
+When using this keyword, the compute must be first active (it must be included
+via a "thermo_style custom"_thermo_style.html command) while the atoms are
+still at their lattice sites (before equilibration).
+
+The temp-ID specified with compute hma command should be same as the fix-ID of Nose-Hoover ("fix nvt"_fix_nh.html) or 
+Berendsen ("fix temp/berendsen"_fix_temp_berendsen.html) thermostat used for the simulation. While using this command, Langevin thermostat 
+("fix langevin"_fix_langevin.html) 
+should be avoided as its extra forces interfere with the HMA implementation. 
+
+ 
+
+NOTE: Compute hma command should be used right after the energy minimization, when the atoms are at their lattice sites. 
+The simulation should not be started before this command has been used in the input script.
+
+
+The following example illustrates the placement of this command in the input script:
+
+
+min_style cg 
+minimize 1e-35 1e-15 50000 500000 
+compute 1 all hma thermostatid u
+fix thermostatid all nvt temp 600.0 600.0 100.0 :pre  
+
+
+
+NOTE: Compute hma should be used when the atoms of the solid do not diffuse. Diffusion will reduce the precision in the potential energy computation.
+
+ 
+NOTE: The "fix_modify energy yes"_fix_modify.html command must also be specified if a fix is to contribute potential energy to this command.
+
+An example input script that uses this compute is included in
+examples/USER/hma/ along with corresponding LAMMPS output showing that the HMA
+properties fluctuate less than the corresponding conventional properties.
+
+[Output info:]
+
+This compute calculates a global vector that includes the n properties
+requested as arguments to the command (the potential energy, pressure and/or heat
+capacity).  The elements of the vector can be accessed by indices 1-n by any
+command that uses global vector values as input.  See the "Howto
+output"_Howto_output.html doc page for an overview of LAMMPS output options.
+
+The vector values calculated by this compute are "extensive".  The
+scalar value will be in energy "units"_units.html.
+
+[Restrictions:]
+
+This compute is part of the USER-MISC package.  It is enabled only
+if LAMMPS was built with that package.  See the "Build
+package"_Build_package.html doc page for more info.
+
+Usage restricted to canonical (NVT) ensemble simulation only.
+
+[Related commands:]
+
+"compute pe"_compute_pe.html, "compute pressure"_compute_pressure.html
+
+"dynamical matrix"_dynamical_matrix.html provides a finite difference
+formulation of the hessian provided by Pair's single_hessian, which is used by
+this compute.
+
+[Default:] none
+
+:line
+
+:link(hma-Moustafa)
+[(Moustafa)] Sabry G. Moustafa, Andrew J. Schultz, and David A. Kofke, {Very fast averaging of thermal properties of crystals by molecular simulation}, 
+"Phys. Rev. E \[92\], 043303 (2015)"_https://link.aps.org/doi/10.1103/PhysRevE.92.043303

doc/src/computes.txt

@@ -47,6 +47,7 @@ Computes :h1
    compute_gyration_shape
    compute_heat_flux
    compute_hexorder_atom
+   compute_hma
    compute_improper
    compute_improper_local
    compute_inertia_chunk

doc/src/dynamical_matrix.txt

@@ -47,6 +47,9 @@ package"_Build_package.html doc page for more info.
 
 "fix phonon"_fix_phonon.html
 
+"compute hma"_compute_hma.html uses an analytic formulation of the hessian
+provided by Pair's single_hessian.
+
 [Default:]
 
 The default settings are file = "dynmat.dyn", binary = no  

doc/src/fix.txt

@@ -244,7 +244,7 @@ accelerated styles exist.
 "mscg"_fix_mscg.html - apply MSCG method for force-matching to generate coarse grain models
 "msst"_fix_msst.html - multi-scale shock technique (MSST) integration
 "mvv/dpd"_fix_mvv_dpd.html - DPD using the modified velocity-Verlet integration algorithm
-"mvv/edpd"_fix_mvv_dpd.html - constant energy DPD using the modified velocity-Verlet algrithm
+"mvv/edpd"_fix_mvv_dpd.html - constant energy DPD using the modified velocity-Verlet algorithm
 "mvv/tdpd"_fix_mvv_dpd.html - constant temperature DPD using the modified velocity-Verlet algorithm
 "neb"_fix_neb.html - nudged elastic band (NEB) spring forces
 "nph"_fix_nh.html - constant NPH time integration via Nose/Hoover
@@ -371,6 +371,7 @@ accelerated styles exist.
 "wall/lj1043"_fix_wall.html - Lennard-Jones 10-4-3 wall
 "wall/lj126"_fix_wall.html - Lennard-Jones 12-6 wall
 "wall/lj93"_fix_wall.html - Lennard-Jones 9-3 wall
+"wall/morse"_fix_wall.html - Morse potential wall
 "wall/piston"_fix_wall_piston.html - moving reflective piston wall
 "wall/reflect"_fix_wall_reflect.html - reflecting wall(s)
 "wall/region"_fix_wall_region.html - use region surface as wall

doc/src/fix_wall.txt

@@ -12,16 +12,18 @@ fix wall/lj126 command :h3
 fix wall/lj1043 command :h3
 fix wall/colloid command :h3
 fix wall/harmonic command :h3
+fix wall/morse command :h3
 
 [Syntax:]
 
 fix ID group-ID style face args ... keyword value ... :pre
 
 ID, group-ID are documented in "fix"_fix.html command :ulb,l
-style = {wall/lj93} or {wall/lj126} or {wall/lj1043} or {wall/colloid} or {wall/harmonic} :l
+style = {wall/lj93} or {wall/lj126} or {wall/lj1043} or {wall/colloid} or {wall/harmonic} or {wall/morse} :l
 one or more face/arg pairs may be appended :l
 face = {xlo} or {xhi} or {ylo} or {yhi} or {zlo} or {zhi} :l
-  args = coord epsilon sigma cutoff
+  args for styles {lj93} or {lj126} or {lj1043} or {colloid} or {harmonic} :l
+    args = coord epsilon sigma cutoff
     coord = position of wall = EDGE or constant or variable
       EDGE = current lo or hi edge of simulation box
       constant = number like 0.0 or -30.0 (distance units)
@@ -31,6 +33,19 @@ face = {xlo} or {xhi} or {ylo} or {yhi} or {zlo} or {zhi} :l
     sigma = size factor for wall-particle interaction (distance units)
       sigma can be a variable (see below)
     cutoff = distance from wall at which wall-particle interaction is cut off (distance units) :pre
+  args for style {morse} :l
+    args = coord D_0 alpha r_0 cutoff
+    coord = position of wall = EDGE or constant or variable
+      EDGE = current lo or hi edge of simulation box
+      constant = number like 0.0 or -30.0 (distance units)
+      variable = "equal-style variable"_variable.html like v_x or v_wiggle
+    D_0 = depth of the potential (energy units)
+      D_0 can be a variable (see below)
+    alpha = width factor for wall-particle interaction (1/distance units)
+      alpha can be a variable (see below)
+    r_0 = distance of the potential minimum from the face of region (distance units)
+      r_0 can be a variable (see below)
+    cutoff = distance from wall at which wall-particle interaction is cut off (distance units) :pre
 zero or more keyword/value pairs may be appended :l
 keyword = {units} or {fld} :l
   {units} value = {lattice} or {box}
@@ -48,6 +63,7 @@ keyword = {units} or {fld} :l
 
 fix wallhi all wall/lj93 xlo -1.0 1.0 1.0 2.5 units box
 fix wallhi all wall/lj93 xhi EDGE 1.0 1.0 2.5
+fix wallhi all wall/morse xhi EDGE 1.0 1.0 1.0 2.5 units box
 fix wallhi all wall/lj126 v_wiggle 23.2 1.0 1.0 2.5
 fix zwalls all wall/colloid zlo 0.0 1.0 1.0 0.858 zhi 40.0 1.0 1.0 0.858 :pre
 
@@ -80,6 +96,10 @@ potential:
 
 :c,image(Eqs/fix_wall_harmonic.jpg)
 
+For style {wall/morse}, the energy E is given by a Morse potential:
+
+:c,image(Eqs/pair_morse.jpg)
+
 In all cases, {r} is the distance from the particle to the wall at
 position {coord}, and Rc is the {cutoff} distance at which the
 particle and wall no longer interact.  The energy of the wall
@@ -147,7 +167,13 @@ constant K, and has units (energy/distance^2).  The input parameter
 spring is at the {cutoff}.  This is a repulsive-only spring since the
 interaction is truncated at the {cutoff}
 
-For any wall, the {epsilon} and/or {sigma} parameter can be specified
+For the {wall/morse} style, the three parameters are in this order:
+{D_0} the depth of the potential, {alpha} the width parameter, and
+{r_0} the location of the minimum.  {D_0} has energy units, {alpha}
+inverse distance units, and {r_0} distance units.
+
+For any wall, the {epsilon} and/or {sigma} and/or {alpha} parameter can
+be specified
 as an "equal-style variable"_variable.html, in which case it should be
 specified as v_name, where name is the variable name.  As with a
 variable wall position, the variable is evaluated each timestep and

doc/src/fix_wall_region.txt

@@ -10,19 +10,27 @@ fix wall/region command :h3
 
 [Syntax:]
 
-fix ID group-ID wall/region region-ID style epsilon sigma cutoff :pre
+fix ID group-ID wall/region region-ID style args ... cutoff :pre
 
-ID, group-ID are documented in "fix"_fix.html command
+ID, group-ID are documented in "fix"_fix.html command :ulb,l
-wall/region = style name of this fix command
+wall/region = style name of this fix command :l
-region-ID = region whose boundary will act as wall
+region-ID = region whose boundary will act as wall :l
-style = {lj93} or {lj126} or {lj1043} or {colloid} or {harmonic}
+style = {lj93} or {lj126} or {lj1043} or {colloid} or {harmonic} or {morse} :l
-epsilon = strength factor for wall-particle interaction (energy or energy/distance^2 units)
+args for styles {lj93} or {lj126} or {lj1043} or {colloid} or {harmonic} = :l
-sigma = size factor for wall-particle interaction (distance units)
+   epsilon = strength factor for wall-particle interaction (energy or energy/distance^2 units)
-cutoff = distance from wall at which wall-particle interaction is cut off (distance units) :ul
+   sigma = size factor for wall-particle interaction (distance units) :pre
+args for style {morse} = :l
+   D_0 = depth of the potential (energy units)
+   alpha = width parameter (1/distance units)
+   r_0 = distance of the potential minimum from wall position (distance units) :pre
+cutoff = distance from wall at which wall-particle interaction is cut off (distance units) :l
+:ule
 
 [Examples:]
 
-fix wall all wall/region mySphere lj93 1.0 1.0 2.5 :pre
+fix wall all wall/region mySphere lj93 1.0 1.0 2.5
+fix wall all wall/region mySphere harmonic 1.0 0.0 2.5
+fix wall all wall/region box_top morse 1.0 1.0 1.5 3.0 :pre
 
 [Description:]
 
@@ -122,15 +130,22 @@ the "pair_style colloid"_pair_colloid.html potential:
 :c,image(Eqs/fix_wall_colloid.jpg)
 
 For style {wall/harmonic}, the energy E is given by a harmonic spring
-potential:
+potential (the distance parameter is ignored):
 
 :c,image(Eqs/fix_wall_harmonic.jpg)
 
+For style {wall/morse}, the energy E is given by the Morse potential:
+
+:c,image(Eqs/pair_morse.jpg)
+
+Unlike other styles, this requires three parameters ({D_0}, {alpha}, {r_0}
+in this order) instead of two like for the other wall styles.
+
 In all cases, {r} is the distance from the particle to the region
 surface, and Rc is the {cutoff} distance at which the particle and
-surface no longer interact.  The energy of the wall potential is
+surface no longer interact.  The cutoff is always the last argument.
-shifted so that the wall-particle interaction energy is 0.0 at the
+The energy of the wall potential is shifted so that the wall-particle
-cutoff distance.
+interaction energy is 0.0 at the cutoff distance.
 
 For a full description of these wall styles, see fix_style
 "wall"_fix_wall.html
@@ -179,7 +194,9 @@ option for this fix.
 
 "fix wall/lj93"_fix_wall.html,
 "fix wall/lj126"_fix_wall.html,
+"fix wall/lj1043"_fix_wall.html,
 "fix wall/colloid"_fix_wall.html,
+"fix wall/harmonic"_fix_wall.html,
 "fix wall/gran"_fix_wall_gran.html
 
 [Default:] none

doc/src/lammps.book

@@ -569,6 +569,7 @@ pair_charmm.html
 pair_class2.html
 pair_colloid.html
 pair_comb.html
+pair_cosine_squared.html
 pair_coul.html
 pair_coul_diel.html
 pair_coul_shield.html

doc/src/pair_cosine_squared.txt

@@ -0,0 +1,108 @@
+"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
+
+pair_style cosine/squared command :h3
+
+[Syntax:]
+
+pair_style cosine/squared cutoff :pre
+
+cutoff = global cutoff for cosine-squared interactions (distance units) :ul
+
+pair_coeff i j eps sigma
+pair_coeff i j eps sigma cutoff
+pair_coeff i j eps sigma wca
+pair_coeff i j eps sigma cutoff wca :pre
+
+i,j = a particle type
+eps = interaction strength, i.e. the depth of the potential minimum (energy units)
+sigma = distance of the potential minimum from 0
+cutoff = the cutoff distance for this pair type, if different from global (distance units)
+wca = if specified a Weeks-Chandler-Andersen potential (with eps strength and minimum at sigma) is added, otherwise not :ul
+
+[Examples:]
+
+pair_style cosine/squared 3.0
+pair_coeff * * 1.0 1.3
+pair_coeff 1 3 1.0 1.3 2.0
+pair_coeff 1 3 1.0 1.3 wca
+pair_coeff 1 3 1.0 1.3 2.0 wca :pre
+
+[Description:]
+
+Style {cosine/squared} computes a potential of the form
+
+:c,image(Eqs/pair_cosine_squared.jpg)
+
+between two point particles, where (sigma, -epsilon) is the location of
+the (rightmost) minimum of the potential, as explained in the syntax
+section above.
+
+This potential was first used in (Cooke)_#CKD for a coarse-grained lipid
+membrane model.  It is generally very useful as a non-specific
+interaction potential because it is fully adjustable in depth and width
+while joining the minimum at (sigma, -epsilon) and zero at (cutoff, 0)
+smoothly, requiring no shifting and causing no related artifacts, tail
+energy calculations etc. This evidently requires {cutoff} to be larger
+than {sigma}.
+
+If the {wca} option is used then a Weeks-Chandler-Andersen potential
+(Weeks)_#WCA is added to the above specified cosine-squared potential,
+specifically the following:
+
+:c,image(Eqs/pair_cosine_squared_wca.jpg)
+
+In this case, and this case only, the {sigma} parameter can be equal to
+{cutoff} (sigma = cutoff) which will result in ONLY the WCA potential
+being used (and print a warning), so the minimum will be attained at
+(sigma, 0). This is a convenience feature that enables a purely
+repulsive potential to be used without a need to define an additional
+pair style and use the hybrid styles.
+
+The energy and force of this pair style for parameters epsilon = 1.0,
+sigma = 1.0, cutoff = 2.5, with and without the WCA potential, are shown
+in the graphs below:
+
+:c,image(JPG/pair_cosine_squared_graphs.jpg)
+
+:line
+
+[Mixing, shift, table, tail correction, restart, rRESPA info]:
+
+Mixing is not supported for this style.
+
+The {shift}, {table} and {tail} options are not relevant for this style.
+
+This pair style writes its information to "binary restart
+files"_restart.html, so pair_style and pair_coeff commands do not need
+to be specified in an input script that reads a restart file.
+
+These pair styles can only be used via the {pair} keyword of the
+"run_style respa"_run_style.html command.  They do not support the
+{inner}, {middle}, {outer} keywords.
+
+:line
+
+[Restrictions:]
+
+The {cosine/squared} style is part of the "USER-MISC" package. It is only
+enabled if LAMMPS is build with that package.  See the "Build
+package"_Build_package.html doc page for more info.
+
+[Related commands:]
+
+"pair_coeff"_pair_coeff.html,
+"pair_style lj/cut"_pair_lj.html
+
+[Default:] none
+
+:link(CKD)
+[(Cooke)] "Cooke, Kremer and Deserno, Phys. Rev. E, 72, 011506 (2005)"
+
+:link(WCA)
+[(Weeks)] "Weeks, Chandler and Andersen, J. Chem. Phys., 54, 5237 (1971)"

doc/src/pair_style.txt

@@ -129,6 +129,7 @@ accelerated styles exist.
 "colloid"_pair_colloid.html - integrated colloidal potential
 "comb"_pair_comb.html - charge-optimized many-body (COMB) potential
 "comb3"_pair_comb.html - charge-optimized many-body (COMB3) potential
+"cosine/squared"_pair_cosine_squared.html - Cooke-Kremer-Deserno membrane model potential
 "coul/cut"_pair_coul.html - cutoff Coulombic potential
 "coul/cut/soft"_pair_fep_soft.html - Coulombic potential with a soft core
 "coul/debye"_pair_coul.html - cutoff Coulombic potential with Debye screening

doc/src/pairs.txt

@@ -24,6 +24,7 @@ Pair Styles :h1
    pair_class2
    pair_colloid
    pair_comb
+   pair_cosine_squared
    pair_coul
    pair_coul_diel
    pair_coul_shield

doc/utils/sphinx-config/false_positives.txt

@@ -45,6 +45,7 @@ Aidan
 aij
 airebo
 Aj
+ajs
 ajaramil
 akohlmey
 Aktulga
@@ -99,6 +100,7 @@ antisymmetry
 anton
 Antonelli
 api
+Apoorva
 Appl
 Apu
 arccos
@@ -356,6 +358,7 @@ Cii
 Cij
 cis
 civ
+CKD
 Clang
 clearstore
 Cleary
@@ -421,6 +424,7 @@ coreshell
 cornflowerblue
 cornsilk
 corotate
+corotation
 corotational
 correlator
 cosineshifted
@@ -521,6 +525,7 @@ Dcut
 de
 dE
 De
+decorrelation
 debye
 Debye
 Decius
@@ -1073,6 +1078,7 @@ Hilger
 histo
 histogrammed
 histogramming
+hma
 hmaktulga
 hoc
 Hochbruck
@@ -1329,6 +1335,8 @@ kmax
 Kmax
 Knizhnik
 knl
+Kofke
+kofke
 Kohlmeyer
 Kohn
 kokkos
@@ -1716,6 +1724,7 @@ Morteza
 Mosayebi
 Moseler
 Moskalev
+Moustafa
 mov
 mpi
 MPI
@@ -2215,6 +2224,7 @@ PTM
 ptr
 pu
 purdue
+Purohit
 pushstore
 pvar
 pw
@@ -2432,6 +2442,7 @@ Ryckaert
 Rycroft
 Rydbergs
 Rz
+Sabry
 saddlebrown
 Sadigh
 saed
@@ -2991,6 +3002,7 @@ Waltham
 wavepacket
 wB
 Wbody
+wca
 webpage
 Weckner
 WeinanE

examples/COUPLE/python/example.py

@@ -0,0 +1,36 @@
+# this example requires the LAMMPS Python package (lammps.py) to be installed
+# and LAMMPS to be loadable as shared library in LD_LIBRARY_PATH
+
+import lammps
+
+def callback(caller, ntimestep, nlocal, tag, x, fext):
+    """
+    This callback receives a caller object that was setup when registering the callback
+
+    In addition to timestep and number of local atoms, the tag and x arrays are passed as
+    NumPy arrays. The fext array is a force array allocated for fix external, which
+    can be used to apply forces to all atoms. Simply update the value in the array,
+    it will be directly written into the LAMMPS C arrays
+    """
+    print("Data passed by caller (optional)", caller)
+    print("Timestep:", ntimestep)
+    print("Number of Atoms:", nlocal)
+    print("Atom Tags:", tag)
+    print("Atom Positions:", x)
+    print("Force Additions:", fext)
+    fext.fill(1.0)
+    print("Force additions after update:", fext)
+    print("="*40)
+
+L = lammps.lammps()
+L.file("in.fix_external")
+
+# you can pass an arbitrary Python object to the callback every time it is called
+# this can be useful if you need more state information such as the LAMMPS ptr to
+# make additional library calls
+custom_object = ["Some data", L]
+
+L.set_fix_external_callback("2", callback, custom_object)
+L.command("run 100")
+
+

examples/COUPLE/python/in.fix_external

@@ -0,0 +1,23 @@
+# LAMMPS input for coupling LAMMPS with Python via fix external
+
+units		metal
+dimension	3
+atom_style	atomic
+atom_modify	sort 0 0.0
+
+lattice		diamond 5.43
+region		box block 0 1 0 1 0 1
+create_box	1 box
+create_atoms	1 box
+mass		1 28.08
+
+velocity	all create 300.0 87293 loop geom
+
+fix		1 all nve
+fix		2 all external pf/callback 1 1
+
+#dump		2 all image 25 image.*.jpg type type &
+#		axes yes 0.8 0.02 view 60 -30
+#dump_modify	2 pad 3
+
+thermo		1

examples/USER/misc/hma/README

@@ -0,0 +1,22 @@
+The example input script sets up a simple FCC crystal using the Lennard-Jones
+potential. The script sets up the HMA compute to calculate the energy, pressure
+and heat capacity.  The output columns are:
+
+1: timestep
+2: measured temperature
+3: potential energy
+4: HMA potential energy
+5: pressure
+6: HMA pressure
+7: HMA heat capacity contribution
+
+Averages of the potential energy (#3 and #4) agree although #4 (HMA) is more precise.
+
+Averages of the pressure (#5 and #6) agree once the ideal gas
+contribution is included; #6 (HMA) is more precise.
+
+The heat capacity can be computed from colume #3 (convential) as
+Cv = Var(#3)/(k T^2)
+
+With HMA, the heat capacity can be computed from column #4 and #7 as
+Cv = #7 + Var(#4)/(k T^2)

examples/USER/misc/hma/hma.in

@@ -0,0 +1,36 @@
+# Harmonically mapped average example
+
+units lj 
+dimension 3 
+boundary p p p 
+atom_style atomic 
+atom_modify map array
+# ---------- Create Atoms ----------------------------
+lattice 	fcc 1.0
+region	box block 0 4 0 4 0 4 units lattice
+create_box	1 box
+lattice	fcc 1.0 orient x 1 0 0 orient y 0 1 0 orient z 0 0 1  
+create_atoms 1 region box
+# ---------- Define Interatomic Potential --------------------- 
+pair_style lj/smooth/linear 3
+pair_coeff * * 1.0 1.0 
+mass 1 1.0
+
+atom_modify sort 0 1
+velocity all create 0.1 45678 dist gaussian
+
+compute u all pe 
+  
+compute p all pressure NULL pair  
+
+compute hma all hma settemp u p 9.579586686264458 cv
+
+timestep 0.005
+
+fix settemp all nvt temp 1.0 1.0 0.5
+thermo_style custom elapsed temp c_u c_hma[1] c_p c_hma[2] c_hma[3]
+thermo_modify format float '%22.15e'
+thermo 500
+run 20000
+thermo 20
+run 200000

examples/kim/log.7Aug19.in.kim-ex.melt.clang.1

@@ -0,0 +1,107 @@
+LAMMPS (7 Aug 2019)
+# 3d Lennard-Jones melt
+#
+# This example requires that the example models provided with
+# the kim-api package are installed.  see the ./lib/kim/README or
+# ./lib/kim/Install.py files for details on how to install these
+# example models.
+#
+
+variable	x index 1
+variable	y index 1
+variable	z index 1
+
+variable	xx equal 20*$x
+variable	xx equal 20*1
+variable	yy equal 20*$y
+variable	yy equal 20*1
+variable	zz equal 20*$z
+variable	zz equal 20*1
+
+kim_init	LennardJones_Ar real
+#=== BEGIN kim-init ==========================================
+units real
+#=== END kim-init ============================================
+
+
+lattice		fcc 4.4300
+Lattice spacing in x,y,z = 4.43 4.43 4.43
+region		box block 0 ${xx} 0 ${yy} 0 ${zz}
+region		box block 0 20 0 ${yy} 0 ${zz}
+region		box block 0 20 0 20 0 ${zz}
+region		box block 0 20 0 20 0 20
+create_box	1 box
+Created orthogonal box = (0 0 0) to (88.6 88.6 88.6)
+  1 by 1 by 1 MPI processor grid
+create_atoms	1 box
+Created 32000 atoms
+  create_atoms CPU = 0.004321 secs
+
+kim_interactions Ar
+#=== BEGIN kim_interactions ==================================
+pair_style kim LennardJones_Ar
+WARNING: KIM Model does not provide `partialParticleEnergy'; energy per atom will be zero (../pair_kim.cpp:974)
+WARNING: KIM Model does not provide `partialParticleVirial'; virial per atom will be zero (../pair_kim.cpp:979)
+pair_coeff * * Ar 
+#=== END kim_interactions ====================================
+
+
+mass		1 39.95
+velocity	all create 200.0 232345 loop geom
+
+neighbor	0.3 bin
+neigh_modify	delay 0 every 1 check yes
+
+fix		1 all nve
+#fix		1 all npt temp 1.0 1.0 1.0 iso 1.0 1.0 3.0
+
+run 		100
+Neighbor list info ...
+  update every 1 steps, delay 0 steps, check yes
+  max neighbors/atom: 2000, page size: 100000
+  master list distance cutoff = 8.45
+  ghost atom cutoff = 8.45
+  binsize = 4.225, bins = 21 21 21
+  1 neighbor lists, perpetual/occasional/extra = 1 0 0
+  (1) pair kim, perpetual
+      attributes: full, newton off, cut 8.45
+      pair build: full/bin/atomonly
+      stencil: full/bin/3d
+      bin: standard
+Setting up Verlet run ...
+  Unit style    : real
+  Current step  : 0
+  Time step     : 1
+Per MPI rank memory allocation (min/avg/max) = 28.12 | 28.12 | 28.12 Mbytes
+Step Temp E_pair E_mol TotEng Press 
+       0          200    145069.63            0    164146.22    128015.94 
+     100    95.179703    154939.42            0    164017.94    131602.75 
+Loop time of 3.48256 on 1 procs for 100 steps with 32000 atoms
+
+Performance: 2.481 ns/day, 9.674 hours/ns, 28.715 timesteps/s
+98.3% CPU use with 1 MPI tasks x no OpenMP threads
+
+MPI task timing breakdown:
+Section |  min time  |  avg time  |  max time  |%varavg| %total
+---------------------------------------------------------------
+Pair    | 3.0502     | 3.0502     | 3.0502     |   0.0 | 87.59
+Neigh   | 0.3646     | 0.3646     | 0.3646     |   0.0 | 10.47
+Comm    | 0.01783    | 0.01783    | 0.01783    |   0.0 |  0.51
+Output  | 6.8e-05    | 6.8e-05    | 6.8e-05    |   0.0 |  0.00
+Modify  | 0.034349   | 0.034349   | 0.034349   |   0.0 |  0.99
+Other   |            | 0.01547    |            |       |  0.44
+
+Nlocal:    32000 ave 32000 max 32000 min
+Histogram: 1 0 0 0 0 0 0 0 0 0
+Nghost:    19911 ave 19911 max 19911 min
+Histogram: 1 0 0 0 0 0 0 0 0 0
+Neighs:    0 ave 0 max 0 min
+Histogram: 1 0 0 0 0 0 0 0 0 0
+FullNghs:  4.25375e+06 ave 4.25375e+06 max 4.25375e+06 min
+Histogram: 1 0 0 0 0 0 0 0 0 0
+
+Total # of neighbors = 4253750
+Ave neighs/atom = 132.93
+Neighbor list builds = 3
+Dangerous builds = 0
+Total wall time: 0:00:03

examples/kim/log.7Aug19.in.kim-ex.melt.clang.4

@@ -0,0 +1,113 @@
+LAMMPS (7 Aug 2019)
+# 3d Lennard-Jones melt
+#
+# This example requires that the example models provided with
+# the kim-api package are installed.  see the ./lib/kim/README or
+# ./lib/kim/Install.py files for details on how to install these
+# example models.
+#
+
+variable	x index 1
+variable	y index 1
+variable	z index 1
+
+variable	xx equal 20*$x
+variable	xx equal 20*1
+variable	yy equal 20*$y
+variable	yy equal 20*1
+variable	zz equal 20*$z
+variable	zz equal 20*1
+
+kim_init	LennardJones_Ar real
+#=== BEGIN kim-init ==========================================
+units real
+#=== END kim-init ============================================
+
+
+lattice		fcc 4.4300
+Lattice spacing in x,y,z = 4.43 4.43 4.43
+region		box block 0 ${xx} 0 ${yy} 0 ${zz}
+region		box block 0 20 0 ${yy} 0 ${zz}
+region		box block 0 20 0 20 0 ${zz}
+region		box block 0 20 0 20 0 20
+create_box	1 box
+Created orthogonal box = (0 0 0) to (88.6 88.6 88.6)
+  1 by 2 by 2 MPI processor grid
+create_atoms	1 box
+Created 32000 atoms
+  create_atoms CPU = 0.000989 secs
+
+kim_interactions Ar
+#=== BEGIN kim_interactions ==================================
+pair_style kim LennardJones_Ar
+WARNING: KIM Model does not provide `partialParticleEnergy'; energy per atom will be zero (../pair_kim.cpp:974)
+WARNING: KIM Model does not provide `partialParticleVirial'; virial per atom will be zero (../pair_kim.cpp:979)
+pair_coeff * * Ar 
+WARNING: KIM Model does not provide `partialParticleEnergy'; energy per atom will be zero (../pair_kim.cpp:974)
+WARNING: KIM Model does not provide `partialParticleVirial'; virial per atom will be zero (../pair_kim.cpp:979)
+#=== END kim_interactions ====================================
+
+
+mass		1 39.95
+velocity	all create 200.0 232345 loop geom
+WARNING: KIM Model does not provide `partialParticleEnergy'; energy per atom will be zero (../pair_kim.cpp:974)
+WARNING: KIM Model does not provide `partialParticleVirial'; virial per atom will be zero (../pair_kim.cpp:979)
+WARNING: KIM Model does not provide `partialParticleEnergy'; energy per atom will be zero (../pair_kim.cpp:974)
+WARNING: KIM Model does not provide `partialParticleVirial'; virial per atom will be zero (../pair_kim.cpp:979)
+
+neighbor	0.3 bin
+neigh_modify	delay 0 every 1 check yes
+
+fix		1 all nve
+#fix		1 all npt temp 1.0 1.0 1.0 iso 1.0 1.0 3.0
+
+run 		100
+Neighbor list info ...
+  update every 1 steps, delay 0 steps, check yes
+  max neighbors/atom: 2000, page size: 100000
+  master list distance cutoff = 8.45
+  ghost atom cutoff = 8.45
+  binsize = 4.225, bins = 21 21 21
+  1 neighbor lists, perpetual/occasional/extra = 1 0 0
+  (1) pair kim, perpetual
+      attributes: full, newton off, cut 8.45
+      pair build: full/bin/atomonly
+      stencil: full/bin/3d
+      bin: standard
+Setting up Verlet run ...
+  Unit style    : real
+  Current step  : 0
+  Time step     : 1
+Per MPI rank memory allocation (min/avg/max) = 9.791 | 9.791 | 9.791 Mbytes
+Step Temp E_pair E_mol TotEng Press 
+       0          200    145069.63            0    164146.22    128015.94 
+     100    95.179703    154939.42            0    164017.94    131602.75 
+Loop time of 0.924494 on 4 procs for 100 steps with 32000 atoms
+
+Performance: 9.346 ns/day, 2.568 hours/ns, 108.167 timesteps/s
+99.6% CPU use with 4 MPI tasks x no OpenMP threads
+
+MPI task timing breakdown:
+Section |  min time  |  avg time  |  max time  |%varavg| %total
+---------------------------------------------------------------
+Pair    | 0.76434    | 0.76847    | 0.77207    |   0.3 | 83.12
+Neigh   | 0.09089    | 0.094446   | 0.099911   |   1.1 | 10.22
+Comm    | 0.038599   | 0.044759   | 0.051381   |   2.1 |  4.84
+Output  | 3.5e-05    | 4e-05      | 4.9e-05    |   0.0 |  0.00
+Modify  | 0.009396   | 0.009685   | 0.009941   |   0.2 |  1.05
+Other   |            | 0.00709    |            |       |  0.77
+
+Nlocal:    8000 ave 8018 max 7967 min
+Histogram: 1 0 0 0 0 0 1 0 0 2
+Nghost:    9131 ave 9164 max 9113 min
+Histogram: 2 0 0 1 0 0 0 0 0 1
+Neighs:    0 ave 0 max 0 min
+Histogram: 4 0 0 0 0 0 0 0 0 0
+FullNghs:  1.06344e+06 ave 1.06594e+06 max 1.05881e+06 min
+Histogram: 1 0 0 0 0 0 1 0 0 2
+
+Total # of neighbors = 4253750
+Ave neighs/atom = 132.93
+Neighbor list builds = 3
+Dangerous builds = 0
+Total wall time: 0:00:00

examples/kim/log.7Aug19.in.kim-pm-query.melt.clang.1

@@ -0,0 +1,124 @@
+LAMMPS (7 Aug 2019)
+# 3d Lennard-Jones melt
+#
+# This example requires that the KIM Portable Model (PM)
+# SW_StillingerWeber_1985_Si__MO_405512056662_005
+# is installed.  This can be done with the command
+#   kim-api-collections-management install user SW_StillingerWeber_1985_Si__MO_405512056662_005
+# If this command does not work, you may need to setup your PATH to find the utility.
+# If you installed the kim-api using the LAMMPS CMake build, you can do the following
+# (where the current working directory is assumed to be the LAMMPS build directory)
+#   source ./kim_build-prefix/bin/kim-api-activate
+# If you installed the kim-api using the LAMMPS Make build, you can do the following
+# (where the current working directory is assumed to be the LAMMPS src directory)
+#   source ../lib/kim/installed-kim-api-X.Y.Z/bin/kim-api-activate
+# (where you should relplace X.Y.Z with the appropriate kim-api version number).
+#
+# Or, see https://openkim.org/doc/obtaining-models for alternative options.
+#
+
+variable	x index 1
+variable	y index 1
+variable	z index 1
+
+variable	xx equal 20*$x
+variable	xx equal 20*1
+variable	yy equal 20*$y
+variable	yy equal 20*1
+variable	zz equal 20*$z
+variable	zz equal 20*1
+
+kim_init	SW_StillingerWeber_1985_Si__MO_405512056662_005 real
+#=== BEGIN kim-init ==========================================
+units real
+#=== END kim-init ============================================
+
+kim_query       a0 get_lattice_constant_cubic crystal=["fcc"] species=["Si"] units=["angstrom"]
+#=== BEGIN kim-query =========================================
+variable a0 string 4.146581932902336
+#=== END kim-query ===========================================
+
+
+lattice		fcc ${a0}
+lattice		fcc 4.146581932902336
+Lattice spacing in x,y,z = 4.14658 4.14658 4.14658
+region		box block 0 ${xx} 0 ${yy} 0 ${zz}
+region		box block 0 20 0 ${yy} 0 ${zz}
+region		box block 0 20 0 20 0 ${zz}
+region		box block 0 20 0 20 0 20
+create_box	1 box
+Created orthogonal box = (0 0 0) to (82.9316 82.9316 82.9316)
+  1 by 1 by 1 MPI processor grid
+create_atoms	1 box
+Created 32000 atoms
+  create_atoms CPU = 0.005415 secs
+
+kim_interactions Si
+#=== BEGIN kim_interactions ==================================
+pair_style kim SW_StillingerWeber_1985_Si__MO_405512056662_005
+pair_coeff * * Si 
+#=== END kim_interactions ====================================
+
+
+mass		1 39.95
+velocity	all create 200.0 232345 loop geom
+
+neighbor	0.3 bin
+neigh_modify	delay 0 every 1 check yes
+
+fix		1 all nve
+#fix		1 all npt temp 1.0 1.0 1.0 iso 1.0 1.0 3.0
+
+run 		100
+Neighbor list info ...
+  update every 1 steps, delay 0 steps, check yes
+  max neighbors/atom: 2000, page size: 100000
+  master list distance cutoff = 4.07118
+  ghost atom cutoff = 4.07118
+  binsize = 2.03559, bins = 41 41 41
+  1 neighbor lists, perpetual/occasional/extra = 1 0 0
+  (1) pair kim, perpetual
+      attributes: full, newton off, cut 4.07118
+      pair build: full/bin/atomonly
+      stencil: full/bin/3d
+      bin: standard
+Setting up Verlet run ...
+  Unit style    : real
+  Current step  : 0
+  Time step     : 1
+Per MPI rank memory allocation (min/avg/max) = 10.36 | 10.36 | 10.36 Mbytes
+Step Temp E_pair E_mol TotEng Press 
+       0          200   -126084.25            0   -107007.66    1528.8768 
+     100    94.450495   -116016.03            0   -107007.07    2282.2685 
+Loop time of 74.6055 on 1 procs for 100 steps with 32000 atoms
+
+Performance: 0.116 ns/day, 207.238 hours/ns, 1.340 timesteps/s
+98.6% CPU use with 1 MPI tasks x no OpenMP threads
+
+MPI task timing breakdown:
+Section |  min time  |  avg time  |  max time  |%varavg| %total
+---------------------------------------------------------------
+Pair    | 74.446     | 74.446     | 74.446     |   0.0 | 99.79
+Neigh   | 0.096611   | 0.096611   | 0.096611   |   0.0 |  0.13
+Comm    | 0.014594   | 0.014594   | 0.014594   |   0.0 |  0.02
+Output  | 7.9e-05    | 7.9e-05    | 7.9e-05    |   0.0 |  0.00
+Modify  | 0.03454    | 0.03454    | 0.03454    |   0.0 |  0.05
+Other   |            | 0.01396    |            |       |  0.02
+
+Nlocal:    32000 ave 32000 max 32000 min
+Histogram: 1 0 0 0 0 0 0 0 0 0
+Nghost:    9667 ave 9667 max 9667 min
+Histogram: 1 0 0 0 0 0 0 0 0 0
+Neighs:    0 ave 0 max 0 min
+Histogram: 1 0 0 0 0 0 0 0 0 0
+FullNghs:  450192 ave 450192 max 450192 min
+Histogram: 1 0 0 0 0 0 0 0 0 0
+
+Total # of neighbors = 450192
+Ave neighs/atom = 14.0685
+Neighbor list builds = 3
+Dangerous builds = 0
+
+Please see the log.cite file for references relevant to this simulation
+
+Total wall time: 0:01:16

examples/kim/log.7Aug19.in.kim-pm-query.melt.clang.4

@@ -0,0 +1,124 @@
+LAMMPS (7 Aug 2019)
+# 3d Lennard-Jones melt
+#
+# This example requires that the KIM Portable Model (PM)
+# SW_StillingerWeber_1985_Si__MO_405512056662_005
+# is installed.  This can be done with the command
+#   kim-api-collections-management install user SW_StillingerWeber_1985_Si__MO_405512056662_005
+# If this command does not work, you may need to setup your PATH to find the utility.
+# If you installed the kim-api using the LAMMPS CMake build, you can do the following
+# (where the current working directory is assumed to be the LAMMPS build directory)
+#   source ./kim_build-prefix/bin/kim-api-activate
+# If you installed the kim-api using the LAMMPS Make build, you can do the following
+# (where the current working directory is assumed to be the LAMMPS src directory)
+#   source ../lib/kim/installed-kim-api-X.Y.Z/bin/kim-api-activate
+# (where you should relplace X.Y.Z with the appropriate kim-api version number).
+#
+# Or, see https://openkim.org/doc/obtaining-models for alternative options.
+#
+
+variable	x index 1
+variable	y index 1
+variable	z index 1
+
+variable	xx equal 20*$x
+variable	xx equal 20*1
+variable	yy equal 20*$y
+variable	yy equal 20*1
+variable	zz equal 20*$z
+variable	zz equal 20*1
+
+kim_init	SW_StillingerWeber_1985_Si__MO_405512056662_005 real
+#=== BEGIN kim-init ==========================================
+units real
+#=== END kim-init ============================================
+
+kim_query       a0 get_lattice_constant_cubic crystal=["fcc"] species=["Si"] units=["angstrom"]
+#=== BEGIN kim-query =========================================
+variable a0 string 4.146581932902336
+#=== END kim-query ===========================================
+
+
+lattice		fcc ${a0}
+lattice		fcc 4.146581932902336
+Lattice spacing in x,y,z = 4.14658 4.14658 4.14658
+region		box block 0 ${xx} 0 ${yy} 0 ${zz}
+region		box block 0 20 0 ${yy} 0 ${zz}
+region		box block 0 20 0 20 0 ${zz}
+region		box block 0 20 0 20 0 20
+create_box	1 box
+Created orthogonal box = (0 0 0) to (82.9316 82.9316 82.9316)
+  1 by 2 by 2 MPI processor grid
+create_atoms	1 box
+Created 32000 atoms
+  create_atoms CPU = 0.000946 secs
+
+kim_interactions Si
+#=== BEGIN kim_interactions ==================================
+pair_style kim SW_StillingerWeber_1985_Si__MO_405512056662_005
+pair_coeff * * Si 
+#=== END kim_interactions ====================================
+
+
+mass		1 39.95
+velocity	all create 200.0 232345 loop geom
+
+neighbor	0.3 bin
+neigh_modify	delay 0 every 1 check yes
+
+fix		1 all nve
+#fix		1 all npt temp 1.0 1.0 1.0 iso 1.0 1.0 3.0
+
+run 		100
+Neighbor list info ...
+  update every 1 steps, delay 0 steps, check yes
+  max neighbors/atom: 2000, page size: 100000
+  master list distance cutoff = 4.07118
+  ghost atom cutoff = 4.07118
+  binsize = 2.03559, bins = 41 41 41
+  1 neighbor lists, perpetual/occasional/extra = 1 0 0
+  (1) pair kim, perpetual
+      attributes: full, newton off, cut 4.07118
+      pair build: full/bin/atomonly
+      stencil: full/bin/3d
+      bin: standard
+Setting up Verlet run ...
+  Unit style    : real
+  Current step  : 0
+  Time step     : 1
+Per MPI rank memory allocation (min/avg/max) = 3.489 | 3.489 | 3.489 Mbytes
+Step Temp E_pair E_mol TotEng Press 
+       0          200   -126084.25            0   -107007.66    1528.8768 
+     100    94.450495   -116016.03            0   -107007.07    2282.2685 
+Loop time of 19.0792 on 4 procs for 100 steps with 32000 atoms
+
+Performance: 0.453 ns/day, 52.998 hours/ns, 5.241 timesteps/s
+99.4% CPU use with 4 MPI tasks x no OpenMP threads
+
+MPI task timing breakdown:
+Section |  min time  |  avg time  |  max time  |%varavg| %total
+---------------------------------------------------------------
+Pair    | 18.78      | 18.855     | 18.937     |   1.5 | 98.83
+Neigh   | 0.026047   | 0.026274   | 0.0266     |   0.1 |  0.14
+Comm    | 0.09039    | 0.17196    | 0.24675    |  15.9 |  0.90
+Output  | 3.9e-05    | 4.975e-05  | 6.1e-05    |   0.0 |  0.00
+Modify  | 0.015667   | 0.015819   | 0.016008   |   0.1 |  0.08
+Other   |            | 0.01008    |            |       |  0.05
+
+Nlocal:    8000 ave 8029 max 7968 min
+Histogram: 1 1 0 0 0 0 0 0 0 2
+Nghost:    4259 ave 4303 max 4202 min
+Histogram: 1 0 0 0 0 0 2 0 0 1
+Neighs:    0 ave 0 max 0 min
+Histogram: 4 0 0 0 0 0 0 0 0 0
+FullNghs:  112548 ave 113091 max 111995 min
+Histogram: 1 0 0 1 0 0 0 1 0 1
+
+Total # of neighbors = 450192
+Ave neighs/atom = 14.0685
+Neighbor list builds = 3
+Dangerous builds = 0
+
+Please see the log.cite file for references relevant to this simulation
+
+Total wall time: 0:00:20

examples/kim/log.7Aug19.in.kim-pm.melt.clang.1

@@ -0,0 +1,118 @@
+LAMMPS (7 Aug 2019)
+# 3d Lennard-Jones melt
+#
+# This example requires that the KIM Portable Model (PM)
+# SW_StillingerWeber_1985_Si__MO_405512056662_005
+# is installed.  This can be done with the command
+#   kim-api-collections-management install user SW_StillingerWeber_1985_Si__MO_405512056662_005
+# If this command does not work, you may need to setup your PATH to find the utility.
+# If you installed the kim-api using the LAMMPS CMake build, you can do the following
+# (where the current working directory is assumed to be the LAMMPS build directory)
+#   source ./kim_build-prefix/bin/kim-api-activate
+# If you installed the kim-api using the LAMMPS Make build, you can do the following
+# (where the current working directory is assumed to be the LAMMPS src directory)
+#   source ../lib/kim/installed-kim-api-X.Y.Z/bin/kim-api-activate
+# (where you should relplace X.Y.Z with the appropriate kim-api version number).
+#
+# Or, see https://openkim.org/doc/obtaining-models for alternative options.
+#
+
+variable	x index 1
+variable	y index 1
+variable	z index 1
+
+variable	xx equal 20*$x
+variable	xx equal 20*1
+variable	yy equal 20*$y
+variable	yy equal 20*1
+variable	zz equal 20*$z
+variable	zz equal 20*1
+
+kim_init	SW_StillingerWeber_1985_Si__MO_405512056662_005 real
+#=== BEGIN kim-init ==========================================
+units real
+#=== END kim-init ============================================
+
+
+lattice		fcc 4.4300
+Lattice spacing in x,y,z = 4.43 4.43 4.43
+region		box block 0 ${xx} 0 ${yy} 0 ${zz}
+region		box block 0 20 0 ${yy} 0 ${zz}
+region		box block 0 20 0 20 0 ${zz}
+region		box block 0 20 0 20 0 20
+create_box	1 box
+Created orthogonal box = (0 0 0) to (88.6 88.6 88.6)
+  1 by 1 by 1 MPI processor grid
+create_atoms	1 box
+Created 32000 atoms
+  create_atoms CPU = 0.003591 secs
+
+kim_interactions Si
+#=== BEGIN kim_interactions ==================================
+pair_style kim SW_StillingerWeber_1985_Si__MO_405512056662_005
+pair_coeff * * Si 
+#=== END kim_interactions ====================================
+
+
+mass		1 39.95
+velocity	all create 200.0 232345 loop geom
+
+neighbor	0.3 bin
+neigh_modify	delay 0 every 1 check yes
+
+fix		1 all nve
+#fix		1 all npt temp 1.0 1.0 1.0 iso 1.0 1.0 3.0
+
+run 		100
+Neighbor list info ...
+  update every 1 steps, delay 0 steps, check yes
+  max neighbors/atom: 2000, page size: 100000
+  master list distance cutoff = 4.07118
+  ghost atom cutoff = 4.07118
+  binsize = 2.03559, bins = 44 44 44
+  1 neighbor lists, perpetual/occasional/extra = 1 0 0
+  (1) pair kim, perpetual
+      attributes: full, newton off, cut 4.07118
+      pair build: full/bin/atomonly
+      stencil: full/bin/3d
+      bin: standard
+Setting up Verlet run ...
+  Unit style    : real
+  Current step  : 0
+  Time step     : 1
+Per MPI rank memory allocation (min/avg/max) = 10.44 | 10.44 | 10.44 Mbytes
+Step Temp E_pair E_mol TotEng Press 
+       0          200   -85249.847            0   -66173.259   -33302.387 
+     100    253.43357    -90346.68            0   -66173.441   -14888.698 
+Loop time of 74.248 on 1 procs for 100 steps with 32000 atoms
+
+Performance: 0.116 ns/day, 206.244 hours/ns, 1.347 timesteps/s
+98.8% CPU use with 1 MPI tasks x no OpenMP threads
+
+MPI task timing breakdown:
+Section |  min time  |  avg time  |  max time  |%varavg| %total
+---------------------------------------------------------------
+Pair    | 74.118     | 74.118     | 74.118     |   0.0 | 99.83
+Neigh   | 0.069623   | 0.069623   | 0.069623   |   0.0 |  0.09
+Comm    | 0.0137     | 0.0137     | 0.0137     |   0.0 |  0.02
+Output  | 7.6e-05    | 7.6e-05    | 7.6e-05    |   0.0 |  0.00
+Modify  | 0.031883   | 0.031883   | 0.031883   |   0.0 |  0.04
+Other   |            | 0.01433    |            |       |  0.02
+
+Nlocal:    32000 ave 32000 max 32000 min
+Histogram: 1 0 0 0 0 0 0 0 0 0
+Nghost:    7760 ave 7760 max 7760 min
+Histogram: 1 0 0 0 0 0 0 0 0 0
+Neighs:    0 ave 0 max 0 min
+Histogram: 1 0 0 0 0 0 0 0 0 0
+FullNghs:  402352 ave 402352 max 402352 min
+Histogram: 1 0 0 0 0 0 0 0 0 0
+
+Total # of neighbors = 402352
+Ave neighs/atom = 12.5735
+Neighbor list builds = 4
+Dangerous builds = 0
+
+Please see the log.cite file for references relevant to this simulation
+
+Total wall time: 0:01:14

examples/kim/log.7Aug19.in.kim-pm.melt.clang.4

@@ -0,0 +1,118 @@
+LAMMPS (7 Aug 2019)
+# 3d Lennard-Jones melt
+#
+# This example requires that the KIM Portable Model (PM)
+# SW_StillingerWeber_1985_Si__MO_405512056662_005
+# is installed.  This can be done with the command
+#   kim-api-collections-management install user SW_StillingerWeber_1985_Si__MO_405512056662_005
+# If this command does not work, you may need to setup your PATH to find the utility.
+# If you installed the kim-api using the LAMMPS CMake build, you can do the following
+# (where the current working directory is assumed to be the LAMMPS build directory)
+#   source ./kim_build-prefix/bin/kim-api-activate
+# If you installed the kim-api using the LAMMPS Make build, you can do the following
+# (where the current working directory is assumed to be the LAMMPS src directory)
+#   source ../lib/kim/installed-kim-api-X.Y.Z/bin/kim-api-activate
+# (where you should relplace X.Y.Z with the appropriate kim-api version number).
+#
+# Or, see https://openkim.org/doc/obtaining-models for alternative options.
+#
+
+variable	x index 1
+variable	y index 1
+variable	z index 1
+
+variable	xx equal 20*$x
+variable	xx equal 20*1
+variable	yy equal 20*$y
+variable	yy equal 20*1
+variable	zz equal 20*$z
+variable	zz equal 20*1
+
+kim_init	SW_StillingerWeber_1985_Si__MO_405512056662_005 real
+#=== BEGIN kim-init ==========================================
+units real
+#=== END kim-init ============================================
+
+
+lattice		fcc 4.4300
+Lattice spacing in x,y,z = 4.43 4.43 4.43
+region		box block 0 ${xx} 0 ${yy} 0 ${zz}
+region		box block 0 20 0 ${yy} 0 ${zz}
+region		box block 0 20 0 20 0 ${zz}
+region		box block 0 20 0 20 0 20
+create_box	1 box
+Created orthogonal box = (0 0 0) to (88.6 88.6 88.6)
+  1 by 2 by 2 MPI processor grid
+create_atoms	1 box
+Created 32000 atoms
+  create_atoms CPU = 0.000997 secs
+
+kim_interactions Si
+#=== BEGIN kim_interactions ==================================
+pair_style kim SW_StillingerWeber_1985_Si__MO_405512056662_005
+pair_coeff * * Si 
+#=== END kim_interactions ====================================
+
+
+mass		1 39.95
+velocity	all create 200.0 232345 loop geom
+
+neighbor	0.3 bin
+neigh_modify	delay 0 every 1 check yes
+
+fix		1 all nve
+#fix		1 all npt temp 1.0 1.0 1.0 iso 1.0 1.0 3.0
+
+run 		100
+Neighbor list info ...
+  update every 1 steps, delay 0 steps, check yes
+  max neighbors/atom: 2000, page size: 100000
+  master list distance cutoff = 4.07118
+  ghost atom cutoff = 4.07118
+  binsize = 2.03559, bins = 44 44 44
+  1 neighbor lists, perpetual/occasional/extra = 1 0 0
+  (1) pair kim, perpetual
+      attributes: full, newton off, cut 4.07118
+      pair build: full/bin/atomonly
+      stencil: full/bin/3d
+      bin: standard
+Setting up Verlet run ...
+  Unit style    : real
+  Current step  : 0
+  Time step     : 1
+Per MPI rank memory allocation (min/avg/max) = 3.517 | 3.517 | 3.517 Mbytes
+Step Temp E_pair E_mol TotEng Press 
+       0          200   -85249.847            0   -66173.259   -33302.387 
+     100    253.43357    -90346.68            0   -66173.441   -14888.698 
+Loop time of 19.0287 on 4 procs for 100 steps with 32000 atoms
+
+Performance: 0.454 ns/day, 52.857 hours/ns, 5.255 timesteps/s
+99.1% CPU use with 4 MPI tasks x no OpenMP threads
+
+MPI task timing breakdown:
+Section |  min time  |  avg time  |  max time  |%varavg| %total
+---------------------------------------------------------------
+Pair    | 18.81      | 18.838     | 18.883     |   0.6 | 99.00
+Neigh   | 0.018598   | 0.01914    | 0.020732   |   0.7 |  0.10
+Comm    | 0.10341    | 0.1475     | 0.17393    |   7.1 |  0.78
+Output  | 6e-05      | 6.225e-05  | 6.7e-05    |   0.0 |  0.00
+Modify  | 0.014839   | 0.014925   | 0.015047   |   0.1 |  0.08
+Other   |            | 0.008997   |            |       |  0.05
+
+Nlocal:    8000 ave 8014 max 7988 min
+Histogram: 1 1 0 0 0 0 1 0 0 1
+Nghost:    3374.75 ave 3389 max 3361 min
+Histogram: 1 0 1 0 0 0 0 1 0 1
+Neighs:    0 ave 0 max 0 min
+Histogram: 4 0 0 0 0 0 0 0 0 0
+FullNghs:  100588 ave 100856 max 100392 min
+Histogram: 1 0 1 0 1 0 0 0 0 1
+
+Total # of neighbors = 402352
+Ave neighs/atom = 12.5735
+Neighbor list builds = 4
+Dangerous builds = 0
+
+Please see the log.cite file for references relevant to this simulation
+
+Total wall time: 0:00:19

examples/kim/log.7Aug19.in.kim-sm.melt.clang.1

@@ -0,0 +1,71 @@
+LAMMPS (7 Aug 2019)
+# 3d Lennard-Jones melt
+#
+# This example requires that the KIM Simulator Model (PM)
+# Sim_LAMMPS_ReaxFF_StrachanVanDuinChakraborty_2003_CHNO__SM_107643900657_000
+# is installed.  This can be done with the command
+#   kim-api-collections-management install user Sim_LAMMPS_ReaxFF_StrachanVanDuinChakraborty_2003_CHNO__SM_107643900657_000
+# If this command does not work, you may need to setup your PATH to find the utility.
+# If you installed the kim-api using the LAMMPS CMake build, you can do the following
+# (where the current working directory is assumed to be the LAMMPS build directory)
+#   source ./kim_build-prefix/bin/kim-api-activate
+# If you installed the kim-api using the LAMMPS Make build, you can do the following
+# (where the current working directory is assumed to be the LAMMPS src directory)
+#   source ../lib/kim/installed-kim-api-X.Y.Z/bin/kim-api-activate
+# (where you should relplace X.Y.Z with the appropriate kim-api version number).
+#
+# See https://openkim.org/doc/obtaining-models for alternative options.
+#
+
+variable	x index 1
+variable	y index 1
+variable	z index 1
+
+variable	xx equal 20*$x
+variable	xx equal 20*1
+variable	yy equal 20*$y
+variable	yy equal 20*1
+variable	zz equal 20*$z
+variable	zz equal 20*1
+
+kim_init	Sim_LAMMPS_ReaxFF_StrachanVanDuinChakraborty_2003_CHNO__SM_107643900657_000 real
+#=== BEGIN kim-init ==========================================
+# Using KIM Simulator Model : Sim_LAMMPS_ReaxFF_StrachanVanDuinChakraborty_2003_CHNO__SM_107643900657_000
+# For Simulator             : LAMMPS 28 Feb 2019
+# Running on                : LAMMPS 7 Aug 2019
+#
+units real
+atom_style charge
+neigh_modify one 4000
+#=== END kim-init ============================================
+
+
+lattice		fcc 4.4300
+Lattice spacing in x,y,z = 4.43 4.43 4.43
+region		box block 0 ${xx} 0 ${yy} 0 ${zz}
+region		box block 0 20 0 ${yy} 0 ${zz}
+region		box block 0 20 0 20 0 ${zz}
+region		box block 0 20 0 20 0 20
+create_box	1 box
+Created orthogonal box = (0 0 0) to (88.6 88.6 88.6)
+  1 by 1 by 1 MPI processor grid
+create_atoms	1 box
+Created 32000 atoms
+  create_atoms CPU = 0.003447 secs
+
+kim_interactions O
+#=== BEGIN kim_interactions ==================================
+pair_style reax/c /var/tmp/kim-simulator-model-parameter-file-directory-6Acs1QDbXgBx/lmp_control safezone 2.0 mincap 100
+ERROR: Unrecognized pair style 'reax/c' is part of the USER-REAXC package which is not enabled in this LAMMPS binary. (../force.cpp:262)
+Last command: pair_style reax/c /var/tmp/kim-simulator-model-parameter-file-directory-6Acs1QDbXgBx/lmp_control safezone 2.0 mincap 100
+--------------------------------------------------------------------------
+Primary job  terminated normally, but 1 process returned
+a non-zero exit code. Per user-direction, the job has been aborted.
+--------------------------------------------------------------------------
+--------------------------------------------------------------------------
+mpirun detected that one or more processes exited with non-zero status, thus causing
+the job to be terminated. The first process to do so was:
+
+  Process name: [[33054,1],0]
+  Exit code:    1
+--------------------------------------------------------------------------

examples/kim/log.7Aug19.in.kim-sm.melt.clang.4

@@ -0,0 +1,60 @@
+LAMMPS (7 Aug 2019)
+# 3d Lennard-Jones melt
+#
+# This example requires that the KIM Simulator Model (PM)
+# Sim_LAMMPS_ReaxFF_StrachanVanDuinChakraborty_2003_CHNO__SM_107643900657_000
+# is installed.  This can be done with the command
+#   kim-api-collections-management install user Sim_LAMMPS_ReaxFF_StrachanVanDuinChakraborty_2003_CHNO__SM_107643900657_000
+# If this command does not work, you may need to setup your PATH to find the utility.
+# If you installed the kim-api using the LAMMPS CMake build, you can do the following
+# (where the current working directory is assumed to be the LAMMPS build directory)
+#   source ./kim_build-prefix/bin/kim-api-activate
+# If you installed the kim-api using the LAMMPS Make build, you can do the following
+# (where the current working directory is assumed to be the LAMMPS src directory)
+#   source ../lib/kim/installed-kim-api-X.Y.Z/bin/kim-api-activate
+# (where you should relplace X.Y.Z with the appropriate kim-api version number).
+#
+# See https://openkim.org/doc/obtaining-models for alternative options.
+#
+
+variable	x index 1
+variable	y index 1
+variable	z index 1
+
+variable	xx equal 20*$x
+variable	xx equal 20*1
+variable	yy equal 20*$y
+variable	yy equal 20*1
+variable	zz equal 20*$z
+variable	zz equal 20*1
+
+kim_init	Sim_LAMMPS_ReaxFF_StrachanVanDuinChakraborty_2003_CHNO__SM_107643900657_000 real
+#=== BEGIN kim-init ==========================================
+# Using KIM Simulator Model : Sim_LAMMPS_ReaxFF_StrachanVanDuinChakraborty_2003_CHNO__SM_107643900657_000
+# For Simulator             : LAMMPS 28 Feb 2019
+# Running on                : LAMMPS 7 Aug 2019
+#
+units real
+atom_style charge
+neigh_modify one 4000
+#=== END kim-init ============================================
+
+
+lattice		fcc 4.4300
+Lattice spacing in x,y,z = 4.43 4.43 4.43
+region		box block 0 ${xx} 0 ${yy} 0 ${zz}
+region		box block 0 20 0 ${yy} 0 ${zz}
+region		box block 0 20 0 20 0 ${zz}
+region		box block 0 20 0 20 0 20
+create_box	1 box
+Created orthogonal box = (0 0 0) to (88.6 88.6 88.6)
+  1 by 2 by 2 MPI processor grid
+create_atoms	1 box
+Created 32000 atoms
+  create_atoms CPU = 0.001307 secs
+
+kim_interactions O
+#=== BEGIN kim_interactions ==================================
+pair_style reax/c /var/tmp/kim-simulator-model-parameter-file-directory-6tmKtZEXzhgv/lmp_control safezone 2.0 mincap 100
+ERROR: Unrecognized pair style 'reax/c' is part of the USER-REAXC package which is not enabled in this LAMMPS binary. (../force.cpp:262)
+Last command: pair_style reax/c /var/tmp/kim-simulator-model-parameter-file-directory-6tmKtZEXzhgv/lmp_control safezone 2.0 mincap 100

examples/kim/log.7Aug19.in.lammps.melt.clang.1

@@ -0,0 +1,92 @@
+LAMMPS (7 Aug 2019)
+# 3d Lennard-Jones melt
+
+variable	x index 1
+variable	y index 1
+variable	z index 1
+
+variable	xx equal 20*$x
+variable	xx equal 20*1
+variable	yy equal 20*$y
+variable	yy equal 20*1
+variable	zz equal 20*$z
+variable	zz equal 20*1
+
+units		real
+
+lattice		fcc 4.4300
+Lattice spacing in x,y,z = 4.43 4.43 4.43
+region		box block 0 ${xx} 0 ${yy} 0 ${zz}
+region		box block 0 20 0 ${yy} 0 ${zz}
+region		box block 0 20 0 20 0 ${zz}
+region		box block 0 20 0 20 0 20
+create_box	1 box
+Created orthogonal box = (0 0 0) to (88.6 88.6 88.6)
+  1 by 1 by 1 MPI processor grid
+create_atoms	1 box
+Created 32000 atoms
+  create_atoms CPU = 0.003037 secs
+
+pair_style	lj/cut 8.1500
+pair_coeff	1 1 0.0104 3.4000
+
+#pair_style      kim LennardJones_Ar
+#pair_coeff      * * Ar
+
+mass		1 39.95
+velocity	all create 200.0 232345 loop geom
+
+neighbor	0.3 bin
+neigh_modify	delay 0 every 1 check yes
+
+fix		1 all nve
+#fix		1 all npt temp 1.0 1.0 1.0 iso 1.0 1.0 3.0
+
+run 		100
+Neighbor list info ...
+  update every 1 steps, delay 0 steps, check yes
+  max neighbors/atom: 2000, page size: 100000
+  master list distance cutoff = 8.45
+  ghost atom cutoff = 8.45
+  binsize = 4.225, bins = 21 21 21
+  1 neighbor lists, perpetual/occasional/extra = 1 0 0
+  (1) pair lj/cut, perpetual
+      attributes: half, newton on
+      pair build: half/bin/atomonly/newton
+      stencil: half/bin/3d/newton
+      bin: standard
+Setting up Verlet run ...
+  Unit style    : real
+  Current step  : 0
+  Time step     : 1
+Per MPI rank memory allocation (min/avg/max) = 19.23 | 19.23 | 19.23 Mbytes
+Step Temp E_pair E_mol TotEng Press 
+       0          200    6290.8194            0    25367.408    6750.7421 
+     100    98.747096    15900.676            0    25319.465    10184.453 
+Loop time of 2.43768 on 1 procs for 100 steps with 32000 atoms
+
+Performance: 3.544 ns/day, 6.771 hours/ns, 41.023 timesteps/s
+97.8% CPU use with 1 MPI tasks x no OpenMP threads
+
+MPI task timing breakdown:
+Section |  min time  |  avg time  |  max time  |%varavg| %total
+---------------------------------------------------------------
+Pair    | 2.1895     | 2.1895     | 2.1895     |   0.0 | 89.82
+Neigh   | 0.17546    | 0.17546    | 0.17546    |   0.0 |  7.20
+Comm    | 0.021001   | 0.021001   | 0.021001   |   0.0 |  0.86
+Output  | 7.9e-05    | 7.9e-05    | 7.9e-05    |   0.0 |  0.00
+Modify  | 0.034253   | 0.034253   | 0.034253   |   0.0 |  1.41
+Other   |            | 0.01735    |            |       |  0.71
+
+Nlocal:    32000 ave 32000 max 32000 min
+Histogram: 1 0 0 0 0 0 0 0 0 0
+Nghost:    19911 ave 19911 max 19911 min
+Histogram: 1 0 0 0 0 0 0 0 0 0
+Neighs:    1.96027e+06 ave 1.96027e+06 max 1.96027e+06 min
+Histogram: 1 0 0 0 0 0 0 0 0 0
+
+Total # of neighbors = 1960266
+Ave neighs/atom = 61.2583
+Neighbor list builds = 3
+Dangerous builds = 0
+Total wall time: 0:00:02

examples/kim/log.7Aug19.in.lammps.melt.clang.4

@@ -0,0 +1,92 @@
+LAMMPS (7 Aug 2019)
+# 3d Lennard-Jones melt
+
+variable	x index 1
+variable	y index 1
+variable	z index 1
+
+variable	xx equal 20*$x
+variable	xx equal 20*1
+variable	yy equal 20*$y
+variable	yy equal 20*1
+variable	zz equal 20*$z
+variable	zz equal 20*1
+
+units		real
+
+lattice		fcc 4.4300
+Lattice spacing in x,y,z = 4.43 4.43 4.43
+region		box block 0 ${xx} 0 ${yy} 0 ${zz}
+region		box block 0 20 0 ${yy} 0 ${zz}
+region		box block 0 20 0 20 0 ${zz}
+region		box block 0 20 0 20 0 20
+create_box	1 box
+Created orthogonal box = (0 0 0) to (88.6 88.6 88.6)
+  1 by 2 by 2 MPI processor grid
+create_atoms	1 box
+Created 32000 atoms
+  create_atoms CPU = 0.001194 secs
+
+pair_style	lj/cut 8.1500
+pair_coeff	1 1 0.0104 3.4000
+
+#pair_style      kim LennardJones_Ar
+#pair_coeff      * * Ar
+
+mass		1 39.95
+velocity	all create 200.0 232345 loop geom
+
+neighbor	0.3 bin
+neigh_modify	delay 0 every 1 check yes
+
+fix		1 all nve
+#fix		1 all npt temp 1.0 1.0 1.0 iso 1.0 1.0 3.0
+
+run 		100
+Neighbor list info ...
+  update every 1 steps, delay 0 steps, check yes
+  max neighbors/atom: 2000, page size: 100000
+  master list distance cutoff = 8.45
+  ghost atom cutoff = 8.45
+  binsize = 4.225, bins = 21 21 21
+  1 neighbor lists, perpetual/occasional/extra = 1 0 0
+  (1) pair lj/cut, perpetual
+      attributes: half, newton on
+      pair build: half/bin/atomonly/newton
+      stencil: half/bin/3d/newton
+      bin: standard
+Setting up Verlet run ...
+  Unit style    : real
+  Current step  : 0
+  Time step     : 1
+Per MPI rank memory allocation (min/avg/max) = 7.633 | 7.633 | 7.633 Mbytes
+Step Temp E_pair E_mol TotEng Press 
+       0          200    6290.8194            0    25367.408    6750.7421 
+     100    98.747096    15900.676            0    25319.465    10184.453 
+Loop time of 0.726239 on 4 procs for 100 steps with 32000 atoms
+
+Performance: 11.897 ns/day, 2.017 hours/ns, 137.696 timesteps/s
+98.7% CPU use with 4 MPI tasks x no OpenMP threads
+
+MPI task timing breakdown:
+Section |  min time  |  avg time  |  max time  |%varavg| %total
+---------------------------------------------------------------
+Pair    | 0.57617    | 0.5835     | 0.59084    |   0.9 | 80.34
+Neigh   | 0.046682   | 0.047783   | 0.048641   |   0.3 |  6.58
+Comm    | 0.065469   | 0.071509   | 0.07899    |   2.3 |  9.85
+Output  | 3.9e-05    | 4.6e-05    | 6.1e-05    |   0.0 |  0.01
+Modify  | 0.013205   | 0.01363    | 0.014044   |   0.3 |  1.88
+Other   |            | 0.009775   |            |       |  1.35
+
+Nlocal:    8000 ave 8012 max 7989 min
+Histogram: 1 0 0 0 2 0 0 0 0 1
+Nghost:    9131 ave 9142 max 9119 min
+Histogram: 1 0 0 0 0 2 0 0 0 1
+Neighs:    490066 ave 491443 max 489273 min
+Histogram: 2 0 0 0 1 0 0 0 0 1
+
+Total # of neighbors = 1960266
+Ave neighs/atom = 61.2583
+Neighbor list builds = 3
+Dangerous builds = 0
+Total wall time: 0:00:00

lib/kim/Install.py

@@ -10,7 +10,7 @@ import sys, os, subprocess, shutil
 from argparse import ArgumentParser
 
 sys.path.append('..')
-from install_helpers import fullpath, geturl
+from install_helpers import fullpath, geturl, checkmd5sum
 
 parser = ArgumentParser(prog='Install.py',
                         description="LAMMPS library build wrapper script")
@@ -18,7 +18,15 @@ parser = ArgumentParser(prog='Install.py',
 # settings
 
 thisdir = fullpath('.')
-version = "kim-api-2.1.2"
+version = "2.1.3"
+
+# known checksums for different KIM-API versions. used to validate the download.
+checksums = { \
+        '2.1.2' : '6ac52e14ef52967fc7858220b208cba5', \
+        '2.1.3' : '6ee829a1bbba5f8b9874c88c4c4ebff8', \
+        }
+
+
 
 # help message
 
@@ -50,7 +58,7 @@ pgroup.add_argument("-n", "--nobuild", action="store_true",
                     help="use the previously downloaded and compiled base KIM API.")
 pgroup.add_argument("-p", "--path",
                     help="specify location of existing KIM API installation.")
-parser.add_argument("-v", "--version", default=version,
+parser.add_argument("-v", "--version", default=version, choices=checksums.keys(),
                     help="set version of KIM API library to download and build (default: %s)" % version)
 parser.add_argument("-a", "--add",
                     help="add single KIM model or model driver. If adding 'everything', then all available OpenKIM models are added (may take a long time)")
@@ -73,6 +81,7 @@ if addflag and addmodelname == "everything":
   everythingflag = True
   buildflag = True
 verboseflag = args.verbose
+version = args.version
 
 if pathflag:
   buildflag = False
@@ -81,7 +90,7 @@ if pathflag:
     sys.exit("KIM API path %s does not exist" % kimdir)
   kimdir = fullpath(kimdir)
 
-url = "https://s3.openkim.org/kim-api/%s.txz" % version
+url = "https://s3.openkim.org/kim-api/kim-api-%s.txz" % version
 
 # set KIM API directory
 
@@ -115,21 +124,28 @@ if buildflag:
   # download entire kim-api tarball
 
   print("Downloading kim-api tarball ...")
-  geturl(url, "%s/%s.txz" % (thisdir, version))
+  filename = "kim-api-%s.txz" % version
+  geturl(url, "%s/%s" % (thisdir, filename))
+
+  # verify downloaded archive integrity via md5 checksum, if known.
+  if version in checksums:
+    if not checkmd5sum(checksums[version], filename):
+      sys.exit("Checksum for KIM-API library does not match")
+
   print("Unpacking kim-api tarball ...")
-  cmd = 'cd "%s"; rm -rf "%s"; tar -xJvf %s.txz' % (thisdir, version, version)
+  cmd = 'cd "%s"; rm -rf "kim-api-%s"; tar -xJvf %s' % (thisdir, version, filename)
   subprocess.check_output(cmd, stderr=subprocess.STDOUT, shell=True)
 
   # configure kim-api
 
   print("Configuring kim-api ...")
-  cmd = 'cd "%s/%s" && mkdir build && cd build && cmake .. -DCMAKE_INSTALL_PREFIX="%s" -DCMAKE_BUILD_TYPE=Release' % (thisdir,version,kimdir)
+  cmd = 'cd "%s/kim-api-%s" && mkdir build && cd build && cmake .. -DCMAKE_INSTALL_PREFIX="%s" -DCMAKE_BUILD_TYPE=Release' % (thisdir,version,kimdir)
   txt = subprocess.check_output(cmd,stderr=subprocess.STDOUT,shell=True)
   if verboseflag: print(txt.decode("UTF-8"))
 
   # build kim-api
   print("Building kim-api ...")
-  cmd = 'cd "%s/%s/build" && make -j2' % (thisdir, version)
+  cmd = 'cd "%s/kim-api-%s/build" && make -j2' % (thisdir, version)
   txt = subprocess.check_output(cmd, stderr=subprocess.STDOUT, shell=True)
   if verboseflag:
     print(txt.decode("UTF-8"))
@@ -137,7 +153,7 @@ if buildflag:
   # install kim-api
 
   print("Installing kim-api ...")
-  cmd = 'cd "%s/%s/build" && make -j2 install' % (thisdir, version)
+  cmd = 'cd "%s/kim-api-%s/build" && make -j2 install' % (thisdir, version)
   txt = subprocess.check_output(cmd, stderr=subprocess.STDOUT, shell=True)
   if verboseflag:
     print(txt.decode("UTF-8"))
@@ -145,7 +161,7 @@ if buildflag:
   # remove source files
 
   print("Removing kim-api source and build files ...")
-  cmd = 'cd "%s"; rm -rf %s; rm -rf %s.txz' % (thisdir, version, version)
+  cmd = 'cd "%s"; rm -rf kim-api-%s; rm -rf kim-api-%s.txz' % (thisdir, version, version)
   subprocess.check_output(cmd, stderr=subprocess.STDOUT, shell=True)
 
   # add all OpenKIM models, if desired

python/lammps.py

@@ -219,6 +219,12 @@ class lammps(object):
     self.c_imageint = get_ctypes_int(self.extract_setting("imageint"))
     self._installed_packages = None
 
+    # add way to insert Python callback for fix external
+    self.callback = {}
+    self.FIX_EXTERNAL_CALLBACK_FUNC = CFUNCTYPE(None, c_void_p, self.c_bigint, c_int, POINTER(self.c_tagint), POINTER(POINTER(c_double)), POINTER(POINTER(c_double)))
+    self.lib.lammps_set_fix_external_callback.argtypes = [c_void_p, c_char_p, self.FIX_EXTERNAL_CALLBACK_FUNC, c_void_p]
+    self.lib.lammps_set_fix_external_callback.restype = None
+
   # shut-down LAMMPS instance
 
   def __del__(self):
@@ -602,6 +608,42 @@ class lammps(object):
         self._installed_packages.append(sb.value.decode())
     return self._installed_packages
 
+  def set_fix_external_callback(self, fix_name, callback, caller=None):
+    import numpy as np
+    def _ctype_to_numpy_int(ctype_int):
+          if ctype_int == c_int32:
+            return np.int32
+          elif ctype_int == c_int64:
+            return np.int64
+          return np.intc
+
+    def callback_wrapper(caller_ptr, ntimestep, nlocal, tag_ptr, x_ptr, fext_ptr):
+      if cast(caller_ptr,POINTER(py_object)).contents:
+        pyCallerObj = cast(caller_ptr,POINTER(py_object)).contents.value
+      else:
+        pyCallerObj = None
+
+      tptr = cast(tag_ptr, POINTER(self.c_tagint * nlocal))
+      tag = np.frombuffer(tptr.contents, dtype=_ctype_to_numpy_int(self.c_tagint))
+      tag.shape = (nlocal)
+
+      xptr = cast(x_ptr[0], POINTER(c_double * nlocal * 3))
+      x = np.frombuffer(xptr.contents)
+      x.shape = (nlocal, 3)
+
+      fptr = cast(fext_ptr[0], POINTER(c_double * nlocal * 3))
+      f = np.frombuffer(fptr.contents)
+      f.shape = (nlocal, 3)
+
+      callback(pyCallerObj, ntimestep, nlocal, tag, x, f)
+
+    cFunc   = self.FIX_EXTERNAL_CALLBACK_FUNC(callback_wrapper)
+    cCaller = cast(pointer(py_object(caller)), c_void_p)
+
+    self.callback[fix_name] = { 'function': cFunc, 'caller': caller }
+
+    self.lib.lammps_set_fix_external_callback(self.lmp, fix_name.encode(), cFunc, cCaller)
+
 # -------------------------------------------------------------------------
 # -------------------------------------------------------------------------
 # -------------------------------------------------------------------------
@@ -872,8 +914,8 @@ class PyLammps(object):
     output = self.__getattr__('run')(*args, **kwargs)
 
     if(lammps.has_mpi4py):
-      output = self.lmp.comm.bcast(output, root=0) 
+      output = self.lmp.comm.bcast(output, root=0)
-    
+
     self.runs += get_thermo_data(output)
     return output
 

src/.gitignore

@@ -801,6 +801,8 @@
 /pair_comb3.h
 /pair_colloid.cpp
 /pair_colloid.h
+/pair_cosine_squared.cpp
+/pair_cosine_squared.h
 /pair_coul_diel.cpp
 /pair_coul_diel.h
 /pair_coul_long.cpp

src/GRANULAR/fix_wall_gran.cpp

@@ -1147,11 +1147,11 @@ void FixWallGran::granular(double rsq, double dx, double dy, double dz,
   else{
     knfac = E; //Hooke
     a = sqrt(dR);
+    Fne = knfac*delta;
     if (normal_model != HOOKE) {
       Fne *= a;
       knfac *= a;
     }
-    Fne = knfac*delta;
     if (normal_model == DMT)
       Fne -= 4*MY_PI*normal_coeffs[3]*Reff;
   }
@@ -1292,10 +1292,12 @@ void FixWallGran::granular(double rsq, double dx, double dy, double dz,
   // rolling resistance
   //****************************************
 
-  if (roll_model != ROLL_NONE) {
+  if (roll_model != ROLL_NONE || twist_model != NONE) {
     relrot1 = omega[0];
     relrot2 = omega[1];
     relrot3 = omega[2];
+  }
+  if (roll_model != ROLL_NONE){
 
     // rolling velocity, see eq. 31 of Wang et al, Particuology v 23, p 49 (2015)
     // This is different from the Marshall papers,

src/KIM/kim_query.cpp

@@ -64,6 +64,8 @@
 #include "input.h"
 #include "modify.h"
 #include "variable.h"
+#include "version.h"
+#include "info.h"
 #include "fix_store_kim.h"
 
 #if defined(LMP_KIM_CURL)
@@ -242,6 +244,10 @@ char *do_query(char *qfunction, char * model_name, int narg, char **arg,
       curl_easy_setopt(handle, CURLOPT_SSL_VERIFYPEER, 0L);
       curl_easy_setopt(handle, CURLOPT_SSL_VERIFYHOST, 0L);
 #endif
+      std::string user_agent = std::string("kim_query--LAMMPS/")
+                               + LAMMPS_VERSION
+                               + " (" + Info::get_os_info() + ")";
+      curl_easy_setopt(handle, CURLOPT_USERAGENT, user_agent.c_str());
 
       curl_easy_setopt(handle, CURLOPT_URL, url.c_str());
       curl_easy_setopt(handle, CURLOPT_FOLLOWLOCATION, 1L);

src/MAKE/MACHINES/Makefile.astra_arm

@@ -0,0 +1,119 @@
+# astra_arm - ThunderX2 ARM/OpenMPI Based, use Arm-PL for FFTW
+
+# need to load the following modules:
+# 1) arm-based developer pack
+
+SHELL = /bin/sh
+
+# ---------------------------------------------------------------------
+# compiler/linker settings
+# specify flags and libraries needed for your compiler
+
+CC =            mpicxx
+CCFLAGS =       -O3 -funroll-loops -fopenmp -mcpu=thunderx2t99 -mtune=thunderx2t99
+DEPFLAGS =      -M
+LINK =          mpicxx
+LINKFLAGS =     -O3 -fopenmp -mcpu=thunderx2t99 -mtune=thunderx2t99
+LIB =           -lstdc++
+ARCHIVE =       ar
+ARFLAGS =       -rcsv
+SIZE =          size
+
+# ---------------------------------------------------------------------
+# LAMMPS-specific settings, all OPTIONAL
+# specify settings for LAMMPS features you will use
+# if you change any -D setting, do full re-compile after "make clean"
+
+# LAMMPS ifdef settings
+# see possible settings in Section 2.2 (step 4) of manual
+
+LMP_INC =       -DLAMMPS_GZIP
+
+# MPI library
+# see discussion in Section 2.2 (step 5) of manual
+# MPI wrapper compiler/linker can provide this info
+# can point to dummy MPI library in src/STUBS as in Makefile.serial
+# use -D MPICH and OMPI settings in INC to avoid C++ lib conflicts
+# INC = path for mpi.h, MPI compiler settings
+# PATH = path for MPI library
+# LIB = name of MPI library
+
+MPI_INC =       
+MPI_PATH = 
+MPI_LIB =
+
+# FFT library
+# see discussion in Section 2.2 (step 6) of manaul
+# can be left blank to use provided KISS FFT library
+# INC = -DFFT setting, e.g. -DFFT_FFTW, FFT compiler settings
+# PATH = path for FFT library
+# LIB = name of FFT library
+
+FFT_INC =    -DFFT_FFTW2 -armpl
+FFT_PATH =   
+FFT_LIB =     
+
+# JPEG and/or PNG library
+# see discussion in Section 2.2 (step 7) of manual
+# only needed if -DLAMMPS_JPEG or -DLAMMPS_PNG listed with LMP_INC
+# INC = path(s) for jpeglib.h and/or png.h
+# PATH = path(s) for JPEG library and/or PNG library
+# LIB = name(s) of JPEG library and/or PNG library
+
+JPG_INC =       
+JPG_PATH = 
+JPG_LIB =
+
+# ---------------------------------------------------------------------
+# build rules and dependencies
+# do not edit this section
+
+include	Makefile.package.settings
+include	Makefile.package
+
+EXTRA_INC = $(LMP_INC) $(PKG_INC) $(MPI_INC) $(FFT_INC) $(JPG_INC) $(PKG_SYSINC)
+EXTRA_PATH = $(PKG_PATH) $(MPI_PATH) $(FFT_PATH) $(JPG_PATH) $(PKG_SYSPATH)
+EXTRA_LIB = $(PKG_LIB) $(MPI_LIB) $(FFT_LIB) $(JPG_LIB) $(PKG_SYSLIB)
+EXTRA_CPP_DEPENDS = $(PKG_CPP_DEPENDS)
+EXTRA_LINK_DEPENDS = $(PKG_LINK_DEPENDS)
+
+# Path to src files
+
+vpath %.cpp ..
+vpath %.h ..
+
+# Link target
+
+$(EXE):	$(OBJ) $(EXTRA_LINK_DEPENDS)
+	$(LINK) $(LINKFLAGS) $(EXTRA_PATH) $(OBJ) $(EXTRA_LIB) $(LIB) -o $(EXE)
+	$(SIZE) $(EXE)
+
+# Library targets
+
+lib:	$(OBJ) $(EXTRA_LINK_DEPENDS)
+	$(ARCHIVE) $(ARFLAGS) $(EXE) $(OBJ)
+
+shlib:	$(OBJ) $(EXTRA_LINK_DEPENDS)
+	$(CC) $(CCFLAGS) $(SHFLAGS) $(SHLIBFLAGS) $(EXTRA_PATH) -o $(EXE) \
+        $(OBJ) $(EXTRA_LIB) $(LIB)
+
+# Compilation rules
+
+%.o:%.cpp $(EXTRA_CPP_DEPENDS)
+	$(CC) $(CCFLAGS) $(SHFLAGS) $(EXTRA_INC) -c $<
+
+%.d:%.cpp $(EXTRA_CPP_DEPENDS)
+	$(CC) $(CCFLAGS) $(EXTRA_INC) $(DEPFLAGS) $< > $@
+
+%.o:%.cu $(EXTRA_CPP_DEPENDS)
+	$(CC) $(CCFLAGS) $(SHFLAGS) $(EXTRA_INC) -c $<
+
+# Individual dependencies
+
+depend : fastdep.exe $(SRC)
+	@./fastdep.exe $(EXTRA_INC) -- $^ > .depend || exit 1
+
+fastdep.exe: ../DEPEND/fastdep.c
+	cc -O -o $@ $<
+
+sinclude .depend

src/MAKE/MACHINES/Makefile.astra_gcc

@@ -0,0 +1,119 @@
+# astra-gcc - ThunderX2 GCC/OpenMPI Based, FFTW
+
+# need to load the following modules:
+# 1) arm-based developer pack
+
+SHELL = /bin/sh
+
+# ---------------------------------------------------------------------
+# compiler/linker settings
+# specify flags and libraries needed for your compiler
+
+CC =            mpicxx
+CCFLAGS =       -O3 -funroll-loops -fopenmp -mcpu=thunderx2t99 -mtune=thunderx2t99
+DEPFLAGS =      -M
+LINK =          mpicxx
+LINKFLAGS =     -O3 -fopenmp -mcpu=thunderx2t99 -mtune=thunderx2t99
+LIB =           -lstdc++
+ARCHIVE =       ar
+ARFLAGS =       -rcsv
+SIZE =          size
+
+# ---------------------------------------------------------------------
+# LAMMPS-specific settings, all OPTIONAL
+# specify settings for LAMMPS features you will use
+# if you change any -D setting, do full re-compile after "make clean"
+
+# LAMMPS ifdef settings
+# see possible settings in Section 2.2 (step 4) of manual
+
+LMP_INC =       -DLAMMPS_GZIP
+
+# MPI library
+# see discussion in Section 2.2 (step 5) of manual
+# MPI wrapper compiler/linker can provide this info
+# can point to dummy MPI library in src/STUBS as in Makefile.serial
+# use -D MPICH and OMPI settings in INC to avoid C++ lib conflicts
+# INC = path for mpi.h, MPI compiler settings
+# PATH = path for MPI library
+# LIB = name of MPI library
+
+MPI_INC =       
+MPI_PATH = 
+MPI_LIB =
+
+# FFT library
+# see discussion in Section 2.2 (step 6) of manaul
+# can be left blank to use provided KISS FFT library
+# INC = -DFFT setting, e.g. -DFFT_FFTW, FFT compiler settings
+# PATH = path for FFT library
+# LIB = name of FFT library
+
+FFT_INC =    -DFFT_FFTW2 -I${FFTW_INC}
+FFT_PATH =   -I${FFTW_INC}
+FFT_LIB =    -L${FFTW_LIB} -lfftw3
+
+# JPEG and/or PNG library
+# see discussion in Section 2.2 (step 7) of manual
+# only needed if -DLAMMPS_JPEG or -DLAMMPS_PNG listed with LMP_INC
+# INC = path(s) for jpeglib.h and/or png.h
+# PATH = path(s) for JPEG library and/or PNG library
+# LIB = name(s) of JPEG library and/or PNG library
+
+JPG_INC =       
+JPG_PATH = 
+JPG_LIB =
+
+# ---------------------------------------------------------------------
+# build rules and dependencies
+# do not edit this section
+
+include	Makefile.package.settings
+include	Makefile.package
+
+EXTRA_INC = $(LMP_INC) $(PKG_INC) $(MPI_INC) $(FFT_INC) $(JPG_INC) $(PKG_SYSINC)
+EXTRA_PATH = $(PKG_PATH) $(MPI_PATH) $(FFT_PATH) $(JPG_PATH) $(PKG_SYSPATH)
+EXTRA_LIB = $(PKG_LIB) $(MPI_LIB) $(FFT_LIB) $(JPG_LIB) $(PKG_SYSLIB)
+EXTRA_CPP_DEPENDS = $(PKG_CPP_DEPENDS)
+EXTRA_LINK_DEPENDS = $(PKG_LINK_DEPENDS)
+
+# Path to src files
+
+vpath %.cpp ..
+vpath %.h ..
+
+# Link target
+
+$(EXE):	$(OBJ) $(EXTRA_LINK_DEPENDS)
+	$(LINK) $(LINKFLAGS) $(EXTRA_PATH) $(OBJ) $(EXTRA_LIB) $(LIB) -o $(EXE)
+	$(SIZE) $(EXE)
+
+# Library targets
+
+lib:	$(OBJ) $(EXTRA_LINK_DEPENDS)
+	$(ARCHIVE) $(ARFLAGS) $(EXE) $(OBJ)
+
+shlib:	$(OBJ) $(EXTRA_LINK_DEPENDS)
+	$(CC) $(CCFLAGS) $(SHFLAGS) $(SHLIBFLAGS) $(EXTRA_PATH) -o $(EXE) \
+        $(OBJ) $(EXTRA_LIB) $(LIB)
+
+# Compilation rules
+
+%.o:%.cpp $(EXTRA_CPP_DEPENDS)
+	$(CC) $(CCFLAGS) $(SHFLAGS) $(EXTRA_INC) -c $<
+
+%.d:%.cpp $(EXTRA_CPP_DEPENDS)
+	$(CC) $(CCFLAGS) $(EXTRA_INC) $(DEPFLAGS) $< > $@
+
+%.o:%.cu $(EXTRA_CPP_DEPENDS)
+	$(CC) $(CCFLAGS) $(SHFLAGS) $(EXTRA_INC) -c $<
+
+# Individual dependencies
+
+depend : fastdep.exe $(SRC)
+	@./fastdep.exe $(EXTRA_INC) -- $^ > .depend || exit 1
+
+fastdep.exe: ../DEPEND/fastdep.c
+	cc -O -o $@ $<
+
+sinclude .depend

src/MAKE/MACHINES/Makefile.theta

@@ -1,7 +1,5 @@
-# knl = Flags for Knights Landing Xeon Phi Processor,Intel Compiler/MPI,MKL FFT
+# theta = Flags for Knights Landing Xeon Phi Processor, Intel Compiler, Cray MPI, MKL FFT
-# module load perftools-base perftools
 # make theta -j 8
-# pat_build -g mpi -u ./lmp_theta
 
 SHELL = /bin/sh
 
@@ -10,24 +8,16 @@ SHELL = /bin/sh
 # specify flags and libraries needed for your compiler
 
 CC =        CC -mkl
-#OPTFLAGS =       -O0
+OPTFLAGS =  -xMIC-AVX512 -O3 -fp-model fast=2 -no-prec-div -qoverride-limits
-OPTFLAGS =      -xMIC-AVX512 -O3 -fp-model fast=2 -no-prec-div -qoverride-limits
 CCFLAGS =   -g -qopenmp -DLAMMPS_MEMALIGN=64 -qno-offload \
                 -fno-alias -ansi-alias -restrict $(OPTFLAGS)
-#CCFLAGS +=      -DLMP_INTEL_NO_TBB
+CCFLAGS +=  -std=c++11
-#CCFLAGS +=      -DLAMMPS_BIGBIG
-#CCFLAGS +=      -D_USE_PAPI
-#CCFLAGS +=      -D_USE_CRAYPAT_API
 SHFLAGS =   -fPIC
 DEPFLAGS =  -M
 
 LINK =      $(CC)
-LINKFLAGS = -g -qopenmp $(OPTFLAGS)
+LINKFLAGS = -g -qopenmp $(OPTFLAGS) -dynamic
-LINKFLAGS += -dynamic
+LIB =       -ltbbmalloc
-LIB =
-#LIB +=           -L${TBBROOT}/lib/intel64/gcc4.7 -ltbbmalloc
-LIB +=           -ltbbmalloc
-#LIB +=          /soft/debuggers/forge-7.0-2017-02-16/lib/64/libdmallocthcxx.a -zmuldefs
 SIZE =      size
 
 ARCHIVE =   ar
@@ -99,35 +89,35 @@ vpath %.h ..
 # Link target
 
 $(EXE): $(OBJ)
-    $(LINK) $(LINKFLAGS) $(EXTRA_PATH) $(OBJ) $(EXTRA_LIB) $(LIB) -o $(EXE)
+	$(LINK) $(LINKFLAGS) $(EXTRA_PATH) $(OBJ) $(EXTRA_LIB) $(LIB) -o $(EXE)
-    $(SIZE) $(EXE)
+	$(SIZE) $(EXE)
 
 # Library targets
 
 lib:    $(OBJ)
-    $(ARCHIVE) $(ARFLAGS) $(EXE) $(OBJ)
+	$(ARCHIVE) $(ARFLAGS) $(EXE) $(OBJ)
 
 shlib:  $(OBJ)
-    $(CC) $(CCFLAGS) $(SHFLAGS) $(SHLIBFLAGS) $(EXTRA_PATH) -o $(EXE) \
+	$(CC) $(CCFLAGS) $(SHFLAGS) $(SHLIBFLAGS) $(EXTRA_PATH) -o $(EXE) \
-        $(OBJ) $(EXTRA_LIB) $(LIB)
+	$(OBJ) $(EXTRA_LIB) $(LIB)
 
 # Compilation rules
 
 %.o:%.cpp
-    $(CC) $(CCFLAGS) $(SHFLAGS) $(EXTRA_INC) -c $<
+	$(CC) $(CCFLAGS) $(SHFLAGS) $(EXTRA_INC) -c $<
 
 %.d:%.cpp
-    $(CC) $(CCFLAGS) $(EXTRA_INC) $(DEPFLAGS) $< > $@
+	$(CC) $(CCFLAGS) $(EXTRA_INC) $(DEPFLAGS) $< > $@
 
 %.o:%.cu
-    $(CC) $(CCFLAGS) $(SHFLAGS) $(EXTRA_INC) -c $<
+	$(CC) $(CCFLAGS) $(SHFLAGS) $(EXTRA_INC) -c $<
 
 # Individual dependencies
 
 depend : fastdep.exe $(SRC)
-    @./fastdep.exe $(EXTRA_INC) -- $^ > .depend || exit 1
+	@./fastdep.exe $(EXTRA_INC) -- $^ > .depend || exit 1
 
 fastdep.exe: ../DEPEND/fastdep.c
-    icc -O -o $@ $<
+	icc -O -o $@ $<
 
 sinclude .depend

src/MAKE/OPTIONS/Makefile.intel_cpu

@@ -10,7 +10,8 @@ CC =		mpiicpc
 OPTFLAGS =      -xHost -O2 -fp-model fast=2 -no-prec-div -qoverride-limits \
                 -qopt-zmm-usage=high
 CCFLAGS =	-qopenmp -qno-offload -fno-alias -ansi-alias -restrict \
-                -DLMP_INTEL_USELRT -DLMP_USE_MKL_RNG $(OPTFLAGS)
+                -DLMP_INTEL_USELRT -DLMP_USE_MKL_RNG $(OPTFLAGS) \
+                -I$(MKLROOT)/include
 SHFLAGS =	-fPIC
 DEPFLAGS =	-M
 

src/MAKE/OPTIONS/Makefile.intel_cpu_intelmpi

@@ -10,7 +10,8 @@ CC =		mpiicpc
 OPTFLAGS =      -xHost -O2 -fp-model fast=2 -no-prec-div -qoverride-limits \
                 -qopt-zmm-usage=high
 CCFLAGS =	-qopenmp -qno-offload -fno-alias -ansi-alias -restrict \
-                -DLMP_INTEL_USELRT -DLMP_USE_MKL_RNG $(OPTFLAGS)
+                -DLMP_INTEL_USELRT -DLMP_USE_MKL_RNG $(OPTFLAGS) \
+                -I$(MKLROOT)/include
 SHFLAGS =	-fPIC
 DEPFLAGS =	-M
 

src/MAKE/OPTIONS/Makefile.intel_cpu_mpich

@@ -10,7 +10,8 @@ CC =		mpicxx -cxx=icc
 OPTFLAGS =      -xHost -O2 -fp-model fast=2 -no-prec-div -qoverride-limits \
                 -qopt-zmm-usage=high
 CCFLAGS =	-qopenmp -qno-offload -fno-alias -ansi-alias -restrict \
-                -DLMP_INTEL_USELRT -DLMP_USE_MKL_RNG $(OPTFLAGS)
+                -DLMP_INTEL_USELRT -DLMP_USE_MKL_RNG $(OPTFLAGS) \
+                -I$(MKLROOT)/include
 SHFLAGS =	-fPIC
 DEPFLAGS =	-M
 

src/MAKE/OPTIONS/Makefile.intel_cpu_openmpi

@@ -11,7 +11,8 @@ CC =		mpicxx
 OPTFLAGS =      -xHost -O2 -fp-model fast=2 -no-prec-div -qoverride-limits \
                 -qopt-zmm-usage=high
 CCFLAGS =	-qopenmp -qno-offload -fno-alias -ansi-alias -restrict \
-                -DLMP_INTEL_USELRT -DLMP_USE_MKL_RNG $(OPTFLAGS)
+                -DLMP_INTEL_USELRT -DLMP_USE_MKL_RNG $(OPTFLAGS) \
+                -I$(MKLROOT)/include
 SHFLAGS =	-fPIC
 DEPFLAGS =	-M
 

src/MAKE/OPTIONS/Makefile.knl

@@ -9,7 +9,8 @@ SHELL = /bin/sh
 CC =		mpiicpc
 OPTFLAGS =      -xMIC-AVX512 -O2 -fp-model fast=2 -no-prec-div -qoverride-limits
 CCFLAGS =	-qopenmp -qno-offload -fno-alias -ansi-alias -restrict \
-                -DLMP_INTEL_USELRT -DLMP_USE_MKL_RNG $(OPTFLAGS)
+                -DLMP_INTEL_USELRT -DLMP_USE_MKL_RNG $(OPTFLAGS) \
+                -I$(MKLROOT)/include
 SHFLAGS =	-fPIC
 DEPFLAGS =	-M
 

src/MANYBODY/pair_tersoff.cpp

@@ -30,9 +30,11 @@
 #include "error.h"
 
 #include "math_const.h"
+#include "math_special.h"
 
 using namespace LAMMPS_NS;
 using namespace MathConst;
+using namespace MathSpecial;
 
 #define MAXLINE 1024
 #define DELTA 4
@@ -604,7 +606,7 @@ double PairTersoff::zeta(Param *param, double rsqij, double rsqik,
   costheta = (delrij[0]*delrik[0] + delrij[1]*delrik[1] +
               delrij[2]*delrik[2]) / (rij*rik);
 
-  if (param->powermint == 3) arg = pow(param->lam3 * (rij-rik),3.0);
+  if (param->powermint == 3) arg = cube(param->lam3 * (rij-rik));
   else arg = param->lam3 * (rij-rik);
 
   if (arg > 69.0776) ex_delr = 1.e30;
@@ -744,7 +746,7 @@ void PairTersoff::ters_zetaterm_d(double prefactor,
 
   fc = ters_fc(rik,param);
   dfc = ters_fc_d(rik,param);
-  if (param->powermint == 3) tmp = pow(param->lam3 * (rij-rik),3.0);
+  if (param->powermint == 3) tmp = cube(param->lam3 * (rij-rik));
   else tmp = param->lam3 * (rij-rik);
 
   if (tmp > 69.0776) ex_delr = 1.e30;
@@ -752,7 +754,7 @@ void PairTersoff::ters_zetaterm_d(double prefactor,
   else ex_delr = exp(tmp);
 
   if (param->powermint == 3)
-    ex_delr_d = 3.0*pow(param->lam3,3.0) * pow(rij-rik,2.0)*ex_delr;
+    ex_delr_d = 3.0*cube(param->lam3) * square(rij-rik)*ex_delr;
   else ex_delr_d = param->lam3 * ex_delr;
 
   cos_theta = vec3_dot(rij_hat,rik_hat);

src/MANYBODY/pair_tersoff_mod.cpp

@@ -25,11 +25,13 @@
 #include "force.h"
 #include "comm.h"
 #include "math_const.h"
+#include "math_special.h"
 #include "memory.h"
 #include "error.h"
 
 using namespace LAMMPS_NS;
 using namespace MathConst;
+using namespace MathSpecial;
 
 #define MAXLINE 1024
 #define DELTA 4
@@ -241,7 +243,7 @@ double PairTersoffMOD::zeta(Param *param, double rsqij, double rsqik,
   costheta = (delrij[0]*delrik[0] + delrij[1]*delrik[1] +
               delrij[2]*delrik[2]) / (rij*rik);
 
-  if (param->powermint == 3) arg = pow(param->lam3 * (rij-rik),3.0);
+  if (param->powermint == 3) arg = cube(param->lam3 * (rij-rik));
   else arg = param->lam3 * (rij-rik);
 
   if (arg > 69.0776) ex_delr = 1.e30;
@@ -314,7 +316,7 @@ void PairTersoffMOD::ters_zetaterm_d(double prefactor,
 
   fc = ters_fc(rik,param);
   dfc = ters_fc_d(rik,param);
-  if (param->powermint == 3) tmp = pow(param->lam3 * (rij-rik),3.0);
+  if (param->powermint == 3) tmp = cube(param->lam3 * (rij-rik));
   else tmp = param->lam3 * (rij-rik);
 
   if (tmp > 69.0776) ex_delr = 1.e30;
@@ -322,7 +324,7 @@ void PairTersoffMOD::ters_zetaterm_d(double prefactor,
   else ex_delr = exp(tmp);
 
   if (param->powermint == 3)
-    ex_delr_d = 3.0*pow(param->lam3,3.0) * pow(rij-rik,2.0)*ex_delr;
+    ex_delr_d = 3.0*cube(param->lam3) * square(rij-rik)*ex_delr;
   else ex_delr_d = param->lam3 * ex_delr;
 
   cos_theta = vec3_dot(rij_hat,rik_hat);

src/MANYBODY/pair_tersoff_zbl.cpp

@@ -28,9 +28,11 @@
 #include "memory.h"
 #include "error.h"
 #include "math_const.h"
+#include "math_special.h"
 
 using namespace LAMMPS_NS;
 using namespace MathConst;
+using namespace MathSpecial;
 
 #define MAXLINE 1024
 #define DELTA 4
@@ -226,7 +228,7 @@ void PairTersoffZBL::repulsive(Param *param, double rsq, double &fforce,
 
   // ZBL repulsive portion
 
-  double esq = pow(global_e,2.0);
+  double esq = square(global_e);
   double a_ij = (0.8854*global_a_0) /
     (pow(param->Z_i,0.23) + pow(param->Z_j,0.23));
   double premult = (param->Z_i * param->Z_j * esq)/(4.0*MY_PI*global_epsilon_0);
@@ -286,5 +288,5 @@ double PairTersoffZBL::F_fermi(double r, Param *param)
 double PairTersoffZBL::F_fermi_d(double r, Param *param)
 {
   return param->ZBLexpscale*exp(-param->ZBLexpscale*(r-param->ZBLcut)) /
-    pow(1.0 + exp(-param->ZBLexpscale*(r-param->ZBLcut)),2.0);
+    square(1.0 + exp(-param->ZBLexpscale*(r-param->ZBLcut)));
 }

src/USER-INTEL/TEST/run_benchmarks.sh

@@ -40,7 +40,7 @@ export LMP_THREAD_LIST="2"     # -- For 2 threads per core w/ HT enabled
 # End settings for your system
 #########################################################################
 
-export WORKLOADS="lj rhodo rhodo_lrt lc sw water eam"
+export WORKLOADS="lj rhodo lc sw water eam airebo dpd tersoff"
 export LMP_ARGS="-pk intel 0 -sf intel -screen none -v d 1"
 export RLMP_ARGS="-pk intel 0 lrt yes -sf intel -screen none -v d 1"
 
@@ -51,7 +51,7 @@ export LOG_DIR=$LOG_DIR_HOST"_"$LOG_DIR_HEADER"_"$DATE_STRING
 mkdir $LOG_DIR
 
 export I_MPI_PIN_DOMAIN=core
-export I_MPI_FABRICS=shm
+#export I_MPI_FABRICS=shm
 export KMP_BLOCKTIME=0
 
 echo -n "Creating restart file...."
@@ -59,13 +59,12 @@ $MPI -np $LMP_CORES $LMP_BIN -in in.lc_generate_restart -log none $LMP_ARGS
 echo "Done."
 for threads in $LMP_THREAD_LIST
 do
-  export OMP_NUM_THREADS=$threads
   for workload in $WORKLOADS
   do
     export LOGFILE=$LOG_DIR/$workload.$LMP_CORES"c"$threads"t".log
     echo "Running $LOGFILE"
     cmd="$MPI -np $LMP_CORES $LMP_BIN -in in.intel.$workload -log $LOGFILE $LMP_ARGS";
-    rthreads=$threads
+    export OMP_NUM_THREADS=$threads
     unset KMP_AFFINITY
     $cmd
 

src/USER-MISC/README

@@ -31,6 +31,7 @@ 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 gyration/shape, Evangelos Voyiatzis, evoyiatzis at gmail.com, 25 July 2019
+compute hma, Andrew Schultz & David Kofke (UB), ajs42 at buffalo.edu & kofke at buffalo.edu, 1 Jul 2019
 compute pressure/cylinder, Cody K. Addington (NCSU), , 2 Oct 2018
 compute momentum, Rupert Nash (University of Edinburgh), r.nash at epcc.ed.ac.uk, 28 June 2019
 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
@@ -69,6 +70,7 @@ improper_style ring, Georgios Vogiatzis, gvog at chemeng.ntua.gr, 25 May 12
 improper_style distance, Paolo Raiteri, p.raiteri at curtin.edu.au, 2 Dec 15
 pair_style agni, Axel Kohlmeyer, akohlmey at gmail.com, 9 Nov 16
 pair_style buck/mdf, Paolo Raiteri, p.raiteri at curtin.edu.au, 2 Dec 15
+pair_style cosine/squared, Eugen Rozic, eugen.rozic.17 at ucl.ac.uk, 9 Aug 19
 pair_style coul/diel, Axel Kohlmeyer, akohlmey at gmail.com, 1 Dec 11
 pair_style coul/shield, Wengen Ouyang (Tel Aviv University), w.g.ouyang at gmail dot com, 30 Mar 18
 pair_style dipole/sf, Mario Orsi, orsimario at gmail.com, 8 Aug 11

src/USER-MISC/compute_hma.cpp

@@ -0,0 +1,516 @@
+/* ----------------------------------------------------------------------
+   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.
+------------------------------------------------------------------------- */
+
+/* ----------------------------------------------------------------------
+This compute implements harmonically-mapped averaging for crystalline solids.
+The current implementation handles atomic crystals.
+
+Computing the heat capacity relies on being able to compute the second
+derivative of the energy with respect to atom positions.  This capability is
+provided by the single2 method in Pair, but is currently only implemented for
+the shifted-force LJ potential (lj/smooth/linear).  Pair::single2 takes a single
+pair and (like Pair::single) returns the energy and sets the force as an out
+parameter, but also sets the elements of 6-element double array out parameter,
+which are the unique components of the atomic Hessian tensor for the pair.  A
+helper method exists (Pair::pairTensor), which will compute the tensor from
+linear derivatives and the vector between the positions.  HMA Heat capacity can
+be computed for other models by implementing single2 in those Pair classes.
+
+More information about HMA is available in these publications:
+
+A. J. Schultz, D. A. Kofke, “Comprehensive high-precision high-accuracy
+equation of state and coexistence properties for classical Lennard-Jones
+crystals and low-temperature fluid phases”, J. Chem. Phys. 149, 204508 (2018)
+https://dx.doi.org/10.1063/1.5053714
+
+S. G. Moustafa, A. J. Schultz, D. A. Kofke, “Harmonically Assisted Methods for
+Computing the Free Energy of Classical Crystals by Molecular Simulation: A
+Comparative Study”, J. Chem. Theory Comput. 13, 825-834 (2017)
+https://dx.doi.org/10.1021/acs.jctc.6b01082
+
+S. G. Moustafa, A. J. Schultz, D. A. Kofke, “Very fast averaging of thermal
+properties of crystals by molecular simulation”, Phys. Rev. E 92, 043303 (2015)
+https://dx.doi.org/10.1103/PhysRevE.92.043303
+------------------------------------------------------------------------- */
+
+#include <cmath>
+#include <cstring>
+#include <mpi.h>
+#include "compute_hma.h"
+#include "atom.h"
+#include "update.h"
+#include "force.h"
+#include "pair.h"
+#include "bond.h"
+#include "angle.h"
+#include "dihedral.h"
+#include "improper.h"
+#include "kspace.h"
+#include "group.h"
+#include "domain.h"
+#include "modify.h"
+#include "fix.h"
+#include "fix_store.h"
+#include "memory.h"
+#include "error.h"
+#include "comm.h"
+#include "neighbor.h"
+#include "neigh_request.h"
+#include "neigh_list.h"
+
+#include <vector>
+
+
+using namespace LAMMPS_NS;
+
+/* ---------------------------------------------------------------------- */
+
+ComputeHMA::ComputeHMA(LAMMPS *lmp, int narg, char **arg) :
+  Compute(lmp, narg, arg), id_temp(NULL), deltaR(NULL)
+{
+  if (narg < 4) error->all(FLERR,"Illegal compute hma command");                    
+  if (igroup) error->all(FLERR,"Compute hma must use group all");
+  if (strcmp(arg[3],"NULL") == 0) {error->all(FLERR,"fix ID specifying the set temperature of canonical simulation is required");}       
+  else {
+    int n = strlen(arg[3]) + 1;                  
+    id_temp = new char[n];   
+    strcpy(id_temp,arg[3]);                         
+  }
+ 
+  create_attribute = 1;                  
+  extscalar = 1;                         
+  timeflag = 1;                         
+
+  // (from compute displace/atom) create a new fix STORE style  
+  // our new fix's id (id_fix)= compute-ID + COMPUTE_STORE 
+  // our new fix's group = same as compute group
+
+  int n = strlen(id) + strlen("_COMPUTE_STORE") + 1;
+  id_fix = new char[n];
+  strcpy(id_fix,id);
+  strcat(id_fix,"_COMPUTE_STORE");
+
+  char **newarg = new char*[6];     
+  newarg[0] = id_fix;
+  newarg[1] = group->names[igroup];
+  newarg[2] = (char *) "STORE";              
+  newarg[3] = (char *) "peratom";
+  newarg[4] = (char *) "1";
+  newarg[5] = (char *) "3";
+  modify->add_fix(6,newarg);        
+  fix = (FixStore *) modify->fix[modify->nfix-1];    
+ 
+  delete [] newarg;                               
+
+  // calculate xu,yu,zu for fix store array
+  // skip if reset from restart file
+
+  if (fix->restart_reset) fix->restart_reset = 0;   
+  else {
+    double **xoriginal = fix->astore;            
+    double **x = atom->x;
+    imageint *image = atom->image;
+    int nlocal = atom->nlocal;
+
+    for (int i = 0; i < nlocal; i++)
+      domain->unmap(x[i],image[i],xoriginal[i]);          
+  }
+
+  vector_flag = 1;
+  extvector = -1;
+  comm_forward = 0; // 3 if 2nd derivative needed
+
+  computeU = computeP = computeCv = -1;
+  returnAnharmonic = 0;
+  size_vector = 0;
+  memory->create(extlist, 3, "hma:extlist");
+  for (int iarg=4; iarg<narg; iarg++) {
+    if (!strcmp(arg[iarg], "u")) {
+      if (computeU>-1) continue;
+      computeU = size_vector;
+      extlist[size_vector] = 1;
+      size_vector++;
+    }
+    else if (!strcmp(arg[iarg], "p")) {
+      if (iarg+2 > narg) error->all(FLERR,"Illegal compute hma command");
+      if (computeP>-1) continue;
+      computeP = size_vector;
+      deltaPcap = force->numeric(FLERR, arg[iarg+1]);
+      extlist[size_vector] = 0;
+      size_vector++;
+      iarg++;
+    }
+    else if (!strcmp(arg[iarg], "cv")) {
+      if (computeCv>-1) continue;
+      computeCv = size_vector;
+      comm_forward = 3;
+      extlist[size_vector] = 1;
+      size_vector++;
+    }
+    else if (!strcmp(arg[iarg], "anharmonic")) {
+      // the first time we're called, we'll grab lattice pressure and energy
+      returnAnharmonic = -1;
+    }
+    else {
+      error->all(FLERR,"Illegal compute hma command");
+    }
+  }
+
+  if (size_vector == 0) {
+    error->all(FLERR,"Illegal compute hma command");
+  }
+  if (size_vector<3) {
+    memory->grow(extlist, size_vector, "hma:extlist");
+  }
+  memory->create(vector, size_vector, "hma:vector");
+
+  if (computeU>-1 || computeCv>-1) {
+    peflag = 1;                             
+  }
+  if (computeP>-1) {
+    pressflag = 1;
+  }
+
+  nmax = 0;
+}
+
+/* ---------------------------------------------------------------------- */
+
+ComputeHMA::~ComputeHMA()
+{
+  // check nfix in case all fixes have already been deleted
+  if (modify->nfix) modify->delete_fix(id_fix);
+
+  delete [] id_fix;
+  delete [] id_temp;
+  memory->destroy(extlist);
+  memory->destroy(vector);
+  memory->destroy(deltaR);
+}
+
+/* ---------------------------------------------------------------------- */
+
+void ComputeHMA::init() {
+  if (computeCv>-1) {
+    if (force->pair == NULL)
+      error->all(FLERR,"No pair style is defined for compute hma cv");
+    if (force->pair->single_enable == 0)
+      error->all(FLERR,"Pair style does not support compute hma cv");
+  }
+
+  int irequest = neighbor->request(this,instance_me);
+  neighbor->requests[irequest]->pair = 0; 
+  neighbor->requests[irequest]->compute = 1; 
+  neighbor->requests[irequest]->occasional = 1; 
+}
+
+void ComputeHMA::init_list(int id, NeighList *ptr)
+{
+  list = ptr;
+}
+
+void ComputeHMA::setup()
+{
+  int dummy=0;
+  int ifix = modify->find_fix(id_temp);
+  if (ifix < 0) error->all(FLERR,"Could not find compute hma temperature ID");
+  double * temperat = (double *) modify->fix[ifix]->extract("t_target",dummy);      
+  if (temperat==NULL) error->all(FLERR,"Could not find compute hma temperature ID");
+  finaltemp = * temperat;       
+
+  // set fix which stores original atom coords
+
+  int ifix2 = modify->find_fix(id_fix);
+  if (ifix2 < 0) error->all(FLERR,"Could not find hma store fix ID");
+  fix = (FixStore *) modify->fix[ifix2];             
+}
+
+/* ---------------------------------------------------------------------- */
+
+void ComputeHMA::compute_vector()
+{
+  invoked_vector = update->ntimestep;         
+
+  // grow deltaR array if necessary
+  if (comm_forward>0 && atom->nmax > nmax) {
+    memory->destroy(deltaR);
+    nmax = atom->nmax;
+    memory->create(deltaR,nmax,3,"hma:deltaR");
+  }
+
+  double **xoriginal = fix->astore;
+  double fdr = 0.0;
+  double **x = atom->x;
+  double **f = atom->f;
+
+  imageint *image = atom->image;
+  int nlocal = atom->nlocal;
+
+  double *h = domain->h;
+  double xprd = domain->xprd;   
+  double yprd = domain->yprd;
+  double zprd = domain->zprd;
+
+  double u = 0.0;
+  if (computeU>-1 || computeCv>-1) {
+    if (force->pair) u += force->pair->eng_vdwl + force->pair->eng_coul;
+    if (force->bond) u += force->bond->energy;
+    if (force->angle) u += force->angle->energy;
+    if (force->dihedral) u += force->dihedral->energy;
+    if (force->improper) u += force->improper->energy;
+  }
+
+  int dimension = domain->dimension;
+  double p = 0, vol = 0;
+  if (computeP>-1) {
+    p = virial_compute(dimension);
+    vol = xprd * yprd;
+    if (dimension == 3) vol *= zprd;
+    p *= force->nktv2p / (dimension*vol);
+    if (returnAnharmonic == -1) {
+      pLat = p;
+    }
+  }
+
+  if (domain->triclinic == 0) {
+    for (int i = 0; i < nlocal; i++) {
+      int xbox = (image[i] & IMGMASK) - IMGMAX;
+      int ybox = (image[i] >> IMGBITS & IMGMASK) - IMGMAX;
+      int zbox = (image[i] >> IMG2BITS) - IMGMAX;
+      double dx = x[i][0] + xbox*xprd - xoriginal[i][0];
+      double dy = x[i][1] + ybox*yprd - xoriginal[i][1];
+      double dz = x[i][2] + zbox*zprd - xoriginal[i][2];
+      if (comm_forward>0) {
+        deltaR[i][0] = dx;
+        deltaR[i][1] = dy;
+        deltaR[i][2] = dz;
+      }
+      fdr += dx*f[i][0] + dy*f[i][1] + dz*f[i][2];
+    }
+  }
+  else {
+    for (int i = 0; i < nlocal; i++) {
+      int xbox = (image[i] & IMGMASK) - IMGMAX;
+      int ybox = (image[i] >> IMGBITS & IMGMASK) - IMGMAX;
+      int zbox = (image[i] >> IMG2BITS) - IMGMAX;
+      double dx = x[i][0] + h[0]*xbox + h[5]*ybox + h[4]*zbox - xoriginal[i][0];
+      double dy = x[i][1] + h[1]*ybox + h[3]*zbox - xoriginal[i][1];
+      double dz = x[i][2] + h[2]*zbox - xoriginal[i][2];
+      if (comm_forward>0) {
+        deltaR[i][0] = dx;
+        deltaR[i][1] = dy;
+        deltaR[i][2] = dz;
+      }
+      fdr += ((dx*f[i][0])+(dy*f[i][1])+(dz*f[i][2]));
+    }
+  }
+
+  double phiSum = 0.0;
+  if (computeCv>-1) {
+    comm->forward_comm_compute(this);
+    int *type = atom->type;
+    double** cutsq = force->pair->cutsq;
+    if (force->pair) {
+      double **x = atom->x;
+      double **f = atom->f;
+      int *type = atom->type;
+      int nlocal = atom->nlocal;
+      double *special_lj = force->special_lj;
+      double *special_coul = force->special_coul;
+      int newton_pair = force->newton_pair;
+
+      if (update->firststep == update->ntimestep) neighbor->build_one(list,1);
+      else neighbor->build_one(list);
+      int inum = list->inum;
+      int *ilist = list->ilist;
+      int *numneigh = list->numneigh;
+      int **firstneigh = list->firstneigh;
+
+      for (int ii = 0; ii < inum; ii++) {
+        int i = ilist[ii];
+        double fac = (newton_pair || i < nlocal) ? 1.0 : 0.5;
+        double* ix = x[i];
+        int itype = type[i];
+        int *jlist = firstneigh[i];
+        int jnum = numneigh[i];
+        double *idr = deltaR[i];
+        for (int jj = 0; jj < jnum; jj++) {
+          int j = jlist[jj];
+          if (!newton_pair && j>=nlocal) fac -= 0.5;
+          double factor_lj = special_lj[sbmask(j)];
+          double factor_coul = special_coul[sbmask(j)];
+          j &= NEIGHMASK;
+          double* jx = x[j];
+          double delr[3];
+          delr[0] = ix[0] - jx[0];
+          delr[1] = ix[1] - jx[1];
+          delr[2] = ix[2] - jx[2];
+          double rsq = delr[0]*delr[0] + delr[1]*delr[1] + delr[2]*delr[2];
+          int jtype = type[j];
+          if (rsq < cutsq[itype][jtype]) {
+            double* jdr = deltaR[j];
+            double fforce, d2u[6];
+            force->pair->single_hessian(i, j, itype, jtype, rsq, delr, factor_coul, factor_lj, fforce, d2u);
+            int m = 0;
+            for (int k=0; k<3; k++) {
+              double a = fac;
+              for (int l=k; l<3; l++) {
+                phiSum += a*(idr[k]*jdr[l]+jdr[k]*idr[l])*d2u[m];
+                phiSum -= a*(idr[k]*idr[l]*d2u[m] + jdr[k]*jdr[l]*d2u[m]);
+                m++;
+                if (k==l) a *= 2;
+              }
+            }
+          }
+        }
+      }
+    }
+  }
+
+  // compute and sum up properties across processors
+
+  double fdrTotal;
+  MPI_Allreduce(&fdr,&fdrTotal,1,MPI_DOUBLE,MPI_SUM,world);
+  double uTotal;
+  if (computeU>-1 || computeCv>-1) {
+    MPI_Allreduce(&u,&uTotal,1,MPI_DOUBLE,MPI_SUM,world);
+    if (returnAnharmonic == -1) {
+      uLat = uTotal;
+    }
+    if (computeU>-1) {
+      if (returnAnharmonic) {
+        vector[computeU] = uTotal - uLat + 0.5*fdrTotal;
+      }
+      else {
+        vector[computeU] = uTotal + 0.5*fdrTotal + 0.5*dimension*(atom->natoms - 1)*force->boltz*finaltemp;
+      }
+    }
+  }
+
+  if (computeP>-1) {
+    double fv = ((deltaPcap)-(force->boltz*finaltemp*force->nktv2p*atom->natoms/vol))/(force->boltz*finaltemp*dimension*(atom->natoms - 1));
+    if (returnAnharmonic) {
+      vector[computeP] = p - pLat + (fv*fdrTotal);
+    }
+    else {
+      vector[computeP] = p + (fv*fdrTotal) + deltaPcap;
+    }
+  }
+
+  if (computeCv>-1) {
+    if (computeU==-1) MPI_Allreduce(&u,&uTotal,1,MPI_DOUBLE,MPI_SUM,world);
+    double buTot;
+    if (returnAnharmonic) {
+      buTot = (uTotal - uLat + 0.5*fdrTotal)/finaltemp;
+    }
+    else {
+      buTot = (uTotal + 0.5*fdrTotal)/finaltemp + 0.5*dimension*(atom->natoms - 1)*force->boltz;
+    }
+    double one = -0.25*(fdr + phiSum)/finaltemp;
+    double Cv;
+    MPI_Allreduce(&one,&Cv,1,MPI_DOUBLE,MPI_SUM,world);
+    vector[computeCv] = Cv + buTot*buTot;
+    if (!returnAnharmonic) {
+      vector[computeCv] += 0.5*dimension*(atom->natoms-1);
+    }
+  }
+  if (returnAnharmonic == -1) {
+    // we have our lattice properties
+    returnAnharmonic = 1;
+  }
+}
+
+double ComputeHMA::virial_compute(int n)
+{
+  double v = 0;
+
+  // sum contributions to virial from forces and fixes
+
+  if (force->pair) v += sumVirial(n, force->pair->virial);
+  if (force->bond) v += sumVirial(n, force->bond->virial);
+  if (force->angle) v += sumVirial(n, force->angle->virial);
+  if (force->dihedral) v += sumVirial(n, force->dihedral->virial);
+  if (force->improper) v += sumVirial(n, force->improper->virial);
+  for (int i = 0; i < modify->nfix; i++)
+    if (modify->fix[i]->thermo_virial) v += sumVirial(n, modify->fix[i]->virial);
+
+  // sum virial across procs
+
+  double virial;
+  MPI_Allreduce(&v,&virial,1,MPI_DOUBLE,MPI_SUM,world);
+
+  // KSpace virial contribution is already summed across procs
+
+  if (force->kspace)
+    for (int i = 0; i < n; i++) virial += force->kspace->virial[i];
+  return virial;
+}
+
+/* ---------------------------------------------------------------------- */
+
+int ComputeHMA::pack_forward_comm(int n, int *list, double *buf,
+                                        int pbc_flag, int *pbc)
+{
+  double **xoriginal = fix->astore;
+  imageint *image = atom->image;
+  double **x = atom->x;
+  double *h = domain->h;
+  double xprd = domain->xprd;   
+  double yprd = domain->yprd;
+  double zprd = domain->zprd;
+
+  int m = 0;
+  for (int ii = 0; ii < n; ii++) {
+    int i = list[ii];
+    buf[m++] = deltaR[i][0];
+    buf[m++] = deltaR[i][1];
+    buf[m++] = deltaR[i][2];
+  }
+  return m;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void ComputeHMA::unpack_forward_comm(int n, int first, double *buf)
+{
+  double **xoriginal = fix->astore;
+  int i,m,last;
+
+  m = 0;
+  last = first + n;
+  for (i = first; i < last; i++) {
+    deltaR[i][0] = buf[m++];
+    deltaR[i][1] = buf[m++];
+    deltaR[i][2] = buf[m++];
+  }
+}
+
+
+/* ----------------------------------------------------------------------
+   initialize one atom's storage values, called when atom is created
+------------------------------------------------------------------------- */
+
+void ComputeHMA::set_arrays(int i)
+{
+  double **xoriginal = fix->astore;
+  double **x = atom->x;
+  xoriginal[i][0] = x[i][0];
+  xoriginal[i][1] = x[i][1];
+  xoriginal[i][2] = x[i][2];
+}
+
+double ComputeHMA::memory_usage()
+{
+  double bytes = nmax * 3 * sizeof(double);
+  return bytes;
+}

src/USER-MISC/compute_hma.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 COMPUTE_CLASS
+
+ComputeStyle(hma,ComputeHMA)
+
+#else
+
+#ifndef LMP_COMPUTE_HMA_H
+#define LMP_COMPUTE_HMA_H
+
+#include "compute.h"
+
+namespace LAMMPS_NS {
+
+class ComputeHMA : public Compute {
+ public:
+  ComputeHMA(class LAMMPS *, int, char **);
+  ~ComputeHMA();
+  void setup();
+  void init();
+  void init_list(int, class NeighList *);
+  void compute_vector();
+  void set_arrays(int);
+  int pack_forward_comm(int, int *, double *, int, int *);
+  void unpack_forward_comm(int, int, double *);
+  double memory_usage();
+
+ private:
+  int nmax;
+  int atomsingroup;
+  char *id_fix;
+  char *id_temp;
+  double finaltemp;
+  class FixStore *fix;
+  double boltz, nktv2p, inv_volume;
+  double deltaPcap;
+  double virial_compute(int);
+  static double sumVirial(int n, double* v) {
+    double x = 0;
+    for (int i=0; i<n; i++) x += v[i];
+    return x;
+  }
+  int computeU, computeP, computeCv;
+  class NeighList *list; // half neighbor list
+  double **deltaR;
+  int returnAnharmonic;
+  double uLat, pLat;
+};
+
+}
+
+#endif
+#endif
+
+ 

src/USER-MISC/pair_cosine_squared.cpp

@@ -0,0 +1,487 @@
+/* ----------------------------------------------------------------------
+   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+   http://lammps.sandia.gov, Sandia National Laboratories
+   Steve Plimpton, sjplimp@sandia.gov
+
+   Copyright (2003) Sandia Corporation.  Under the terms of Contract
+   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+   certain rights in this software.  This software is distributed under
+   the GNU General Public License.
+
+   See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+/* ----------------------------------------------------------------------
+   Contributing authors: Eugen Rozic (University College London)
+------------------------------------------------------------------------- */
+
+#include "pair_cosine_squared.h"
+#include <cmath>
+#include <cstdlib>
+#include <cstring>
+#include "atom.h"
+#include "comm.h"
+#include "force.h"
+#include "neighbor.h"
+#include "neigh_list.h"
+#include "neigh_request.h"
+#include "update.h"
+#include "integrate.h"
+#include "respa.h"
+#include "math_const.h"
+#include "memory.h"
+#include "error.h"
+
+using namespace LAMMPS_NS;
+using namespace MathConst;
+
+/* ---------------------------------------------------------------------- */
+
+PairCosineSquared::PairCosineSquared(LAMMPS *lmp) : Pair(lmp)
+{
+  writedata = 1;
+}
+
+/* ---------------------------------------------------------------------- */
+
+PairCosineSquared::~PairCosineSquared()
+{
+  if (allocated) {
+    memory->destroy(setflag);
+    memory->destroy(cutsq);
+
+    memory->destroy(epsilon);
+    memory->destroy(sigma);
+    memory->destroy(w);
+    memory->destroy(cut);
+    memory->destroy(wcaflag);
+
+    memory->destroy(lj12_e);
+    memory->destroy(lj6_e);
+    memory->destroy(lj12_f);
+    memory->destroy(lj6_f);
+  }
+}
+
+/* ----------------------------------------------------------------------
+   allocate all arrays
+------------------------------------------------------------------------- */
+
+void PairCosineSquared::allocate()
+{
+  allocated = 1;
+  int n = atom->ntypes;
+  memory->create(setflag, n+1, n+1, "pair:setflag");
+  memory->create(cutsq, n+1, n+1, "pair:cutsq");
+
+  memory->create(cut, n+1, n+1, "pair:cut");
+  memory->create(epsilon, n+1, n+1, "pair:epsilon");
+  memory->create(sigma, n+1, n+1, "pair:sigma");
+  memory->create(w, n+1, n+1, "pair:w");
+  memory->create(wcaflag, n+1, n+1, "pair:wcaflag");
+
+  memory->create(lj12_e, n+1, n+1, "pair:lj12_e");
+  memory->create(lj6_e, n+1, n+1, "pair:lj6_e");
+  memory->create(lj12_f, n+1, n+1, "pair:lj12_f");
+  memory->create(lj6_f, n+1, n+1, "pair:lj6_f");
+
+  for (int i = 1; i <= n; i++) {
+    for (int j = i; j <= n; j++) {
+      setflag[i][j] = 0;
+      wcaflag[i][j] = 0;
+    }
+  }
+}
+
+/* ----------------------------------------------------------------------
+   global settings
+------------------------------------------------------------------------- */
+
+void PairCosineSquared::settings(int narg, char **arg)
+{
+  if (narg != 1) {
+    error->all(FLERR, "Illegal pair_style command (wrong number of params)");
+  }
+
+  cut_global = force->numeric(FLERR, arg[0]);
+
+  // reset cutoffs that have been explicitly set
+
+  if (allocated) {
+    int i, j;
+    for (i = 1; i <= atom->ntypes; i++)
+      for (j = i+1; j <= atom->ntypes; j++)
+        if (setflag[i][j])
+          cut[i][j] = cut_global;
+  }
+}
+
+
+/* ----------------------------------------------------------------------
+   set coeffs for one or more type pairs
+------------------------------------------------------------------------- */
+
+void PairCosineSquared::coeff(int narg, char **arg)
+{
+  if (narg < 4 || narg > 6)
+    error->all(FLERR, "Incorrect args for pair coefficients (too few or too many)");
+  
+  if (!allocated)
+    allocate();
+
+  int ilo, ihi, jlo, jhi;
+  force->bounds(FLERR, arg[0], atom->ntypes, ilo, ihi);
+  force->bounds(FLERR, arg[1], atom->ntypes, jlo, jhi);
+
+  double epsilon_one = force->numeric(FLERR, arg[2]);
+  double sigma_one = force->numeric(FLERR, arg[3]);
+
+  double cut_one = cut_global;
+  double wca_one = 0;
+  if (narg == 6) {
+    cut_one = force->numeric(FLERR, arg[4]);
+    if (strcmp(arg[5], "wca") == 0) {
+      wca_one = 1;
+    } else {
+      error->all(FLERR, "Incorrect args for pair coefficients (unknown option)");
+    }
+  } else if (narg == 5) {
+    if (strcmp(arg[4], "wca") == 0) {
+      wca_one = 1;
+    } else {
+      cut_one = force->numeric(FLERR, arg[4]);
+    }
+  }
+
+  if (cut_one < sigma_one) {
+    error->all(FLERR, "Incorrect args for pair coefficients (cutoff < sigma)");
+  } else if (cut_one == sigma_one) {
+    if (wca_one == 0) {
+      error->all(FLERR, "Incorrect args for pair coefficients (cutoff = sigma w/o wca)");
+    } else {
+      error->warning(FLERR, "Cosine/squared set to WCA only (cutoff = sigma)");
+    }
+  }
+
+  int count = 0;
+  for (int i = ilo; i <= ihi; i++) {
+    for (int j = MAX(jlo,i); j <= jhi; j++) {
+      epsilon[i][j] = epsilon_one;
+      sigma[i][j] = sigma_one;
+      cut[i][j] = cut_one;
+      wcaflag[i][j] = wca_one;
+      setflag[i][j] = 1;
+      count++;
+    }
+  }
+
+  if (count == 0)
+    error->all(FLERR, "Incorrect args for pair coefficients (none set)");
+}
+
+/* ----------------------------------------------------------------------
+   init specific to this pair style (unneccesary)
+------------------------------------------------------------------------- */
+
+/*
+void PairCosineSquared::init_style()
+{
+  neighbor->request(this,instance_me);
+}
+*/
+
+/* ----------------------------------------------------------------------
+   init for one type pair i,j and corresponding j,i
+------------------------------------------------------------------------- */
+
+double PairCosineSquared::init_one(int i, int j)
+{
+  if (setflag[i][j] == 0)
+    error->all(FLERR, "Mixing not supported in pair_style cosine/squared");
+
+  epsilon[j][i] = epsilon[i][j];
+  sigma[j][i] = sigma[i][j];
+  cut[j][i] = cut[i][j];
+  wcaflag[j][i] = wcaflag[i][j];
+
+  w[j][i] = w[i][j] = cut[i][j] - sigma[i][j];
+
+  if (wcaflag[i][j]) {
+    lj12_e[j][i] = lj12_e[i][j] = epsilon[i][j] * pow(sigma[i][j], 12.0);
+    lj6_e[j][i] = lj6_e[i][j] = 2.0 * epsilon[i][j] * pow(sigma[i][j], 6.0);
+    lj12_f[j][i] = lj12_f[i][j] = 12.0 * epsilon[i][j] * pow(sigma[i][j], 12.0);
+    lj6_f[j][i] = lj6_f[i][j] = 12.0 * epsilon[i][j] * pow(sigma[i][j], 6.0);
+  }
+
+  // Note: cutsq is set in pair.cpp
+
+  return cut[i][j];
+}
+
+/* ----------------------------------------------------------------------
+   this is here to throw errors & warnings for given options
+------------------------------------------------------------------------- */
+
+void PairCosineSquared::modify_params(int narg, char **arg)
+{
+  Pair::modify_params(narg, arg);
+
+  int iarg = 0;
+  while (iarg < narg) {
+    if (strcmp(arg[iarg], "mix") == 0) {
+      error->all(FLERR, "pair_modify mix not supported for pair_style cosine/squared");
+    } else if (strcmp(arg[iarg], "shift") == 0) {
+      error->warning(FLERR, "pair_modify shift has no effect on pair_style cosine/squared");
+      offset_flag = 0;
+    } else if (strcmp(arg[iarg], "tail") == 0) {
+      error->warning(FLERR, "pair_modify tail has no effect on pair_style cosine/squared");
+      tail_flag = 0;
+    }
+    iarg++;
+  }
+}
+
+/* ----------------------------------------------------------------------
+   proc 0 writes to restart file
+------------------------------------------------------------------------- */
+
+void PairCosineSquared::write_restart(FILE *fp)
+{
+  write_restart_settings(fp);
+
+  int i, j;
+  for (i = 1; i <= atom->ntypes; i++)
+    for (j = i; j <= atom->ntypes; j++) {
+      fwrite(&setflag[i][j], sizeof(int), 1, fp);
+      if (setflag[i][j]) {
+        fwrite(&epsilon[i][j], sizeof(double), 1, fp);
+        fwrite(&sigma[i][j], sizeof(double), 1, fp);
+        fwrite(&cut[i][j], sizeof(double), 1, fp);
+        fwrite(&wcaflag[i][j], sizeof(int), 1, fp);
+      }
+    }
+}
+
+/* ----------------------------------------------------------------------
+   proc 0 reads from restart file, bcasts
+------------------------------------------------------------------------- */
+
+void PairCosineSquared::read_restart(FILE *fp)
+{
+  read_restart_settings(fp);
+  allocate();
+
+  int i,j;
+  int me = comm->me;
+  for (i = 1; i <= atom->ntypes; i++) {
+    for (j = i; j <= atom->ntypes; j++) {
+      if (me == 0)
+        fread(&setflag[i][j], sizeof(int), 1, fp);
+      MPI_Bcast(&setflag[i][j], 1, MPI_INT, 0, world);
+      if (setflag[i][j]) {
+        if (me == 0) {
+          fread(&epsilon[i][j], sizeof(double), 1, fp);
+          fread(&sigma[i][j], sizeof(double), 1, fp);
+          fread(&cut[i][j], sizeof(double), 1, fp);
+          fread(&wcaflag[i][j], sizeof(int), 1, fp);
+        }
+        MPI_Bcast(&epsilon[i][j], 1, MPI_DOUBLE, 0, world);
+        MPI_Bcast(&sigma[i][j], 1, MPI_DOUBLE, 0, world);
+        MPI_Bcast(&cut[i][j], 1, MPI_DOUBLE, 0, world);
+        MPI_Bcast(&wcaflag[i][j], 1, MPI_INT, 0, world);
+      }
+    }
+  }
+}
+
+/* ----------------------------------------------------------------------
+   proc 0 writes to restart file
+------------------------------------------------------------------------- */
+
+void PairCosineSquared::write_restart_settings(FILE *fp)
+{
+  fwrite(&cut_global, sizeof(double), 1, fp);
+}
+
+/* ----------------------------------------------------------------------
+   proc 0 reads from restart file, bcasts
+------------------------------------------------------------------------- */
+
+void PairCosineSquared::read_restart_settings(FILE *fp)
+{
+  int me = comm->me;
+  if (me == 0) {
+    fread(&cut_global, sizeof(double), 1, fp);
+  }
+  MPI_Bcast(&cut_global, 1, MPI_DOUBLE, 0, world);
+}
+
+/* ----------------------------------------------------------------------
+   proc 0 writes to data file
+------------------------------------------------------------------------- */
+
+void PairCosineSquared::write_data(FILE *fp)
+{
+  for (int i = 1; i <= atom->ntypes; i++)
+    fprintf(fp, "%d %g %g %g %d\n", i, epsilon[i][i], sigma[i][i],
+            cut[i][i], wcaflag[i][i]);
+}
+
+/* ----------------------------------------------------------------------
+   proc 0 writes all pairs to data file
+------------------------------------------------------------------------- */
+
+void PairCosineSquared::write_data_all(FILE *fp)
+{
+  for (int i = 1; i <= atom->ntypes; i++)
+    for (int j = i; j <= atom->ntypes; j++)
+      fprintf(fp, "%d %d %g %g %g %d\n", i, j, epsilon[i][j], sigma[i][j],
+              cut[i][j], wcaflag[i][j]);
+}
+
+/* ---------------------------------------------------------------------- */
+
+void PairCosineSquared::compute(int eflag, int vflag)
+{
+  int i, j, ii, jj, inum, jnum, itype, jtype;
+  int *ilist, *jlist, *numneigh, **firstneigh;
+  double xtmp, ytmp, ztmp, delx, dely, delz, evdwl, fpair;
+  double r, rsq, r2inv, r6inv;
+  double factor_lj, force_lj, force_cos, cosone;
+
+  evdwl = 0.0;
+  if (eflag || vflag)
+    ev_setup(eflag, vflag);
+  else
+    evflag = vflag_fdotr = 0;
+
+  double **x = atom->x;
+  double **f = atom->f;
+  int *type = atom->type;
+  int nlocal = atom->nlocal;
+  double *special_lj = force->special_lj;
+  int newton_pair = force->newton_pair;
+
+  inum = list->inum;
+  ilist = list->ilist;
+  numneigh = list->numneigh;
+  firstneigh = list->firstneigh;
+
+  // loop over neighbors of my atoms
+
+  for (ii = 0; ii < inum; ii++) {
+    i = ilist[ii];
+    xtmp = x[i][0];
+    ytmp = x[i][1];
+    ztmp = x[i][2];
+    itype = type[i];
+    jlist = firstneigh[i];
+    jnum = numneigh[i];
+
+    for (jj = 0; jj < jnum; jj++) {
+      j = jlist[jj];
+      factor_lj = special_lj[sbmask(j)];
+      j &= NEIGHMASK;
+
+      delx = xtmp - x[j][0];
+      dely = ytmp - x[j][1];
+      delz = ztmp - x[j][2];
+      rsq = delx*delx + dely*dely + delz*delz;
+      jtype = type[j];
+
+      if (rsq < cutsq[itype][jtype]) {
+
+        /*
+          This is exactly what the "single" method does, in fact it could be called
+          here instead of repeating the code but here energy calculation is optional
+          so a little bit of calculation is possibly saved
+        */
+
+        r = sqrt(rsq);
+
+        if (r <= sigma[itype][jtype]) {
+          if (wcaflag[itype][jtype]) {
+            r2inv = 1.0/rsq;
+            r6inv = r2inv*r2inv*r2inv;
+            force_lj = r6inv*(lj12_f[itype][jtype]*r6inv - lj6_f[itype][jtype]);
+            fpair = factor_lj*force_lj*r2inv;
+            if (eflag) {
+              evdwl = factor_lj*r6inv*
+                      (lj12_e[itype][jtype]*r6inv - lj6_e[itype][jtype]);
+              if (sigma[itype][jtype] == cut[itype][jtype]) {
+                // this is the WCA-only case (it requires this shift by definition)
+                evdwl += factor_lj*epsilon[itype][jtype];
+              }
+            }
+          } else {
+            fpair = 0.0;
+            if (eflag) {
+              evdwl = -factor_lj*epsilon[itype][jtype];
+            }
+          }
+        } else {
+          force_cos = -(MY_PI*epsilon[itype][jtype] / (2.0*w[itype][jtype])) *
+                      sin(MY_PI*(r-sigma[itype][jtype]) / w[itype][jtype]);
+          fpair = factor_lj*force_cos / r;
+          if (eflag) {
+            cosone = cos(MY_PI*(r-sigma[itype][jtype]) / (2.0*w[itype][jtype]));
+            evdwl = -factor_lj*epsilon[itype][jtype]*cosone*cosone;
+          }
+        }
+
+        f[i][0] += delx*fpair;
+        f[i][1] += dely*fpair;
+        f[i][2] += delz*fpair;
+
+        if (newton_pair || j < nlocal) {
+          f[j][0] -= delx*fpair;
+          f[j][1] -= dely*fpair;
+          f[j][2] -= delz*fpair;
+        }
+
+        if (evflag)
+          ev_tally(i, j, nlocal, newton_pair, evdwl, 0.0, fpair, delx, dely, delz);
+      }
+    }
+  }
+
+  if (vflag_fdotr)
+    virial_fdotr_compute();
+}
+
+/* ----------------------------------------------------------------------
+   This is used be pair_write;
+   it is called only if rsq < cutsq[itype][jtype], no need to check that
+------------------------------------------------------------------------- */
+
+double PairCosineSquared::single(int /* i */, int /* j */, int itype, int jtype, double rsq,
+                         double /* factor_coul */, double factor_lj,
+                         double &fforce)
+{
+  double r, r2inv, r6inv, cosone, force, energy;
+  
+  r = sqrt(rsq);
+
+  if (r <= sigma[itype][jtype]) {
+    if (wcaflag[itype][jtype]) {
+      r2inv = 1.0/rsq;
+      r6inv = r2inv*r2inv*r2inv;
+      force = r6inv*(lj12_f[itype][jtype]*r6inv - lj6_f[itype][jtype])*r2inv;
+      energy = r6inv*(lj12_e[itype][jtype]*r6inv - lj6_e[itype][jtype]);
+      if (sigma[itype][jtype] == cut[itype][jtype]) {
+        // this is the WCA-only case (it requires this shift by definition)
+        energy += epsilon[itype][jtype];
+      }
+    } else {
+      force = 0.0;
+      energy = -epsilon[itype][jtype];
+    }
+  } else {
+    cosone = cos(MY_PI*(r-sigma[itype][jtype]) / (2.0*w[itype][jtype]));
+    force = -(MY_PI*epsilon[itype][jtype] / (2.0*w[itype][jtype])) * 
+                 sin(MY_PI*(r-sigma[itype][jtype]) / w[itype][jtype]) / r;
+    energy = -epsilon[itype][jtype]*cosone*cosone;
+  }
+  fforce = factor_lj*force;
+  return factor_lj*energy;
+}
+

src/USER-MISC/pair_cosine_squared.h

@@ -0,0 +1,100 @@
+/* -*- c++ -*- ----------------------------------------------------------
+   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+   http://lammps.sandia.gov, Sandia National Laboratories
+   Steve Plimpton, sjplimp@sandia.gov
+
+   Copyright (2003) Sandia Corporation.  Under the terms of Contract
+   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+   certain rights in this software.  This software is distributed under
+   the GNU General Public License.
+
+   See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+/* ----------------------------------------------------------------------
+   Contributing authors: Eugen Rozic (University College London)
+------------------------------------------------------------------------- */
+
+#ifdef PAIR_CLASS
+
+PairStyle(cosine/squared, PairCosineSquared)
+
+#else
+
+#ifndef LMP_PAIR_LJ_COS_SQ_H
+#define LMP_PAIR_LJ_COS_SQ_H
+
+#include "pair.h"
+
+namespace LAMMPS_NS {
+
+class PairCosineSquared : public Pair {
+ public:
+  PairCosineSquared(class LAMMPS *);
+  virtual ~PairCosineSquared();
+  void settings(int, char **);
+  void coeff(int, char **);
+  // void init_style();
+  double init_one(int, int);
+  void modify_params(int, char **);
+  void write_restart(FILE *);
+  void read_restart(FILE *);
+  void write_restart_settings(FILE *);
+  void read_restart_settings(FILE *);
+  void write_data(FILE *);
+  void write_data_all(FILE *);
+  virtual void compute(int, int);
+  double single(int, int, int, int, double, double, double, double &);
+  // void *extract(const char *, int &);
+
+/* RESPA stuff not implemented...
+  void compute_inner();
+  void compute_middle();
+  void compute_outer(int, int);
+*/
+
+ protected:
+  double cut_global;
+  double **epsilon, **sigma, **w, **cut;
+  int **wcaflag;
+  double **lj12_e, **lj6_e, **lj12_f, **lj6_f;
+
+  virtual void allocate();
+};
+
+}  // namespace LAMMPS_NS
+
+#endif
+#endif
+
+/* ERROR/WARNING messages:
+
+E: Illegal ... command
+
+Self-explanatory.  Check the input script syntax and compare to the
+documentation for the command.  You can use -echo screen as a
+command-line option when running LAMMPS to see the offending line.
+
+E: Incorrect args for pair coefficients
+
+Self-explanatory.  Check the input script or data file.
+
+E: Mixing not supported in pair_style cosine/squared
+
+Self-explanatory. All coefficients need to be specified explicitly.
+
+E: pair_modify mix not supported for pair_style cosine/squared
+
+Same as above, only when calling "pair_modify" command
+
+W: pair_modify shift/tail is meaningless for pair_style cosine/squared
+
+This style by definition gets to zero at cutoff distance, so there is nothing
+to shift and there is no tail contribution
+
+W: Cosine/squared set to WCA only (cutoff = sigma)
+
+If cutoff is equal to sigma (minimum) then this pair style basically
+degenerates/reverts to only WCA. This is for convenience.
+
+*/
+

src/comm_tiled.h

@@ -14,7 +14,6 @@
 #ifndef LMP_COMM_TILED_H
 #define LMP_COMM_TILED_H
 
-#include <mpi.h>
 #include "comm.h"
 
 namespace LAMMPS_NS {

src/fix_wall.cpp

@@ -23,6 +23,7 @@
 #include "respa.h"
 #include "error.h"
 #include "force.h"
+#include "utils.h"
 
 using namespace LAMMPS_NS;
 using namespace FixConst;
@@ -99,6 +100,19 @@ FixWall::FixWall(LAMMPS *lmp, int narg, char **arg) :
         estyle[nwall] = CONSTANT;
       }
 
+      if (utils::strmatch(style,"^wall/morse")) {
+        if (strstr(arg[iarg+3],"v_") == arg[iarg+3]) {
+          int n = strlen(&arg[iarg+3][2]) + 1;
+          astr[nwall] = new char[n];
+          strcpy(astr[nwall],&arg[iarg+3][2]);
+          astyle[nwall] = VARIABLE;
+        } else {
+          alpha[nwall] = force->numeric(FLERR,arg[iarg+3]);
+          astyle[nwall] = CONSTANT;
+        }
+        ++iarg;
+      }
+
       if (strstr(arg[iarg+3],"v_") == arg[iarg+3]) {
         int n = strlen(&arg[iarg+3][2]) + 1;
         sstr[nwall] = new char[n];

src/fix_wall.h

@@ -45,12 +45,12 @@ class FixWall : public Fix {
   virtual void wall_particle(int, int, double) = 0;
 
  protected:
-  double epsilon[6],sigma[6],cutoff[6];
+  double epsilon[6],sigma[6],alpha[6],cutoff[6];
   double ewall[7],ewall_all[7];
   double xscale,yscale,zscale;
-  int estyle[6],sstyle[6],wstyle[6];
+  int estyle[6],sstyle[6],astyle[6],wstyle[6];
   int eindex[6],sindex[6];
-  char *estr[6],*sstr[6];
+  char *estr[6],*sstr[6],*astr[6];
   int varflag;                // 1 if any wall position,epsilon,sigma is a var
   int eflag;                  // per-wall flag for energy summation
   int ilevel_respa;

src/fix_wall_morse.cpp

@@ -0,0 +1,87 @@
+/* ----------------------------------------------------------------------
+   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 "fix_wall_morse.h"
+#include <cmath>
+#include "atom.h"
+#include "error.h"
+
+using namespace LAMMPS_NS;
+using namespace FixConst;
+
+/* ---------------------------------------------------------------------- */
+
+FixWallMorse::FixWallMorse(LAMMPS *lmp, int narg, char **arg) :
+  FixWall(lmp, narg, arg)
+{
+  dynamic_group_allow = 1;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixWallMorse::precompute(int m)
+{
+  coeff1[m] = 2.0 * epsilon[m] * alpha[m];
+  const double alpha_dr = -alpha[m] * (cutoff[m] - sigma[m]);
+  offset[m] = epsilon[m] * (exp(2.0*alpha_dr) - 2.0*exp(alpha_dr));
+}
+
+/* ----------------------------------------------------------------------
+   interaction of all particles in group with a wall
+   m = index of wall coeffs
+   which = xlo,xhi,ylo,yhi,zlo,zhi
+   error if any particle is on or behind wall
+------------------------------------------------------------------------- */
+
+void FixWallMorse::wall_particle(int m, int which, double coord)
+{
+  double delta,fwall;
+  double vn;
+
+  double **x = atom->x;
+  double **f = atom->f;
+  int *mask = atom->mask;
+  int nlocal = atom->nlocal;
+
+  int dim = which / 2;
+  int side = which % 2;
+  if (side == 0) side = -1;
+
+  int onflag = 0;
+
+  for (int i = 0; i < nlocal; i++) {
+    if (mask[i] & groupbit) {
+      if (side < 0) delta = x[i][dim] - coord;
+      else delta = coord - x[i][dim];
+      if (delta >= cutoff[m]) continue;
+      if (delta <= 0.0) {
+        onflag = 1;
+        continue;
+      }
+      double dr = delta - sigma[m];
+      double dexp = exp(-alpha[m] * dr);
+      fwall = side * coeff1[m] * (dexp*dexp - dexp) / delta;
+      ewall[0] += epsilon[m] * (dexp*dexp - 2.0*dexp) - offset[m];
+      f[i][dim] -= fwall;
+      ewall[m+1] += fwall;
+
+      if (evflag) {
+        if (side < 0) vn = -fwall*delta;
+        else vn = fwall*delta;
+        v_tally(dim, i, vn);
+      }
+    }
+  }
+
+  if (onflag) error->one(FLERR,"Particle on or inside fix wall surface");
+}

src/fix_wall_morse.h

@@ -0,0 +1,49 @@
+/* -*- 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(wall/morse,FixWallMorse)
+
+#else
+
+#ifndef LMP_FIX_WALL_MORSE_H
+#define LMP_FIX_WALL_MORSE_H
+
+#include "fix_wall.h"
+
+namespace LAMMPS_NS {
+
+class FixWallMorse : public FixWall {
+ public:
+  FixWallMorse(class LAMMPS *, int, char **);
+  void precompute(int);
+  void wall_particle(int, int, double);
+
+ private:
+  double coeff1[6],offset[6];
+};
+
+}
+
+#endif
+#endif
+
+/* ERROR/WARNING messages:
+
+E: Particle on or inside fix wall surface
+
+Particles must be "exterior" to the wall in order for energy/force to
+be calculated.
+
+*/

src/fix_wall_region.cpp

@@ -28,7 +28,7 @@ using namespace LAMMPS_NS;
 using namespace FixConst;
 using namespace MathConst;
 
-enum{LJ93,LJ126,LJ1043,COLLOID,HARMONIC};
+enum{LJ93,LJ126,LJ1043,COLLOID,HARMONIC,MORSE};
 
 /* ---------------------------------------------------------------------- */
 
@@ -36,7 +36,7 @@ FixWallRegion::FixWallRegion(LAMMPS *lmp, int narg, char **arg) :
   Fix(lmp, narg, arg),
   idregion(NULL)
 {
-  if (narg != 8) error->all(FLERR,"Illegal fix wall/region command");
+  if (narg < 8) error->all(FLERR,"Illegal fix wall/region command");
 
   scalar_flag = 1;
   vector_flag = 1;
@@ -62,13 +62,28 @@ FixWallRegion::FixWallRegion(LAMMPS *lmp, int narg, char **arg) :
   else if (strcmp(arg[4],"lj1043") == 0) style = LJ1043;
   else if (strcmp(arg[4],"colloid") == 0) style = COLLOID;
   else if (strcmp(arg[4],"harmonic") == 0) style = HARMONIC;
+  else if (strcmp(arg[4],"morse") == 0) style = MORSE;
   else error->all(FLERR,"Illegal fix wall/region command");
 
   if (style != COLLOID) dynamic_group_allow = 1;
 
-  epsilon = force->numeric(FLERR,arg[5]);
+  if (style == MORSE) {
-  sigma = force->numeric(FLERR,arg[6]);
+    if (narg != 9)
-  cutoff = force->numeric(FLERR,arg[7]);
+      error->all(FLERR,"Illegal fix wall/region command");
+
+    epsilon = force->numeric(FLERR,arg[5]);
+    alpha = force->numeric(FLERR,arg[6]);
+    sigma = force->numeric(FLERR,arg[7]);
+    cutoff = force->numeric(FLERR,arg[8]);
+
+  } else {
+    if (narg != 8)
+      error->all(FLERR,"Illegal fix wall/region command");
+
+    epsilon = force->numeric(FLERR,arg[5]);
+    sigma = force->numeric(FLERR,arg[6]);
+    cutoff = force->numeric(FLERR,arg[7]);
+  }
 
   if (cutoff <= 0.0) error->all(FLERR,"Fix wall/region cutoff <= 0.0");
 
@@ -154,12 +169,15 @@ void FixWallRegion::init()
     coeff5 = coeff1 * 10.0;
     coeff6 = coeff2 * 4.0;
     coeff7 = coeff3 * 3.0;
-
     double rinv = 1.0/cutoff;
     double r2inv = rinv*rinv;
     double r4inv = r2inv*r2inv;
     offset = coeff1*r4inv*r4inv*r2inv - coeff2*r4inv -
       coeff3*pow(cutoff+coeff4,-3.0);
+  } else if (style == MORSE) {
+    coeff1 = 2 * epsilon * alpha;
+    double alpha_dr = -alpha * (cutoff - sigma);
+    offset = epsilon * (exp(2.0*alpha_dr) - 2.0*exp(alpha_dr));
   } else if (style == COLLOID) {
     coeff1 = -4.0/315.0 * epsilon * pow(sigma,6.0);
     coeff2 = -2.0/3.0 * epsilon;
@@ -250,6 +268,7 @@ void FixWallRegion::post_force(int vflag)
         if (style == LJ93) lj93(region->contact[m].r);
         else if (style == LJ126) lj126(region->contact[m].r);
         else if (style == LJ1043) lj1043(region->contact[m].r);
+        else if (style == MORSE) morse(region->contact[m].r);
         else if (style == COLLOID) colloid(region->contact[m].r,radius[i]);
         else harmonic(region->contact[m].r);
 
@@ -284,7 +303,7 @@ void FixWallRegion::post_force(int vflag)
 
 /* ---------------------------------------------------------------------- */
 
-void FixWallRegion::post_force_respa(int vflag, int ilevel, int /*iloop*/)
+void FixWallRegion::post_force_respa(int vflag, int ilevel, int /* iloop */)
 {
   if (ilevel == ilevel_respa) post_force(vflag);
 }
@@ -372,6 +391,19 @@ void FixWallRegion::lj1043(double r)
     coeff3*pow(r+coeff4,-3.0) - offset;
 }
 
+/* ----------------------------------------------------------------------
+   Morse interaction for particle with wall
+   compute eng and fwall = magnitude of wall force
+------------------------------------------------------------------------- */
+
+void FixWallRegion::morse(double r)
+{
+  double dr = r - sigma;
+  double dexp = exp(-alpha * dr);
+  fwall = coeff1 * (dexp*dexp - dexp) / r;
+  eng = epsilon * (dexp*dexp - 2.0*dexp) - offset;
+}
+
 /* ----------------------------------------------------------------------
    colloid interaction for finite-size particle of rad with wall
    compute eng and fwall = magnitude of wall force

src/fix_wall_region.h

@@ -41,6 +41,7 @@ class FixWallRegion : public Fix {
  private:
   int style,iregion;
   double epsilon,sigma,cutoff;
+  double alpha;
   int eflag;
   double ewall[4],ewall_all[4];
   int ilevel_respa;
@@ -53,6 +54,7 @@ class FixWallRegion : public Fix {
   void lj93(double);
   void lj126(double);
   void lj1043(double);
+  void morse(double);
   void colloid(double, double);
   void harmonic(double);
 };

src/irregular.h

@@ -14,7 +14,6 @@
 #ifndef LMP_IRREGULAR_H
 #define LMP_IRREGULAR_H
 
-#include <mpi.h>
 #include "pointers.h"
 
 namespace LAMMPS_NS {

src/library.cpp

@@ -37,6 +37,7 @@
 #include "error.h"
 #include "force.h"
 #include "info.h"
+#include "fix_external.h"
 
 #if defined(LAMMPS_EXCEPTIONS)
 #include "exceptions.h"
@@ -1595,7 +1596,7 @@ void lammps_create_atoms(void *ptr, int n, tagint *id, int *type,
 
     if (lmp->atom->natoms != natoms_prev + n) {
       char str[128];
-      sprintf(str,"Library warning in lammps_create_atoms, "
+      snprintf(str, 128, "Library warning in lammps_create_atoms, "
               "invalid total atoms " BIGINT_FORMAT " " BIGINT_FORMAT,
               lmp->atom->natoms,natoms_prev+n);
       if (lmp->comm->me == 0)
@@ -1605,6 +1606,40 @@ void lammps_create_atoms(void *ptr, int n, tagint *id, int *type,
   END_CAPTURE
 }
 
+/* ----------------------------------------------------------------------
+   find fix external with given ID and set the callback function
+   and caller pointer
+------------------------------------------------------------------------- */
+
+void lammps_set_fix_external_callback(void *ptr, char *id, FixExternalFnPtr callback_ptr, void * caller)
+{
+  LAMMPS *lmp = (LAMMPS *) ptr;
+  FixExternal::FnPtr callback = (FixExternal::FnPtr) callback_ptr;
+
+  BEGIN_CAPTURE
+  {
+    int ifix = lmp->modify->find_fix(id);
+    if (ifix < 0) {
+      char str[128];
+      snprintf(str, 128, "Can not find fix with ID '%s'!", id);
+      lmp->error->all(FLERR,str);
+    }
+
+    Fix *fix = lmp->modify->fix[ifix];
+
+    if (strcmp("external",fix->style) != 0){
+      char str[128];
+      snprintf(str, 128, "Fix '%s' is not of style external!", id);
+      lmp->error->all(FLERR,str);
+    }
+
+    FixExternal * fext = (FixExternal*) fix;
+    fext->set_callback(callback, caller);
+  }
+  END_CAPTURE
+}
+
+
 // ----------------------------------------------------------------------
 // library API functions for accessing LAMMPS configuration
 // ----------------------------------------------------------------------

src/library.h

@@ -58,6 +58,14 @@ void lammps_gather_atoms_subset(void *, char *, int, int, int, int *, void *);
 void lammps_scatter_atoms(void *, char *, int, int, void *);
 void lammps_scatter_atoms_subset(void *, char *, int, int, int, int *, void *);
 
+#ifdef LAMMPS_BIGBIG
+typedef void (*FixExternalFnPtr)(void *, int64_t, int, int64_t *, double **, double **);
+void lammps_set_fix_external_callback(void *, char *, FixExternalFnPtr, void*);
+#else
+typedef void (*FixExternalFnPtr)(void *, int, int, int *, double **, double **);
+void lammps_set_fix_external_callback(void *, char *, FixExternalFnPtr, void*);
+#endif
+
 int lammps_config_has_package(char * package_name);
 int lammps_config_package_count();
 int lammps_config_package_name(int index, char * buffer, int max_size);

src/pair.cpp

@@ -54,6 +54,7 @@ Pair::Pair(LAMMPS *lmp) : Pointers(lmp)
   comm_forward = comm_reverse = comm_reverse_off = 0;
 
   single_enable = 1;
+  single_hessian_enable = 0;
   restartinfo = 1;
   respa_enable = 0;
   one_coeff = 0;
@@ -1743,6 +1744,18 @@ void Pair::init_bitmap(double inner, double outer, int ntablebits,
 
 /* ---------------------------------------------------------------------- */
 
+void Pair::hessian_twobody(double fforce, double dfac, double delr[3], double phiTensor[6]) {
+  int m = 0;
+  for (int k=0; k<3; k++) {
+    phiTensor[m] = fforce;
+    for (int l=k; l<3; l++) {
+      if (l>k) phiTensor[m] = 0;
+      phiTensor[m++] += delr[k]*delr[l] * dfac;
+    }
+  }
+}
+/* ---------------------------------------------------------------------- */
+
 double Pair::memory_usage()
 {
   double bytes = comm->nthreads*maxeatom * sizeof(double);

src/pair.h

@@ -46,6 +46,7 @@ class Pair : protected Pointers {
   int comm_reverse_off;          // size of reverse comm even if newton off
 
   int single_enable;             // 1 if single() routine exists
+  int single_hessian_enable;     // 1 if single_hessian() routine exists
   int restartinfo;               // 1 if pair style writes restart info
   int respa_enable;              // 1 if inner/middle/outer rRESPA routines
   int one_coeff;                 // 1 if allows only one coeff * * call
@@ -148,6 +149,16 @@ class Pair : protected Pointers {
     return 0.0;
   }
 
+  void hessian_twobody(double fforce, double dfac, double delr[3], double phiTensor[6]);
+
+  virtual double single_hessian(int, int, int, int,
+                        double, double[3], double, double,
+                        double& fforce, double d2u[6]) {
+    fforce = 0.0;
+    for (int i=0; i<6; i++) d2u[i] = 0;
+    return 0.0;
+  }
+
   virtual void settings(int, char **) = 0;
   virtual void coeff(int, char **) = 0;
 

src/pair_lj_smooth_linear.cpp

@@ -29,7 +29,9 @@ using namespace LAMMPS_NS;
 
 /* ---------------------------------------------------------------------- */
 
-PairLJSmoothLinear::PairLJSmoothLinear(LAMMPS *lmp) : Pair(lmp) {}
+PairLJSmoothLinear::PairLJSmoothLinear(LAMMPS *lmp) : Pair(lmp) {
+  single_hessian_enable = 1;
+}
 
 /* ---------------------------------------------------------------------- */
 
@@ -345,3 +347,26 @@ double PairLJSmoothLinear::single(int /*i*/, int /*j*/, int itype, int jtype,
 
   return factor_lj*philj;
 }
+
+double PairLJSmoothLinear::single_hessian(int /*i*/, int /*j*/, int itype, int jtype, double rsq,
+                         double delr[3], double /*factor_coul*/, double factor_lj,
+                         double &fforce, double d2u[6])
+{
+  double r2inv,r6inv,forcelj,philj,r,rinv;
+
+  r2inv = 1.0/rsq;
+  r6inv = r2inv*r2inv*r2inv;
+  rinv  = sqrt(r2inv);
+  r     = sqrt(rsq);
+  forcelj = r6inv*(lj1[itype][jtype]*r6inv-lj2[itype][jtype]);
+  forcelj = rinv*forcelj - dljcut[itype][jtype];
+  fforce = factor_lj*forcelj*rinv;
+
+  philj = r6inv*(lj3[itype][jtype]*r6inv-lj4[itype][jtype]);
+  philj = philj - ljcut[itype][jtype]
+                + (r-cut[itype][jtype])*dljcut[itype][jtype];
+
+  double d2r = factor_lj * r6inv * (13.0*lj1[itype][jtype]*r6inv - 7.0*lj2[itype][jtype])/rsq;
+  hessian_twobody(fforce, -(fforce + d2r) / rsq, delr, d2u);
+  return factor_lj*philj;
+}

src/pair_lj_smooth_linear.h

@@ -38,6 +38,7 @@ class PairLJSmoothLinear : public Pair {
   void write_restart_settings(FILE *);
   void read_restart_settings(FILE *);
   double single(int, int, int, int, double, double, double, double &);
+  double single_hessian(int, int, int, int, double, double[3], double, double, double&, double[6]);
 
  protected:
   double cut_global;

src/rcb.h

@@ -14,7 +14,6 @@
 #ifndef LAMMPS_RCB_H
 #define LAMMPS_RCB_H
 
-#include <mpi.h>
 #include "pointers.h"
 
 namespace LAMMPS_NS {

src/read_dump.h

@@ -22,7 +22,6 @@ CommandStyle(read_dump,ReadDump)
 #ifndef LMP_READ_DUMP_H
 #define LMP_READ_DUMP_H
 
-#include <mpi.h>
 #include "pointers.h"
 
 namespace LAMMPS_NS {

src/variable.cpp

@@ -18,6 +18,7 @@
 #include <cstring>
 #include <cctype>
 #include <unistd.h>
+#include <string>
 #include "universe.h"
 #include "atom.h"
 #include "update.h"
@@ -1302,8 +1303,12 @@ double Variable::evaluate(char *str, Tree **tree, int ivar)
         if (word[0] == 'C') lowercase = 0;
 
         int icompute = modify->find_compute(word+2);
-        if (icompute < 0)
+        if (icompute < 0) {
-          print_var_error(FLERR,"Invalid compute ID in variable formula",ivar);
+          std::string mesg = "Invalid compute ID '";
+          mesg += (word+2);
+          mesg += "' in variable formula";
+          print_var_error(FLERR,mesg.c_str(),ivar);
+        }
         Compute *compute = modify->compute[icompute];
 
         // parse zero or one or two trailing brackets
@@ -1604,9 +1609,10 @@ double Variable::evaluate(char *str, Tree **tree, int ivar)
 
         int ifix = modify->find_fix(word+2);
         if (ifix < 0) {
-          char msg[128];
+          std::string mesg = "Invalid fix ID '";
-          snprintf(msg,128,"Invalid fix ID '%s' in variable formula",word+2);
+          mesg += (word+2);
-          print_var_error(FLERR,msg,ivar);
+          mesg += "' in variable formula";
+          print_var_error(FLERR,mesg.c_str(),ivar);
         }
         Fix *fix = modify->fix[ifix];
 
@@ -3792,8 +3798,12 @@ int Variable::group_function(char *word, char *contents, Tree **tree,
   // group to operate on
 
   int igroup = group->find(args[0]);
-  if (igroup == -1)
+  if (igroup == -1) {
-    print_var_error(FLERR,"Group ID in variable formula does not exist",ivar);
+    std::string mesg = "Group ID '";
+    mesg += args[0];
+    mesg += "' in variable formula does not exist";
+    print_var_error(FLERR,mesg.c_str(),ivar);
+  }
 
   // match word to group function
 
@@ -4001,8 +4011,12 @@ int Variable::group_function(char *word, char *contents, Tree **tree,
 int Variable::region_function(char *id, int ivar)
 {
   int iregion = domain->find_region(id);
-  if (iregion == -1)
+  if (iregion == -1) {
-    print_var_error(FLERR,"Region ID in variable formula does not exist",ivar);
+    std::string mesg = "Region ID '";
+    mesg += id;
+    mesg += "' in variable formula does not exist";
+    print_var_error(FLERR,mesg.c_str(),ivar);
+  }
 
   // init region in case sub-regions have been deleted
 
@@ -4080,9 +4094,10 @@ int Variable::special_function(char *word, char *contents, Tree **tree,
 
       int icompute = modify->find_compute(&args[0][2]);
       if (icompute < 0) {
-        char msg[128];
+        std::string mesg = "Invalid compute ID '";
-        snprintf(msg,128,"Invalid compute ID '%s' in variable formula",word+2);
+        mesg += (args[0]+2);
-        print_var_error(FLERR,msg,ivar);
+        mesg += "' in variable formula";
+        print_var_error(FLERR,mesg.c_str(),ivar);
       }
       compute = modify->compute[icompute];
       if (index == 0 && compute->vector_flag) {
@@ -4123,13 +4138,20 @@ int Variable::special_function(char *word, char *contents, Tree **tree,
       } else index = 0;
 
       int ifix = modify->find_fix(&args[0][2]);
-      if (ifix < 0)
+      if (ifix < 0) {
-        print_var_error(FLERR,"Invalid fix ID in variable formula",ivar);
+        std::string mesg = "Invalid fix ID '";
+        mesg += (args[0]+2);
+        mesg += "' in variable formula";
+        print_var_error(FLERR,mesg.c_str(),ivar);
+      }
       fix = modify->fix[ifix];
       if (index == 0 && fix->vector_flag) {
-        if (update->whichflag > 0 && update->ntimestep % fix->global_freq)
+        if (update->whichflag > 0 && update->ntimestep % fix->global_freq) {
-          print_var_error(FLERR,"Fix in variable not computed at "
+          std::string mesg = "Fix with ID '";
-                          "compatible time",ivar);
+          mesg += (args[0]+2);
+          mesg += "' in variable formula not computed at compatible time";
+          print_var_error(FLERR,mesg.c_str(),ivar);
+        }
         nvec = fix->size_vector;
         nstride = 1;
       } else if (index && fix->array_flag) {
@@ -4336,9 +4358,12 @@ int Variable::special_function(char *word, char *contents, Tree **tree,
       print_var_error(FLERR,"Invalid special function in variable formula",ivar);
 
     int ivar = find(args[0]);
-    if (ivar < 0)
+    if (ivar < 0) {
-      print_var_error(FLERR,"Variable ID in "
+      std::string mesg = "Variable ID '";
-                      "variable formula does not exist",ivar);
+      mesg += args[0];
+      mesg += "' in variable formula does not exist";
+      print_var_error(FLERR,mesg.c_str(),ivar);
+    }
 
     // SCALARFILE has single current value, read next one
     // save value in tree or on argstack