diff --git a/cmake/CMakeLists.txt b/cmake/CMakeLists.txt
index d7fb090e4ab977cc907d3b624dcbed62ae3fbfeb..7ce7ca23e16423be98be34edf51e4b4321931817 100644
--- a/cmake/CMakeLists.txt
+++ b/cmake/CMakeLists.txt
@@ -56,6 +56,9 @@ option(BUILD_SHARED_LIBS "Build shared libs" OFF)
 if(BUILD_SHARED_LIBS) # for all pkg libs, mpi_stubs and linalg
   set(CMAKE_POSITION_INDEPENDENT_CODE ON)
 endif()
+option(DEVELOPER_MODE "Enable developer mode" OFF)
+mark_as_advanced(DEVELOPER_MODE)
+option(CMAKE_VERBOSE_MAKEFILE "Generate verbose Makefiles" OFF)
 include(GNUInstallDirs)
 
 set(LAMMPS_LINK_LIBS)
@@ -182,7 +185,7 @@ if(PKG_KSPACE)
   if(NOT FFT STREQUAL "KISSFFT")
     find_package(${FFT} REQUIRED)
     if(NOT FFT STREQUAL "FFTW3F")
-	    add_definitions(-DFFT_FFTW)
+      add_definitions(-DFFT_FFTW)
     else()
       add_definitions(-DFFT_${FFT})
     endif()
@@ -370,13 +373,12 @@ if(PKG_USER-VTK)
 endif()
 
 if(PKG_KIM)
-  find_package(KIM QUIET)
-  if(NOT KIM_FOUND)
-    message(STATUS "KIM not found - we will build our own")
+  option(DOWNLOAD_KIM "Download kim-api (instead of using the system's one)" OFF)
+  if(DOWNLOAD_KIM)
     include(ExternalProject)
     ExternalProject_Add(kim_build
-      URL https://github.com/openkim/kim-api/archive/v1.9.4.tar.gz
-      URL_MD5 f4d35a1705eed46d64c7c0ab448ff3e0
+      URL https://github.com/openkim/kim-api/archive/v1.9.5.tar.gz
+      URL_MD5 9f66efc128da33039e30659f36fc6d00
       BUILD_IN_SOURCE 1
       CONFIGURE_COMMAND <SOURCE_DIR>/configure --prefix=<INSTALL_DIR>
       )
@@ -384,6 +386,11 @@ if(PKG_KIM)
     set(KIM_INCLUDE_DIRS ${INSTALL_DIR}/include/kim-api-v1)
     set(KIM_LIBRARIES ${INSTALL_DIR}/lib/libkim-api-v1.so)
     list(APPEND LAMMPS_DEPS kim_build)
+  else()
+    find_package(KIM)
+    if(NOT KIM_FOUND)
+      message(FATAL_ERROR "KIM not found, help CMake to find it by setting KIM_LIBRARY and KIM_INCLUDE_DIR, or set DOWNLOAD_KIM=ON to download it")
+    endif()
   endif()
   list(APPEND LAMMPS_LINK_LIBS ${KIM_LIBRARIES})
   include_directories(${KIM_INCLUDE_DIRS})
@@ -603,6 +610,29 @@ if(PKG_OPT)
 endif()
 
 if(PKG_USER-INTEL)
+    if(NOT DEVELOPER_MODE)
+      if(NOT CMAKE_CXX_COMPILER_ID STREQUAL "Intel")
+        message(FATAL_ERROR "USER-INTEL is only useful together with intel compiler")
+      endif()
+      if(CMAKE_CXX_COMPILER_VERSION VERSION_LESS 16)
+        message(FATAL_ERROR "USER-INTEL is needed at least 2016 intel compiler, found ${CMAKE_CXX_COMPILER_VERSION}")
+      endif()
+    endif()
+    option(INJECT_INTEL_FLAG "Inject OMG fast flags for USER-INTEL" ON)
+    if(INJECT_INTEL_FLAG AND CMAKE_CXX_COMPILER_ID STREQUAL "Intel")
+      if(CMAKE_CXX_COMPILER_VERSION VERSION_EQUAL 17.3 OR CMAKE_CXX_COMPILER_VERSION VERSION_EQUAL 17.4)
+        set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -xCOMMON-AVX512")
+      else()
+        set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -xHost")
+      endif()
+      include(CheckCXXCompilerFlag)
+      foreach(_FLAG -qopenmp -qno-offload -fno-alias -ansi-alias -restrict -DLMP_INTEL_USELRT -DLMP_USE_MKL_RNG -O2 "-fp-model fast=2" -no-prec-div -qoverride-limits -qopt-zmm-usage=high)
+        check_cxx_compiler_flag("${__FLAG}" COMPILER_SUPPORTS${_FLAG})
+        if(COMPILER_SUPPORTS${_FLAG})
+          set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} ${_FLAG}")
+        endif()
+      endforeach()
+    endif()
     set(USER-INTEL_SOURCES_DIR ${LAMMPS_SOURCE_DIR}/USER-INTEL)
     set(USER-INTEL_SOURCES ${USER-INTEL_SOURCES_DIR}/intel_preprocess.h
                            ${USER-INTEL_SOURCES_DIR}/intel_buffers.h
@@ -844,5 +874,5 @@ if(ENABLED_GPU)
   endif()
 endif()
 if(PKG_KSPACE)
-	message(STATUS "Using ${FFT} as FFT")
+  message(STATUS "Using ${FFT} as FFT")
 endif()
diff --git a/doc/src/fix_bond_react.txt b/doc/src/fix_bond_react.txt
index 4a7c2428c16aa050a24f71d685960d9d952b65a4..f85ef9bc1aa6ea034364e5358ac1eafa2824f0d6 100644
--- a/doc/src/fix_bond_react.txt
+++ b/doc/src/fix_bond_react.txt
@@ -20,14 +20,15 @@ ID, group-ID are documented in "fix"_fix.html command. Group-ID is ignored. :ulb
 bond/react = style name of this fix command :l
 zero or more common keyword/value pairs may be appended directly after 'bond/react' :l
 these apply to all reaction specifications (below) :l
-common_keyword = {stabilization}
-  {stabilization} values = group-ID xmax
-    group-ID = user-assigned ID of an internally-created dynamic group that excludes reacting atoms, and can be used by a subsequent time integration fix such as nvt, npt, or nve (cannot be 'all')
-  {xmax} value = distance
-    distance = xmax value that is used by an internally created "nve/limit"_fix_nve_limit.html integrator
-react = mandatory argument indicating new reaction specification
-  react-ID = user-assigned name for the reaction
-  react-group-ID = only atoms in this group are available for the reaction
+common_keyword = {stabilization} :l
+  {stabilization} values = {no} or {yes} {group-ID} {xmax}
+    {no} = no reaction site stabilization
+    {yes} = perform reaction site stabilization
+      {group-ID} = user-assigned ID for all non-reacting atoms (group created internally)
+      {xmax} = xmax value that is used by an internally created "nve/limit"_fix_nve_limit.html integrator :pre
+react = mandatory argument indicating new reaction specification :l
+  react-ID = user-assigned name for the reaction :l
+  react-group-ID = only atoms in this group are available for the reaction :l
   Nevery = attempt reaction every this many steps :l
   Rmin = bonding pair atoms must be separated by more than Rmin to initiate reaction (distance units) :l
   Rmax = bonding pair atoms must be separated by less than Rmax to initiate reaction (distance units) :l
@@ -47,7 +48,7 @@ react = mandatory argument indicating new reaction specification
 
 molecule mol1 pre_reacted_topology.txt
 molecule mol2 post_reacted_topology.txt
-fix 5 all bond/react stabilization no react myrxn1 all 1 0 3.25 mol1 mol2 map_file.txt
+fix 5 all bond/react stabilization no react myrxn1 all 1 0 3.25 mol1 mol2 map_file.txt :pre
 
 molecule mol1 pre_reacted_rxn1.txt
 molecule mol2 post_reacted_rxn1.txt
@@ -56,12 +57,12 @@ molecule mol4 post_reacted_rxn2.txt
 fix 5 all bond/react stabilization yes nvt_grp .03 &
   react myrxn1 all 1 0 3.25 mol1 mol2 map_file_rxn1.txt prob 0.50 12345 &
   react myrxn2 all 1 0 2.75 mol3 mol4 map_file_rxn2.txt prob 0.25 12345
-fix 6 nvt_grp nvt temp 300 300 100 # system-wide thermostat must be defined after bond/react :pre
+fix 6 nvt_grp nvt temp 300 300 100 # set thermostat after bond/react :pre
 
 [Description:]
 
 Initiate complex covalent bonding (topology) changes. These topology
-changes will be referred to as "reactions" throughout this
+changes will be referred to as 'reactions' throughout this
 documentation. Topology changes are defined in pre- and post-reaction
 molecule templates and can include creation and deletion of bonds,
 angles, dihedrals, impropers, bond-types, angle-types, dihedral-types,
@@ -81,10 +82,10 @@ occurred 3) build a molecule template of the reaction site after the
 reaction has occurred 4) create a map that relates the
 template-atom-IDs of each atom between pre- and post-reaction molecule
 templates 5) fill a simulation box with molecules and run a simulation
-with fix/bond react.
+with fix bond/react.
 
 Only one 'fix bond/react' command can be used at a time. Multiple
-reactions can be simultaneously applied by specifying multiple 'react'
+reactions can be simultaneously applied by specifying multiple {react}
 arguments to a single 'fix bond/react' command. This syntax is
 necessary because the 'common keywords' are applied to all reactions.
 
@@ -99,10 +100,11 @@ typically be set to the maximum distance that non-reacting atoms move
 during the simulation.
 
 The group-ID set using the {stabilization} keyword should be a
-previously unused group-ID. The fix bond/react command creates a
-"dynamic group"_group.html of this name that excludes reacting atoms.
-This dynamic group-ID should then be used by a subsequent system-wide
-time integrator, as shown in the second example above. It is currently
+previously unused group-ID. It cannot be specified as 'all'. The fix
+bond/react command creates a "dynamic group"_group.html of this name
+that includes all non-reacting atoms. This dynamic group-ID should
+then be used by a subsequent system-wide time integrator such as nvt,
+npt, or nve, as shown in the second example above. It is currently
 necessary to place the time integration command after the fix
 bond/react command due to the internal dynamic grouping performed by
 fix bond/react.
@@ -111,9 +113,9 @@ NOTE: The internally created group currently applies to all atoms in
 the system, i.e. you should generally not have a separate thermostat
 which acts on the 'all' group.
 
-The following comments pertain to each 'react' argument:
+The following comments pertain to each {react} argument:
 
-A check for possible new reaction sites is performed every Nevery
+A check for possible new reaction sites is performed every {Nevery}
 timesteps.
 
 Two conditions must be met for a reaction to occur. First a bonding
@@ -124,20 +126,20 @@ modified to match the post-reaction template.
 
 A bonding atom pair will be identified if several conditions are met.
 First, a pair of atoms within the specified react-group-ID of type
-typei and typej must separated by a distance between Rmin and Rmax. It
-is possible that multiple bonding atom pairs are identified: if the
-bonding atoms in the pre-reacted template are not 1-2, 1-3, or 1-4
-neighbors, the closest bonding atom partner is set as its bonding
-partner; otherwise, the farthest potential partner is chosen. Then, if
-both an atomi and atomj have each other as their nearest bonding
-partners, these two atoms are identified as the bonding atom pair of
-the reaction site. Once this unique bonding atom pair is identified
-for each reaction, there could two or more reactions that involve a
-given atom on the same timestep. If this is the case, only one such
-reaction is permitted to occur. This reaction is chosen randomly from
-all potential reactions. This capability allows e.g. for different
-reaction pathways to proceed from identical reaction sites with
-user-specified probabilities.
+typei and typej must separated by a distance between {Rmin} and
+{Rmax}. It is possible that multiple bonding atom pairs are
+identified: if the bonding atoms in the pre-reacted template are not
+1-2, 1-3, or 1-4 neighbors, the closest bonding atom partner is set as
+its bonding partner; otherwise, the farthest potential partner is
+chosen. Then, if both an atomi and atomj have each other as their
+nearest bonding partners, these two atoms are identified as the
+bonding atom pair of the reaction site. Once this unique bonding atom
+pair is identified for each reaction, there could two or more
+reactions that involve a given atom on the same timestep. If this is
+the case, only one such reaction is permitted to occur. This reaction
+is chosen randomly from all potential reactions. This capability
+allows e.g. for different reaction pathways to proceed from identical
+reaction sites with user-specified probabilities.
 
 The pre-reacted molecule template is specified by a molecule command.
 This molecule template file contains a sample reaction site and its
@@ -175,77 +177,43 @@ A discussion of correctly handling this is also provided on the
 
 The map file is a text document with the following format:
 
-Format of the map file
-
-A map file has a header and a body. The header appears first. The
-first line of the header is always skipped; it typically contains a
-description of the file.  Lines can have a trailing comment starting
-with '#' that is ignored. If the line is blank (only whitespace after
-comment is deleted), it is skipped. If the line contains a header
-keyword, the corresponding value(s) is read from the line. If it
-doesn't contain a header keyword, the line begins the body of the
-file.
-
-The header contains one mandatory keyword and one optional keyword.
-The mandatory keyword is 'equivalences' and the optional keyword is
-'edgeIDs.' These specify the number of atoms in the pre- and
-post-reacted templates and the number of edge atoms in pre-reacted
-template, respectively.
-
-The body contains two mandatory sections and one optional section. The
-first section begins with the keyword 'BondingIDs' and lists the atom
-IDs of the bonding atom pair in the pre-reacted molecule template. The
-second mandatory section begins with the keyword 'Equivalences' and
-lists a one-to-one correspondence between atom IDs of the pre- and
-post-reacted templates. The optional section begins with the keyword
-'EdgeIDs' and list the atom IDs of edge atoms in the pre-reacted
+A map file has a header and a body. The header of map file the
+contains one mandatory keyword and one optional keyword. The mandatory
+keyword is 'equivalences' and the optional keyword is 'edgeIDs':
+
+N {equivalences} = # of atoms N in the reaction molecule templates
+N {edgeIDs} = # of edge atoms N in the pre-reacted molecule template :pre
+
+The body of the map file contains two mandatory sections and one
+optional section. The first mandatory section begins with the keyword
+'BondingIDs' and lists the atom IDs of the bonding atom pair in the
+pre-reacted molecule template. The second mandatory section begins
+with the keyword 'Equivalences' and lists a one-to-one correspondence
+between atom IDs of the pre- and post-reacted templates. The first
+column is an atom ID of the pre-reacted molecule template, and the
+second column is the corresponding atom ID of the post-reacted
+molecule template. The optional section begins with the keyword
+'EdgeIDs' and lists the atom IDs of edge atoms in the pre-reacted
 molecule template.
 
-Format of the header of the map file
-
-These are the recognized header keywords. Header lines can come in any
-order. The value(s) are read from the beginning of the line. Thus the
-keyword 'equivalences' should be in a line like "25 equivalences."
-
-equivalences = # of atoms in the pre- and post-reacted molecule
-templates edgeIDs = # of edge atoms in the pre-reacted molecule template :pre
-
-The edgeIDs keyword is optional.
-
-Format of the body of the map file
-
-These are the section keywords for the body of the file.
-
-BondingIDs, EdgeIDs = list of atom IDs of bonding and edge atoms in
-the pre-reacted molecule template
-
-Equivalences = a two column list where the first column is an atom ID
-of the pre-reacted molecule template, and the second column is the
-corresponding atom ID of the post-reacted molecule template
-
-The bondingIDs section will always contain two atom IDs, corresponding
-to the bonding atom pairs of the pre-reacted map file. The
-Equivalences section will contain as many rows as there are atoms in
-the pre- and post-reacted molecule templates. The edgeIDs section is
-optional, but would contain an atom ID for each edge atom in the
-pre-reacted molecule template.
-
 A sample map file is given below:
 
 :line
 
-# This is a map file :pre
+# this is a map file :pre
 
 2 edgeIDs
 7 equivalences :pre
 
 BondingIDs :pre
 
-3 5 :pre
+3
+5 :pre
 
 EdgeIDs :pre
 
-1 7 :pre
+1
+7 :pre
 
 Equivalences :pre
 
@@ -264,13 +232,13 @@ within LAMMPS that store bond topology are updated to reflect the
 post-reacted molecule template. All force fields with fixed bonds,
 angles, dihedrals or impropers are supported.
 
-A few capabilities to note: 1) You may specify as many 'react'
+A few capabilities to note: 1) You may specify as many {react}
 arguments as desired. For example, you could break down a complicated
 reaction mechanism into several reaction steps, each defined by its
-own 'react' argument. 2) While typically a bond is formed or removed
+own {react} argument. 2) While typically a bond is formed or removed
 between the bonding atom pairs specified in the pre-reacted molecule
 template, this is not required. 3) By reversing the order of the pre-
-and post- reacted molecule templates in another 'react' argument, you
+and post- reacted molecule templates in another {react} argument, you
 can allow for the possibility of one or more reverse reactions.
 
 The optional keywords deal with the probability of a given reaction
@@ -304,7 +272,7 @@ you can use the internally-created dynamic group named
 would thermostat the group of all atoms currently involved in a
 reaction:
 
-fix 1 bond_react_MASTER_group temp/rescale 1 300 300 10 1
+fix 1 bond_react_MASTER_group temp/rescale 1 300 300 10 1 :pre
 
 NOTE: This command must be added after the fix bond/react command, and
 will apply to all reactions.
@@ -324,10 +292,11 @@ local command.
 [Restart, fix_modify, output, run start/stop, minimize info:]
 
 No information about this fix is written to "binary restart
-files"_restart.html.  None of the "fix_modify"_fix_modify.html options
-are relevant to this fix.
+files"_restart.html, aside from internally-created per-atom
+properties. None of the "fix_modify"_fix_modify.html options are
+relevant to this fix.
 
-This fix computes one statistic for each 'react' argument that it
+This fix computes one statistic for each {react} argument that it
 stores in a global vector, of length 'number of react arguments', that
 can be accessed by various "output
 commands"_Section_howto.html#howto_15. The vector values calculated by
@@ -359,5 +328,5 @@ The option defaults are stabilization = no, stabilize_steps = 60
 
 :line
 
-:link(Gissinger) [(Gissinger)] Gissinger, Jensen and Wise, Polymer,
-128, 211 (2017).
+:link(Gissinger)
+[(Gissinger)] Gissinger, Jensen and Wise, Polymer, 128, 211 (2017).
diff --git a/doc/src/pair_gw.txt b/doc/src/pair_gw.txt
index fcf63b1bc4e0da31deef339a15da096b245bd4f1..8f1251cf1f3a96aa50e8f9339a69f20b893336bb 100644
--- a/doc/src/pair_gw.txt
+++ b/doc/src/pair_gw.txt
@@ -95,9 +95,9 @@ This pair style can only be used via the {pair} keyword of the
 
 [Restrictions:]
 
-This pair style is part of the USER-MISC package. It is only enabled
-if LAMMPS was built with that package.  See
-the "Making LAMMPS"_Section_start.html#start_3 section for more info.
+This pair style is part of the MANYBODY package. It is only enabled if
+LAMMPS was built with that package.  See the "Making
+LAMMPS"_Section_start.html#start_3 section for more info.
 
 This pair style requires the "newton"_newton.html setting to be "on"
 for pair interactions.
@@ -117,4 +117,5 @@ appropriate units if your simulation doesn't use "metal" units.
 :line
 
 :link(Gao)
-[(Gao)] Gao and Weber, Nuclear Instruments and Methods in Physics Research B 191 (2012) 504.
+[(Gao)] Gao and Weber, Nuclear Instruments and Methods in Physics
+Research B 191 (2012) 504.
diff --git a/doc/src/pair_reaxc.txt b/doc/src/pair_reaxc.txt
index b9dc6e0ed861dccef018578664c70ba510927222..39759b3111825077e3ee0e4dffe54b2eeecd4915 100644
--- a/doc/src/pair_reaxc.txt
+++ b/doc/src/pair_reaxc.txt
@@ -47,13 +47,14 @@ the "(Aktulga)"_#Aktulga paper. The {reax/c} style was initially
 implemented as a stand-alone C code and is now integrated into LAMMPS
 as a package.
 
-The {reax/c/kk} style is a Kokkos version of the ReaxFF potential that is
-derived from the {reax/c} style. The Kokkos version can run on GPUs and
-can also use OpenMP multithreading. For more information about the Kokkos package,
-see "Section 4"_Section_packages.html#kokkos and "Section 5.3.3"_accelerate_kokkos.html.
-One important consideration when using the {reax/c/kk} style is the choice of either
-half or full neighbor lists. This setting can be changed using the Kokkos "package"_package.html
-command.
+The {reax/c/kk} style is a Kokkos version of the ReaxFF potential that
+is derived from the {reax/c} style. The Kokkos version can run on GPUs
+and can also use OpenMP multithreading. For more information about the
+Kokkos package, see "Section 4"_Section_packages.html#kokkos and
+"Section 5.3.3"_accelerate_kokkos.html.  One important consideration
+when using the {reax/c/kk} style is the choice of either half or full
+neighbor lists. This setting can be changed using the Kokkos
+"package"_package.html command.
 
 The {reax/c} style differs from the "pair_style reax"_pair_reax.html
 command in the lo-level implementation details.  The {reax} style is a
@@ -80,9 +81,8 @@ parameterizations for different classes of materials.  You can submit
 a contact request at the Materials Computation Center (MCC) website
 "https://www.mri.psu.edu/materials-computation-center/connect-mcc"_https://www.mri.psu.edu/materials-computation-center/connect-mcc,
 describing the material(s) you are interested in modeling with ReaxFF.
-They can tell
-you what is currently available or what it would take to create a
-suitable ReaxFF parameterization.
+They can tell you what is currently available or what it would take to
+create a suitable ReaxFF parameterization.
 
 The {cfile} setting can be specified as NULL, in which case default
 settings are used. A control file can be specified which defines
@@ -120,28 +120,31 @@ assign to each atom will be used for computing the electrostatic
 interactions in the system.
 See the "fix qeq/reax"_fix_qeq_reax.html command for details.
 
-Using the optional keyword {lgvdw} with the value {yes} turns on
-the low-gradient correction of the ReaxFF/C for long-range
-London Dispersion, as described in the "(Liu)"_#Liu_2011 paper. Force field
+Using the optional keyword {lgvdw} with the value {yes} turns on the
+low-gradient correction of the ReaxFF/C for long-range London
+Dispersion, as described in the "(Liu)"_#Liu_2011 paper. Force field
 file {ffield.reax.lg} is designed for this correction, and is trained
 for several energetic materials (see "Liu"). When using lg-correction,
 recommended value for parameter {thb} is 0.01, which can be set in the
 control file.  Note: Force field files are different for the original
-or lg corrected pair styles, using wrong ffield file generates an error message.
+or lg corrected pair styles, using wrong ffield file generates an
+error message.
 
 Using the optional keyword {enobonds} with the value {yes}, the energy
 of atoms with no bonds (i.e. isolated atoms) is included in the total
 potential energy and the per-atom energy of that atom.  If the value
-{no} is specified then the energy of atoms with no bonds is set to zero.
-The latter behavior is usual not desired, as it causes discontinuities
-in the potential energy when the bonding of an atom drops to zero.
+{no} is specified then the energy of atoms with no bonds is set to
+zero.  The latter behavior is usual not desired, as it causes
+discontinuities in the potential energy when the bonding of an atom
+drops to zero.
 
 Optional keywords {safezone} and {mincap} are used for allocating
-reax/c arrays.  Increasing these values can avoid memory problems, such
-as segmentation faults and bondchk failed errors, that could occur under
-certain conditions. These keywords aren't used by the Kokkos version, which
-instead uses a more robust memory allocation scheme that checks if the sizes of
-the arrays have been exceeded and automatically allocates more memory.
+reax/c arrays.  Increasing these values can avoid memory problems,
+such as segmentation faults and bondchk failed errors, that could
+occur under certain conditions. These keywords aren't used by the
+Kokkos version, which instead uses a more robust memory allocation
+scheme that checks if the sizes of the arrays have been exceeded and
+automatically allocates more memory.
 
 The thermo variable {evdwl} stores the sum of all the ReaxFF potential
 energy contributions, with the exception of the Coulombic and charge
@@ -153,7 +156,8 @@ This pair style tallies a breakdown of the total ReaxFF potential
 energy into sub-categories, which can be accessed via the "compute
 pair"_compute_pair.html command as a vector of values of length 14.
 The 14 values correspond to the following sub-categories (the variable
-names in italics match those used in the original FORTRAN ReaxFF code):
+names in italics match those used in the original FORTRAN ReaxFF
+code):
 
 {eb} = bond energy
 {ea} = atom energy
@@ -340,8 +344,8 @@ reax"_pair_reax.html
 
 [Default:]
 
-The keyword defaults are checkqeq = yes, enobonds = yes, lgvdw = no, safezone = 1.2,
-mincap = 50.
+The keyword defaults are checkqeq = yes, enobonds = yes, lgvdw = no,
+safezone = 1.2, mincap = 50.
 
 :line
 
diff --git a/doc/src/pair_sw.txt b/doc/src/pair_sw.txt
index 6ed8f00236e2e4a4f441401c39823586802818a3..4932fe55d335bb101db35e478f0681ffe416a69d 100644
--- a/doc/src/pair_sw.txt
+++ b/doc/src/pair_sw.txt
@@ -192,8 +192,8 @@ This pair style can only be used via the {pair} keyword of the
 [Restrictions:]
 
 This pair style is part of the MANYBODY package.  It is only enabled
-if LAMMPS was built with that package.  See
-the "Making LAMMPS"_Section_start.html#start_3 section for more info.
+if LAMMPS was built with that package.  See the "Making
+LAMMPS"_Section_start.html#start_3 section for more info.
 
 This pair style requires the "newton"_newton.html setting to be "on"
 for pair interactions.
diff --git a/examples/COUPLE/lammps_quest/lmpqst.cpp b/examples/COUPLE/lammps_quest/lmpqst.cpp
index e7f2d92c1dfeb0d512ab64b24d05ab34f2b42e3a..76a89f1d064b289be0e68b6252dd262a7d8ec549 100644
--- a/examples/COUPLE/lammps_quest/lmpqst.cpp
+++ b/examples/COUPLE/lammps_quest/lmpqst.cpp
@@ -6,10 +6,10 @@
 //         in.lammps = LAMMPS input script
 //         in.quest = Quest input script
 
-#include "mpi.h"
-#include "stdio.h"
-#include "stdlib.h"
-#include "string.h"
+#include <mpi.h>
+#include <cstdio>
+#include <cstdlib>
+#include <cstring>
 #include "stdint.h"
 
 #include "many2one.h"
diff --git a/examples/COUPLE/lammps_spparks/lmpspk.cpp b/examples/COUPLE/lammps_spparks/lmpspk.cpp
index da6d2523696bebe87c759b56a874ca41f6bd4061..2e234c89f1fadb6bbbcdde2a0c579de1a8c97b54 100644
--- a/examples/COUPLE/lammps_spparks/lmpspk.cpp
+++ b/examples/COUPLE/lammps_spparks/lmpspk.cpp
@@ -7,10 +7,10 @@
 //         Sfactor = multiplier on strain effect
 //         in.spparks = SPPARKS input script
 
-#include "mpi.h"
-#include "stdio.h"
-#include "stdlib.h"
-#include "string.h"
+#include <mpi.h>
+#include <cstdio>
+#include <cstdlib>
+#include <cstring>
 
 #include "lammps_data_write.h"
 #include "many2many.h"
diff --git a/examples/COUPLE/library/error.cpp b/examples/COUPLE/library/error.cpp
index f538dce23f103035181620a99f4d0e3a27f31f42..b7a59b8c72344159a77fcd7326f0d606035f7110 100644
--- a/examples/COUPLE/library/error.cpp
+++ b/examples/COUPLE/library/error.cpp
@@ -1,6 +1,6 @@
 #include <mpi.h>
-#include <stdio.h>
-#include <stdlib.h>
+#include <cstdio>
+#include <cstdlib>
 #include "error.h"
 
 /* ---------------------------------------------------------------------- */
diff --git a/examples/COUPLE/library/files.cpp b/examples/COUPLE/library/files.cpp
index 0d98e69a9e0bb80d4594fef1b2c39ce50f45a777..d4d707a6c91fdfcee12e9c446515ba147019bfb9 100644
--- a/examples/COUPLE/library/files.cpp
+++ b/examples/COUPLE/library/files.cpp
@@ -1,5 +1,5 @@
-#include <stdio.h>
-#include <string.h>
+#include <cstdio>
+#include <cstring>
 #include "files.h"
 
 #define MAXLINE 256
diff --git a/examples/COUPLE/library/irregular.cpp b/examples/COUPLE/library/irregular.cpp
index aea9b8332ca4b60ab523c44e7f19d1fcd0308c3a..b822fc859d4dbb25de9cfff9eb4924c3d032e471 100644
--- a/examples/COUPLE/library/irregular.cpp
+++ b/examples/COUPLE/library/irregular.cpp
@@ -1,6 +1,6 @@
-#include "stdio.h"
-#include "stdlib.h"
-#include "string.h"
+#include <cstdio>
+#include <cstdlib>
+#include <cstring>
 #include "irregular.h"
 #include "memory.h"
 #include "error.h"
diff --git a/examples/COUPLE/library/lammps_data_write.cpp b/examples/COUPLE/library/lammps_data_write.cpp
index cb525871e20f054e44dafc94dedca9ef6f18f7d2..9d1039ff4732e79b2a6823fe752a894a74cc3173 100644
--- a/examples/COUPLE/library/lammps_data_write.cpp
+++ b/examples/COUPLE/library/lammps_data_write.cpp
@@ -1,6 +1,6 @@
 #include <mpi.h>
-#include <stdlib.h>
-#include <string.h>
+#include <cstdlib>
+#include <cstring>
 #include "lammps_data_write.h"
 #include "memory.h"
 #include "error.h"
diff --git a/examples/COUPLE/library/lammps_data_write.h b/examples/COUPLE/library/lammps_data_write.h
index 9a6de2b63c3406142c719d00873078a83dc9b8aa..0192c17e2984356c330022b55c7f3124038d1bb1 100644
--- a/examples/COUPLE/library/lammps_data_write.h
+++ b/examples/COUPLE/library/lammps_data_write.h
@@ -1,7 +1,7 @@
 #ifndef LAMMPS_DATA_WRITE_H
 #define LAMMPS_DATA_WRITE_H
 
-#include <stdio.h>
+#include <cstdio>
 #include "send2one.h"
 
 class LAMMPSDataWrite : public Send2One {
diff --git a/examples/COUPLE/library/many2many.cpp b/examples/COUPLE/library/many2many.cpp
index 3ca9ca63bbcb5f36e9193d4eb514d3f05aa580e2..cbf5da8e6ef57dc307940d2894c2e0ef6e9f286f 100644
--- a/examples/COUPLE/library/many2many.cpp
+++ b/examples/COUPLE/library/many2many.cpp
@@ -1,6 +1,6 @@
 #include <mpi.h>
-#include <stdlib.h>
-#include <stdio.h>
+#include <cstdlib>
+#include <cstdio>
 #include "many2many.h"
 #include "irregular.h"
 #include "memory.h"
diff --git a/examples/COUPLE/library/many2one.cpp b/examples/COUPLE/library/many2one.cpp
index 5b5467aa61d325abf2afb324e35320a49b4718ad..ba227bf4719a162559350888ec0a548f53f79b4e 100644
--- a/examples/COUPLE/library/many2one.cpp
+++ b/examples/COUPLE/library/many2one.cpp
@@ -1,6 +1,6 @@
-#include "mpi.h"
-#include "stdio.h"
-#include "stdlib.h"
+#include <mpi.h>
+#include <cstdio>
+#include <cstdlib>
 #include "many2one.h"
 #include "memory.h"
 
diff --git a/examples/COUPLE/library/memory.cpp b/examples/COUPLE/library/memory.cpp
index 9d299a8376e27c7126d9a62d881cd3eae73584e1..98ef8bafc02cbad582fb9b97dcc703a0bcdcb84b 100644
--- a/examples/COUPLE/library/memory.cpp
+++ b/examples/COUPLE/library/memory.cpp
@@ -1,6 +1,6 @@
 #include <mpi.h>
-#include <stdlib.h>
-#include <stdio.h>
+#include <cstdlib>
+#include <cstdio>
 #include "memory.h"
 #include "error.h"
 
diff --git a/examples/COUPLE/library/one2many.cpp b/examples/COUPLE/library/one2many.cpp
index 5c4a1799ba8ac3ee78271bcef27ed921755d5462..aa4a03720a54729c2937abf089642317b1285f34 100644
--- a/examples/COUPLE/library/one2many.cpp
+++ b/examples/COUPLE/library/one2many.cpp
@@ -1,5 +1,5 @@
 #include <mpi.h>
-#include <stdlib.h>
+#include <cstdlib>
 #include "one2many.h"
 #include "memory.h"
 
diff --git a/examples/COUPLE/library/send2one.cpp b/examples/COUPLE/library/send2one.cpp
index d969b4678d3563413d8d0011197c2aed121d66bb..d35cf19a4ecee245b740ae1959ae1c2c0a3b06c8 100644
--- a/examples/COUPLE/library/send2one.cpp
+++ b/examples/COUPLE/library/send2one.cpp
@@ -1,6 +1,6 @@
-#include "mpi.h"
-#include "stdlib.h"
-#include "stdio.h"
+#include <mpi.h>
+#include <cstdlib>
+#include <cstdio>
 #include "send2one.h"
 #include "memory.h"
 #include "error.h"
diff --git a/examples/COUPLE/multiple/multiple.cpp b/examples/COUPLE/multiple/multiple.cpp
index f2de959bac43c37cdc8472f3d6b02e221a1ec6cd..c58c3c78a97624b8eb0ff9237eec4d134549da8d 100644
--- a/examples/COUPLE/multiple/multiple.cpp
+++ b/examples/COUPLE/multiple/multiple.cpp
@@ -23,10 +23,10 @@
 //         Tdelta = incremental temperature for each of N runs
 // See README for compilation instructions
 
-#include "stdio.h"
-#include "stdlib.h"
-#include "string.h"
-#include "mpi.h"
+#include <mpi.h>
+#include <cstdio>
+#include <cstdlib>
+#include <cstring>
 
 #include "lammps.h"         // these are LAMMPS include files
 #include "input.h"
diff --git a/examples/COUPLE/simple/simple.cpp b/examples/COUPLE/simple/simple.cpp
index 912f4b868954e913e92010354bc6b3a296068297..a76ed8a6e9f4785ee19305279a1eb614d33c4621 100644
--- a/examples/COUPLE/simple/simple.cpp
+++ b/examples/COUPLE/simple/simple.cpp
@@ -19,15 +19,16 @@
 //         in.lammps = LAMMPS input script
 // See README for compilation instructions
 
-#include "stdio.h"
-#include "stdlib.h"
-#include "string.h"
-#include "mpi.h"
-
-#include "lammps.h"         // these are LAMMPS include files
-#include "input.h"
-#include "atom.h"
-#include "library.h"
+#include <cstdio>
+#include <cstdlib>
+#include <cstring>
+#include <mpi.h>
+
+// these are LAMMPS include files
+#include <lammps/lammps.h>
+#include <lammps/input.h>
+#include <lammps/atom.h>
+#include <lammps/library.h>
 
 using namespace LAMMPS_NS;
 
diff --git a/lib/atc/Function.cpp b/lib/atc/Function.cpp
index 6004937a018ffe58fcc6f4c0ace6cf916b881f83..12396937bd75a149285d92ae3435964dcc399738 100644
--- a/lib/atc/Function.cpp
+++ b/lib/atc/Function.cpp
@@ -63,11 +63,11 @@ namespace ATC {
     int narg = nargs -1;
 #ifdef _WIN32
     double *dargs = (double *) _alloca(sizeof(double) * narg);
-#endif
+#else
     double *dargs = (double *) alloca(sizeof(double) * narg);
 #endif
     for (int i = 0; i < narg; ++i) dargs[i] = atof(args[i+1]);
-  
+
     return function(type, narg, dargs);
   }
 
diff --git a/lib/kim/Install.py b/lib/kim/Install.py
index 3f1d9fb19141e67be56a1fdeb898ce6853aba2f2..d098250906bfabd19e74db7ad2a657d92a4c4ead 100644
--- a/lib/kim/Install.py
+++ b/lib/kim/Install.py
@@ -21,7 +21,7 @@ Syntax from lib dir: python Install.py -b -v version  -a kim-name
 specify one or more options, order does not matter
 
   -v = version of KIM API library to use
-       default = kim-api-v1.9.4 (current as of Apr 2018)
+       default = kim-api-v1.9.5 (current as of May 2018)
   -b = download and build base KIM API library with example Models
        this will delete any previous installation in the current folder
   -n = do NOT download and build base KIM API library.
@@ -109,7 +109,7 @@ nargs = len(args)
 if nargs == 0: error()
 
 thisdir = os.environ['PWD']
-version = "kim-api-v1.9.4"
+version = "kim-api-v1.9.5"
 
 buildflag = False
 everythingflag = False
diff --git a/lib/kokkos/core/src/impl/Kokkos_Atomic_Fetch_Sub.hpp b/lib/kokkos/core/src/impl/Kokkos_Atomic_Fetch_Sub.hpp
index 28aca0aeed25090ec2c6bb7655c9ead3757b841d..3f58c55396d5ca2f841b066f043ab8d0e6b9e660 100644
--- a/lib/kokkos/core/src/impl/Kokkos_Atomic_Fetch_Sub.hpp
+++ b/lib/kokkos/core/src/impl/Kokkos_Atomic_Fetch_Sub.hpp
@@ -70,6 +70,20 @@ __inline__ __device__
 unsigned int atomic_fetch_sub( volatile unsigned int * const dest , const unsigned int val )
 { return atomicSub((unsigned int*)dest,val); }
 
+__inline__ __device__
+unsigned int atomic_fetch_sub( volatile int64_t * const dest , const int64_t val )
+{ return atomic_fetch_add(dest,-val); }
+
+__inline__ __device__
+unsigned int atomic_fetch_sub( volatile float * const dest , const float val )
+{ return atomicAdd((float*)dest,-val); }
+
+#if ( 600 <= __CUDA_ARCH__ )
+__inline__ __device__
+unsigned int atomic_fetch_sub( volatile double * const dest , const double val )
+{ return atomicAdd((double*)dest,-val); }
+#endif
+
 template < typename T >
 __inline__ __device__
 T atomic_fetch_sub( volatile T * const dest ,
diff --git a/src/DIPOLE/pair_lj_long_dipole_long.cpp b/src/DIPOLE/pair_lj_long_dipole_long.cpp
index 1635bd4f636f33d3b4f678aef45a4d2b4727fe2f..21fc035854d7a2fabeba3b6fdf54eb24055c4aab 100644
--- a/src/DIPOLE/pair_lj_long_dipole_long.cpp
+++ b/src/DIPOLE/pair_lj_long_dipole_long.cpp
@@ -452,7 +452,7 @@ void PairLJLongDipoleLong::compute(int eflag, int vflag)
       ni = sbmask(j);                                   // special index
       j &= NEIGHMASK;
 
-      { register double *xj = x0+(j+(j<<1));
+      { double *xj = x0+(j+(j<<1));
         d[0] = xi[0] - xj[0];                           // pair vector
         d[1] = xi[1] - xj[1];
         d[2] = xi[2] - xj[2]; }
@@ -463,9 +463,9 @@ void PairLJLongDipoleLong::compute(int eflag, int vflag)
       if (order3 && (rsq < cut_coulsq)) {               // dipole
         memcpy(muj, jmu = mu0+(j<<2), sizeof(vector));
         {                                               // series real space
-          register double r = sqrt(rsq);
-          register double x = g_ewald*r;
-          register double f = exp(-x*x)*qqrd2e;
+          double r = sqrt(rsq);
+          double x = g_ewald*r;
+          double f = exp(-x*x)*qqrd2e;
 
           B0 = 1.0/(1.0+EWALD_P*x);                     // eqn 2.8
           B0 *= ((((A5*B0+A4)*B0+A3)*B0+A2)*B0+A1)*f/r;
@@ -524,8 +524,8 @@ void PairLJLongDipoleLong::compute(int eflag, int vflag)
 
       if (rsq < cut_ljsqi[typej]) {                     // lj
         if (order6) {                                   // long-range lj
-          register double rn = r2inv*r2inv*r2inv;
-          register double x2 = g2*rsq, a2 = 1.0/x2;
+          double rn = r2inv*r2inv*r2inv;
+          double x2 = g2*rsq, a2 = 1.0/x2;
           x2 = a2*exp(-x2)*lj4i[typej];
           if (ni < 0) {
             force_lj =
@@ -533,7 +533,7 @@ void PairLJLongDipoleLong::compute(int eflag, int vflag)
             if (eflag) evdwl = rn*lj3i[typej]-g6*((a2+1.0)*a2+0.5)*x2;
           }
           else {                                        // special case
-            register double f = special_lj[ni], t = rn*(1.0-f);
+            double f = special_lj[ni], t = rn*(1.0-f);
             force_lj = f*(rn *= rn)*lj1i[typej]-
               g8*(((6.0*a2+6.0)*a2+3.0)*a2+1.0)*x2*rsq+t*lj2i[typej];
             if (eflag) evdwl =
@@ -541,13 +541,13 @@ void PairLJLongDipoleLong::compute(int eflag, int vflag)
           }
         }
         else {                                          // cut lj
-          register double rn = r2inv*r2inv*r2inv;
+          double rn = r2inv*r2inv*r2inv;
           if (ni < 0) {
             force_lj = rn*(rn*lj1i[typej]-lj2i[typej]);
             if (eflag) evdwl = rn*(rn*lj3i[typej]-lj4i[typej])-offseti[typej];
           }
           else {                                        // special case
-            register double f = special_lj[ni];
+            double f = special_lj[ni];
             force_lj = f*rn*(rn*lj1i[typej]-lj2i[typej]);
             if (eflag) evdwl = f*(
                 rn*(rn*lj3i[typej]-lj4i[typej])-offseti[typej]);
@@ -559,14 +559,14 @@ void PairLJLongDipoleLong::compute(int eflag, int vflag)
 
       fpair = force_coul+force_lj;                      // force
       if (newton_pair || j < nlocal) {
-        register double *fj = f0+(j+(j<<1));
+        double *fj = f0+(j+(j<<1));
         fi[0] += fx = d[0]*fpair+force_d[0]; fj[0] -= fx;
         fi[1] += fy = d[1]*fpair+force_d[1]; fj[1] -= fy;
         fi[2] += fz = d[2]*fpair+force_d[2]; fj[2] -= fz;
         tqi[0] += mui[1]*ti[2]-mui[2]*ti[1];            // torque
         tqi[1] += mui[2]*ti[0]-mui[0]*ti[2];
         tqi[2] += mui[0]*ti[1]-mui[1]*ti[0];
-        register double *tqj = tq0+(j+(j<<1));
+        double *tqj = tq0+(j+(j<<1));
         tqj[0] += muj[1]*tj[2]-muj[2]*tj[1];
         tqj[1] += muj[2]*tj[0]-muj[0]*tj[2];
         tqj[2] += muj[0]*tj[1]-muj[1]*tj[0];
@@ -608,9 +608,9 @@ double PairLJLongDipoleLong::single(int i, int j, int itype, int jtype,
     double G0, G1, G2, B0, B1, B2, B3, mudi, mudj, muij;
     vector d = {xi[0]-xj[0], xi[1]-xj[1], xi[2]-xj[2]};
     {                                                   // series real space
-      register double r = sqrt(rsq);
-      register double x = g_ewald*r;
-      register double f = exp(-x*x)*qqrd2e;
+      double r = sqrt(rsq);
+      double x = g_ewald*r;
+      double f = exp(-x*x)*qqrd2e;
 
       B0 = 1.0/(1.0+EWALD_P*x);                 // eqn 2.8
       B0 *= ((((A5*B0+A4)*B0+A3)*B0+A2)*B0+A1)*f/r;
@@ -644,7 +644,7 @@ double PairLJLongDipoleLong::single(int i, int j, int itype, int jtype,
   if (rsq < cut_ljsq[itype][jtype]) {                   // lennard-jones
     r6inv = r2inv*r2inv*r2inv;
     if (ewald_order&0x40) {                             // long-range
-      register double x2 = g2*rsq, a2 = 1.0/x2, t = r6inv*(1.0-factor_lj);
+      double x2 = g2*rsq, a2 = 1.0/x2, t = r6inv*(1.0-factor_lj);
       x2 = a2*exp(-x2)*lj4[itype][jtype];
       force_lj = factor_lj*(r6inv *= r6inv)*lj1[itype][jtype]-
         g8*(((6.0*a2+6.0)*a2+3.0)*a2+a2)*x2*rsq+t*lj2[itype][jtype];
diff --git a/src/KOKKOS/atom_vec_dpd_kokkos.cpp b/src/KOKKOS/atom_vec_dpd_kokkos.cpp
index 539f8e673698b67b570bad160ff45648c01c2567..24b23c6ccf667f076c5fe6e51157a7fc4a555dc4 100644
--- a/src/KOKKOS/atom_vec_dpd_kokkos.cpp
+++ b/src/KOKKOS/atom_vec_dpd_kokkos.cpp
@@ -158,6 +158,7 @@ void AtomVecDPDKokkos::copy(int i, int j, int delflag)
 {
   sync(Host,X_MASK | V_MASK | TAG_MASK | TYPE_MASK |
             MASK_MASK | IMAGE_MASK | DPDTHETA_MASK |
+            UCG_MASK | UCGNEW_MASK |
             UCOND_MASK | UMECH_MASK | UCHEM_MASK | DVECTOR_MASK);
 
   h_tag[j] = h_tag[i];
@@ -183,6 +184,7 @@ void AtomVecDPDKokkos::copy(int i, int j, int delflag)
 
   modified(Host,X_MASK | V_MASK | TAG_MASK | TYPE_MASK |
                 MASK_MASK | IMAGE_MASK | DPDTHETA_MASK |
+                UCG_MASK | UCGNEW_MASK |
                 UCOND_MASK | UMECH_MASK | UCHEM_MASK | DVECTOR_MASK);
 }
 
@@ -1029,7 +1031,7 @@ int AtomVecDPDKokkos::pack_comm_hybrid(int n, int *list, double *buf)
   int i,j,m;
 
   sync(Host,DPDTHETA_MASK | UCOND_MASK |
-            UMECH_MASK | UCHEM_MASK | UCG_MASK | UCGNEW_MASK);
+            UMECH_MASK | UCHEM_MASK);
 
   m = 0;
   for (i = 0; i < n; i++) {
@@ -1234,7 +1236,7 @@ int AtomVecDPDKokkos::unpack_comm_hybrid(int n, int first, double *buf)
   }
 
   modified(Host,DPDTHETA_MASK | UCOND_MASK |
-                UMECH_MASK | UCHEM_MASK | UCG_MASK | UCGNEW_MASK);
+                UMECH_MASK | UCHEM_MASK );
 
   return m;
 }
@@ -1645,6 +1647,8 @@ int AtomVecDPDKokkos::unpack_restart(double *buf)
   h_uCond[nlocal] = buf[m++];
   h_uMech[nlocal] = buf[m++];
   h_uChem[nlocal] = buf[m++];
+  h_uCG[nlocal] = 0.0;
+  h_uCGnew[nlocal] = 0.0;
 
   double **extra = atom->extra;
   if (atom->nextra_store) {
@@ -1654,6 +1658,7 @@ int AtomVecDPDKokkos::unpack_restart(double *buf)
 
   modified(Host,X_MASK | V_MASK | TAG_MASK | TYPE_MASK |
                 MASK_MASK | IMAGE_MASK | DPDTHETA_MASK |
+                UCG_MASK | UCGNEW_MASK |
                 UCOND_MASK | UMECH_MASK | UCHEM_MASK | DVECTOR_MASK);
 
   atom->nlocal++;
diff --git a/src/KSPACE/ewald_disp.cpp b/src/KSPACE/ewald_disp.cpp
index 8cdaa044aa737917ebfa9c04ff26391cbb2eb3eb..78d7b38e4d424a990b687dca472cb294b318b64a 100644
--- a/src/KSPACE/ewald_disp.cpp
+++ b/src/KSPACE/ewald_disp.cpp
@@ -778,7 +778,7 @@ void EwaldDisp::compute_ek()
         cek->re += zxyz.re*ci[i]; (cek++)->im += zxyz.im*ci[i];
       }
       if (func[3]) {
-        register double muk = mui[0]*h->x+mui[1]*h->y+mui[2]*h->z; ++h;
+        double muk = mui[0]*h->x+mui[1]*h->y+mui[2]*h->z; ++h;
         cek->re += zxyz.re*muk; (cek++)->im += zxyz.im*muk;
       }
     }
@@ -819,7 +819,7 @@ void EwaldDisp::compute_force()
     cek = cek_global;
     memset(sum, 0, EWALD_MAX_NSUMS*sizeof(vector));
     if (func[3]) {
-      register double di = c[3];
+      double di = c[3];
       mui[0] = di*(mu++)[0]; mui[1] = di*(mu++)[0]; mui[2] = di*(mu++)[0];
       mu++;
     }
@@ -830,33 +830,33 @@ void EwaldDisp::compute_force()
       }
       C_CRMULT(zc, z[k->z].z, zxy);
       if (func[0]) {                                        // 1/r
-        register double im = *(ke++)*(zc.im*cek->re+cek->im*zc.re);
+        double im = *(ke++)*(zc.im*cek->re+cek->im*zc.re);
         if (func[3]) cek_coul = cek;
         ++cek;
         sum[0][0] += h->x*im; sum[0][1] += h->y*im; sum[0][2] += h->z*im;
       }
       if (func[1]) {                                        // geometric 1/r^6
-        register double im = *(ke++)*(zc.im*cek->re+cek->im*zc.re); ++cek;
+        double im = *(ke++)*(zc.im*cek->re+cek->im*zc.re); ++cek;
         sum[1][0] += h->x*im; sum[1][1] += h->y*im; sum[1][2] += h->z*im;
       }
       if (func[2]) {                                        // arithmetic 1/r^6
-        register double im, c = *(ke++);
+        double im, c = *(ke++);
         for (i=2; i<9; ++i) {
           im = c*(zc.im*cek->re+cek->im*zc.re); ++cek;
           sum[i][0] += h->x*im; sum[i][1] += h->y*im; sum[i][2] += h->z*im;
         }
       }
       if (func[3]) {                                        // dipole
-        register double im = *(ke)*(zc.im*cek->re+
+        double im = *(ke)*(zc.im*cek->re+
             cek->im*zc.re)*(mui[0]*h->x+mui[1]*h->y+mui[2]*h->z);
-        register double im2 = *(ke)*(zc.re*cek->re-
+        double im2 = *(ke)*(zc.re*cek->re-
             cek->im*zc.im);
         sum[9][0] += h->x*im; sum[9][1] += h->y*im; sum[9][2] += h->z*im;
         t[0] += -mui[1]*h->z*im2 + mui[2]*h->y*im2;        // torque
         t[1] += -mui[2]*h->x*im2 + mui[0]*h->z*im2;
         t[2] += -mui[0]*h->y*im2 + mui[1]*h->x*im2;
         if (func[0]) {                                      // charge-dipole
-          register double qi = *(q)*c[0];
+          double qi = *(q)*c[0];
           im = - *(ke)*(zc.re*cek_coul->re -
               cek_coul->im*zc.im)*(mui[0]*h->x+mui[1]*h->y+mui[2]*h->z);
           im += *(ke)*(zc.re*cek->re - cek->im*zc.im)*qi;
@@ -873,17 +873,17 @@ void EwaldDisp::compute_force()
       }
     }
     if (func[0]) {                                        // 1/r
-      register double qi = *(q++)*c[0];
+      double qi = *(q++)*c[0];
       f[0] -= sum[0][0]*qi; f[1] -= sum[0][1]*qi; f[2] -= sum[0][2]*qi;
     }
     if (func[1]) {                                        // geometric 1/r^6
-      register double bi = B[*type]*c[1];
+      double bi = B[*type]*c[1];
       f[0] -= sum[1][0]*bi; f[1] -= sum[1][1]*bi; f[2] -= sum[1][2]*bi;
     }
     if (func[2]) {                                        // arithmetic 1/r^6
-      register double *bi = B+7*type[0]+7;
+      double *bi = B+7*type[0]+7;
       for (i=2; i<9; ++i) {
-        register double c2 = (--bi)[0]*c[2];
+        double c2 = (--bi)[0]*c[2];
         f[0] -= sum[i][0]*c2; f[1] -= sum[i][1]*c2; f[2] -= sum[i][2]*c2;
       }
     }
@@ -962,7 +962,7 @@ void EwaldDisp::compute_energy()
     if (func[1]) {                                        // geometric 1/r^6
       sum[1] += *(ke++)*(cek->re*cek->re+cek->im*cek->im); ++cek; }
     if (func[2]) {                                        // arithmetic 1/r^6
-      register double r =
+      double r =
             (cek[0].re*cek[6].re+cek[0].im*cek[6].im)+
             (cek[1].re*cek[5].re+cek[1].im*cek[5].im)+
             (cek[2].re*cek[4].re+cek[2].im*cek[4].im)+
@@ -1013,7 +1013,7 @@ void EwaldDisp::compute_energy_peratom()
     cek = cek_global;
     memset(sum, 0, EWALD_MAX_NSUMS*sizeof(double));
     if (func[3]) {
-      register double di = c[3];
+      double di = c[3];
       mui[0] = di*(mu++)[0]; mui[1] = di*(mu++)[0]; mui[2] = di*(mu++)[0];
       mu++;
     }
@@ -1031,7 +1031,7 @@ void EwaldDisp::compute_energy_peratom()
       if (func[1]) {                                        // geometric 1/r^6
         sum[1] += *(ke++)*(cek->re*zc.re - cek->im*zc.im); ++cek; }
       if (func[2]) {                                        // arithmetic 1/r^6
-        register double im, c = *(ke++);
+        double im, c = *(ke++);
         for (i=2; i<9; ++i) {
           im = c*(cek->re*zc.re - cek->im*zc.im); ++cek;
           sum[i] += im;
@@ -1041,7 +1041,7 @@ void EwaldDisp::compute_energy_peratom()
         double muk = (mui[0]*h->x+mui[1]*h->y+mui[2]*h->z);
         sum[9] += *(ke)*(cek->re*zc.re - cek->im*zc.im)*muk;
         if (func[0]) {                                      // charge-dipole
-          register double qj = *(q)*c[0];
+          double qj = *(q)*c[0];
           sum[9] += *(ke)*(cek_coul->im*zc.re + cek_coul->re*zc.im)*muk;
           sum[9] -= *(ke)*(cek->re*zc.im + cek->im*zc.re)*qj;
         }
@@ -1051,17 +1051,17 @@ void EwaldDisp::compute_energy_peratom()
     }
 
     if (func[0]) {                                        // 1/r
-      register double qj = *(q++)*c[0];
+      double qj = *(q++)*c[0];
       *eatomj += sum[0]*qj - energy_self_peratom[j][0];
     }
     if (func[1]) {                                        // geometric 1/r^6
-      register double bj = B[*type]*c[1];
+      double bj = B[*type]*c[1];
       *eatomj += sum[1]*bj - energy_self_peratom[j][1];
     }
     if (func[2]) {                                        // arithmetic 1/r^6
-      register double *bj = B+7*type[0]+7;
+      double *bj = B+7*type[0]+7;
       for (i=2; i<9; ++i) {
-        register double c2 = (--bj)[0]*c[2];
+        double c2 = (--bj)[0]*c[2];
         *eatomj += 0.5*sum[i]*c2;
       }
       *eatomj -= energy_self_peratom[j][2];
@@ -1097,19 +1097,19 @@ void EwaldDisp::compute_virial()
   memset(sum, 0, EWALD_NFUNCS*sizeof(shape));
   for (int k=0; k<nkvec; ++k) {                      // sum over k vectors
     if (func[0]) {                                         // 1/r
-      register double r = cek->re*cek->re+cek->im*cek->im;
+      double r = cek->re*cek->re+cek->im*cek->im;
       if (func[3]) cek_coul = cek;
       ++cek;
       sum[0][0] += *(kv++)*r; sum[0][1] += *(kv++)*r; sum[0][2] += *(kv++)*r;
       sum[0][3] += *(kv++)*r; sum[0][4] += *(kv++)*r; sum[0][5] += *(kv++)*r;
     }
     if (func[1]) {                                        // geometric 1/r^6
-      register double r = cek->re*cek->re+cek->im*cek->im; ++cek;
+      double r = cek->re*cek->re+cek->im*cek->im; ++cek;
       sum[1][0] += *(kv++)*r; sum[1][1] += *(kv++)*r; sum[1][2] += *(kv++)*r;
       sum[1][3] += *(kv++)*r; sum[1][4] += *(kv++)*r; sum[1][5] += *(kv++)*r;
     }
     if (func[2]) {                                        // arithmetic 1/r^6
-      register double r =
+      double r =
             (cek[0].re*cek[6].re+cek[0].im*cek[6].im)+
             (cek[1].re*cek[5].re+cek[1].im*cek[5].im)+
             (cek[2].re*cek[4].re+cek[2].im*cek[4].im)+
@@ -1118,12 +1118,12 @@ void EwaldDisp::compute_virial()
       sum[2][3] += *(kv++)*r; sum[2][4] += *(kv++)*r; sum[2][5] += *(kv++)*r;
     }
     if (func[3]) {
-      register double r = cek->re*cek->re+cek->im*cek->im;
+      double r = cek->re*cek->re+cek->im*cek->im;
       sum[3][0] += *(kv++)*r; sum[3][1] += *(kv++)*r; sum[3][2] += *(kv++)*r;
       sum[3][3] += *(kv++)*r; sum[3][4] += *(kv++)*r; sum[3][5] += *(kv++)*r;
       if (func[0]) {                                      // charge-dipole
         kv -= 6;
-        register double r = 2.0*(cek->re*cek_coul->im - cek->im*cek_coul->re);
+        double r = 2.0*(cek->re*cek_coul->im - cek->im*cek_coul->re);
         sum[3][0] += *(kv++)*r; sum[3][1] += *(kv++)*r; sum[3][2] += *(kv++)*r;
         sum[3][3] += *(kv++)*r; sum[3][4] += *(kv++)*r; sum[3][5] += *(kv++)*r;
       }
@@ -1173,7 +1173,7 @@ void EwaldDisp::compute_virial_dipole()
     cek = cek_global;
     memset(&sum[0], 0, 6*sizeof(double));
     if (func[3]) {
-      register double di = c[3];
+      double di = c[3];
       mui[0] = di*(mu++)[0]; mui[1] = di*(mu++)[0]; mui[2] = di*(mu++)[0];
       mu++;
     }
@@ -1267,7 +1267,7 @@ void EwaldDisp::compute_virial_peratom()
     cek = cek_global;
     memset(sum, 0, EWALD_MAX_NSUMS*sizeof(shape));
     if (func[3]) {
-      register double di = c[3];
+      double di = c[3];
       mui[0] = di*(mu++)[0]; mui[1] = di*(mu++)[0]; mui[2] = di*(mu++)[0];
       mu++;
     }
@@ -1279,7 +1279,7 @@ void EwaldDisp::compute_virial_peratom()
       C_CRMULT(zc, z[k->z].z, zxy);
       if (func[0]) {                                        // 1/r
           if (func[3]) cek_coul = cek;
-          register double r = cek->re*zc.re - cek->im*zc.im; ++cek;
+          double r = cek->re*zc.re - cek->im*zc.im; ++cek;
           sum[0][0] += *(kv++)*r;
           sum[0][1] += *(kv++)*r;
           sum[0][2] += *(kv++)*r;
@@ -1288,7 +1288,7 @@ void EwaldDisp::compute_virial_peratom()
           sum[0][5] += *(kv++)*r;
       }
       if (func[1]) {                                        // geometric 1/r^6
-          register double r = cek->re*zc.re - cek->im*zc.im; ++cek;
+          double r = cek->re*zc.re - cek->im*zc.im; ++cek;
           sum[1][0] += *(kv++)*r;
           sum[1][1] += *(kv++)*r;
           sum[1][2] += *(kv++)*r;
@@ -1297,7 +1297,7 @@ void EwaldDisp::compute_virial_peratom()
           sum[1][5] += *(kv++)*r;
       }
       if (func[2]) {                                        // arithmetic 1/r^6
-        register double r;
+        double r;
         for (i=2; i<9; ++i) {
           r = cek->re*zc.re - cek->im*zc.im; ++cek;
           sum[i][0] += *(kv++)*r;
@@ -1312,7 +1312,7 @@ void EwaldDisp::compute_virial_peratom()
       }
       if (func[3]) {                                        // dipole
          double muk = (mui[0]*h->x+mui[1]*h->y+mui[2]*h->z);
-         register double
+         double
            r = (cek->re*zc.re - cek->im*zc.im)*muk;
          sum[9][0] += *(kv++)*r;
          sum[9][1] += *(kv++)*r;
@@ -1322,7 +1322,7 @@ void EwaldDisp::compute_virial_peratom()
          sum[9][5] += *(kv++)*r;
          if (func[0]) {                                      // charge-dipole
            kv -= 6;
-           register double qj = *(q)*c[0];
+           double qj = *(q)*c[0];
            r = (cek_coul->im*zc.re + cek_coul->re*zc.im)*muk;
            r += -(cek->re*zc.im + cek->im*zc.re)*qj;
            sum[9][0] += *(kv++)*r; sum[9][1] += *(kv++)*r; sum[9][2] += *(kv++)*r;
@@ -1333,17 +1333,17 @@ void EwaldDisp::compute_virial_peratom()
     }
 
     if (func[0]) {                                        // 1/r
-      register double qi = *(q++)*c[0];
+      double qi = *(q++)*c[0];
       for (int n = 0; n < 6; n++) vatomj[n] += sum[0][n]*qi;
     }
     if (func[1]) {                                        // geometric 1/r^6
-      register double bi = B[*type]*c[1];
+      double bi = B[*type]*c[1];
       for (int n = 0; n < 6; n++) vatomj[n] += sum[1][n]*bi;
     }
     if (func[2]) {                                        // arithmetic 1/r^6
-      register double *bj = B+7*type[0]+7;
+      double *bj = B+7*type[0]+7;
       for (i=2; i<9; ++i) {
-        register double c2 = (--bj)[0]*c[2];
+        double c2 = (--bj)[0]*c[2];
         for (int n = 0; n < 6; n++) vatomj[n] += 0.5*sum[i][n]*c2;
       }
     }
diff --git a/src/KSPACE/pair_buck_long_coul_long.cpp b/src/KSPACE/pair_buck_long_coul_long.cpp
index 55825b56571c492eca8c3262a75916e29f82ae60..d5bb7b6c5b6106686e2bd1f641b2a17de7118c05 100644
--- a/src/KSPACE/pair_buck_long_coul_long.cpp
+++ b/src/KSPACE/pair_buck_long_coul_long.cpp
@@ -483,7 +483,7 @@ void PairBuckLongCoulLong::compute(int eflag, int vflag)
       ni = sbmask(j);
       j &= NEIGHMASK;
 
-      { register double *xj = x0+(j+(j<<1));
+      { double *xj = x0+(j+(j<<1));
         d[0] = xi[0] - xj[0];                                // pair vector
         d[1] = xi[1] - xj[1];
         d[2] = xi[2] - xj[2]; }
@@ -494,25 +494,25 @@ void PairBuckLongCoulLong::compute(int eflag, int vflag)
 
       if (order1 && (rsq < cut_coulsq)) {                // coulombic
         if (!ncoultablebits || rsq <= tabinnersq) {        // series real space
-          register double x = g_ewald*r;
-          register double s = qri*q[j], t = 1.0/(1.0+EWALD_P*x);
+          double x = g_ewald*r;
+          double s = qri*q[j], t = 1.0/(1.0+EWALD_P*x);
           if (ni == 0) {
             s *= g_ewald*exp(-x*x);
             force_coul = (t *= ((((t*A5+A4)*t+A3)*t+A2)*t+A1)*s/x)+EWALD_F*s;
             if (eflag) ecoul = t;
           }
           else {                                        // special case
-            register double f = s*(1.0-special_coul[ni])/r;
+            double f = s*(1.0-special_coul[ni])/r;
             s *= g_ewald*exp(-x*x);
             force_coul = (t *= ((((t*A5+A4)*t+A3)*t+A2)*t+A1)*s/x)+EWALD_F*s-f;
             if (eflag) ecoul = t-f;
           }
         }                                                // table real space
         else {
-          register union_int_float_t t;
+          union_int_float_t t;
           t.f = rsq;
-          register const int k = (t.i & ncoulmask) >> ncoulshiftbits;
-          register double f = (rsq-rtable[k])*drtable[k], qiqj = qi*q[j];
+          const int k = (t.i & ncoulmask) >> ncoulshiftbits;
+          double f = (rsq-rtable[k])*drtable[k], qiqj = qi*q[j];
           if (ni == 0) {
             force_coul = qiqj*(ftable[k]+f*dftable[k]);
             if (eflag) ecoul = qiqj*(etable[k]+f*detable[k]);
@@ -527,11 +527,11 @@ void PairBuckLongCoulLong::compute(int eflag, int vflag)
       else force_coul = ecoul = 0.0;
 
       if (rsq < cut_bucksqi[typej]) {                        // buckingham
-        register double rn = r2inv*r2inv*r2inv,
+        double rn = r2inv*r2inv*r2inv,
                         expr = exp(-r*rhoinvi[typej]);
         if (order6) {                                        // long-range
           if (!ndisptablebits || rsq <= tabinnerdispsq) {
-            register double x2 = g2*rsq, a2 = 1.0/x2;
+            double x2 = g2*rsq, a2 = 1.0/x2;
             x2 = a2*exp(-x2)*buckci[typej];
             if (ni == 0) {
               force_buck =
@@ -539,7 +539,7 @@ void PairBuckLongCoulLong::compute(int eflag, int vflag)
               if (eflag) evdwl = expr*buckai[typej]-g6*((a2+1.0)*a2+0.5)*x2;
             }
             else {                                        // special case
-              register double f = special_lj[ni], t = rn*(1.0-f);
+              double f = special_lj[ni], t = rn*(1.0-f);
               force_buck = f*r*expr*buck1i[typej]-
                 g8*(((6.0*a2+6.0)*a2+3.0)*a2+1.0)*x2*rsq+t*buck2i[typej];
               if (eflag) evdwl = f*expr*buckai[typej] -
@@ -547,16 +547,16 @@ void PairBuckLongCoulLong::compute(int eflag, int vflag)
             }
           }
           else {                                              //table real space
-            register union_int_float_t disp_t;
+            union_int_float_t disp_t;
             disp_t.f = rsq;
-            register const int disp_k = (disp_t.i & ndispmask)>>ndispshiftbits;
-            register double f_disp = (rsq-rdisptable[disp_k])*drdisptable[disp_k];
+            const int disp_k = (disp_t.i & ndispmask)>>ndispshiftbits;
+            double f_disp = (rsq-rdisptable[disp_k])*drdisptable[disp_k];
             if (ni == 0) {
               force_buck = r*expr*buck1i[typej]-(fdisptable[disp_k]+f_disp*dfdisptable[disp_k])*buckci[typej];
               if (eflag) evdwl = expr*buckai[typej]-(edisptable[disp_k]+f_disp*dedisptable[disp_k])*buckci[typej];
             }
             else {                                             //speial case
-              register double f = special_lj[ni], t = rn*(1.0-f);
+              double f = special_lj[ni], t = rn*(1.0-f);
               force_buck = f*r*expr*buck1i[typej] -(fdisptable[disp_k]+f_disp*dfdisptable[disp_k])*buckci[typej] +t*buck2i[typej];
               if (eflag) evdwl = f*expr*buckai[typej] -(edisptable[disp_k]+f_disp*dedisptable[disp_k])*buckci[typej]+t*buckci[typej];
             }
@@ -569,7 +569,7 @@ void PairBuckLongCoulLong::compute(int eflag, int vflag)
                          rn*buckci[typej]-offseti[typej];
           }
           else {                                        // special case
-            register double f = special_lj[ni];
+            double f = special_lj[ni];
             force_buck = f*(r*expr*buck1i[typej]-rn*buck2i[typej]);
             if (eflag)
               evdwl = f*(expr*buckai[typej]-rn*buckci[typej]-offseti[typej]);
@@ -581,7 +581,7 @@ void PairBuckLongCoulLong::compute(int eflag, int vflag)
       fpair = (force_coul+force_buck)*r2inv;
 
       if (newton_pair || j < nlocal) {
-        register double *fj = f0+(j+(j<<1)), f;
+        double *fj = f0+(j+(j<<1)), f;
         fi[0] += f = d[0]*fpair; fj[0] -= f;
         fi[1] += f = d[1]*fpair; fj[1] -= f;
         fi[2] += f = d[2]*fpair; fj[2] -= f;
@@ -640,7 +640,7 @@ void PairBuckLongCoulLong::compute_inner()
       ni = sbmask(j);
       j &= NEIGHMASK;
 
-      { register double *xj = x0+(j+(j<<1));
+      { double *xj = x0+(j+(j<<1));
         d[0] = xi[0] - xj[0];                                // pair vector
         d[1] = xi[1] - xj[1];
         d[2] = xi[2] - xj[2]; }
@@ -654,7 +654,7 @@ void PairBuckLongCoulLong::compute_inner()
           qri*q[j]/r : qri*q[j]/r*special_coul[ni];
 
       if (rsq < cut_bucksqi[typej = type[j]]) {                // buckingham
-        register double rn = r2inv*r2inv*r2inv,
+        double rn = r2inv*r2inv*r2inv,
                         expr = exp(-r*rhoinvi[typej]);
         force_buck = ni == 0 ?
           (r*expr*buck1i[typej]-rn*buck2i[typej]) :
@@ -665,12 +665,12 @@ void PairBuckLongCoulLong::compute_inner()
       fpair = (force_coul + force_buck) * r2inv;
 
       if (rsq > cut_out_on_sq) {                        // switching
-        register double rsw = (sqrt(rsq) - cut_out_on)/cut_out_diff;
+        double rsw = (sqrt(rsq) - cut_out_on)/cut_out_diff;
         fpair  *= 1.0 + rsw*rsw*(2.0*rsw-3.0);
       }
 
       if (newton_pair || j < nlocal) {                        // force update
-        register double *fj = f0+(j+(j<<1)), f;
+        double *fj = f0+(j+(j<<1)), f;
         fi[0] += f = d[0]*fpair; fj[0] -= f;
         fi[1] += f = d[1]*fpair; fj[1] -= f;
         fi[2] += f = d[2]*fpair; fj[2] -= f;
@@ -730,7 +730,7 @@ void PairBuckLongCoulLong::compute_middle()
       ni = sbmask(j);
       j &= NEIGHMASK;
 
-      { register double *xj = x0+(j+(j<<1));
+      { double *xj = x0+(j+(j<<1));
         d[0] = xi[0] - xj[0];                                // pair vector
         d[1] = xi[1] - xj[1];
         d[2] = xi[2] - xj[2]; }
@@ -745,7 +745,7 @@ void PairBuckLongCoulLong::compute_middle()
           qri*q[j]/r : qri*q[j]/r*special_coul[ni];
 
       if (rsq < cut_bucksqi[typej = type[j]]) {                // buckingham
-        register double rn = r2inv*r2inv*r2inv,
+        double rn = r2inv*r2inv*r2inv,
                         expr = exp(-r*rhoinvi[typej]);
         force_buck = ni == 0 ?
           (r*expr*buck1i[typej]-rn*buck2i[typej]) :
@@ -756,16 +756,16 @@ void PairBuckLongCoulLong::compute_middle()
       fpair = (force_coul + force_buck) * r2inv;
 
       if (rsq < cut_in_on_sq) {                                // switching
-        register double rsw = (sqrt(rsq) - cut_in_off)/cut_in_diff;
+        double rsw = (sqrt(rsq) - cut_in_off)/cut_in_diff;
         fpair  *= rsw*rsw*(3.0 - 2.0*rsw);
       }
       if (rsq > cut_out_on_sq) {
-        register double rsw = (sqrt(rsq) - cut_out_on)/cut_out_diff;
+        double rsw = (sqrt(rsq) - cut_out_on)/cut_out_diff;
         fpair  *= 1.0 + rsw*rsw*(2.0*rsw-3.0);
       }
 
       if (newton_pair || j < nlocal) {                        // force update
-        register double *fj = f0+(j+(j<<1)), f;
+        double *fj = f0+(j+(j<<1)), f;
         fi[0] += f = d[0]*fpair; fj[0] -= f;
         fi[1] += f = d[1]*fpair; fj[1] -= f;
         fi[2] += f = d[2]*fpair; fj[2] -= f;
@@ -831,7 +831,7 @@ void PairBuckLongCoulLong::compute_outer(int eflag, int vflag)
       ni = sbmask(j);
       j &= NEIGHMASK;
 
-      { register double *xj = x0+(j+(j<<1));
+      { double *xj = x0+(j+(j<<1));
         d[0] = xi[0] - xj[0];                                // pair vector
         d[1] = xi[1] - xj[1];
         d[2] = xi[2] - xj[2]; }
@@ -845,36 +845,36 @@ void PairBuckLongCoulLong::compute_outer(int eflag, int vflag)
       respa_buck = 0.0;
       respa_flag = rsq < cut_in_on_sq ? 1 : 0;
       if (respa_flag && (rsq > cut_in_off_sq)) {
-        register double rsw = (r-cut_in_off)/cut_in_diff;
+        double rsw = (r-cut_in_off)/cut_in_diff;
         frespa = 1-rsw*rsw*(3.0-2.0*rsw);
       }
 
       if (order1 && (rsq < cut_coulsq)) {                // coulombic
         if (!ncoultablebits || rsq <= tabinnersq) {        // series real space
-          register double s = qri*q[j];
+          double s = qri*q[j];
           if (respa_flag)                                // correct for respa
             respa_coul = ni == 0 ? frespa*s/r : frespa*s/r*special_coul[ni];
-          register double x = g_ewald*r, t = 1.0/(1.0+EWALD_P*x);
+          double x = g_ewald*r, t = 1.0/(1.0+EWALD_P*x);
           if (ni == 0) {
             s *= g_ewald*exp(-x*x);
             force_coul = (t *= ((((t*A5+A4)*t+A3)*t+A2)*t+A1)*s/x)+EWALD_F*s-respa_coul;
             if (eflag) ecoul = t;
           }
           else {                                        // correct for special
-            register double ri = s*(1.0-special_coul[ni])/r; s *= g_ewald*exp(-x*x);
+            double ri = s*(1.0-special_coul[ni])/r; s *= g_ewald*exp(-x*x);
             force_coul = (t *= ((((t*A5+A4)*t+A3)*t+A2)*t+A1)*s/x)+EWALD_F*s-ri-respa_coul;
             if (eflag) ecoul = t-ri;
           }
         }                                                // table real space
         else {
           if (respa_flag) {
-            register double s = qri*q[j];
+            double s = qri*q[j];
             respa_coul = ni == 0 ? frespa*s/r : frespa*s/r*special_coul[ni];
           }
-          register union_int_float_t t;
+          union_int_float_t t;
           t.f = rsq;
-          register const int k = (t.i & ncoulmask) >> ncoulshiftbits;
-          register double f = (rsq-rtable[k])*drtable[k], qiqj = qi*q[j];
+          const int k = (t.i & ncoulmask) >> ncoulshiftbits;
+          double f = (rsq-rtable[k])*drtable[k], qiqj = qi*q[j];
           if (ni == 0) {
             force_coul = qiqj*(ftable[k]+f*dftable[k]);
             if (eflag) ecoul = qiqj*(etable[k]+f*detable[k]);
@@ -892,14 +892,14 @@ void PairBuckLongCoulLong::compute_outer(int eflag, int vflag)
       else force_coul = respa_coul = ecoul = 0.0;
 
       if (rsq < cut_bucksqi[typej]) {                        // buckingham
-        register double rn = r2inv*r2inv*r2inv,
+        double rn = r2inv*r2inv*r2inv,
                         expr = exp(-r*rhoinvi[typej]);
         if (respa_flag) respa_buck = ni == 0 ?                 // correct for respa
             frespa*(r*expr*buck1i[typej]-rn*buck2i[typej]) :
             frespa*(r*expr*buck1i[typej]-rn*buck2i[typej])*special_lj[ni];
         if (order6) {                                        // long-range form
           if (!ndisptablebits || rsq <= tabinnerdispsq) {
-            register double x2 = g2*rsq, a2 = 1.0/x2;
+            double x2 = g2*rsq, a2 = 1.0/x2;
             x2 = a2*exp(-x2)*buckci[typej];
             if (ni == 0) {
               force_buck =
@@ -907,7 +907,7 @@ void PairBuckLongCoulLong::compute_outer(int eflag, int vflag)
               if (eflag) evdwl = expr*buckai[typej]-g6*((a2+1.0)*a2+0.5)*x2;
             }
             else {                                        // correct for special
-              register double f = special_lj[ni], t = rn*(1.0-f);
+              double f = special_lj[ni], t = rn*(1.0-f);
               force_buck = f*r*expr*buck1i[typej]-
                 g8*(((6.0*a2+6.0)*a2+3.0)*a2+1.0)*x2*rsq+t*buck2i[typej]-respa_buck;
               if (eflag) evdwl = f*expr*buckai[typej] -
@@ -915,17 +915,17 @@ void PairBuckLongCoulLong::compute_outer(int eflag, int vflag)
             }
           }
           else {          // table real space
-            register union_int_float_t disp_t;
+            union_int_float_t disp_t;
             disp_t.f = rsq;
-            register const int disp_k = (disp_t.i & ndispmask)>>ndispshiftbits;
-            register double f_disp = (rsq-rdisptable[disp_k])*drdisptable[disp_k];
-            register double rn = r2inv*r2inv*r2inv;
+            const int disp_k = (disp_t.i & ndispmask)>>ndispshiftbits;
+            double f_disp = (rsq-rdisptable[disp_k])*drdisptable[disp_k];
+            double rn = r2inv*r2inv*r2inv;
             if (ni == 0) {
               force_buck = r*expr*buck1i[typej]-(fdisptable[disp_k]+f_disp*dfdisptable[disp_k])*buckci[typej]-respa_buck;
               if (eflag) evdwl =  expr*buckai[typej]-(edisptable[disp_k]+f_disp*dedisptable[disp_k])*buckci[typej];
             }
             else {                             //special case
-              register double f = special_lj[ni], t = rn*(1.0-f);
+              double f = special_lj[ni], t = rn*(1.0-f);
               force_buck = f*r*expr*buck1i[typej]-(fdisptable[disp_k]+f_disp*dfdisptable[disp_k])*buckci[typej]+t*buck2i[typej]-respa_buck;
               if (eflag) evdwl = f*expr*buckai[typej]-(edisptable[disp_k]+f_disp*dedisptable[disp_k])*buckci[typej]+t*buckci[typej];
             }
@@ -938,7 +938,7 @@ void PairBuckLongCoulLong::compute_outer(int eflag, int vflag)
               evdwl = expr*buckai[typej]-rn*buckci[typej]-offseti[typej];
           }
           else {                                        // correct for special
-            register double f = special_lj[ni];
+            double f = special_lj[ni];
             force_buck = f*(r*expr*buck1i[typej]-rn*buck2i[typej])-respa_buck;
             if (eflag)
               evdwl = f*(expr*buckai[typej]-rn*buckci[typej]-offseti[typej]);
@@ -950,7 +950,7 @@ void PairBuckLongCoulLong::compute_outer(int eflag, int vflag)
       fpair = (force_coul+force_buck)*r2inv;
 
       if (newton_pair || j < nlocal) {
-        register double *fj = f0+(j+(j<<1)), f;
+        double *fj = f0+(j+(j<<1)), f;
         fi[0] += f = d[0]*fpair; fj[0] -= f;
         fi[1] += f = d[1]*fpair; fj[1] -= f;
         fi[2] += f = d[2]*fpair; fj[2] -= f;
@@ -985,17 +985,17 @@ double PairBuckLongCoulLong::single(int i, int j, int itype, int jtype,
 
   if ((ewald_order&2) && (rsq < cut_coulsq)) {                // coulombic
     if (!ncoultablebits || rsq <= tabinnersq) {                // series real space
-      register double x = g_ewald*r;
-      register double s = force->qqrd2e*q[i]*q[j], t = 1.0/(1.0+EWALD_P*x);
+      double x = g_ewald*r;
+      double s = force->qqrd2e*q[i]*q[j], t = 1.0/(1.0+EWALD_P*x);
       f = s*(1.0-factor_coul)/r; s *= g_ewald*exp(-x*x);
       force_coul = (t *= ((((t*A5+A4)*t+A3)*t+A2)*t+A1)*s/x)+EWALD_F*s-f;
       eng += t-f;
     }
     else {                                                // table real space
-      register union_int_float_t t;
+      union_int_float_t t;
       t.f = rsq;
-      register const int k = (t.i & ncoulmask) >> ncoulshiftbits;
-      register double f = (rsq-rtable[k])*drtable[k], qiqj = q[i]*q[j];
+      const int k = (t.i & ncoulmask) >> ncoulshiftbits;
+      double f = (rsq-rtable[k])*drtable[k], qiqj = q[i]*q[j];
       t.f = (1.0-factor_coul)*(ctable[k]+f*dctable[k]);
       force_coul = qiqj*(ftable[k]+f*dftable[k]-t.f);
       eng += qiqj*(etable[k]+f*detable[k]-t.f);
@@ -1003,10 +1003,10 @@ double PairBuckLongCoulLong::single(int i, int j, int itype, int jtype,
   } else force_coul = 0.0;
 
   if (rsq < cut_bucksq[itype][jtype]) {                        // buckingham
-    register double expr = factor_buck*exp(-sqrt(rsq)*rhoinv[itype][jtype]);
+    double expr = factor_buck*exp(-sqrt(rsq)*rhoinv[itype][jtype]);
     r6inv = r2inv*r2inv*r2inv;
     if (ewald_order&64) {                                // long-range
-      register double x2 = g2*rsq, a2 = 1.0/x2, t = r6inv*(1.0-factor_buck);
+      double x2 = g2*rsq, a2 = 1.0/x2, t = r6inv*(1.0-factor_buck);
       x2 = a2*exp(-x2)*buck_c[itype][jtype];
       force_buck = buck1[itype][jtype]*r*expr-
                g8*(((6.0*a2+6.0)*a2+3.0)*a2+a2)*x2*rsq+t*buck2[itype][jtype];
diff --git a/src/KSPACE/pair_lj_long_coul_long.cpp b/src/KSPACE/pair_lj_long_coul_long.cpp
index 32d6d9e0a3e82e602d838beade4b656496e689b3..2de527ab396e8f0a0a27bcdda23f62707b563181 100644
--- a/src/KSPACE/pair_lj_long_coul_long.cpp
+++ b/src/KSPACE/pair_lj_long_coul_long.cpp
@@ -478,7 +478,7 @@ void PairLJLongCoulLong::compute(int eflag, int vflag)
       ni = sbmask(j);
       j &= NEIGHMASK;
 
-      { register double *xj = x0+(j+(j<<1));
+      { double *xj = x0+(j+(j<<1));
         d[0] = xi[0] - xj[0];                                // pair vector
         d[1] = xi[1] - xj[1];
         d[2] = xi[2] - xj[2]; }
@@ -488,8 +488,8 @@ void PairLJLongCoulLong::compute(int eflag, int vflag)
 
       if (order1 && (rsq < cut_coulsq)) {                // coulombic
         if (!ncoultablebits || rsq <= tabinnersq) {        // series real space
-          register double r = sqrt(rsq), x = g_ewald*r;
-          register double s = qri*q[j], t = 1.0/(1.0+EWALD_P*x);
+          double r = sqrt(rsq), x = g_ewald*r;
+          double s = qri*q[j], t = 1.0/(1.0+EWALD_P*x);
           if (ni == 0) {
             s *= g_ewald*exp(-x*x);
             force_coul = (t *= ((((t*A5+A4)*t+A3)*t+A2)*t+A1)*s/x)+EWALD_F*s;
@@ -502,10 +502,10 @@ void PairLJLongCoulLong::compute(int eflag, int vflag)
           }
         }                                                // table real space
         else {
-          register union_int_float_t t;
+          union_int_float_t t;
           t.f = rsq;
-          register const int k = (t.i & ncoulmask)>>ncoulshiftbits;
-          register double f = (rsq-rtable[k])*drtable[k], qiqj = qi*q[j];
+          const int k = (t.i & ncoulmask)>>ncoulshiftbits;
+          double f = (rsq-rtable[k])*drtable[k], qiqj = qi*q[j];
           if (ni == 0) {
             force_coul = qiqj*(ftable[k]+f*dftable[k]);
             if (eflag) ecoul = qiqj*(etable[k]+f*detable[k]);
@@ -522,8 +522,8 @@ void PairLJLongCoulLong::compute(int eflag, int vflag)
       if (rsq < cut_ljsqi[typej]) {                        // lj
         if (order6) {                                        // long-range lj
           if(!ndisptablebits || rsq <= tabinnerdispsq) {                // series real space
-            register double rn = r2inv*r2inv*r2inv;
-            register double x2 = g2*rsq, a2 = 1.0/x2;
+            double rn = r2inv*r2inv*r2inv;
+            double x2 = g2*rsq, a2 = 1.0/x2;
             x2 = a2*exp(-x2)*lj4i[typej];
             if (ni == 0) {
               force_lj =
@@ -532,7 +532,7 @@ void PairLJLongCoulLong::compute(int eflag, int vflag)
                 evdwl = rn*lj3i[typej]-g6*((a2+1.0)*a2+0.5)*x2;
             }
             else {                                        // special case
-              register double f = special_lj[ni], t = rn*(1.0-f);
+              double f = special_lj[ni], t = rn*(1.0-f);
               force_lj = f*(rn *= rn)*lj1i[typej]-
               g8*(((6.0*a2+6.0)*a2+3.0)*a2+1.0)*x2*rsq+t*lj2i[typej];
               if (eflag)
@@ -540,30 +540,30 @@ void PairLJLongCoulLong::compute(int eflag, int vflag)
             }
           }
           else {                        // table real space
-            register union_int_float_t disp_t;
+            union_int_float_t disp_t;
             disp_t.f = rsq;
-            register const int disp_k = (disp_t.i & ndispmask)>>ndispshiftbits;
-            register double f_disp = (rsq-rdisptable[disp_k])*drdisptable[disp_k];
-            register double rn = r2inv*r2inv*r2inv;
+            const int disp_k = (disp_t.i & ndispmask)>>ndispshiftbits;
+            double f_disp = (rsq-rdisptable[disp_k])*drdisptable[disp_k];
+            double rn = r2inv*r2inv*r2inv;
             if (ni == 0) {
               force_lj = (rn*=rn)*lj1i[typej]-(fdisptable[disp_k]+f_disp*dfdisptable[disp_k])*lj4i[typej];
               if (eflag) evdwl = rn*lj3i[typej]-(edisptable[disp_k]+f_disp*dedisptable[disp_k])*lj4i[typej];
             }
             else {                  // special case
-              register double f = special_lj[ni], t = rn*(1.0-f);
+              double f = special_lj[ni], t = rn*(1.0-f);
               force_lj = f*(rn *= rn)*lj1i[typej]-(fdisptable[disp_k]+f_disp*dfdisptable[disp_k])*lj4i[typej]+t*lj2i[typej];
               if (eflag) evdwl = f*rn*lj3i[typej]-(edisptable[disp_k]+f_disp*dedisptable[disp_k])*lj4i[typej]+t*lj4i[typej];
             }
           }
         }
         else {                                                // cut lj
-          register double rn = r2inv*r2inv*r2inv;
+          double rn = r2inv*r2inv*r2inv;
           if (ni == 0) {
             force_lj = rn*(rn*lj1i[typej]-lj2i[typej]);
             if (eflag) evdwl = rn*(rn*lj3i[typej]-lj4i[typej])-offseti[typej];
           }
           else {                                        // special case
-            register double f = special_lj[ni];
+            double f = special_lj[ni];
             force_lj = f*rn*(rn*lj1i[typej]-lj2i[typej]);
             if (eflag)
               evdwl = f * (rn*(rn*lj3i[typej]-lj4i[typej])-offseti[typej]);
@@ -576,7 +576,7 @@ void PairLJLongCoulLong::compute(int eflag, int vflag)
       fpair = (force_coul+force_lj)*r2inv;
 
       if (newton_pair || j < nlocal) {
-        register double *fj = f0+(j+(j<<1)), f;
+        double *fj = f0+(j+(j<<1)), f;
         fi[0] += f = d[0]*fpair; fj[0] -= f;
         fi[1] += f = d[1]*fpair; fj[1] -= f;
         fi[2] += f = d[2]*fpair; fj[2] -= f;
@@ -635,7 +635,7 @@ void PairLJLongCoulLong::compute_inner()
       ni = sbmask(j);
       j &= NEIGHMASK;
 
-      { register double *xj = x0+(j+(j<<1));
+      { double *xj = x0+(j+(j<<1));
         d[0] = xi[0] - xj[0];                                // pair vector
         d[1] = xi[1] - xj[1];
         d[2] = xi[2] - xj[2]; }
@@ -650,7 +650,7 @@ void PairLJLongCoulLong::compute_inner()
       }
 
       if (rsq < cut_ljsqi[typej = type[j]]) {                // lennard-jones
-        register double rn = r2inv*r2inv*r2inv;
+        double rn = r2inv*r2inv*r2inv;
         force_lj = ni == 0 ?
           rn*(rn*lj1i[typej]-lj2i[typej]) :
           rn*(rn*lj1i[typej]-lj2i[typej])*special_lj[ni];
@@ -660,12 +660,12 @@ void PairLJLongCoulLong::compute_inner()
       fpair = (force_coul + force_lj) * r2inv;
 
       if (rsq > cut_out_on_sq) {                        // switching
-        register double rsw = (sqrt(rsq) - cut_out_on)/cut_out_diff;
+        double rsw = (sqrt(rsq) - cut_out_on)/cut_out_diff;
         fpair  *= 1.0 + rsw*rsw*(2.0*rsw-3.0);
       }
 
       if (newton_pair || j < nlocal) {                        // force update
-        register double *fj = f0+(j+(j<<1)), f;
+        double *fj = f0+(j+(j<<1)), f;
         fi[0] += f = d[0]*fpair; fj[0] -= f;
         fi[1] += f = d[1]*fpair; fj[1] -= f;
         fi[2] += f = d[2]*fpair; fj[2] -= f;
@@ -725,7 +725,7 @@ void PairLJLongCoulLong::compute_middle()
       ni = sbmask(j);
       j &= NEIGHMASK;
 
-      { register double *xj = x0+(j+(j<<1));
+      { double *xj = x0+(j+(j<<1));
         d[0] = xi[0] - xj[0];                                // pair vector
         d[1] = xi[1] - xj[1];
         d[2] = xi[2] - xj[2]; }
@@ -739,7 +739,7 @@ void PairLJLongCoulLong::compute_middle()
           qri*q[j]*sqrt(r2inv) : qri*q[j]*sqrt(r2inv)*special_coul[ni];
 
       if (rsq < cut_ljsqi[typej = type[j]]) {                // lennard-jones
-        register double rn = r2inv*r2inv*r2inv;
+        double rn = r2inv*r2inv*r2inv;
         force_lj = ni == 0 ?
           rn*(rn*lj1i[typej]-lj2i[typej]) :
           rn*(rn*lj1i[typej]-lj2i[typej])*special_lj[ni];
@@ -749,16 +749,16 @@ void PairLJLongCoulLong::compute_middle()
       fpair = (force_coul + force_lj) * r2inv;
 
       if (rsq < cut_in_on_sq) {                                // switching
-        register double rsw = (sqrt(rsq) - cut_in_off)/cut_in_diff;
+        double rsw = (sqrt(rsq) - cut_in_off)/cut_in_diff;
         fpair  *= rsw*rsw*(3.0 - 2.0*rsw);
       }
       if (rsq > cut_out_on_sq) {
-        register double rsw = (sqrt(rsq) - cut_out_on)/cut_out_diff;
+        double rsw = (sqrt(rsq) - cut_out_on)/cut_out_diff;
         fpair  *= 1.0 + rsw*rsw*(2.0*rsw-3.0);
       }
 
       if (newton_pair || j < nlocal) {                        // force update
-        register double *fj = f0+(j+(j<<1)), f;
+        double *fj = f0+(j+(j<<1)), f;
         fi[0] += f = d[0]*fpair; fj[0] -= f;
         fi[1] += f = d[1]*fpair; fj[1] -= f;
         fi[2] += f = d[2]*fpair; fj[2] -= f;
@@ -823,7 +823,7 @@ void PairLJLongCoulLong::compute_outer(int eflag, int vflag)
       ni = sbmask(j);
       j &= NEIGHMASK;
 
-      { register double *xj = x0+(j+(j<<1));
+      { double *xj = x0+(j+(j<<1));
         d[0] = xi[0] - xj[0];                                // pair vector
         d[1] = xi[1] - xj[1];
         d[2] = xi[2] - xj[2]; }
@@ -836,16 +836,16 @@ void PairLJLongCoulLong::compute_outer(int eflag, int vflag)
       respa_lj = 0;
       respa_flag = rsq < cut_in_on_sq ? 1 : 0;
       if (respa_flag && (rsq > cut_in_off_sq)) {
-        register double rsw = (sqrt(rsq)-cut_in_off)/cut_in_diff;
+        double rsw = (sqrt(rsq)-cut_in_off)/cut_in_diff;
         frespa = 1-rsw*rsw*(3.0-2.0*rsw);
       }
 
       if (order1 && (rsq < cut_coulsq)) {                // coulombic
         if (!ncoultablebits || rsq <= tabinnersq) {        // series real space
-          register double r = sqrt(rsq), s = qri*q[j];
+          double r = sqrt(rsq), s = qri*q[j];
           if (respa_flag)                                // correct for respa
             respa_coul = ni == 0 ? frespa*s/r : frespa*s/r*special_coul[ni];
-          register double x = g_ewald*r, t = 1.0/(1.0+EWALD_P*x);
+          double x = g_ewald*r, t = 1.0/(1.0+EWALD_P*x);
           if (ni == 0) {
             s *= g_ewald*exp(-x*x);
             force_coul = (t *= ((((t*A5+A4)*t+A3)*t+A2)*t+A1)*s/x)+EWALD_F*s-respa_coul;
@@ -859,13 +859,13 @@ void PairLJLongCoulLong::compute_outer(int eflag, int vflag)
         }                                                // table real space
         else {
           if (respa_flag) {
-            register double r = sqrt(rsq), s = qri*q[j];
+            double r = sqrt(rsq), s = qri*q[j];
             respa_coul = ni == 0 ? frespa*s/r : frespa*s/r*special_coul[ni];
           }
-          register union_int_float_t t;
+          union_int_float_t t;
           t.f = rsq;
-          register const int k = (t.i & ncoulmask) >> ncoulshiftbits;
-          register double f = (rsq-rtable[k])*drtable[k], qiqj = qi*q[j];
+          const int k = (t.i & ncoulmask) >> ncoulshiftbits;
+          double f = (rsq-rtable[k])*drtable[k], qiqj = qi*q[j];
           if (ni == 0) {
             force_coul = qiqj*(ftable[k]+f*dftable[k]);
             if (eflag) ecoul = qiqj*(etable[k]+f*detable[k]);
@@ -884,13 +884,13 @@ void PairLJLongCoulLong::compute_outer(int eflag, int vflag)
       else force_coul = respa_coul = ecoul = 0.0;
 
       if (rsq < cut_ljsqi[typej]) {                        // lennard-jones
-        register double rn = r2inv*r2inv*r2inv;
+        double rn = r2inv*r2inv*r2inv;
         if (respa_flag) respa_lj = ni == 0 ?                 // correct for respa
             frespa*rn*(rn*lj1i[typej]-lj2i[typej]) :
             frespa*rn*(rn*lj1i[typej]-lj2i[typej])*special_lj[ni];
         if (order6) {                                        // long-range form
           if (!ndisptablebits || rsq <= tabinnerdispsq) {
-            register double x2 = g2*rsq, a2 = 1.0/x2;
+            double x2 = g2*rsq, a2 = 1.0/x2;
             x2 = a2*exp(-x2)*lj4i[typej];
             if (ni == 0) {
               force_lj =
@@ -898,7 +898,7 @@ void PairLJLongCoulLong::compute_outer(int eflag, int vflag)
               if (eflag) evdwl = rn*lj3i[typej]-g6*((a2+1.0)*a2+0.5)*x2;
             }
             else {                                        // correct for special
-              register double f = special_lj[ni], t = rn*(1.0-f);
+              double f = special_lj[ni], t = rn*(1.0-f);
               force_lj = f*(rn *= rn)*lj1i[typej]-
                 g8*(((6.0*a2+6.0)*a2+3.0)*a2+1.0)*x2*rsq+t*lj2i[typej]-respa_lj;
               if (eflag)
@@ -906,17 +906,17 @@ void PairLJLongCoulLong::compute_outer(int eflag, int vflag)
             }
           }
           else {                        // table real space
-            register union_int_float_t disp_t;
+            union_int_float_t disp_t;
             disp_t.f = rsq;
-            register const int disp_k = (disp_t.i & ndispmask)>>ndispshiftbits;
-            register double f_disp = (rsq-rdisptable[disp_k])*drdisptable[disp_k];
-            register double rn = r2inv*r2inv*r2inv;
+            const int disp_k = (disp_t.i & ndispmask)>>ndispshiftbits;
+            double f_disp = (rsq-rdisptable[disp_k])*drdisptable[disp_k];
+            double rn = r2inv*r2inv*r2inv;
             if (ni == 0) {
               force_lj = (rn*=rn)*lj1i[typej]-(fdisptable[disp_k]+f_disp*dfdisptable[disp_k])*lj4i[typej]-respa_lj;
               if (eflag) evdwl = rn*lj3i[typej]-(edisptable[disp_k]+f_disp*dedisptable[disp_k])*lj4i[typej];
             }
             else {                  // special case
-              register double f = special_lj[ni], t = rn*(1.0-f);
+              double f = special_lj[ni], t = rn*(1.0-f);
               force_lj = f*(rn *= rn)*lj1i[typej]-(fdisptable[disp_k]+f_disp*dfdisptable[disp_k])*lj4i[typej]+t*lj2i[typej]-respa_lj;
               if (eflag) evdwl = f*rn*lj3i[typej]-(edisptable[disp_k]+f_disp*dedisptable[disp_k])*lj4i[typej]+t*lj4i[typej];
             }
@@ -928,7 +928,7 @@ void PairLJLongCoulLong::compute_outer(int eflag, int vflag)
             if (eflag) evdwl = rn*(rn*lj3i[typej]-lj4i[typej])-offseti[typej];
           }
           else {                                        // correct for special
-            register double f = special_lj[ni];
+            double f = special_lj[ni];
             force_lj = f*rn*(rn*lj1i[typej]-lj2i[typej])-respa_lj;
             if (eflag)
               evdwl = f*(rn*(rn*lj3i[typej]-lj4i[typej])-offseti[typej]);
@@ -940,7 +940,7 @@ void PairLJLongCoulLong::compute_outer(int eflag, int vflag)
       fpair = (force_coul+force_lj)*r2inv;
 
       if (newton_pair || j < nlocal) {
-        register double *fj = f0+(j+(j<<1)), f;
+        double *fj = f0+(j+(j<<1)), f;
         fi[0] += f = d[0]*fpair; fj[0] -= f;
         fi[1] += f = d[1]*fpair; fj[1] -= f;
         fi[2] += f = d[2]*fpair; fj[2] -= f;
@@ -974,17 +974,17 @@ double PairLJLongCoulLong::single(int i, int j, int itype, int jtype,
   r2inv = 1.0/rsq;
   if ((ewald_order&2) && (rsq < cut_coulsq)) {                // coulombic
     if (!ncoultablebits || rsq <= tabinnersq) {                // series real space
-      register double r = sqrt(rsq), x = g_ewald*r;
-      register double s = force->qqrd2e*q[i]*q[j], t = 1.0/(1.0+EWALD_P*x);
+      double r = sqrt(rsq), x = g_ewald*r;
+      double s = force->qqrd2e*q[i]*q[j], t = 1.0/(1.0+EWALD_P*x);
       r = s*(1.0-factor_coul)/r; s *= g_ewald*exp(-x*x);
       force_coul = (t *= ((((t*A5+A4)*t+A3)*t+A2)*t+A1)*s/x)+EWALD_F*s-r;
       eng += t-r;
     }
     else {                                                // table real space
-      register union_int_float_t t;
+      union_int_float_t t;
       t.f = rsq;
-      register const int k = (t.i & ncoulmask) >> ncoulshiftbits;
-      register double f = (rsq-rtable[k])*drtable[k], qiqj = q[i]*q[j];
+      const int k = (t.i & ncoulmask) >> ncoulshiftbits;
+      double f = (rsq-rtable[k])*drtable[k], qiqj = q[i]*q[j];
       t.f = (1.0-factor_coul)*(ctable[k]+f*dctable[k]);
       force_coul = qiqj*(ftable[k]+f*dftable[k]-t.f);
       eng += qiqj*(etable[k]+f*detable[k]-t.f);
@@ -994,7 +994,7 @@ double PairLJLongCoulLong::single(int i, int j, int itype, int jtype,
   if (rsq < cut_ljsq[itype][jtype]) {                        // lennard-jones
     r6inv = r2inv*r2inv*r2inv;
     if (ewald_order&64) {                                // long-range
-      register double x2 = g2*rsq, a2 = 1.0/x2, t = r6inv*(1.0-factor_lj);
+      double x2 = g2*rsq, a2 = 1.0/x2, t = r6inv*(1.0-factor_lj);
       x2 = a2*exp(-x2)*lj4[itype][jtype];
       force_lj = factor_lj*(r6inv *= r6inv)*lj1[itype][jtype]-
                g8*(((6.0*a2+6.0)*a2+3.0)*a2+a2)*x2*rsq+t*lj2[itype][jtype];
diff --git a/src/KSPACE/pair_lj_long_tip4p_long.cpp b/src/KSPACE/pair_lj_long_tip4p_long.cpp
index 78316b9daad147c3fb7b236760c13a99a08c89dd..a571e30fc377b6bafb48126e149576882aa4d27d 100644
--- a/src/KSPACE/pair_lj_long_tip4p_long.cpp
+++ b/src/KSPACE/pair_lj_long_tip4p_long.cpp
@@ -190,8 +190,8 @@ void PairLJLongTIP4PLong::compute(int eflag, int vflag)
         r2inv = 1.0/rsq;
         if (order6) {                   // long-range lj
           if (!ndisptablebits || rsq <= tabinnerdispsq) {
-            register double rn = r2inv*r2inv*r2inv;
-            register double x2 = g2*rsq, a2 = 1.0/x2;
+            double rn = r2inv*r2inv*r2inv;
+            double x2 = g2*rsq, a2 = 1.0/x2;
             x2 = a2*exp(-x2)*lj4i[jtype];
             if (ni == 0) {
               forcelj =
@@ -200,7 +200,7 @@ void PairLJLongTIP4PLong::compute(int eflag, int vflag)
                 evdwl = rn*lj3i[jtype]-g6*((a2+1.0)*a2+0.5)*x2;
             }
             else {                  // special case
-              register double f = special_lj[ni], t = rn*(1.0-f);
+              double f = special_lj[ni], t = rn*(1.0-f);
               forcelj = f*(rn *= rn)*lj1i[jtype]-
                 g8*(((6.0*a2+6.0)*a2+3.0)*a2+1.0)*x2*rsq+t*lj2i[jtype];
               if (eflag)
@@ -208,30 +208,30 @@ void PairLJLongTIP4PLong::compute(int eflag, int vflag)
             }
           }
           else {                                        // table real space
-            register union_int_float_t disp_t;
+            union_int_float_t disp_t;
             disp_t.f = rsq;
-            register const int disp_k = (disp_t.i & ndispmask)>>ndispshiftbits;
-            register double f_disp = (rsq-rdisptable[disp_k])*drdisptable[disp_k];
-            register double rn = r2inv*r2inv*r2inv;
+            const int disp_k = (disp_t.i & ndispmask)>>ndispshiftbits;
+            double f_disp = (rsq-rdisptable[disp_k])*drdisptable[disp_k];
+            double rn = r2inv*r2inv*r2inv;
             if (ni == 0) {
               forcelj = (rn*=rn)*lj1i[jtype]-(fdisptable[disp_k]+f_disp*dfdisptable[disp_k])*lj4i[jtype];
               if (eflag) evdwl = rn*lj3i[jtype]-(edisptable[disp_k]+f_disp*dedisptable[disp_k])*lj4i[jtype];
             }
             else {                  // special case
-              register double f = special_lj[ni], t = rn*(1.0-f);
+              double f = special_lj[ni], t = rn*(1.0-f);
               forcelj = f*(rn *= rn)*lj1i[jtype]-(fdisptable[disp_k]+f_disp*dfdisptable[disp_k])*lj4i[jtype]+t*lj2i[jtype];
               if (eflag) evdwl = f*rn*lj3i[jtype]-(edisptable[disp_k]+f_disp*dedisptable[disp_k])*lj4i[jtype]+t*lj4i[jtype];
             }
           }
         }
         else {                      // cut lj
-          register double rn = r2inv*r2inv*r2inv;
+          double rn = r2inv*r2inv*r2inv;
           if (ni == 0) {
             forcelj = rn*(rn*lj1i[jtype]-lj2i[jtype]);
             if (eflag) evdwl = rn*(rn*lj3i[jtype]-lj4i[jtype])-offseti[jtype];
           }
           else {                    // special case
-            register double f = special_lj[ni];
+            double f = special_lj[ni];
             forcelj = f*rn*(rn*lj1i[jtype]-lj2i[jtype]);
             if (eflag)
               evdwl = f * (rn*(rn*lj3i[jtype]-lj4i[jtype])-offseti[jtype]);
@@ -573,15 +573,15 @@ void PairLJLongTIP4PLong::compute_inner()
 
       if (rsq < cut_ljsq[itype][jtype] && rsq < cut_out_off_sq ) {  // lj
         r2inv = 1.0/rsq;
-        register double rn = r2inv*r2inv*r2inv;
+        double rn = r2inv*r2inv*r2inv;
         if (ni == 0) forcelj = rn*(rn*lj1i[jtype]-lj2i[jtype]);
         else {                  // special case
-          register double f = special_lj[ni];
+          double f = special_lj[ni];
           forcelj = f*rn*(rn*lj1i[jtype]-lj2i[jtype]);
         }
 
         if (rsq > cut_out_on_sq) {                        // switching
-          register double rsw = (sqrt(rsq) - cut_out_on)/cut_out_diff;
+          double rsw = (sqrt(rsq) - cut_out_on)/cut_out_diff;
           forcelj  *= 1.0 + rsw*rsw*(2.0*rsw-3.0);
         }
 
@@ -642,7 +642,7 @@ void PairLJLongTIP4PLong::compute_inner()
           }
 
           if (rsq > cut_out_on_sq) {                        // switching
-            register double rsw = (sqrt(rsq) - cut_out_on)/cut_out_diff;
+            double rsw = (sqrt(rsq) - cut_out_on)/cut_out_diff;
             forcecoul  *= 1.0 + rsw*rsw*(2.0*rsw-3.0);
           }
 
@@ -826,19 +826,19 @@ void PairLJLongTIP4PLong::compute_middle()
 
       if (rsq < cut_ljsq[itype][jtype] && rsq >= cut_in_off_sq && rsq <= cut_out_off_sq ) {  // lj
         r2inv = 1.0/rsq;
-        register double rn = r2inv*r2inv*r2inv;
+        double rn = r2inv*r2inv*r2inv;
         if (ni == 0) forcelj = rn*(rn*lj1i[jtype]-lj2i[jtype]);
         else {                  // special case
-          register double f = special_lj[ni];
+          double f = special_lj[ni];
           forcelj = f*rn*(rn*lj1i[jtype]-lj2i[jtype]);
         }
 
         if (rsq < cut_in_on_sq) {                                // switching
-          register double rsw = (sqrt(rsq) - cut_in_off)/cut_in_diff;
+          double rsw = (sqrt(rsq) - cut_in_off)/cut_in_diff;
           forcelj  *= rsw*rsw*(3.0 - 2.0*rsw);
         }
         if (rsq > cut_out_on_sq) {
-          register double rsw = (sqrt(rsq) - cut_out_on)/cut_out_diff;
+          double rsw = (sqrt(rsq) - cut_out_on)/cut_out_diff;
           forcelj  *= 1.0 + rsw*rsw*(2.0*rsw-3.0);
         }
 
@@ -899,11 +899,11 @@ void PairLJLongTIP4PLong::compute_middle()
           }
 
           if (rsq < cut_in_on_sq) {                                // switching
-            register double rsw = (sqrt(rsq) - cut_in_off)/cut_in_diff;
+            double rsw = (sqrt(rsq) - cut_in_off)/cut_in_diff;
             forcecoul  *= rsw*rsw*(3.0 - 2.0*rsw);
           }
           if (rsq > cut_out_on_sq) {
-            register double rsw = (sqrt(rsq) - cut_out_on)/cut_out_diff;
+            double rsw = (sqrt(rsq) - cut_out_on)/cut_out_diff;
             forcecoul  *= 1.0 + rsw*rsw*(2.0*rsw-3.0);
           }
 
@@ -1112,18 +1112,18 @@ void PairLJLongTIP4PLong::compute_outer(int eflag, int vflag)
         frespa = 1.0;                                       // check whether and how to compute respa corrections
         respa_flag = rsq < cut_in_on_sq ? 1 : 0;
         if (respa_flag && (rsq > cut_in_off_sq)) {
-          register double rsw = (sqrt(rsq)-cut_in_off)/cut_in_diff;
+          double rsw = (sqrt(rsq)-cut_in_off)/cut_in_diff;
           frespa = 1-rsw*rsw*(3.0-2.0*rsw);
         }
 
         r2inv = 1.0/rsq;
-        register double rn = r2inv*r2inv*r2inv;
+        double rn = r2inv*r2inv*r2inv;
         if (respa_flag) respa_lj = ni == 0 ?                 // correct for respa
                           frespa*rn*(rn*lj1i[jtype]-lj2i[jtype]) :
                           frespa*rn*(rn*lj1i[jtype]-lj2i[jtype])*special_lj[ni];
         if (order6) {                                        // long-range form
           if (!ndisptablebits || rsq <= tabinnerdispsq) {
-            register double x2 = g2*rsq, a2 = 1.0/x2;
+            double x2 = g2*rsq, a2 = 1.0/x2;
             x2 = a2*exp(-x2)*lj4i[jtype];
             if (ni == 0) {
               forcelj =
@@ -1131,7 +1131,7 @@ void PairLJLongTIP4PLong::compute_outer(int eflag, int vflag)
               if (eflag) evdwl = rn*lj3i[jtype]-g6*((a2+1.0)*a2+0.5)*x2;
             }
             else {                                        // correct for special
-              register double f = special_lj[ni], t = rn*(1.0-f);
+              double f = special_lj[ni], t = rn*(1.0-f);
               forcelj = f*(rn *= rn)*lj1i[jtype]-
                 g8*(((6.0*a2+6.0)*a2+3.0)*a2+1.0)*x2*rsq+t*lj2i[jtype]-respa_lj;
               if (eflag)
@@ -1139,16 +1139,16 @@ void PairLJLongTIP4PLong::compute_outer(int eflag, int vflag)
             }
           }
           else {                        // table real space
-            register union_int_float_t disp_t;
+            union_int_float_t disp_t;
             disp_t.f = rsq;
-            register const int disp_k = (disp_t.i & ndispmask)>>ndispshiftbits;
-            register double f_disp = (rsq-rdisptable[disp_k])*drdisptable[disp_k];
+            const int disp_k = (disp_t.i & ndispmask)>>ndispshiftbits;
+            double f_disp = (rsq-rdisptable[disp_k])*drdisptable[disp_k];
             if (ni == 0) {
               forcelj = (rn*=rn)*lj1i[jtype]-(fdisptable[disp_k]+f_disp*dfdisptable[disp_k])*lj4i[jtype]-respa_lj;
               if (eflag) evdwl = rn*lj3i[jtype]-(edisptable[disp_k]+f_disp*dedisptable[disp_k])*lj4i[jtype];
             }
             else {                  // special case
-              register double f = special_lj[ni], t = rn*(1.0-f);
+              double f = special_lj[ni], t = rn*(1.0-f);
               forcelj = f*(rn *= rn)*lj1i[jtype]-(fdisptable[disp_k]+f_disp*dfdisptable[disp_k])*lj4i[jtype]+t*lj2i[jtype]-respa_lj;
               if (eflag) evdwl = f*rn*lj3i[jtype]-(edisptable[disp_k]+f_disp*dedisptable[disp_k])*lj4i[jtype]+t*lj4i[jtype];
             }
@@ -1160,7 +1160,7 @@ void PairLJLongTIP4PLong::compute_outer(int eflag, int vflag)
             if (eflag) evdwl = rn*(rn*lj3i[jtype]-lj4i[jtype])-offseti[jtype];
           }
           else {                                        // correct for special
-            register double f = special_lj[ni];
+            double f = special_lj[ni];
             forcelj = f*rn*(rn*lj1i[jtype]-lj2i[jtype])-respa_lj;
             if (eflag)
               evdwl = f*(rn*(rn*lj3i[jtype]-lj4i[jtype])-offseti[jtype]);
@@ -1225,16 +1225,16 @@ void PairLJLongTIP4PLong::compute_outer(int eflag, int vflag)
           frespa = 1.0;                                       // check whether and how to compute respa corrections
           respa_flag = rsq < cut_in_on_sq ? 1 : 0;
           if (respa_flag && (rsq > cut_in_off_sq)) {
-            register double rsw = (sqrt(rsq)-cut_in_off)/cut_in_diff;
+            double rsw = (sqrt(rsq)-cut_in_off)/cut_in_diff;
             frespa = 1-rsw*rsw*(3.0-2.0*rsw);
           }
 
           r2inv = 1.0 / rsq;
           if (!ncoultablebits || rsq <= tabinnersq) {        // series real space
-            register double r = sqrt(rsq), s = qri*q[j];
+            double r = sqrt(rsq), s = qri*q[j];
             if (respa_flag)                                // correct for respa
               respa_coul = ni == 0 ? frespa*s/r : frespa*s/r*special_coul[ni];
-            register double x = g_ewald*r, t = 1.0/(1.0+EWALD_P*x);
+            double x = g_ewald*r, t = 1.0/(1.0+EWALD_P*x);
             if (ni == 0) {
               s *= g_ewald*exp(-x*x);
               forcecoul = (t *= ((((t*A5+A4)*t+A3)*t+A2)*t+A1)*s/x)+EWALD_F*s-respa_coul;
@@ -1248,13 +1248,13 @@ void PairLJLongTIP4PLong::compute_outer(int eflag, int vflag)
           }                                                // table real space
           else {
             if (respa_flag) {
-              register double r = sqrt(rsq), s = qri*q[j];
+              double r = sqrt(rsq), s = qri*q[j];
               respa_coul = ni == 0 ? frespa*s/r : frespa*s/r*special_coul[ni];
             }
-            register union_int_float_t t;
+            union_int_float_t t;
             t.f = rsq;
-            register const int k = (t.i & ncoulmask) >> ncoulshiftbits;
-            register double f = (t.f-rtable[k])*drtable[k], qiqj = qtmp*q[j];
+            const int k = (t.i & ncoulmask) >> ncoulshiftbits;
+            double f = (t.f-rtable[k])*drtable[k], qiqj = qtmp*q[j];
             if (ni == 0) {
               forcecoul = qiqj*(ftable[k]+f*dftable[k]);
               if (eflag) ecoul = qiqj*(etable[k]+f*detable[k]);
diff --git a/src/OPT/pair_lj_long_coul_long_opt.cpp b/src/OPT/pair_lj_long_coul_long_opt.cpp
index 91ffc2c01ed20e96f317cf82cd48d778cb1233a0..6bd1e8cb9682411d4b3e2b7034beace732b95c2b 100644
--- a/src/OPT/pair_lj_long_coul_long_opt.cpp
+++ b/src/OPT/pair_lj_long_coul_long_opt.cpp
@@ -574,7 +574,7 @@ void PairLJLongCoulLongOpt::eval()
       ni = sbmask(j);
       j &= NEIGHMASK;
 
-      { register double *xj = x0+(j+(j<<1));
+      { double *xj = x0+(j+(j<<1));
         d[0] = xi[0] - xj[0];                                // pair vector
         d[1] = xi[1] - xj[1];
         d[2] = xi[2] - xj[2]; }
@@ -584,8 +584,8 @@ void PairLJLongCoulLongOpt::eval()
 
       if (ORDER1 && (rsq < cut_coulsq)) {                // coulombic
         if (!CTABLE || rsq <= tabinnersq) {        // series real space
-          register double r = sqrt(rsq), x = g_ewald*r;
-          register double s = qri*q[j], t = 1.0/(1.0+EWALD_P*x);
+          double r = sqrt(rsq), x = g_ewald*r;
+          double s = qri*q[j], t = 1.0/(1.0+EWALD_P*x);
           if (ni == 0) {
             s *= g_ewald*exp(-x*x);
             force_coul = (t *= ((((t*A5+A4)*t+A3)*t+A2)*t+A1)*s/x)+EWALD_F*s;
@@ -598,10 +598,10 @@ void PairLJLongCoulLongOpt::eval()
           }
         }                                                // table real space
         else {
-          register union_int_float_t t;
+          union_int_float_t t;
           t.f = rsq;
-          register const int k = (t.i & ncoulmask)>>ncoulshiftbits;
-          register double f = (rsq-rtable[k])*drtable[k], qiqj = qi*q[j];
+          const int k = (t.i & ncoulmask)>>ncoulshiftbits;
+          double f = (rsq-rtable[k])*drtable[k], qiqj = qi*q[j];
           if (ni == 0) {
             force_coul = qiqj*(ftable[k]+f*dftable[k]);
             if (EFLAG) ecoul = qiqj*(etable[k]+f*detable[k]);
@@ -618,8 +618,8 @@ void PairLJLongCoulLongOpt::eval()
       if (rsq < cut_ljsqi[typej]) {                        // lj
         if (ORDER6) {                                        // long-range lj
           if(!LJTABLE || rsq <= tabinnerdispsq) {               // series real space
-            register double rn = r2inv*r2inv*r2inv;
-            register double x2 = g2*rsq, a2 = 1.0/x2;
+            double rn = r2inv*r2inv*r2inv;
+            double x2 = g2*rsq, a2 = 1.0/x2;
             x2 = a2*exp(-x2)*lj4i[typej];
             if (ni == 0) {
               force_lj =
@@ -628,7 +628,7 @@ void PairLJLongCoulLongOpt::eval()
                 evdwl = rn*lj3i[typej]-g6*((a2+1.0)*a2+0.5)*x2;
             }
             else {                                        // special case
-              register double f = special_lj[ni], t = rn*(1.0-f);
+              double f = special_lj[ni], t = rn*(1.0-f);
               force_lj = f*(rn *= rn)*lj1i[typej]-
               g8*(((6.0*a2+6.0)*a2+3.0)*a2+1.0)*x2*rsq+t*lj2i[typej];
               if (EFLAG)
@@ -636,30 +636,30 @@ void PairLJLongCoulLongOpt::eval()
             }
           }
           else {                        // table real space
-            register union_int_float_t disp_t;
+            union_int_float_t disp_t;
             disp_t.f = rsq;
-            register const int disp_k = (disp_t.i & ndispmask)>>ndispshiftbits;
-            register double f_disp = (rsq-rdisptable[disp_k])*drdisptable[disp_k];
-            register double rn = r2inv*r2inv*r2inv;
+            const int disp_k = (disp_t.i & ndispmask)>>ndispshiftbits;
+            double f_disp = (rsq-rdisptable[disp_k])*drdisptable[disp_k];
+            double rn = r2inv*r2inv*r2inv;
             if (ni == 0) {
               force_lj = (rn*=rn)*lj1i[typej]-(fdisptable[disp_k]+f_disp*dfdisptable[disp_k])*lj4i[typej];
               if (EFLAG) evdwl = rn*lj3i[typej]-(edisptable[disp_k]+f_disp*dedisptable[disp_k])*lj4i[typej];
             }
             else {                  // special case
-              register double f = special_lj[ni], t = rn*(1.0-f);
+              double f = special_lj[ni], t = rn*(1.0-f);
               force_lj = f*(rn *= rn)*lj1i[typej]-(fdisptable[disp_k]+f_disp*dfdisptable[disp_k])*lj4i[typej]+t*lj2i[typej];
               if (EFLAG) evdwl = f*rn*lj3i[typej]-(edisptable[disp_k]+f_disp*dedisptable[disp_k])*lj4i[typej]+t*lj4i[typej];
             }
           }
         }
         else {                                                // cut lj
-          register double rn = r2inv*r2inv*r2inv;
+          double rn = r2inv*r2inv*r2inv;
           if (ni == 0) {
             force_lj = rn*(rn*lj1i[typej]-lj2i[typej]);
             if (EFLAG) evdwl = rn*(rn*lj3i[typej]-lj4i[typej])-offseti[typej];
           }
           else {                                        // special case
-            register double f = special_lj[ni];
+            double f = special_lj[ni];
             force_lj = f*rn*(rn*lj1i[typej]-lj2i[typej]);
             if (EFLAG)
               evdwl = f * (rn*(rn*lj3i[typej]-lj4i[typej])-offseti[typej]);
@@ -671,7 +671,7 @@ void PairLJLongCoulLongOpt::eval()
       fpair = (force_coul+force_lj)*r2inv;
 
       if (NEWTON_PAIR || j < nlocal) {
-        register double *fj = f0+(j+(j<<1)), f;
+        double *fj = f0+(j+(j<<1)), f;
         fi[0] += f = d[0]*fpair; fj[0] -= f;
         fi[1] += f = d[1]*fpair; fj[1] -= f;
         fi[2] += f = d[2]*fpair; fj[2] -= f;
@@ -742,7 +742,7 @@ void PairLJLongCoulLongOpt::eval_outer()
       ni = sbmask(j);
       j &= NEIGHMASK;
 
-      { register double *xj = x0+(j+(j<<1));
+      { double *xj = x0+(j+(j<<1));
         d[0] = xi[0] - xj[0];                                // pair vector
         d[1] = xi[1] - xj[1];
         d[2] = xi[2] - xj[2]; }
@@ -755,16 +755,16 @@ void PairLJLongCoulLongOpt::eval_outer()
       respa_lj = 0;
       respa_flag = rsq < cut_in_on_sq ? 1 : 0;
       if (respa_flag && (rsq > cut_in_off_sq)) {
-        register double rsw = (sqrt(rsq)-cut_in_off)/cut_in_diff;
+        double rsw = (sqrt(rsq)-cut_in_off)/cut_in_diff;
         frespa = 1-rsw*rsw*(3.0-2.0*rsw);
       }
 
       if (ORDER1 && (rsq < cut_coulsq)) {                // coulombic
         if (!CTABLE || rsq <= tabinnersq) {        // series real space
-          register double r = sqrt(rsq), s = qri*q[j];
+          double r = sqrt(rsq), s = qri*q[j];
           if (respa_flag)                                // correct for respa
             respa_coul = ni == 0 ? frespa*s/r : frespa*s/r*special_coul[ni];
-          register double x = g_ewald*r, t = 1.0/(1.0+EWALD_P*x);
+          double x = g_ewald*r, t = 1.0/(1.0+EWALD_P*x);
           if (ni == 0) {
             s *= g_ewald*exp(-x*x);
             force_coul = (t *= ((((t*A5+A4)*t+A3)*t+A2)*t+A1)*s/x)+EWALD_F*s-respa_coul;
@@ -778,13 +778,13 @@ void PairLJLongCoulLongOpt::eval_outer()
         }                                                // table real space
         else {
           if (respa_flag) {
-            register double r = sqrt(rsq), s = qri*q[j];
+            double r = sqrt(rsq), s = qri*q[j];
             respa_coul = ni == 0 ? frespa*s/r : frespa*s/r*special_coul[ni];
           }
-          register union_int_float_t t;
+          union_int_float_t t;
           t.f = rsq;
-          register const int k = (t.i & ncoulmask) >> ncoulshiftbits;
-          register double f = (rsq-rtable[k])*drtable[k], qiqj = qi*q[j];
+          const int k = (t.i & ncoulmask) >> ncoulshiftbits;
+          double f = (rsq-rtable[k])*drtable[k], qiqj = qi*q[j];
           if (ni == 0) {
             force_coul = qiqj*(ftable[k]+f*dftable[k]);
             if (EFLAG) ecoul = qiqj*(etable[k]+f*detable[k]);
@@ -803,13 +803,13 @@ void PairLJLongCoulLongOpt::eval_outer()
       else force_coul = respa_coul = ecoul = 0.0;
 
       if (rsq < cut_ljsqi[typej]) {                        // lennard-jones
-        register double rn = r2inv*r2inv*r2inv;
+        double rn = r2inv*r2inv*r2inv;
         if (respa_flag) respa_lj = ni == 0 ?                 // correct for respa
             frespa*rn*(rn*lj1i[typej]-lj2i[typej]) :
             frespa*rn*(rn*lj1i[typej]-lj2i[typej])*special_lj[ni];
         if (ORDER6) {                                        // long-range form
           if (!LJTABLE || rsq <= tabinnerdispsq) {
-            register double x2 = g2*rsq, a2 = 1.0/x2;
+            double x2 = g2*rsq, a2 = 1.0/x2;
             x2 = a2*exp(-x2)*lj4i[typej];
             if (ni == 0) {
               force_lj =
@@ -817,7 +817,7 @@ void PairLJLongCoulLongOpt::eval_outer()
               if (EFLAG) evdwl = rn*lj3i[typej]-g6*((a2+1.0)*a2+0.5)*x2;
             }
             else {                                        // correct for special
-              register double f = special_lj[ni], t = rn*(1.0-f);
+              double f = special_lj[ni], t = rn*(1.0-f);
               force_lj = f*(rn *= rn)*lj1i[typej]-
                 g8*(((6.0*a2+6.0)*a2+3.0)*a2+1.0)*x2*rsq+t*lj2i[typej]-respa_lj;
               if (EFLAG)
@@ -825,17 +825,17 @@ void PairLJLongCoulLongOpt::eval_outer()
             }
           }
           else {                        // table real space
-            register union_int_float_t disp_t;
+            union_int_float_t disp_t;
             disp_t.f = rsq;
-            register const int disp_k = (disp_t.i & ndispmask)>>ndispshiftbits;
-            register double f_disp = (rsq-rdisptable[disp_k])*drdisptable[disp_k];
-            register double rn = r2inv*r2inv*r2inv;
+            const int disp_k = (disp_t.i & ndispmask)>>ndispshiftbits;
+            double f_disp = (rsq-rdisptable[disp_k])*drdisptable[disp_k];
+            double rn = r2inv*r2inv*r2inv;
             if (ni == 0) {
               force_lj = (rn*=rn)*lj1i[typej]-(fdisptable[disp_k]+f_disp*dfdisptable[disp_k])*lj4i[typej]-respa_lj;
               if (EFLAG) evdwl = rn*lj3i[typej]-(edisptable[disp_k]+f_disp*dedisptable[disp_k])*lj4i[typej];
             }
             else {                  // special case
-              register double f = special_lj[ni], t = rn*(1.0-f);
+              double f = special_lj[ni], t = rn*(1.0-f);
               force_lj = f*(rn *= rn)*lj1i[typej]-(fdisptable[disp_k]+f_disp*dfdisptable[disp_k])*lj4i[typej]+t*lj2i[typej]-respa_lj;
               if (EFLAG) evdwl = f*rn*lj3i[typej]-(edisptable[disp_k]+f_disp*dedisptable[disp_k])*lj4i[typej]+t*lj4i[typej];
             }
@@ -847,7 +847,7 @@ void PairLJLongCoulLongOpt::eval_outer()
             if (EFLAG) evdwl = rn*(rn*lj3i[typej]-lj4i[typej])-offseti[typej];
           }
           else {                                        // correct for special
-            register double f = special_lj[ni];
+            double f = special_lj[ni];
             force_lj = f*rn*(rn*lj1i[typej]-lj2i[typej])-respa_lj;
             if (EFLAG)
               evdwl = f*(rn*(rn*lj3i[typej]-lj4i[typej])-offseti[typej]);
@@ -859,7 +859,7 @@ void PairLJLongCoulLongOpt::eval_outer()
       fpair = (force_coul+force_lj)*r2inv;
 
       if (NEWTON_PAIR || j < nlocal) {
-        register double *fj = f0+(j+(j<<1)), f;
+        double *fj = f0+(j+(j<<1)), f;
         fi[0] += f = d[0]*fpair; fj[0] -= f;
         fi[1] += f = d[1]*fpair; fj[1] -= f;
         fi[2] += f = d[2]*fpair; fj[2] -= f;
diff --git a/src/SNAP/pair_snap.cpp b/src/SNAP/pair_snap.cpp
index 2b1816950c28c11c3b23c288aacb4fd33be5348f..eaefb0325db0ad40720492e601b7b9173eecc82f 100644
--- a/src/SNAP/pair_snap.cpp
+++ b/src/SNAP/pair_snap.cpp
@@ -1529,11 +1529,11 @@ void PairSNAP::coeff(int narg, char **arg)
       sna[tid]->grow_rij(nmax);
   }
 
-  if (comm->me == 0)
-    if (ncoeff != sna[0]->ncoeff) {
+  if (ncoeff != sna[0]->ncoeff) {
+    if (comm->me == 0)
       printf("ncoeff = %d snancoeff = %d \n",ncoeff,sna[0]->ncoeff);
-      error->all(FLERR,"Incorrect SNAP parameter file");
-    }
+    error->all(FLERR,"Incorrect SNAP parameter file");
+  }
 
   // Calculate maximum cutoff for all elements
 
diff --git a/src/SNAP/pair_snap.h b/src/SNAP/pair_snap.h
index b60ab3c3e42657fba6fff2f316aa6e93aa5056a5..d39cb0f8d44a73b37ff5c4688b4c3f690e52f232 100644
--- a/src/SNAP/pair_snap.h
+++ b/src/SNAP/pair_snap.h
@@ -37,8 +37,11 @@ public:
   virtual double init_one(int, int);
   virtual double memory_usage();
 
+  double rcutfac, quadraticflag; // declared public to workaround gcc 4.9
+  int ncoeff;                    //  compiler bug, manifest in KOKKOS package
+
 protected:
-  int ncoeff, ncoeffq, ncoeffall;
+  int ncoeffq, ncoeffall;
   double **bvec, ***dbvec;
   class SNA** sna;
   int nmax;
@@ -97,8 +100,8 @@ protected:
   double *wjelem;               // elements weights
   double **coeffelem;           // element bispectrum coefficients
   int *map;                     // mapping from atom types to elements
-  int twojmax, diagonalstyle, switchflag, bzeroflag, quadraticflag;
-  double rcutfac, rfac0, rmin0, wj1, wj2;
+  int twojmax, diagonalstyle, switchflag, bzeroflag;
+  double rfac0, rmin0, wj1, wj2;
   int rcutfacflag, twojmaxflag; // flags for required parameters
 };
 
diff --git a/src/USER-DPD/atom_vec_dpd.cpp b/src/USER-DPD/atom_vec_dpd.cpp
index 6c8e6f91a9b4ea16b778228731be2eb3d2639f94..4cf6356100c0bd4182bbc9f05eb95a0292876830 100644
--- a/src/USER-DPD/atom_vec_dpd.cpp
+++ b/src/USER-DPD/atom_vec_dpd.cpp
@@ -757,6 +757,8 @@ int AtomVecDPD::unpack_restart(double *buf)
   uCond[nlocal] = buf[m++];
   uMech[nlocal] = buf[m++];
   uChem[nlocal] = buf[m++];
+  uCG[nlocal] = 0.0;
+  uCGnew[nlocal] = 0.0;
 
   double **extra = atom->extra;
   if (atom->nextra_store) {
diff --git a/src/USER-MISC/fix_bond_react.cpp b/src/USER-MISC/fix_bond_react.cpp
index 4120ce0df5f8eb715008478af71d77320d440c7c..4d4b536c36e2d53b4049d0fca54d78cf52e7b557 100644
--- a/src/USER-MISC/fix_bond_react.cpp
+++ b/src/USER-MISC/fix_bond_react.cpp
@@ -430,16 +430,6 @@ void FixBondReact::post_constructor()
     delete [] newarg;
   }
 
-  // limit_tags: these are recently reacted atoms being relaxed
-  // per-atom properties already initialized to zero (not in group)
-  // let's do it anyway for clarity
-  int flag;
-  int index = atom->find_custom("limit_tags",flag); //here's where error would happen
-  int *i_limit_tags = atom->ivector[index];
-
-  for (int i = 0; i < atom->nlocal; i++)
-    i_limit_tags[i] = 0;
-
   // create master_group if not already existing
   if (group->find(master_group) == -1) {
     group->find_or_create(master_group);
@@ -456,11 +446,16 @@ void FixBondReact::post_constructor()
 
   // on to statted_tags (system-wide thermostat)
   // intialize per-atom statted_flags to 1
-  index = atom->find_custom("statted_tags",flag);
-  int *i_statted_tags = atom->ivector[index];
-
-  for (int i = 0; i < atom->nlocal; i++)
-    i_statted_tags[i] = 1;
+  // (only if not already initialized by restart)
+  // NOTE: limit_tags and react_tags automaticaly intitialized to zero (unless read from restart)
+  if (fix2->restart_reset != 1) {
+    int flag;
+    int index = atom->find_custom("statted_tags",flag);
+    int *i_statted_tags = atom->ivector[index];
+
+    for (int i = 0; i < atom->nlocal; i++)
+      i_statted_tags[i] = 1;
+  }
 
   if (stabilization_flag == 1) {
     // create exclude_group if not already existing
@@ -497,16 +492,6 @@ void FixBondReact::post_constructor()
 
   }
 
-  //react_tags: this per-atom property is the ID of the 'react' argument which recently caused atom to react
-  //so that atoms which wander between processors may be released to global thermostat at the proper time
-
-  //per-atom values initalized to 0
-  index = atom->find_custom("react_tags",flag);
-  int *i_react_tags = atom->ivector[index];
-
-  for (int i = 0; i < atom->nlocal; i++)
-    i_react_tags[i] = 0;
-
   // currently must redefine dynamic groups so they are updated at proper time
   // -> should double check as to why
 
@@ -553,6 +538,10 @@ void FixBondReact::init()
     if (strcmp(modify->fix[i]->style,"bond/react") == 0) count++;
   if (count > 1 && comm->me == 0) error->warning(FLERR,"More than one fix bond/react");
 
+  if (force->newton_bond == 0 && comm->me == 0) error->warning(FLERR,"Fewer reactions may occur in some cases "
+                                                           "when 'newton off' is set for bonded interactions. "
+                                                           "The reccomended setting is 'newton on'.");
+
   if (strstr(update->integrate_style,"respa"))
     nlevels_respa = ((Respa *) update->integrate)->nlevels;
 
diff --git a/src/USER-MISC/pair_edip.cpp b/src/USER-MISC/pair_edip.cpp
index f556a099c97370ddd350196a613e31f524d9087f..2976d7ad73c994fcb4237e6d05edd642f0f2c749 100644
--- a/src/USER-MISC/pair_edip.cpp
+++ b/src/USER-MISC/pair_edip.cpp
@@ -94,7 +94,7 @@ void PairEDIP::compute(int eflag, int vflag)
   int itype,jtype,ktype,ijparam,ikparam;
   double xtmp,ytmp,ztmp,evdwl;
   int *ilist,*jlist,*numneigh,**firstneigh;
-  register int preForceCoord_counter;
+  int preForceCoord_counter;
 
   double invR_ij;
   double invR_ik;
diff --git a/src/USER-MISC/pair_edip_multi.cpp b/src/USER-MISC/pair_edip_multi.cpp
index 7a4193f55cbebb25944b8266a71444d3b265af0c..6be57eee74c9c32ffb92f28f27391b225e10d6ed 100644
--- a/src/USER-MISC/pair_edip_multi.cpp
+++ b/src/USER-MISC/pair_edip_multi.cpp
@@ -106,7 +106,7 @@ void PairEDIPMulti::compute(int eflag, int vflag)
   int itype,jtype,ktype,ijparam,ikparam,ijkparam;
   double xtmp,ytmp,ztmp,evdwl;
   int *ilist,*jlist,*numneigh,**firstneigh;
-  register int preForceCoord_counter;
+  int preForceCoord_counter;
 
   double zeta_i;
   double dzetair;
diff --git a/src/USER-OMP/pair_buck_long_coul_long_omp.cpp b/src/USER-OMP/pair_buck_long_coul_long_omp.cpp
index 4a67cf9544fda593f82862e976498f1537d3a3c5..acb95ab76cd8c49f63940edb5c2300e9ae910aa1 100644
--- a/src/USER-OMP/pair_buck_long_coul_long_omp.cpp
+++ b/src/USER-OMP/pair_buck_long_coul_long_omp.cpp
@@ -678,7 +678,7 @@ void PairBuckLongCoulLongOMP::eval(int iifrom, int iito, ThrData * const thr)
       ni = sbmask(j);
       j &= NEIGHMASK;
 
-      { register const double *xj = x0+(j+(j<<1));
+      { const double *xj = x0+(j+(j<<1));
         d[0] = xi[0] - xj[0];                                // pair vector
         d[1] = xi[1] - xj[1];
         d[2] = xi[2] - xj[2]; }
@@ -689,25 +689,25 @@ void PairBuckLongCoulLongOMP::eval(int iifrom, int iito, ThrData * const thr)
 
       if (ORDER1 && (rsq < cut_coulsq)) {                // coulombic
         if (!CTABLE || rsq <= tabinnersq) {        // series real space
-          register double x = g_ewald*r;
-          register double s = qri*q[j], t = 1.0/(1.0+EWALD_P*x);
+          double x = g_ewald*r;
+          double s = qri*q[j], t = 1.0/(1.0+EWALD_P*x);
           if (ni == 0) {
             s *= g_ewald*exp(-x*x);
             force_coul = (t *= ((((t*A5+A4)*t+A3)*t+A2)*t+A1)*s/x)+EWALD_F*s;
             if (EFLAG) ecoul = t;
           }
           else {                                        // special case
-            register double f = s*(1.0-special_coul[ni])/r;
+            double f = s*(1.0-special_coul[ni])/r;
             s *= g_ewald*exp(-x*x);
             force_coul = (t *= ((((t*A5+A4)*t+A3)*t+A2)*t+A1)*s/x)+EWALD_F*s-f;
             if (EFLAG) ecoul = t-f;
           }
         }                                                // table real space
         else {
-          register union_int_float_t t;
+          union_int_float_t t;
           t.f = rsq;
-          register const int k = (t.i & ncoulmask) >> ncoulshiftbits;
-          register double f = (rsq-rtable[k])*drtable[k], qiqj = qi*q[j];
+          const int k = (t.i & ncoulmask) >> ncoulshiftbits;
+          double f = (rsq-rtable[k])*drtable[k], qiqj = qi*q[j];
           if (ni == 0) {
             force_coul = qiqj*(ftable[k]+f*dftable[k]);
             if (EFLAG) ecoul = qiqj*(etable[k]+f*detable[k]);
@@ -722,11 +722,11 @@ void PairBuckLongCoulLongOMP::eval(int iifrom, int iito, ThrData * const thr)
       else force_coul = ecoul = 0.0;
 
       if (rsq < cut_bucksqi[typej]) {                        // buckingham
-        register double rn = r2inv*r2inv*r2inv,
+        double rn = r2inv*r2inv*r2inv,
                         expr = exp(-r*rhoinvi[typej]);
         if (ORDER6) {                                        // long-range
           if (!DISPTABLE || rsq <= tabinnerdispsq) {
-            register double x2 = g2*rsq, a2 = 1.0/x2;
+            double x2 = g2*rsq, a2 = 1.0/x2;
             x2 = a2*exp(-x2)*buckci[typej];
             if (ni == 0) {
               force_buck =
@@ -734,7 +734,7 @@ void PairBuckLongCoulLongOMP::eval(int iifrom, int iito, ThrData * const thr)
               if (EFLAG) evdwl = expr*buckai[typej]-g6*((a2+1.0)*a2+0.5)*x2;
             }
             else {                                        // special case
-              register double f = special_lj[ni], t = rn*(1.0-f);
+              double f = special_lj[ni], t = rn*(1.0-f);
               force_buck = f*r*expr*buck1i[typej]-
                 g8*(((6.0*a2+6.0)*a2+3.0)*a2+1.0)*x2*rsq+t*buck2i[typej];
               if (EFLAG) evdwl = f*expr*buckai[typej] -
@@ -742,16 +742,16 @@ void PairBuckLongCoulLongOMP::eval(int iifrom, int iito, ThrData * const thr)
             }
           }
           else {                                              //table real space
-            register union_int_float_t disp_t;
+            union_int_float_t disp_t;
             disp_t.f = rsq;
-            register const int disp_k = (disp_t.i & ndispmask)>>ndispshiftbits;
-            register double f_disp = (rsq-rdisptable[disp_k])*drdisptable[disp_k];
+            const int disp_k = (disp_t.i & ndispmask)>>ndispshiftbits;
+            double f_disp = (rsq-rdisptable[disp_k])*drdisptable[disp_k];
             if (ni == 0) {
               force_buck = r*expr*buck1i[typej]-(fdisptable[disp_k]+f_disp*dfdisptable[disp_k])*buckci[typej];
               if (EFLAG) evdwl = expr*buckai[typej]-(edisptable[disp_k]+f_disp*dedisptable[disp_k])*buckci[typej];
             }
             else {                                             //speial case
-              register double f = special_lj[ni], t = rn*(1.0-f);
+              double f = special_lj[ni], t = rn*(1.0-f);
               force_buck = f*r*expr*buck1i[typej] -(fdisptable[disp_k]+f_disp*dfdisptable[disp_k])*buckci[typej] +t*buck2i[typej];
               if (EFLAG) evdwl = f*expr*buckai[typej] -(edisptable[disp_k]+f_disp*dedisptable[disp_k])*buckci[typej]+t*buckci[typej];
             }
@@ -764,7 +764,7 @@ void PairBuckLongCoulLongOMP::eval(int iifrom, int iito, ThrData * const thr)
                          rn*buckci[typej]-offseti[typej];
           }
           else {                                        // special case
-            register double f = special_lj[ni];
+            double f = special_lj[ni];
             force_buck = f*(r*expr*buck1i[typej]-rn*buck2i[typej]);
             if (EFLAG)
               evdwl = f*(expr*buckai[typej]-rn*buckci[typej]-offseti[typej]);
@@ -776,7 +776,7 @@ void PairBuckLongCoulLongOMP::eval(int iifrom, int iito, ThrData * const thr)
       fpair = (force_coul+force_buck)*r2inv;
 
       if (NEWTON_PAIR || j < nlocal) {
-        register double *fj = f0+(j+(j<<1)), f;
+        double *fj = f0+(j+(j<<1)), f;
         fi[0] += f = d[0]*fpair; fj[0] -= f;
         fi[1] += f = d[1]*fpair; fj[1] -= f;
         fi[2] += f = d[2]*fpair; fj[2] -= f;
@@ -842,7 +842,7 @@ void PairBuckLongCoulLongOMP::eval_inner(int iifrom, int iito, ThrData * const t
       ni = sbmask(j);
       j &= NEIGHMASK;
 
-      { register const double *xj = x0+(j+(j<<1));
+      { const double *xj = x0+(j+(j<<1));
         d[0] = xi[0] - xj[0];                                // pair vector
         d[1] = xi[1] - xj[1];
         d[2] = xi[2] - xj[2]; }
@@ -856,7 +856,7 @@ void PairBuckLongCoulLongOMP::eval_inner(int iifrom, int iito, ThrData * const t
           qri*q[j]/r : qri*q[j]/r*special_coul[ni];
 
       if (rsq < cut_bucksqi[typej = type[j]]) {                // buckingham
-        register double rn = r2inv*r2inv*r2inv,
+        double rn = r2inv*r2inv*r2inv,
                         expr = exp(-r*rhoinvi[typej]);
         force_buck = ni == 0 ?
           (r*expr*buck1i[typej]-rn*buck2i[typej]) :
@@ -867,12 +867,12 @@ void PairBuckLongCoulLongOMP::eval_inner(int iifrom, int iito, ThrData * const t
       fpair = (force_coul + force_buck) * r2inv;
 
       if (rsq > cut_out_on_sq) {                        // switching
-        register double rsw = (sqrt(rsq) - cut_out_on)/cut_out_diff;
+        double rsw = (sqrt(rsq) - cut_out_on)/cut_out_diff;
         fpair  *= 1.0 + rsw*rsw*(2.0*rsw-3.0);
       }
 
       if (newton_pair || j < nlocal) {                        // force update
-        register double *fj = f0+(j+(j<<1)), f;
+        double *fj = f0+(j+(j<<1)), f;
         fi[0] += f = d[0]*fpair; fj[0] -= f;
         fi[1] += f = d[1]*fpair; fj[1] -= f;
         fi[2] += f = d[2]*fpair; fj[2] -= f;
@@ -939,7 +939,7 @@ void PairBuckLongCoulLongOMP::eval_middle(int iifrom, int iito, ThrData * const
       ni = sbmask(j);
       j &= NEIGHMASK;
 
-      { register const double *xj = x0+(j+(j<<1));
+      { const double *xj = x0+(j+(j<<1));
         d[0] = xi[0] - xj[0];                                // pair vector
         d[1] = xi[1] - xj[1];
         d[2] = xi[2] - xj[2]; }
@@ -954,7 +954,7 @@ void PairBuckLongCoulLongOMP::eval_middle(int iifrom, int iito, ThrData * const
           qri*q[j]/r : qri*q[j]/r*special_coul[ni];
 
       if (rsq < cut_bucksqi[typej = type[j]]) {                // buckingham
-        register double rn = r2inv*r2inv*r2inv,
+        double rn = r2inv*r2inv*r2inv,
                         expr = exp(-r*rhoinvi[typej]);
         force_buck = ni == 0 ?
           (r*expr*buck1i[typej]-rn*buck2i[typej]) :
@@ -965,16 +965,16 @@ void PairBuckLongCoulLongOMP::eval_middle(int iifrom, int iito, ThrData * const
       fpair = (force_coul + force_buck) * r2inv;
 
       if (rsq < cut_in_on_sq) {                                // switching
-        register double rsw = (sqrt(rsq) - cut_in_off)/cut_in_diff;
+        double rsw = (sqrt(rsq) - cut_in_off)/cut_in_diff;
         fpair  *= rsw*rsw*(3.0 - 2.0*rsw);
       }
       if (rsq > cut_out_on_sq) {
-        register double rsw = (sqrt(rsq) - cut_out_on)/cut_out_diff;
+        double rsw = (sqrt(rsq) - cut_out_on)/cut_out_diff;
         fpair  *= 1.0 + rsw*rsw*(2.0*rsw-3.0);
       }
 
       if (newton_pair || j < nlocal) {                        // force update
-        register double *fj = f0+(j+(j<<1)), f;
+        double *fj = f0+(j+(j<<1)), f;
         fi[0] += f = d[0]*fpair; fj[0] -= f;
         fi[1] += f = d[1]*fpair; fj[1] -= f;
         fi[2] += f = d[2]*fpair; fj[2] -= f;
@@ -1042,7 +1042,7 @@ void PairBuckLongCoulLongOMP::eval_outer(int iiform, int iito, ThrData * const t
       ni = sbmask(j);
       j &= NEIGHMASK;
 
-      { register const double *xj = x0+(j+(j<<1));
+      { const double *xj = x0+(j+(j<<1));
         d[0] = xi[0] - xj[0];                                // pair vector
         d[1] = xi[1] - xj[1];
         d[2] = xi[2] - xj[2]; }
@@ -1056,36 +1056,36 @@ void PairBuckLongCoulLongOMP::eval_outer(int iiform, int iito, ThrData * const t
       respa_buck = 0.0;
       respa_flag = rsq < cut_in_on_sq ? 1 : 0;
       if (respa_flag && (rsq > cut_in_off_sq)) {
-        register double rsw = (r-cut_in_off)/cut_in_diff;
+        double rsw = (r-cut_in_off)/cut_in_diff;
         frespa = 1-rsw*rsw*(3.0-2.0*rsw);
       }
 
       if (ORDER1 && (rsq < cut_coulsq)) {                // coulombic
         if (!CTABLE || rsq <= tabinnersq) {        // series real space
-          register double s = qri*q[j];
+          double s = qri*q[j];
           if (respa_flag)                                // correct for respa
             respa_coul = ni == 0 ? frespa*s/r : frespa*s/r*special_coul[ni];
-          register double x = g_ewald*r, t = 1.0/(1.0+EWALD_P*x);
+          double x = g_ewald*r, t = 1.0/(1.0+EWALD_P*x);
           if (ni == 0) {
             s *= g_ewald*exp(-x*x);
             force_coul = (t *= ((((t*A5+A4)*t+A3)*t+A2)*t+A1)*s/x)+EWALD_F*s-respa_coul;
             if (EFLAG) ecoul = t;
           }
           else {                                        // correct for special
-            register double ri = s*(1.0-special_coul[ni])/r; s *= g_ewald*exp(-x*x);
+            double ri = s*(1.0-special_coul[ni])/r; s *= g_ewald*exp(-x*x);
             force_coul = (t *= ((((t*A5+A4)*t+A3)*t+A2)*t+A1)*s/x)+EWALD_F*s-ri-respa_coul;
             if (EFLAG) ecoul = t-ri;
           }
         }                                                // table real space
         else {
           if (respa_flag) {
-            register double s = qri*q[j];
+            double s = qri*q[j];
             respa_coul = ni == 0 ? frespa*s/r : frespa*s/r*special_coul[ni];
           }
-          register union_int_float_t t;
+          union_int_float_t t;
           t.f = rsq;
-          register const int k = (t.i & ncoulmask) >> ncoulshiftbits;
-          register double f = (rsq-rtable[k])*drtable[k], qiqj = qi*q[j];
+          const int k = (t.i & ncoulmask) >> ncoulshiftbits;
+          double f = (rsq-rtable[k])*drtable[k], qiqj = qi*q[j];
           if (ni == 0) {
             force_coul = qiqj*(ftable[k]+f*dftable[k]);
             if (EFLAG) ecoul = qiqj*(etable[k]+f*detable[k]);
@@ -1103,14 +1103,14 @@ void PairBuckLongCoulLongOMP::eval_outer(int iiform, int iito, ThrData * const t
       else force_coul = respa_coul = ecoul = 0.0;
 
       if (rsq < cut_bucksqi[typej]) {                        // buckingham
-        register double rn = r2inv*r2inv*r2inv,
+        double rn = r2inv*r2inv*r2inv,
                         expr = exp(-r*rhoinvi[typej]);
         if (respa_flag) respa_buck = ni == 0 ?                 // correct for respa
             frespa*(r*expr*buck1i[typej]-rn*buck2i[typej]) :
             frespa*(r*expr*buck1i[typej]-rn*buck2i[typej])*special_lj[ni];
         if (ORDER6) {                                        // long-range form
           if (!DISPTABLE || rsq <= tabinnerdispsq) {
-            register double x2 = g2*rsq, a2 = 1.0/x2;
+            double x2 = g2*rsq, a2 = 1.0/x2;
             x2 = a2*exp(-x2)*buckci[typej];
             if (ni == 0) {
               force_buck =
@@ -1118,7 +1118,7 @@ void PairBuckLongCoulLongOMP::eval_outer(int iiform, int iito, ThrData * const t
               if (EFLAG) evdwl = expr*buckai[typej]-g6*((a2+1.0)*a2+0.5)*x2;
             }
             else {                                        // correct for special
-              register double f = special_lj[ni], t = rn*(1.0-f);
+              double f = special_lj[ni], t = rn*(1.0-f);
               force_buck = f*r*expr*buck1i[typej]-
                 g8*(((6.0*a2+6.0)*a2+3.0)*a2+1.0)*x2*rsq+t*buck2i[typej]-respa_buck;
               if (EFLAG) evdwl = f*expr*buckai[typej] -
@@ -1126,17 +1126,17 @@ void PairBuckLongCoulLongOMP::eval_outer(int iiform, int iito, ThrData * const t
             }
           }
           else {          // table real space
-            register union_int_float_t disp_t;
+            union_int_float_t disp_t;
             disp_t.f = rsq;
-            register const int disp_k = (disp_t.i & ndispmask)>>ndispshiftbits;
-            register double f_disp = (rsq-rdisptable[disp_k])*drdisptable[disp_k];
-            register double rn = r2inv*r2inv*r2inv;
+            const int disp_k = (disp_t.i & ndispmask)>>ndispshiftbits;
+            double f_disp = (rsq-rdisptable[disp_k])*drdisptable[disp_k];
+            double rn = r2inv*r2inv*r2inv;
             if (ni == 0) {
               force_buck = r*expr*buck1i[typej]-(fdisptable[disp_k]+f_disp*dfdisptable[disp_k])*buckci[typej]-respa_buck;
               if (EFLAG) evdwl =  expr*buckai[typej]-(edisptable[disp_k]+f_disp*dedisptable[disp_k])*buckci[typej];
             }
             else {                             //special case
-              register double f = special_lj[ni], t = rn*(1.0-f);
+              double f = special_lj[ni], t = rn*(1.0-f);
               force_buck = f*r*expr*buck1i[typej]-(fdisptable[disp_k]+f_disp*dfdisptable[disp_k])*buckci[typej]+t*buck2i[typej]-respa_buck;
               if (EFLAG) evdwl = f*expr*buckai[typej]-(edisptable[disp_k]+f_disp*dedisptable[disp_k])*buckci[typej]+t*buckci[typej];
             }
@@ -1149,7 +1149,7 @@ void PairBuckLongCoulLongOMP::eval_outer(int iiform, int iito, ThrData * const t
               evdwl = expr*buckai[typej]-rn*buckci[typej]-offseti[typej];
           }
           else {                                        // correct for special
-            register double f = special_lj[ni];
+            double f = special_lj[ni];
             force_buck = f*(r*expr*buck1i[typej]-rn*buck2i[typej])-respa_buck;
             if (EFLAG)
               evdwl = f*(expr*buckai[typej]-rn*buckci[typej]-offseti[typej]);
@@ -1161,7 +1161,7 @@ void PairBuckLongCoulLongOMP::eval_outer(int iiform, int iito, ThrData * const t
       fpair = (force_coul+force_buck)*r2inv;
 
       if (NEWTON_PAIR || j < nlocal) {
-        register double *fj = f0+(j+(j<<1)), f;
+        double *fj = f0+(j+(j<<1)), f;
         fi[0] += f = d[0]*fpair; fj[0] -= f;
         fi[1] += f = d[1]*fpair; fj[1] -= f;
         fi[2] += f = d[2]*fpair; fj[2] -= f;
diff --git a/src/USER-OMP/pair_edip_omp.cpp b/src/USER-OMP/pair_edip_omp.cpp
index 5ad8c808ba5e39a0ce728f1375065f296b247048..3392160d925f0e504b8a3fb8450f0202f681dd37 100644
--- a/src/USER-OMP/pair_edip_omp.cpp
+++ b/src/USER-OMP/pair_edip_omp.cpp
@@ -84,7 +84,7 @@ void PairEDIPOMP::eval(int iifrom, int iito, ThrData * const thr)
   int itype,jtype,ktype,ijparam,ikparam;
   double xtmp,ytmp,ztmp,evdwl;
   int *ilist,*jlist,*numneigh,**firstneigh;
-  register int preForceCoord_counter;
+  int preForceCoord_counter;
 
   double invR_ij;
   double invR_ik;
diff --git a/src/USER-OMP/pair_lj_long_coul_long_omp.cpp b/src/USER-OMP/pair_lj_long_coul_long_omp.cpp
index 0afc230e92c54c40f1a7578353ee367a80ef82b2..391840285b94e1d54d1d7ab5c98cb16285302954 100644
--- a/src/USER-OMP/pair_lj_long_coul_long_omp.cpp
+++ b/src/USER-OMP/pair_lj_long_coul_long_omp.cpp
@@ -673,7 +673,7 @@ void PairLJLongCoulLongOMP::eval(int iifrom, int iito, ThrData * const thr)
       ni = sbmask(j);
       j &= NEIGHMASK;
 
-      { register const double *xj = x0+(j+(j<<1));
+      { const double *xj = x0+(j+(j<<1));
         d[0] = xi[0] - xj[0];                                // pair vector
         d[1] = xi[1] - xj[1];
         d[2] = xi[2] - xj[2]; }
@@ -683,8 +683,8 @@ void PairLJLongCoulLongOMP::eval(int iifrom, int iito, ThrData * const thr)
 
       if (ORDER1 && (rsq < cut_coulsq)) {                // coulombic
         if (!CTABLE || rsq <= tabinnersq) {        // series real space
-          register double r = sqrt(rsq), x = g_ewald*r;
-          register double s = qri*q[j], t = 1.0/(1.0+EWALD_P*x);
+          double r = sqrt(rsq), x = g_ewald*r;
+          double s = qri*q[j], t = 1.0/(1.0+EWALD_P*x);
           if (ni == 0) {
             s *= g_ewald*exp(-x*x);
             force_coul = (t *= ((((t*A5+A4)*t+A3)*t+A2)*t+A1)*s/x)+EWALD_F*s;
@@ -697,10 +697,10 @@ void PairLJLongCoulLongOMP::eval(int iifrom, int iito, ThrData * const thr)
           }
         }                                                // table real space
         else {
-          register union_int_float_t t;
+          union_int_float_t t;
           t.f = rsq;
-          register const int k = (t.i & ncoulmask)>>ncoulshiftbits;
-          register double f = (rsq-rtable[k])*drtable[k], qiqj = qi*q[j];
+          const int k = (t.i & ncoulmask)>>ncoulshiftbits;
+          double f = (rsq-rtable[k])*drtable[k], qiqj = qi*q[j];
           if (ni == 0) {
             force_coul = qiqj*(ftable[k]+f*dftable[k]);
             if (EFLAG) ecoul = qiqj*(etable[k]+f*detable[k]);
@@ -717,8 +717,8 @@ void PairLJLongCoulLongOMP::eval(int iifrom, int iito, ThrData * const thr)
       if (rsq < cut_ljsqi[typej]) {                        // lj
         if (ORDER6) {                                        // long-range lj
           if(!LJTABLE || rsq <= tabinnerdispsq) {            //series real space
-            register double rn = r2inv*r2inv*r2inv;
-            register double x2 = g2*rsq, a2 = 1.0/x2;
+            double rn = r2inv*r2inv*r2inv;
+            double x2 = g2*rsq, a2 = 1.0/x2;
             x2 = a2*exp(-x2)*lj4i[typej];
             if (ni == 0) {
               force_lj =
@@ -727,7 +727,7 @@ void PairLJLongCoulLongOMP::eval(int iifrom, int iito, ThrData * const thr)
                 evdwl = rn*lj3i[typej]-g6*((a2+1.0)*a2+0.5)*x2;
             }
             else {                                        // special case
-              register double f = special_lj[ni], t = rn*(1.0-f);
+              double f = special_lj[ni], t = rn*(1.0-f);
               force_lj = f*(rn *= rn)*lj1i[typej]-
                 g8*(((6.0*a2+6.0)*a2+3.0)*a2+1.0)*x2*rsq+t*lj2i[typej];
               if (EFLAG)
@@ -735,30 +735,30 @@ void PairLJLongCoulLongOMP::eval(int iifrom, int iito, ThrData * const thr)
             }
           }
           else {                                        // table real space
-            register union_int_float_t disp_t;
+            union_int_float_t disp_t;
             disp_t.f = rsq;
-            register const int disp_k = (disp_t.i & ndispmask)>>ndispshiftbits;
-            register double f_disp = (rsq-rdisptable[disp_k])*drdisptable[disp_k];
-            register double rn = r2inv*r2inv*r2inv;
+            const int disp_k = (disp_t.i & ndispmask)>>ndispshiftbits;
+            double f_disp = (rsq-rdisptable[disp_k])*drdisptable[disp_k];
+            double rn = r2inv*r2inv*r2inv;
             if (ni == 0) {
               force_lj = (rn*=rn)*lj1i[typej]-(fdisptable[disp_k]+f_disp*dfdisptable[disp_k])*lj4i[typej];
               if (EFLAG) evdwl = rn*lj3i[typej]-(edisptable[disp_k]+f_disp*dedisptable[disp_k])*lj4i[typej];
             }
             else {                                      // special case
-              register double f = special_lj[ni], t = rn*(1.0-f);
+              double f = special_lj[ni], t = rn*(1.0-f);
               force_lj = f*(rn *= rn)*lj1i[typej]-(fdisptable[disp_k]+f_disp*dfdisptable[disp_k])*lj4i[typej]+t*lj2i[typej];
               if (EFLAG) evdwl = f*rn*lj3i[typej]-(edisptable[disp_k]+f_disp*dedisptable[disp_k])*lj4i[typej]+t*lj4i[typej];
             }
           }
         }
         else {                                                // cut lj
-          register double rn = r2inv*r2inv*r2inv;
+          double rn = r2inv*r2inv*r2inv;
           if (ni == 0) {
             force_lj = rn*(rn*lj1i[typej]-lj2i[typej]);
             if (EFLAG) evdwl = rn*(rn*lj3i[typej]-lj4i[typej])-offseti[typej];
           }
           else {                                        // special case
-            register double f = special_lj[ni];
+            double f = special_lj[ni];
             force_lj = f*rn*(rn*lj1i[typej]-lj2i[typej]);
             if (EFLAG)
               evdwl = f * (rn*(rn*lj3i[typej]-lj4i[typej])-offseti[typej]);
@@ -770,7 +770,7 @@ void PairLJLongCoulLongOMP::eval(int iifrom, int iito, ThrData * const thr)
       fpair = (force_coul+force_lj)*r2inv;
 
       if (NEWTON_PAIR || j < nlocal) {
-        register double *fj = f0+(j+(j<<1)), f;
+        double *fj = f0+(j+(j<<1)), f;
         fi[0] += f = d[0]*fpair; fj[0] -= f;
         fi[1] += f = d[1]*fpair; fj[1] -= f;
         fi[2] += f = d[2]*fpair; fj[2] -= f;
@@ -834,7 +834,7 @@ void PairLJLongCoulLongOMP::eval_inner(int iifrom, int iito, ThrData * const thr
       ni = sbmask(j);
       j &= NEIGHMASK;
 
-      { register const double *xj = x0+(j+(j<<1));
+      { const double *xj = x0+(j+(j<<1));
         d[0] = xi[0] - xj[0];                                // pair vector
         d[1] = xi[1] - xj[1];
         d[2] = xi[2] - xj[2]; }
@@ -849,7 +849,7 @@ void PairLJLongCoulLongOMP::eval_inner(int iifrom, int iito, ThrData * const thr
       }
 
       if (rsq < cut_ljsqi[typej = type[j]]) {                // lennard-jones
-        register double rn = r2inv*r2inv*r2inv;
+        double rn = r2inv*r2inv*r2inv;
         force_lj = ni == 0 ?
           rn*(rn*lj1i[typej]-lj2i[typej]) :
           rn*(rn*lj1i[typej]-lj2i[typej])*special_lj[ni];
@@ -859,12 +859,12 @@ void PairLJLongCoulLongOMP::eval_inner(int iifrom, int iito, ThrData * const thr
       fpair = (force_coul + force_lj) * r2inv;
 
       if (rsq > cut_out_on_sq) {                        // switching
-        register double rsw = (sqrt(rsq) - cut_out_on)/cut_out_diff;
+        double rsw = (sqrt(rsq) - cut_out_on)/cut_out_diff;
         fpair  *= 1.0 + rsw*rsw*(2.0*rsw-3.0);
       }
 
       if (newton_pair || j < nlocal) {                        // force update
-        register double *fj = f0+(j+(j<<1)), f;
+        double *fj = f0+(j+(j<<1)), f;
         fi[0] += f = d[0]*fpair; fj[0] -= f;
         fi[1] += f = d[1]*fpair; fj[1] -= f;
         fi[2] += f = d[2]*fpair; fj[2] -= f;
@@ -932,7 +932,7 @@ void PairLJLongCoulLongOMP::eval_middle(int iifrom, int iito, ThrData * const th
       ni = sbmask(j);
       j &= NEIGHMASK;
 
-      { register const double *xj = x0+(j+(j<<1));
+      { const double *xj = x0+(j+(j<<1));
         d[0] = xi[0] - xj[0];                                // pair vector
         d[1] = xi[1] - xj[1];
         d[2] = xi[2] - xj[2]; }
@@ -946,7 +946,7 @@ void PairLJLongCoulLongOMP::eval_middle(int iifrom, int iito, ThrData * const th
           qri*q[j]*sqrt(r2inv) : qri*q[j]*sqrt(r2inv)*special_coul[ni];
 
       if (rsq < cut_ljsqi[typej = type[j]]) {                // lennard-jones
-        register double rn = r2inv*r2inv*r2inv;
+        double rn = r2inv*r2inv*r2inv;
         force_lj = ni == 0 ?
           rn*(rn*lj1i[typej]-lj2i[typej]) :
           rn*(rn*lj1i[typej]-lj2i[typej])*special_lj[ni];
@@ -956,16 +956,16 @@ void PairLJLongCoulLongOMP::eval_middle(int iifrom, int iito, ThrData * const th
       fpair = (force_coul + force_lj) * r2inv;
 
       if (rsq < cut_in_on_sq) {                                // switching
-        register double rsw = (sqrt(rsq) - cut_in_off)/cut_in_diff;
+        double rsw = (sqrt(rsq) - cut_in_off)/cut_in_diff;
         fpair  *= rsw*rsw*(3.0 - 2.0*rsw);
       }
       if (rsq > cut_out_on_sq) {
-        register double rsw = (sqrt(rsq) - cut_out_on)/cut_out_diff;
+        double rsw = (sqrt(rsq) - cut_out_on)/cut_out_diff;
         fpair  *= 1.0 + rsw*rsw*(2.0*rsw-3.0);
       }
 
       if (newton_pair || j < nlocal) {                        // force update
-        register double *fj = f0+(j+(j<<1)), f;
+        double *fj = f0+(j+(j<<1)), f;
         fi[0] += f = d[0]*fpair; fj[0] -= f;
         fi[1] += f = d[1]*fpair; fj[1] -= f;
         fi[2] += f = d[2]*fpair; fj[2] -= f;
@@ -1034,7 +1034,7 @@ void PairLJLongCoulLongOMP::eval_outer(int iiform, int iito, ThrData * const thr
       ni = sbmask(j);
       j &= NEIGHMASK;
 
-      { register const double *xj = x0+(j+(j<<1));
+      { const double *xj = x0+(j+(j<<1));
         d[0] = xi[0] - xj[0];                                // pair vector
         d[1] = xi[1] - xj[1];
         d[2] = xi[2] - xj[2]; }
@@ -1047,16 +1047,16 @@ void PairLJLongCoulLongOMP::eval_outer(int iiform, int iito, ThrData * const thr
       respa_lj = 0;
       respa_flag = rsq < cut_in_on_sq ? 1 : 0;
       if (respa_flag && (rsq > cut_in_off_sq)) {
-        register double rsw = (sqrt(rsq)-cut_in_off)/cut_in_diff;
+        double rsw = (sqrt(rsq)-cut_in_off)/cut_in_diff;
         frespa = 1-rsw*rsw*(3.0-2.0*rsw);
       }
 
       if (ORDER1 && (rsq < cut_coulsq)) {                // coulombic
         if (!CTABLE || rsq <= tabinnersq) {        // series real space
-          register double r = sqrt(rsq), s = qri*q[j];
+          double r = sqrt(rsq), s = qri*q[j];
           if (respa_flag)                                // correct for respa
             respa_coul = ni == 0 ? frespa*s/r : frespa*s/r*special_coul[ni];
-          register double x = g_ewald*r, t = 1.0/(1.0+EWALD_P*x);
+          double x = g_ewald*r, t = 1.0/(1.0+EWALD_P*x);
           if (ni == 0) {
             s *= g_ewald*exp(-x*x);
             force_coul = (t *= ((((t*A5+A4)*t+A3)*t+A2)*t+A1)*s/x)+EWALD_F*s-respa_coul;
@@ -1070,13 +1070,13 @@ void PairLJLongCoulLongOMP::eval_outer(int iiform, int iito, ThrData * const thr
         }                                                // table real space
         else {
           if (respa_flag) {
-            register double r = sqrt(rsq), s = qri*q[j];
+            double r = sqrt(rsq), s = qri*q[j];
             respa_coul = ni == 0 ? frespa*s/r : frespa*s/r*special_coul[ni];
           }
-          register union_int_float_t t;
+          union_int_float_t t;
           t.f = rsq;
-          register const int k = (t.i & ncoulmask) >> ncoulshiftbits;
-          register double f = (rsq-rtable[k])*drtable[k], qiqj = qi*q[j];
+          const int k = (t.i & ncoulmask) >> ncoulshiftbits;
+          double f = (rsq-rtable[k])*drtable[k], qiqj = qi*q[j];
           if (ni == 0) {
             force_coul = qiqj*(ftable[k]+f*dftable[k]);
             if (EFLAG) ecoul = qiqj*(etable[k]+f*detable[k]);
@@ -1095,13 +1095,13 @@ void PairLJLongCoulLongOMP::eval_outer(int iiform, int iito, ThrData * const thr
       else force_coul = respa_coul = ecoul = 0.0;
 
       if (rsq < cut_ljsqi[typej]) {                        // lennard-jones
-        register double rn = r2inv*r2inv*r2inv;
+        double rn = r2inv*r2inv*r2inv;
         if (respa_flag) respa_lj = ni == 0 ?                 // correct for respa
             frespa*rn*(rn*lj1i[typej]-lj2i[typej]) :
             frespa*rn*(rn*lj1i[typej]-lj2i[typej])*special_lj[ni];
         if (ORDER6) {                                        // long-range form
           if (!LJTABLE || rsq <= tabinnerdispsq) {
-            register double x2 = g2*rsq, a2 = 1.0/x2;
+            double x2 = g2*rsq, a2 = 1.0/x2;
             x2 = a2*exp(-x2)*lj4i[typej];
             if (ni == 0) {
               force_lj =
@@ -1109,7 +1109,7 @@ void PairLJLongCoulLongOMP::eval_outer(int iiform, int iito, ThrData * const thr
               if (EFLAG) evdwl = rn*lj3i[typej]-g6*((a2+1.0)*a2+0.5)*x2;
             }
             else {                                        // correct for special
-              register double f = special_lj[ni], t = rn*(1.0-f);
+              double f = special_lj[ni], t = rn*(1.0-f);
               force_lj = f*(rn *= rn)*lj1i[typej]-
                 g8*(((6.0*a2+6.0)*a2+3.0)*a2+1.0)*x2*rsq+t*lj2i[typej]-respa_lj;
               if (EFLAG)
@@ -1117,17 +1117,17 @@ void PairLJLongCoulLongOMP::eval_outer(int iiform, int iito, ThrData * const thr
             }
           }
           else {                                                // table real space
-            register union_int_float_t disp_t;
+            union_int_float_t disp_t;
             disp_t.f = rsq;
-            register const int disp_k = (disp_t.i & ndispmask)>>ndispshiftbits;
-            register double f_disp = (rsq-rdisptable[disp_k])*drdisptable[disp_k];
-            register double rn = r2inv*r2inv*r2inv;
+            const int disp_k = (disp_t.i & ndispmask)>>ndispshiftbits;
+            double f_disp = (rsq-rdisptable[disp_k])*drdisptable[disp_k];
+            double rn = r2inv*r2inv*r2inv;
             if (ni == 0) {
               force_lj = (rn*=rn)*lj1i[typej]-(fdisptable[disp_k]+f_disp*dfdisptable[disp_k])*lj4i[typej]-respa_lj;
               if (EFLAG) evdwl = rn*lj3i[typej]-(edisptable[disp_k]+f_disp*dedisptable[disp_k])*lj4i[typej];
             }
             else {                                      // special case
-              register double f = special_lj[ni], t = rn*(1.0-f);
+              double f = special_lj[ni], t = rn*(1.0-f);
               force_lj = f*(rn *= rn)*lj1i[typej]-(fdisptable[disp_k]+f_disp*dfdisptable[disp_k])*lj4i[typej]+t*lj2i[typej]-respa_lj;
               if (EFLAG) evdwl = f*rn*lj3i[typej]-(edisptable[disp_k]+f_disp*dedisptable[disp_k])*lj4i[typej]+t*lj4i[typej];
             }
@@ -1139,7 +1139,7 @@ void PairLJLongCoulLongOMP::eval_outer(int iiform, int iito, ThrData * const thr
             if (EFLAG) evdwl = rn*(rn*lj3i[typej]-lj4i[typej])-offseti[typej];
           }
           else {                                        // correct for special
-            register double f = special_lj[ni];
+            double f = special_lj[ni];
             force_lj = f*rn*(rn*lj1i[typej]-lj2i[typej])-respa_lj;
             if (EFLAG)
               evdwl = f*(rn*(rn*lj3i[typej]-lj4i[typej])-offseti[typej]);
@@ -1151,7 +1151,7 @@ void PairLJLongCoulLongOMP::eval_outer(int iiform, int iito, ThrData * const thr
       fpair = (force_coul+force_lj)*r2inv;
 
       if (NEWTON_PAIR || j < nlocal) {
-        register double *fj = f0+(j+(j<<1)), f;
+        double *fj = f0+(j+(j<<1)), f;
         fi[0] += f = d[0]*fpair; fj[0] -= f;
         fi[1] += f = d[1]*fpair; fj[1] -= f;
         fi[2] += f = d[2]*fpair; fj[2] -= f;
diff --git a/src/USER-OMP/pair_lj_long_tip4p_long_omp.cpp b/src/USER-OMP/pair_lj_long_tip4p_long_omp.cpp
index 8c0d72c74180a5dc61a32f6da8576196c5ee825e..19278d19bb6fbbe7c6fd70da5d94f49d4bffac14 100644
--- a/src/USER-OMP/pair_lj_long_tip4p_long_omp.cpp
+++ b/src/USER-OMP/pair_lj_long_tip4p_long_omp.cpp
@@ -812,8 +812,8 @@ void PairLJLongTIP4PLongOMP::eval(int iifrom, int iito, ThrData * const thr)
         r2inv = 1.0/rsq;
         if (ORDER6) {                   // long-range lj
           if (!LJTABLE || rsq <= tabinnerdispsq) {
-            register double rn = r2inv*r2inv*r2inv;
-            register double x2 = g2*rsq, a2 = 1.0/x2;
+            double rn = r2inv*r2inv*r2inv;
+            double x2 = g2*rsq, a2 = 1.0/x2;
             x2 = a2*exp(-x2)*lj4i[jtype];
             if (ni == 0) {
               forcelj =
@@ -822,7 +822,7 @@ void PairLJLongTIP4PLongOMP::eval(int iifrom, int iito, ThrData * const thr)
                 evdwl = rn*lj3i[jtype]-g6*((a2+1.0)*a2+0.5)*x2;
             }
             else {                  // special case
-              register double f = special_lj[ni], t = rn*(1.0-f);
+              double f = special_lj[ni], t = rn*(1.0-f);
               forcelj = f*(rn *= rn)*lj1i[jtype]-
                 g8*(((6.0*a2+6.0)*a2+3.0)*a2+1.0)*x2*rsq+t*lj2i[jtype];
               if (EFLAG)
@@ -830,30 +830,30 @@ void PairLJLongTIP4PLongOMP::eval(int iifrom, int iito, ThrData * const thr)
             }
           }
           else {                                        // table real space
-            register union_int_float_t disp_t;
+            union_int_float_t disp_t;
             disp_t.f = rsq;
-            register const int disp_k = (disp_t.i & ndispmask)>>ndispshiftbits;
-            register double f_disp = (rsq-rdisptable[disp_k])*drdisptable[disp_k];
-            register double rn = r2inv*r2inv*r2inv;
+            const int disp_k = (disp_t.i & ndispmask)>>ndispshiftbits;
+            double f_disp = (rsq-rdisptable[disp_k])*drdisptable[disp_k];
+            double rn = r2inv*r2inv*r2inv;
             if (ni == 0) {
               forcelj = (rn*=rn)*lj1i[jtype]-(fdisptable[disp_k]+f_disp*dfdisptable[disp_k])*lj4i[jtype];
               if (EFLAG) evdwl = rn*lj3i[jtype]-(edisptable[disp_k]+f_disp*dedisptable[disp_k])*lj4i[jtype];
             }
             else {                  // special case
-              register double f = special_lj[ni], t = rn*(1.0-f);
+              double f = special_lj[ni], t = rn*(1.0-f);
               forcelj = f*(rn *= rn)*lj1i[jtype]-(fdisptable[disp_k]+f_disp*dfdisptable[disp_k])*lj4i[jtype]+t*lj2i[jtype];
               if (EFLAG) evdwl = f*rn*lj3i[jtype]-(edisptable[disp_k]+f_disp*dedisptable[disp_k])*lj4i[jtype]+t*lj4i[jtype];
             }
           }
         }
         else {                      // cut lj
-          register double rn = r2inv*r2inv*r2inv;
+          double rn = r2inv*r2inv*r2inv;
           if (ni == 0) {
             forcelj = rn*(rn*lj1i[jtype]-lj2i[jtype]);
             if (EFLAG) evdwl = rn*(rn*lj3i[jtype]-lj4i[jtype])-offseti[jtype];
           }
           else {                    // special case
-            register double f = special_lj[ni];
+            double f = special_lj[ni];
             forcelj = f*rn*(rn*lj1i[jtype]-lj2i[jtype]);
             if (EFLAG)
               evdwl = f * (rn*(rn*lj3i[jtype]-lj4i[jtype])-offseti[jtype]);
@@ -1191,15 +1191,15 @@ void PairLJLongTIP4PLongOMP::eval_inner(int iifrom, int iito, ThrData * const th
 
       if (rsq < cut_ljsq[itype][jtype] && rsq < cut_out_off_sq ) {  // lj
         r2inv = 1.0/rsq;
-        register double rn = r2inv*r2inv*r2inv;
+        double rn = r2inv*r2inv*r2inv;
         if (ni == 0) forcelj = rn*(rn*lj1i[jtype]-lj2i[jtype]);
         else {                  // special case
-          register double f = special_lj[ni];
+          double f = special_lj[ni];
           forcelj = f*rn*(rn*lj1i[jtype]-lj2i[jtype]);
         }
 
         if (rsq > cut_out_on_sq) {                        // switching
-          register double rsw = (sqrt(rsq) - cut_out_on)/cut_out_diff;
+          double rsw = (sqrt(rsq) - cut_out_on)/cut_out_diff;
           forcelj  *= 1.0 + rsw*rsw*(2.0*rsw-3.0);
         }
 
@@ -1260,7 +1260,7 @@ void PairLJLongTIP4PLongOMP::eval_inner(int iifrom, int iito, ThrData * const th
           }
 
           if (rsq > cut_out_on_sq) {                        // switching
-            register double rsw = (sqrt(rsq) - cut_out_on)/cut_out_diff;
+            double rsw = (sqrt(rsq) - cut_out_on)/cut_out_diff;
             forcecoul  *= 1.0 + rsw*rsw*(2.0*rsw-3.0);
           }
 
@@ -1445,19 +1445,19 @@ void PairLJLongTIP4PLongOMP::eval_middle(int iifrom, int iito, ThrData * const t
 
       if (rsq < cut_ljsq[itype][jtype] && rsq >= cut_in_off_sq && rsq <= cut_out_off_sq ) {  // lj
         r2inv = 1.0/rsq;
-        register double rn = r2inv*r2inv*r2inv;
+        double rn = r2inv*r2inv*r2inv;
         if (ni == 0) forcelj = rn*(rn*lj1i[jtype]-lj2i[jtype]);
         else {                  // special case
-          register double f = special_lj[ni];
+          double f = special_lj[ni];
           forcelj = f*rn*(rn*lj1i[jtype]-lj2i[jtype]);
         }
 
         if (rsq < cut_in_on_sq) {                                // switching
-          register double rsw = (sqrt(rsq) - cut_in_off)/cut_in_diff;
+          double rsw = (sqrt(rsq) - cut_in_off)/cut_in_diff;
           forcelj  *= rsw*rsw*(3.0 - 2.0*rsw);
         }
         if (rsq > cut_out_on_sq) {
-          register double rsw = (sqrt(rsq) - cut_out_on)/cut_out_diff;
+          double rsw = (sqrt(rsq) - cut_out_on)/cut_out_diff;
           forcelj  *= 1.0 + rsw*rsw*(2.0*rsw-3.0);
         }
 
@@ -1518,11 +1518,11 @@ void PairLJLongTIP4PLongOMP::eval_middle(int iifrom, int iito, ThrData * const t
           }
 
           if (rsq < cut_in_on_sq) {                                // switching
-            register double rsw = (sqrt(rsq) - cut_in_off)/cut_in_diff;
+            double rsw = (sqrt(rsq) - cut_in_off)/cut_in_diff;
             forcecoul  *= rsw*rsw*(3.0 - 2.0*rsw);
           }
           if (rsq > cut_out_on_sq) {
-            register double rsw = (sqrt(rsq) - cut_out_on)/cut_out_diff;
+            double rsw = (sqrt(rsq) - cut_out_on)/cut_out_diff;
             forcecoul  *= 1.0 + rsw*rsw*(2.0*rsw-3.0);
           }
 
@@ -1719,18 +1719,18 @@ void PairLJLongTIP4PLongOMP::eval_outer(int iifrom, int iito, ThrData * const th
         frespa = 1.0;                                       // check whether and how to compute respa corrections
         respa_flag = rsq < cut_in_on_sq ? 1 : 0;
         if (respa_flag && (rsq > cut_in_off_sq)) {
-          register double rsw = (sqrt(rsq)-cut_in_off)/cut_in_diff;
+          double rsw = (sqrt(rsq)-cut_in_off)/cut_in_diff;
           frespa = 1-rsw*rsw*(3.0-2.0*rsw);
         }
 
         r2inv = 1.0/rsq;
-        register double rn = r2inv*r2inv*r2inv;
+        double rn = r2inv*r2inv*r2inv;
         if (respa_flag) respa_lj = ni == 0 ?                 // correct for respa
                           frespa*rn*(rn*lj1i[jtype]-lj2i[jtype]) :
                           frespa*rn*(rn*lj1i[jtype]-lj2i[jtype])*special_lj[ni];
         if (ORDER6) {                                        // long-range form
           if (!ndisptablebits || rsq <= tabinnerdispsq) {
-            register double x2 = g2*rsq, a2 = 1.0/x2;
+            double x2 = g2*rsq, a2 = 1.0/x2;
             x2 = a2*exp(-x2)*lj4i[jtype];
             if (ni == 0) {
               forcelj =
@@ -1738,7 +1738,7 @@ void PairLJLongTIP4PLongOMP::eval_outer(int iifrom, int iito, ThrData * const th
               if (EFLAG) evdwl = rn*lj3i[jtype]-g6*((a2+1.0)*a2+0.5)*x2;
             }
             else {                                        // correct for special
-              register double f = special_lj[ni], t = rn*(1.0-f);
+              double f = special_lj[ni], t = rn*(1.0-f);
               forcelj = f*(rn *= rn)*lj1i[jtype]-
                 g8*(((6.0*a2+6.0)*a2+3.0)*a2+1.0)*x2*rsq+t*lj2i[jtype]-respa_lj;
               if (EFLAG)
@@ -1746,16 +1746,16 @@ void PairLJLongTIP4PLongOMP::eval_outer(int iifrom, int iito, ThrData * const th
             }
           }
           else {                        // table real space
-            register union_int_float_t disp_t;
+            union_int_float_t disp_t;
             disp_t.f = rsq;
-            register const int disp_k = (disp_t.i & ndispmask)>>ndispshiftbits;
-            register double f_disp = (rsq-rdisptable[disp_k])*drdisptable[disp_k];
+            const int disp_k = (disp_t.i & ndispmask)>>ndispshiftbits;
+            double f_disp = (rsq-rdisptable[disp_k])*drdisptable[disp_k];
             if (ni == 0) {
               forcelj = (rn*=rn)*lj1i[jtype]-(fdisptable[disp_k]+f_disp*dfdisptable[disp_k])*lj4i[jtype]-respa_lj;
               if (EFLAG) evdwl = rn*lj3i[jtype]-(edisptable[disp_k]+f_disp*dedisptable[disp_k])*lj4i[jtype];
             }
             else {                  // special case
-              register double f = special_lj[ni], t = rn*(1.0-f);
+              double f = special_lj[ni], t = rn*(1.0-f);
               forcelj = f*(rn *= rn)*lj1i[jtype]-(fdisptable[disp_k]+f_disp*dfdisptable[disp_k])*lj4i[jtype]+t*lj2i[jtype]-respa_lj;
               if (EFLAG) evdwl = f*rn*lj3i[jtype]-(edisptable[disp_k]+f_disp*dedisptable[disp_k])*lj4i[jtype]+t*lj4i[jtype];
             }
@@ -1767,7 +1767,7 @@ void PairLJLongTIP4PLongOMP::eval_outer(int iifrom, int iito, ThrData * const th
             if (EFLAG) evdwl = rn*(rn*lj3i[jtype]-lj4i[jtype])-offseti[jtype];
           }
           else {                                        // correct for special
-            register double f = special_lj[ni];
+            double f = special_lj[ni];
             forcelj = f*rn*(rn*lj1i[jtype]-lj2i[jtype])-respa_lj;
             if (EFLAG)
               evdwl = f*(rn*(rn*lj3i[jtype]-lj4i[jtype])-offseti[jtype]);
@@ -1832,16 +1832,16 @@ void PairLJLongTIP4PLongOMP::eval_outer(int iifrom, int iito, ThrData * const th
           frespa = 1.0;                                       // check whether and how to compute respa corrections
           respa_flag = rsq < cut_in_on_sq ? 1 : 0;
           if (respa_flag && (rsq > cut_in_off_sq)) {
-            register double rsw = (sqrt(rsq)-cut_in_off)/cut_in_diff;
+            double rsw = (sqrt(rsq)-cut_in_off)/cut_in_diff;
             frespa = 1-rsw*rsw*(3.0-2.0*rsw);
           }
 
           r2inv = 1.0 / rsq;
           if (!CTABLE || rsq <= tabinnersq) {        // series real space
-            register double r = sqrt(rsq), s = qri*q[j];
+            double r = sqrt(rsq), s = qri*q[j];
             if (respa_flag)                                // correct for respa
               respa_coul = ni == 0 ? frespa*s/r : frespa*s/r*special_coul[ni];
-            register double x = g_ewald*r, t = 1.0/(1.0+EWALD_P*x);
+            double x = g_ewald*r, t = 1.0/(1.0+EWALD_P*x);
             if (ni == 0) {
               s *= g_ewald*exp(-x*x);
               forcecoul = (t *= ((((t*A5+A4)*t+A3)*t+A2)*t+A1)*s/x)+EWALD_F*s-respa_coul;
@@ -1855,13 +1855,13 @@ void PairLJLongTIP4PLongOMP::eval_outer(int iifrom, int iito, ThrData * const th
           }                                                // table real space
           else {
             if (respa_flag) {
-              register double r = sqrt(rsq), s = qri*q[j];
+              double r = sqrt(rsq), s = qri*q[j];
               respa_coul = ni == 0 ? frespa*s/r : frespa*s/r*special_coul[ni];
             }
-            register union_int_float_t t;
+            union_int_float_t t;
             t.f = rsq;
-            register const int k = (t.i & ncoulmask) >> ncoulshiftbits;
-            register double f = (t.f-rtable[k])*drtable[k], qiqj = qtmp*q[j];
+            const int k = (t.i & ncoulmask) >> ncoulshiftbits;
+            double f = (t.f-rtable[k])*drtable[k], qiqj = qtmp*q[j];
             if (ni == 0) {
               forcecoul = qiqj*(ftable[k]+f*dftable[k]);
               if (EFLAG) ecoul = qiqj*(etable[k]+f*detable[k]);
diff --git a/src/compute_rdf.cpp b/src/compute_rdf.cpp
index 263e7f799ba02eb254dcc422fdddd86f2bfa798f..f2635c09485ea9770b4adcd0770d207a5be83330 100644
--- a/src/compute_rdf.cpp
+++ b/src/compute_rdf.cpp
@@ -169,7 +169,7 @@ void ComputeRDF::init()
       cutghost = comm->cutghostuser;
 
     if (mycutneigh > cutghost)
-      error->all(FLERR,"Compure rdf cutoff exceeds ghost atom range - "
+      error->all(FLERR,"Compute rdf cutoff exceeds ghost atom range - "
                  "use comm_modify cutoff command");
     if (force->pair && mycutneigh < force->pair->cutforce + skin)
       if (comm->me == 0)
diff --git a/src/math_complex.h b/src/math_complex.h
index e4d889c4e6fa438a5febe706a1041af9d85ec648..4cf25b6dacaa83cdccfd2c9c732607b1ca4ca84c 100644
--- a/src/math_complex.h
+++ b/src/math_complex.h
@@ -32,12 +32,12 @@ typedef struct complex {
   d.im = x.re*y.im+x.im*y.re; }
 
 #define C_RMULT(d, x, y) { \
-  register complex t = x; \
+  complex t = x; \
   d.re = t.re*y.re-t.im*y.im; \
   d.im = t.re*y.im+t.im*y.re; }
 
 #define C_CRMULT(d, x, y) { \
-  register complex t = x; \
+  complex t = x; \
   d.re = t.re*y.re-t.im*y.im; \
   d.im = -t.re*y.im-t.im*y.re; }
 
@@ -62,7 +62,7 @@ typedef struct complex {
   d.im = y; }
 
 #define C_ANGLE(d, angle) { \
-  register double a = angle; \
+  double a = angle; \
   d.re = cos(a); \
   d.im = sin(a); }
 
diff --git a/src/math_vector.h b/src/math_vector.h
index c2227572729a96367f3ca35d1d11fd12c1391592..ff260b794d19a7ee4c599420e4a88ca278465e84 100644
--- a/src/math_vector.h
+++ b/src/math_vector.h
@@ -62,7 +62,7 @@ inline void vec_neg(vector &dest) {                                // -a
   dest[2] = -dest[2]; }
 
 inline void vec_norm(vector &dest) {                                 // a/|a|
-  register double f = sqrt(vec_dot(dest, dest));
+  double f = sqrt(vec_dot(dest, dest));
   dest[0] /= f;
   dest[1] /= f;
   dest[2] /= f; }
@@ -190,7 +190,7 @@ inline void form_subtr(shape &dest, form &src) {                // m_a-m_b
   dest[3] -= src[3]; dest[4] -= src[4]; dest[5] -= src[5]; }
 
 inline int form_inv(form &m_inv, form &m) {                        // m^-1
-  register double det = form_det(m);
+  double det = form_det(m);
   if (fzero(det)) return 0;
   m_inv[0] = (m[1]*m[2]-m[3]*m[3])/det;
   m_inv[1] = (m[0]*m[2]-m[4]*m[4])/det;
@@ -345,7 +345,7 @@ inline void form4_unit(form4 &dest) {
   dest[0] = dest[1] = dest[2] = dest[3] = 1.0; }
 
 inline double form4_det(form4 &m) {
-  register double f = m[6]*m[7]-m[5]*m[8];
+  double f = m[6]*m[7]-m[5]*m[8];
   return m[0]*(
       m[1]*(m[2]*m[3]-m[4]*m[4])+
       m[5]*(2.0*m[4]*m[7]-m[2]*m[5])-m[3]*m[7]*m[7])+f*f+
@@ -355,7 +355,7 @@ inline double form4_det(form4 &m) {
         m[9]*(m[4]*m[4]-m[2]*m[3])); }
 
 inline int form4_inv(form4 &m_inv, form4 &m) {
-  register double det = form4_det(m);
+  double det = form4_det(m);
   if (fzero(det)) return 0;
   m_inv[0] = (m[1]*(m[2]*m[3]-m[4]*m[4])+
       m[5]*(2.0*m[4]*m[7]-m[2]*m[5])-m[3]*m[7]*m[7])/det;
diff --git a/src/pack.h b/src/pack.h
index 5db35117c55162ed2045ca37f1b2850a1c8922e3..066535f5c9fc576094ef9a059e090a2d10393286 100644
--- a/src/pack.h
+++ b/src/pack.h
@@ -55,8 +55,8 @@ struct pack_plan_3d {
 
 static void pack_3d(PACK_DATA *data, PACK_DATA *buf, struct pack_plan_3d *plan)
 {
-  register int in,out,fast,mid,slow;
-  register int nfast,nmid,nslow,nstride_line,nstride_plane,plane;
+  int in,out,fast,mid,slow;
+  int nfast,nmid,nslow,nstride_line,nstride_plane,plane;
 
   nfast = plan->nfast;
   nmid = plan->nmid;
@@ -81,8 +81,8 @@ static void pack_3d(PACK_DATA *data, PACK_DATA *buf, struct pack_plan_3d *plan)
 
 static void unpack_3d(PACK_DATA *buf, PACK_DATA *data, struct pack_plan_3d *plan)
 {
-  register int in,out,fast,mid,slow;
-  register int nfast,nmid,nslow,nstride_line,nstride_plane,plane;
+  int in,out,fast,mid,slow;
+  int nfast,nmid,nslow,nstride_line,nstride_plane,plane;
 
   nfast = plan->nfast;
   nmid = plan->nmid;
@@ -107,8 +107,8 @@ static void unpack_3d(PACK_DATA *buf, PACK_DATA *data, struct pack_plan_3d *plan
 
 static void unpack_3d_permute1_1(PACK_DATA *buf, PACK_DATA *data, struct pack_plan_3d *plan)
 {
-  register int in,out,fast,mid,slow;
-  register int nfast,nmid,nslow,nstride_line,nstride_plane,plane;
+  int in,out,fast,mid,slow;
+  int nfast,nmid,nslow,nstride_line,nstride_plane,plane;
 
   nfast = plan->nfast;
   nmid = plan->nmid;
@@ -133,8 +133,8 @@ static void unpack_3d_permute1_1(PACK_DATA *buf, PACK_DATA *data, struct pack_pl
 
 static void unpack_3d_permute1_2(PACK_DATA *buf, PACK_DATA *data, struct pack_plan_3d *plan)
 {
-  register int in,out,fast,mid,slow;
-  register int nfast,nmid,nslow,nstride_line,nstride_plane,plane;
+  int in,out,fast,mid,slow;
+  int nfast,nmid,nslow,nstride_line,nstride_plane,plane;
 
   nfast = plan->nfast;
   nmid = plan->nmid;
@@ -162,8 +162,8 @@ static void unpack_3d_permute1_2(PACK_DATA *buf, PACK_DATA *data, struct pack_pl
 static void unpack_3d_permute1_n(PACK_DATA *buf, PACK_DATA *data, struct pack_plan_3d *plan)
 
 {
-  register int in,out,iqty,instart,fast,mid,slow;
-  register int nfast,nmid,nslow,nstride_line,nstride_plane,plane,nqty;
+  int in,out,iqty,instart,fast,mid,slow;
+  int nfast,nmid,nslow,nstride_line,nstride_plane,plane,nqty;
 
   nfast = plan->nfast;
   nmid = plan->nmid;
@@ -192,8 +192,8 @@ static void unpack_3d_permute1_n(PACK_DATA *buf, PACK_DATA *data, struct pack_pl
 static void unpack_3d_permute2_1(PACK_DATA *buf, PACK_DATA *data, struct pack_plan_3d *plan)
 
 {
-  register int in,out,fast,mid,slow;
-  register int nfast,nmid,nslow,nstride_line,nstride_plane;
+  int in,out,fast,mid,slow;
+  int nfast,nmid,nslow,nstride_line,nstride_plane;
 
   nfast = plan->nfast;
   nmid = plan->nmid;
@@ -218,8 +218,8 @@ static void unpack_3d_permute2_1(PACK_DATA *buf, PACK_DATA *data, struct pack_pl
 static void unpack_3d_permute2_2(PACK_DATA *buf, PACK_DATA *data, struct pack_plan_3d *plan)
 
 {
-  register int in,out,fast,mid,slow;
-  register int nfast,nmid,nslow,nstride_line,nstride_plane;
+  int in,out,fast,mid,slow;
+  int nfast,nmid,nslow,nstride_line,nstride_plane;
 
   nfast = plan->nfast;
   nmid = plan->nmid;
@@ -246,8 +246,8 @@ static void unpack_3d_permute2_2(PACK_DATA *buf, PACK_DATA *data, struct pack_pl
 static void unpack_3d_permute2_n(PACK_DATA *buf, PACK_DATA *data, struct pack_plan_3d *plan)
 
 {
-  register int in,out,iqty,instart,fast,mid,slow;
-  register int nfast,nmid,nslow,nstride_line,nstride_plane,nqty;
+  int in,out,iqty,instart,fast,mid,slow;
+  int nfast,nmid,nslow,nstride_line,nstride_plane,nqty;
 
   nfast = plan->nfast;
   nmid = plan->nmid;
@@ -283,9 +283,9 @@ static void unpack_3d_permute2_n(PACK_DATA *buf, PACK_DATA *data, struct pack_pl
 static void pack_3d(PACK_DATA *data, PACK_DATA *buf, struct pack_plan_3d *plan)
 
 {
-  register PACK_DATA *in,*out,*begin,*end;
-  register int mid,slow;
-  register int nfast,nmid,nslow,nstride_line,nstride_plane,plane;
+  PACK_DATA *in,*out,*begin,*end;
+  int mid,slow;
+  int nfast,nmid,nslow,nstride_line,nstride_plane,plane;
 
   nfast = plan->nfast;
   nmid = plan->nmid;
@@ -312,9 +312,9 @@ static void pack_3d(PACK_DATA *data, PACK_DATA *buf, struct pack_plan_3d *plan)
 static void unpack_3d(PACK_DATA *buf, PACK_DATA *data, struct pack_plan_3d *plan)
 
 {
-  register PACK_DATA *in,*out,*begin,*end;
-  register int mid,slow;
-  register int nfast,nmid,nslow,nstride_line,nstride_plane,plane;
+  PACK_DATA *in,*out,*begin,*end;
+  int mid,slow;
+  int nfast,nmid,nslow,nstride_line,nstride_plane,plane;
 
   nfast = plan->nfast;
   nmid = plan->nmid;
@@ -341,9 +341,9 @@ static void unpack_3d(PACK_DATA *buf, PACK_DATA *data, struct pack_plan_3d *plan
 static void unpack_3d_permute1_1(PACK_DATA *buf, PACK_DATA *data, struct pack_plan_3d *plan)
 
 {
-  register PACK_DATA *in,*out,*begin,*end;
-  register int mid,slow;
-  register int nfast,nmid,nslow,nstride_line,nstride_plane,plane;
+  PACK_DATA *in,*out,*begin,*end;
+  int mid,slow;
+  int nfast,nmid,nslow,nstride_line,nstride_plane,plane;
 
   nfast = plan->nfast;
   nmid = plan->nmid;
@@ -370,9 +370,9 @@ static void unpack_3d_permute1_1(PACK_DATA *buf, PACK_DATA *data, struct pack_pl
 static void unpack_3d_permute1_2(PACK_DATA *buf, PACK_DATA *data, struct pack_plan_3d *plan)
 
 {
-  register PACK_DATA *in,*out,*begin,*end;
-  register int mid,slow;
-  register int nfast,nmid,nslow,nstride_line,nstride_plane,plane;
+  PACK_DATA *in,*out,*begin,*end;
+  int mid,slow;
+  int nfast,nmid,nslow,nstride_line,nstride_plane,plane;
 
   nfast = plan->nfast;
   nmid = plan->nmid;
@@ -401,9 +401,9 @@ static void unpack_3d_permute1_2(PACK_DATA *buf, PACK_DATA *data, struct pack_pl
 static void unpack_3d_permute1_n(PACK_DATA *buf, PACK_DATA *data, struct pack_plan_3d *plan)
 
 {
-  register PACK_DATA *in,*out,*instart,*begin,*end;
-  register int iqty,mid,slow;
-  register int nfast,nmid,nslow,nstride_line,nstride_plane,plane,nqty;
+  PACK_DATA *in,*out,*instart,*begin,*end;
+  int iqty,mid,slow;
+  int nfast,nmid,nslow,nstride_line,nstride_plane,plane,nqty;
 
   nfast = plan->nfast;
   nmid = plan->nmid;
@@ -433,9 +433,9 @@ static void unpack_3d_permute1_n(PACK_DATA *buf, PACK_DATA *data, struct pack_pl
 static void unpack_3d_permute2_1(PACK_DATA *buf, PACK_DATA *data, struct pack_plan_3d *plan)
 
 {
-  register PACK_DATA *in,*out,*begin,*end;
-  register int mid,slow;
-  register int nfast,nmid,nslow,nstride_line,nstride_plane;
+  PACK_DATA *in,*out,*begin,*end;
+  int mid,slow;
+  int nfast,nmid,nslow,nstride_line,nstride_plane;
 
   nfast = plan->nfast;
   nmid = plan->nmid;
@@ -461,9 +461,9 @@ static void unpack_3d_permute2_1(PACK_DATA *buf, PACK_DATA *data, struct pack_pl
 static void unpack_3d_permute2_2(PACK_DATA *buf, PACK_DATA *data, struct pack_plan_3d *plan)
 
 {
-  register PACK_DATA *in,*out,*begin,*end;
-  register int mid,slow;
-  register int nfast,nmid,nslow,nstride_line,nstride_plane;
+  PACK_DATA *in,*out,*begin,*end;
+  int mid,slow;
+  int nfast,nmid,nslow,nstride_line,nstride_plane;
 
   nfast = plan->nfast;
   nmid = plan->nmid;
@@ -491,9 +491,9 @@ static void unpack_3d_permute2_2(PACK_DATA *buf, PACK_DATA *data, struct pack_pl
 static void unpack_3d_permute2_n(PACK_DATA *buf, PACK_DATA *data, struct pack_plan_3d *plan)
 
 {
-  register PACK_DATA *in,*out,*instart,*begin,*end;
-  register int iqty,mid,slow;
-  register int nfast,nmid,nslow,nstride_line,nstride_plane,nqty;
+  PACK_DATA *in,*out,*instart,*begin,*end;
+  int iqty,mid,slow;
+  int nfast,nmid,nslow,nstride_line,nstride_plane,nqty;
 
   nfast = plan->nfast;
   nmid = plan->nmid;
@@ -532,9 +532,9 @@ static void unpack_3d_permute2_n(PACK_DATA *buf, PACK_DATA *data, struct pack_pl
 static void pack_3d(PACK_DATA *data, PACK_DATA *buf, struct pack_plan_3d *plan)
 
 {
-  register PACK_DATA *in,*out;
-  register int mid,slow,size;
-  register int nfast,nmid,nslow,nstride_line,nstride_plane,plane,upto;
+  PACK_DATA *in,*out;
+  int mid,slow,size;
+  int nfast,nmid,nslow,nstride_line,nstride_plane,plane,upto;
 
   nfast = plan->nfast;
   nmid = plan->nmid;
@@ -561,9 +561,9 @@ static void pack_3d(PACK_DATA *data, PACK_DATA *buf, struct pack_plan_3d *plan)
 static void unpack_3d(PACK_DATA *buf, PACK_DATA *data, struct pack_plan_3d *plan)
 
 {
-  register PACK_DATA *in,*out;
-  register int mid,slow,size;
-  register int nfast,nmid,nslow,nstride_line,nstride_plane,plane,upto;
+  PACK_DATA *in,*out;
+  int mid,slow,size;
+  int nfast,nmid,nslow,nstride_line,nstride_plane,plane,upto;
 
   nfast = plan->nfast;
   nmid = plan->nmid;
@@ -590,9 +590,9 @@ static void unpack_3d(PACK_DATA *buf, PACK_DATA *data, struct pack_plan_3d *plan
 static void unpack_3d_permute1_1(PACK_DATA *buf, PACK_DATA *data, struct pack_plan_3d *plan)
 
 {
-  register PACK_DATA *in,*out,*begin,*end;
-  register int mid,slow;
-  register int nfast,nmid,nslow,nstride_line,nstride_plane,plane;
+  PACK_DATA *in,*out,*begin,*end;
+  int mid,slow;
+  int nfast,nmid,nslow,nstride_line,nstride_plane,plane;
 
   nfast = plan->nfast;
   nmid = plan->nmid;
@@ -619,9 +619,9 @@ static void unpack_3d_permute1_1(PACK_DATA *buf, PACK_DATA *data, struct pack_pl
 static void unpack_3d_permute1_2(PACK_DATA *buf, PACK_DATA *data, struct pack_plan_3d *plan)
 
 {
-  register PACK_DATA *in,*out,*begin,*end;
-  register int mid,slow;
-  register int nfast,nmid,nslow,nstride_line,nstride_plane,plane;
+  PACK_DATA *in,*out,*begin,*end;
+  int mid,slow;
+  int nfast,nmid,nslow,nstride_line,nstride_plane,plane;
 
   nfast = plan->nfast;
   nmid = plan->nmid;
@@ -650,9 +650,9 @@ static void unpack_3d_permute1_2(PACK_DATA *buf, PACK_DATA *data, struct pack_pl
 static void unpack_3d_permute1_n(PACK_DATA *buf, PACK_DATA *data, struct pack_plan_3d *plan)
 
 {
-  register PACK_DATA *in,*out,*instart,*begin,*end;
-  register int iqty,mid,slow;
-  register int nfast,nmid,nslow,nstride_line,nstride_plane,plane,nqty;
+  PACK_DATA *in,*out,*instart,*begin,*end;
+  int iqty,mid,slow;
+  int nfast,nmid,nslow,nstride_line,nstride_plane,plane,nqty;
 
   nfast = plan->nfast;
   nmid = plan->nmid;
@@ -682,9 +682,9 @@ static void unpack_3d_permute1_n(PACK_DATA *buf, PACK_DATA *data, struct pack_pl
 static void unpack_3d_permute2_1(PACK_DATA *buf, PACK_DATA *data, struct pack_plan_3d *plan)
 
 {
-  register PACK_DATA *in,*out,*begin,*end;
-  register int mid,slow;
-  register int nfast,nmid,nslow,nstride_line,nstride_plane;
+  PACK_DATA *in,*out,*begin,*end;
+  int mid,slow;
+  int nfast,nmid,nslow,nstride_line,nstride_plane;
 
   nfast = plan->nfast;
   nmid = plan->nmid;
@@ -710,9 +710,9 @@ static void unpack_3d_permute2_1(PACK_DATA *buf, PACK_DATA *data, struct pack_pl
 static void unpack_3d_permute2_2(PACK_DATA *buf, PACK_DATA *data, struct pack_plan_3d *plan)
 
 {
-  register PACK_DATA *in,*out,*begin,*end;
-  register int mid,slow;
-  register int nfast,nmid,nslow,nstride_line,nstride_plane;
+  PACK_DATA *in,*out,*begin,*end;
+  int mid,slow;
+  int nfast,nmid,nslow,nstride_line,nstride_plane;
 
   nfast = plan->nfast;
   nmid = plan->nmid;
@@ -740,9 +740,9 @@ static void unpack_3d_permute2_2(PACK_DATA *buf, PACK_DATA *data, struct pack_pl
 static void unpack_3d_permute2_n(PACK_DATA *buf, PACK_DATA *data, struct pack_plan_3d *plan)
 
 {
-  register PACK_DATA *in,*out,*instart,*begin,*end;
-  register int iqty,mid,slow;
-  register int nfast,nmid,nslow,nstride_line,nstride_plane,nqty;
+  PACK_DATA *in,*out,*instart,*begin,*end;
+  int iqty,mid,slow;
+  int nfast,nmid,nslow,nstride_line,nstride_plane,nqty;
 
   nfast = plan->nfast;
   nmid = plan->nmid;
diff --git a/src/pair.cpp b/src/pair.cpp
index 3b839dff497808a033a97f9c551381b80658e267..73dbb9428e37939a7a407ad70d8b9ea516a8a133 100644
--- a/src/pair.cpp
+++ b/src/pair.cpp
@@ -563,7 +563,7 @@ void Pair::init_tables_disp(double cut_lj_global)
       rsq_lookup.i |= maskhi;
     }
     rsq = rsq_lookup.f;
-    register double x2 = g2*rsq, a2 = 1.0/x2;
+    double x2 = g2*rsq, a2 = 1.0/x2;
     x2 = a2*exp(-x2);
 
     rdisptable[i] = rsq_lookup.f;
@@ -609,7 +609,7 @@ void Pair::init_tables_disp(double cut_lj_global)
   if (rsq_lookup.f < (cut_lj_globalsq = cut_lj_global * cut_lj_global)) {
     rsq_lookup.f = cut_lj_globalsq;
 
-    register double x2 = g2*rsq, a2 = 1.0/x2;
+    double x2 = g2*rsq, a2 = 1.0/x2;
     x2 = a2*exp(-x2);
     f_tmp = g8*(((6.0*a2+6.0)*a2+3.0)*a2+1.0)*x2*rsq;
     e_tmp = g6*((a2+1.0)*a2+0.5)*x2;