Add ACKS2 charge equilibration method to REAXFF and support for electric fields in QEq

doc/src/fix_acks2_reaxff.rst

@@ -0,0 +1,97 @@
+.. index:: fix acks2/reaxff
+
+fix acks2/reaxff command
+====================
+
+Syntax
+""""""
+
+.. parsed-literal::
+
+   fix ID group-ID acks2/reaxff Nevery cutlo cuthi tolerance params args
+
+* ID, group-ID are documented in :doc:`fix <fix>` command
+* acks2/reaxff = style name of this fix command
+* Nevery = perform QEq every this many steps
+* cutlo,cuthi = lo and hi cutoff for Taper radius
+* tolerance = precision to which charges will be equilibrated
+* params = reaxff or a filename
+
+Examples
+""""""""
+
+.. code-block:: LAMMPS
+
+   fix 1 all acks2/reaxff 1 0.0 10.0 1.0e-6 reaxff
+   fix 1 all acks2/reaxff 1 0.0 10.0 1.0e-6 param.acks2
+
+Description
+"""""""""""
+
+Perform the charge equilibration (QEq) method as described in :ref:`(Verstraelen) <Verstraelen>`.  It is
+typically used in conjunction with the ReaxFF force field model as
+implemented in the :doc:`pair_style reaxff <pair_reaxff>` command, but
+it can be used with any potential in LAMMPS, so long as it defines and
+uses charges on each atom.  For more technical details about the
+charge equilibration performed by fix acks2/reaxff, see the
+:ref:`(O'Hearn) <O'Hearn>` paper.
+
+The QEq method minimizes the electrostatic energy of the system by
+adjusting the partial charge on individual atoms based on interactions
+with their neighbors.  It requires some parameters for each atom type.
+If the *params* setting above is the word "reaxff", then these are
+extracted from the :doc:`pair_style reaxff <pair_reaxff>` command and
+the ReaxFF force field file it reads in.  If a file name is specified
+for *params*\ , then the parameters are taken from the specified file
+and the file must contain one line for each atom type.  The latter
+form must be used when performing QeQ with a non-ReaxFF potential.
+Each line should be formatted as follows:
+
+.. parsed-literal::
+
+   bond_softness
+   itype chi eta gamma bcut
+
+where the first line is the global parameter *bond_softness*. The remaining
+1 to Ntypes lines include *itype*, the atom type from 1 to Ntypes, *chi*, the
+electronegativity in eV, *eta*, the self-Coulomb
+potential in eV, *gamma*, the valence orbital
+exponent, and *bcut*, the bond cutoff distance.  Note that these 4 quantities are also in the ReaxFF
+potential file, except that eta is defined here as twice the eta value
+in the ReaxFF file. Note that unlike the rest of LAMMPS, the units
+of this fix are hard-coded to be A, eV, and electronic charge.
+
+**Restart, fix_modify, output, run start/stop, minimize info:**
+
+No information about this fix is written to :doc:`binary restart files <restart>`.  No global scalar or vector or per-atom
+quantities are stored by this fix for access by various :doc:`output commands <Howto_output>`.  No parameter of this fix can be used
+with the *start/stop* keywords of the :doc:`run <run>` command.
+
+This fix is invoked during :doc:`energy minimization <minimize>`.
+
+Restrictions
+""""""""""""
+
+This fix is part of the REAXFF package.  It is only enabled if
+LAMMPS was built with that package.  See the :doc:`Build package <Build_package>` doc page for more info.
+
+This fix does not correctly handle interactions
+involving multiple periodic images of the same atom. Hence, it should not
+be used for periodic cell dimensions less than 10 angstroms.
+
+Related commands
+""""""""""""""""
+
+:doc:`pair_style reaxff <pair_reaxff>`
+
+**Default:** none
+
+----------
+
+.. _O'Hearn:
+
+**(O'Hearn)** O'Hearn, Alperen, Aktulga, SIAM J. Sci. Comput., 42(1), C1–C22 (2020).
+
+.. _Verstraelen:
+
+**(Verstraelen)** Verstraelen, Ayers, Speybroeck, Waroquier, J. Chem. Phys. 138, 074108 (2013).

doc/src/pair_reaxff.rst

@@ -20,7 +20,7 @@ Syntax
   .. parsed-literal::
 
      keyword = *checkqeq* or *lgvdw* or *safezone* or *mincap* or *minhbonds*
-       *checkqeq* value = *yes* or *no* = whether or not to require qeq/reaxff fix
+       *checkqeq* value = *yes* or *no* = whether or not to require qeq/reaxff or acks2/reax fix
        *enobonds* value = *yes* or *no* = whether or not to tally energy of atoms with no bonds
        *lgvdw* value = *yes* or *no* = whether or not to use a low gradient vdW correction
        *safezone* = factor used for array allocation
@@ -119,7 +119,8 @@ The ReaxFF parameter files provided were created using a charge
 equilibration (QEq) model for handling the electrostatic interactions.
 Therefore, by default, LAMMPS requires that either the
 :doc:`fix qeq/reaxff <fix_qeq_reaxff>` or the
-:doc:`fix qeq/shielded <fix_qeq>` command be used with
+:doc:`fix qeq/shielded <fix_qeq>` or :doc:`fix acks2/reaxff <fix_acks2>`
+command be used with
 *pair_style reaxff* when simulating a ReaxFF model, to equilibrate
 the charges each timestep.
 
@@ -128,7 +129,8 @@ for the QEq fixes, allowing a simulation to be run without charge
 equilibration. In this case, the static charges you assign to each
 atom will be used for computing the electrostatic interactions in
 the system. See the :doc:`fix qeq/reaxff <fix_qeq_reaxff>` or
-:doc:`fix qeq/shielded <fix_qeq>` command documentation for more details.
+:doc:`fix qeq/shielded <fix_qeq>` or :doc:`fix acks2/reaxff <fix_acks2>`
+command documentation for more details.
 
 Using the optional keyword *lgvdw* with the value *yes* turns on the
 low-gradient correction of ReaxFF for long-range London Dispersion,
@@ -352,7 +354,8 @@ Related commands
 """"""""""""""""
 
 :doc:`pair_coeff <pair_coeff>`, :doc:`fix qeq/reaxff <fix_qeq_reaxff>`,
-:doc:`fix reaxff/bonds <fix_reaxff_bonds>`, :doc:`fix reaxff/species <fix_reaxff_species>`
+:doc:`fix acks2/reax <fix_acks2_reax>`, :doc:`fix reaxff/bonds <fix_reaxff_bonds>`,
+:doc:`fix reaxff/species <fix_reaxff_species>`
 
 Default
 """""""

src/.gitignore

@@ -613,6 +613,8 @@
 /fft3d_wrap.h
 /fix_accelerate_cos.cpp
 /fix_accelerate_cos.h
+/fix_acks2_reaxff.cpp
+/fix_acks2_reaxff.h
 /fix_adapt_fep.cpp
 /fix_adapt_fep.h
 /fix_addtorque.cpp

src/KOKKOS/Install.sh

@@ -109,6 +109,8 @@ action domain_kokkos.h
 action fftdata_kokkos.h fft3d.h
 action fft3d_kokkos.cpp fft3d.cpp
 action fft3d_kokkos.h fft3d.h
+action fix_acks2_reaxff_kokkos.cpp fix_acks2_reaxff.cpp
+action fix_acks2_reaxff_kokkos.h fix_acks2_reaxff.h
 action fix_deform_kokkos.cpp
 action fix_deform_kokkos.h
 action fix_enforce2d_kokkos.cpp

src/KOKKOS/fix_acks2_reaxff_kokkos.h

@@ -0,0 +1,368 @@
+/* -*- c++ -*- ----------------------------------------------------------
+   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+   http://lammps.sandia.gov, Sandia National Laboratories
+   Steve Plimpton, sjplimp@sandia.gov
+
+   Copyright (2003) Sandia Corporation.  Under the terms of Contract
+   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+   certain rights in this software.  This software is distributed under
+   the GNU General Public License.
+
+   See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+#ifdef FIX_CLASS
+
+FixStyle(acks2/reaxff/kk,FixACKS2ReaxFFKokkos<LMPDeviceType>)
+FixStyle(acks2/reaxff/kk/device,FixACKS2ReaxFFKokkos<LMPDeviceType>)
+FixStyle(acks2/reaxff/kk/host,FixACKS2ReaxFFKokkos<LMPHostType>)
+
+#else
+
+#ifndef LMP_FIX_ACKS2_REAXFF_KOKKOS_H
+#define LMP_FIX_ACKS2_REAXFF_KOKKOS_H
+
+#include "fix_acks2_reaxff.h"
+#include "kokkos_type.h"
+#include "neigh_list.h"
+#include "neigh_list_kokkos.h"
+
+namespace LAMMPS_NS {
+
+struct TagACKS2Zero{};
+struct TagACKS2InitMatvec{};
+struct TagACKS2SparseMatvec1{};
+struct TagACKS2SparseMatvec2{};
+
+template<int NEIGHFLAG>
+struct TagACKS2SparseMatvec3_Half{};
+
+struct TagACKS2SparseMatvec3_Full{};
+struct TagACKS2Norm1{};
+struct TagACKS2Norm2{};
+struct TagACKS2Norm3{};
+struct TagACKS2Dot1{};
+struct TagACKS2Dot2{};
+struct TagACKS2Dot3{};
+struct TagACKS2Dot4{};
+struct TagACKS2Dot5{};
+struct TagACKS2Precon1A{};
+struct TagACKS2Precon1B{};
+struct TagACKS2Precon2{};
+struct TagACKS2Add{};
+struct TagACKS2ZeroQGhosts{};
+struct TagACKS2CalculateQ1{};
+struct TagACKS2CalculateQ2{};
+
+template<class DeviceType>
+class FixACKS2ReaxFFKokkos : public FixACKS2ReaxFF {
+ public:
+  typedef DeviceType device_type;
+  typedef double value_type;
+  typedef ArrayTypes<DeviceType> AT;
+  FixACKS2ReaxFFKokkos(class LAMMPS *, int, char **);
+  ~FixACKS2ReaxFFKokkos();
+
+  void cleanup_copy();
+  void init();
+  void setup_pre_force(int);
+  void pre_force(int);
+
+  DAT::tdual_ffloat_1d get_s() {return k_s;}
+
+  KOKKOS_INLINE_FUNCTION
+  void num_neigh_item(int, int&) const;
+
+  KOKKOS_INLINE_FUNCTION
+  void operator()(TagACKS2Zero, const int&) const;
+
+  KOKKOS_INLINE_FUNCTION
+  void operator()(TagACKS2InitMatvec, const int&) const;
+
+  template<int NEIGHFLAG>
+  KOKKOS_INLINE_FUNCTION
+  void compute_h_item(int, int &, const bool &) const;
+
+  template<int NEIGHFLAG>
+  KOKKOS_INLINE_FUNCTION
+  void compute_h_team(const typename Kokkos::TeamPolicy <DeviceType> ::member_type &team, int, int) const;
+
+  template<int NEIGHFLAG>
+  KOKKOS_INLINE_FUNCTION
+  void compute_x_item(int, int &, const bool &) const;
+
+  template<int NEIGHFLAG>
+  KOKKOS_INLINE_FUNCTION
+  void compute_x_team(const typename Kokkos::TeamPolicy <DeviceType> ::member_type &team, int, int) const;
+
+  KOKKOS_INLINE_FUNCTION
+  void operator()(TagACKS2SparseMatvec1, const int&) const;
+
+  KOKKOS_INLINE_FUNCTION
+  void operator()(TagACKS2SparseMatvec2, const int&) const;
+
+  template<int NEIGHFLAG>
+  KOKKOS_INLINE_FUNCTION
+  void operator()(TagACKS2SparseMatvec3_Half<NEIGHFLAG>, const int&) const;
+
+  typedef typename Kokkos::TeamPolicy<DeviceType, TagACKS2SparseMatvec3_Full>::member_type membertype;
+  KOKKOS_INLINE_FUNCTION
+  void operator() (TagACKS2SparseMatvec3_Full, const membertype &team) const;
+
+  KOKKOS_INLINE_FUNCTION
+  void operator()(TagACKS2Norm1, const int&, double&) const;
+
+  KOKKOS_INLINE_FUNCTION
+  void operator()(TagACKS2Norm2, const int&, double&) const;
+
+  KOKKOS_INLINE_FUNCTION
+  void operator()(TagACKS2Norm3, const int&, double&) const;
+
+  KOKKOS_INLINE_FUNCTION
+  void operator()(TagACKS2Dot1, const int&, double&) const;
+
+  KOKKOS_INLINE_FUNCTION
+  void operator()(TagACKS2Dot2, const int&, double&) const;
+
+  KOKKOS_INLINE_FUNCTION
+  void operator()(TagACKS2Dot3, const int&, double&) const;
+
+  KOKKOS_INLINE_FUNCTION
+  void operator()(TagACKS2Dot4, const int&, double&) const;
+
+  KOKKOS_INLINE_FUNCTION
+  void operator()(TagACKS2Dot5, const int&, double&) const;
+
+  KOKKOS_INLINE_FUNCTION
+  void operator()(TagACKS2Precon1A, const int&) const;
+
+  KOKKOS_INLINE_FUNCTION
+  void operator()(TagACKS2Precon1B, const int&) const;
+
+  KOKKOS_INLINE_FUNCTION
+  void operator()(TagACKS2Precon2, const int&) const;
+
+  KOKKOS_INLINE_FUNCTION
+  void operator()(TagACKS2Add, const int&) const;
+
+  KOKKOS_INLINE_FUNCTION
+  void operator()(TagACKS2ZeroQGhosts, const int&) const;
+
+  KOKKOS_INLINE_FUNCTION
+  void operator()(TagACKS2CalculateQ1, const int&) const;
+
+  KOKKOS_INLINE_FUNCTION
+  void operator()(TagACKS2CalculateQ2, const int&) const;
+
+  KOKKOS_INLINE_FUNCTION
+  double calculate_H_k(const F_FLOAT &r, const F_FLOAT &shld) const;
+
+  KOKKOS_INLINE_FUNCTION
+  double calculate_X_k(const F_FLOAT &r, const F_FLOAT &bcut) const;
+
+  struct params_acks2{
+    KOKKOS_INLINE_FUNCTION
+    params_acks2(){chi=0;eta=0;gamma=0;bcut_acks2=0;};
+    KOKKOS_INLINE_FUNCTION
+    params_acks2(int i){chi=0;eta=0;gamma=0;bcut_acks2=0;};
+    F_FLOAT chi, eta, gamma, bcut_acks2;
+  };
+
+ private:
+  int inum;
+  int allocated_flag, last_allocate;
+  int need_dup,comm_me_0_flag;
+
+  typename AT::t_int_scalar d_mfill_offset;
+
+  typedef Kokkos::DualView<int***,DeviceType> tdual_int_1d;
+  Kokkos::DualView<params_acks2*,Kokkos::LayoutRight,DeviceType> k_params;
+  typename Kokkos::DualView<params_acks2*, Kokkos::LayoutRight,DeviceType>::t_dev_const params;
+
+  typename ArrayTypes<DeviceType>::t_x_array x;
+  typename ArrayTypes<DeviceType>::t_v_array v;
+  typename ArrayTypes<DeviceType>::t_f_array_const f;
+  typename ArrayTypes<DeviceType>::t_ffloat_1d_randomread mass;
+  typename ArrayTypes<DeviceType>::t_ffloat_1d q;
+  typename ArrayTypes<DeviceType>::t_int_1d type, mask;
+  typename ArrayTypes<DeviceType>::t_tagint_1d tag;
+
+  typename ArrayTypes<DeviceType>::t_neighbors_2d d_neighbors;
+  typename ArrayTypes<DeviceType>::t_int_1d_randomread d_ilist, d_numneigh;
+
+  DAT::tdual_ffloat_1d k_tap;
+  typename AT::t_ffloat_1d d_tap;
+
+  DAT::tdual_ffloat_2d k_bcut;
+  typename AT::t_ffloat_2d d_bcut;
+
+  typename AT::t_int_1d d_firstnbr;
+  typename AT::t_int_1d d_numnbrs;
+  typename AT::t_int_1d d_jlist;
+  typename AT::t_ffloat_1d d_val;
+
+  typename AT::t_int_1d d_firstnbr_X;
+  typename AT::t_int_1d d_numnbrs_X;
+  typename AT::t_int_1d d_jlist_X;
+  typename AT::t_ffloat_1d d_val_X;
+
+  DAT::tdual_ffloat_1d k_s, k_X_diag, k_chi_field;
+  typename AT::t_ffloat_1d d_Hdia_inv,d_Xdia_inv, d_X_diag, d_chi_field, d_b_s,  d_s;
+  typename AT::t_ffloat_1d_randomread r_b_s, r_s;
+
+  DAT::tdual_ffloat_1d k_d, k_q_hat, k_z, k_y;
+  typename AT::t_ffloat_1d d_p, d_q, d_r, d_d, d_g, d_q_hat, d_r_hat, d_y, d_z, d_bb, d_xx;
+  typename AT::t_ffloat_1d_randomread r_p, r_r, r_d;
+
+  DAT::tdual_ffloat_2d k_shield, k_s_hist, k_s_hist_X, k_s_hist_last;
+  typename AT::t_ffloat_2d d_shield, d_s_hist, d_s_hist_X, d_s_hist_last;
+  typename AT::t_ffloat_2d_randomread r_s_hist, r_s_hist_X, r_s_hist_last;
+
+  Kokkos::Experimental::ScatterView<F_FLOAT*, typename AT::t_ffloat_1d::array_layout, typename KKDevice<DeviceType>::value, Kokkos::Experimental::ScatterSum, Kokkos::Experimental::ScatterDuplicated> dup_X_diag;
+  Kokkos::Experimental::ScatterView<F_FLOAT*, typename AT::t_ffloat_1d::array_layout, typename KKDevice<DeviceType>::value, Kokkos::Experimental::ScatterSum, Kokkos::Experimental::ScatterNonDuplicated> ndup_X_diag;
+
+  Kokkos::Experimental::ScatterView<F_FLOAT*, typename AT::t_ffloat_1d::array_layout, typename KKDevice<DeviceType>::value, Kokkos::Experimental::ScatterSum, Kokkos::Experimental::ScatterDuplicated> dup_bb;
+  Kokkos::Experimental::ScatterView<F_FLOAT*, typename AT::t_ffloat_1d::array_layout, typename KKDevice<DeviceType>::value, Kokkos::Experimental::ScatterSum, Kokkos::Experimental::ScatterNonDuplicated> ndup_bb;
+
+  void init_shielding_k();
+  void init_hist();
+  void allocate_matrix();
+  void allocate_array();
+  void deallocate_array();
+  int bicgstab_solve();
+  void calculate_Q();
+
+  int neighflag;
+  int nlocal,nall,nmax,newton_pair;
+  int count, isuccess;
+  double alpha, beta, omega, cutsq;
+
+  int iswap;
+  int first;
+  typename AT::t_int_2d d_sendlist;
+  typename AT::t_xfloat_1d_um v_buf;
+
+  void grow_arrays(int);
+  void copy_arrays(int, int, int);
+  int pack_exchange(int, double *);
+  int unpack_exchange(int, double *);
+  void get_chi_field();
+  double memory_usage();
+
+  void sparse_matvec_acks2(typename AT::t_ffloat_1d &, typename AT::t_ffloat_1d &);
+};
+
+template <class DeviceType>
+struct FixACKS2ReaxFFKokkosNumNeighFunctor  {
+  typedef DeviceType device_type;
+  typedef int value_type;
+  FixACKS2ReaxFFKokkos<DeviceType> c;
+  FixACKS2ReaxFFKokkosNumNeighFunctor(FixACKS2ReaxFFKokkos<DeviceType>* c_ptr):c(*c_ptr) {
+    c.cleanup_copy();
+  };
+  KOKKOS_INLINE_FUNCTION
+  void operator()(const int ii, int &maxneigh) const {
+    c.num_neigh_item(ii, maxneigh);
+  }
+};
+
+template <class DeviceType, int NEIGHFLAG>
+struct FixACKS2ReaxFFKokkosComputeHFunctor {
+  int atoms_per_team, vector_length;
+  typedef int value_type;
+  typedef Kokkos::ScratchMemorySpace<DeviceType> scratch_space;
+  FixACKS2ReaxFFKokkos<DeviceType> c;
+
+  FixACKS2ReaxFFKokkosComputeHFunctor(FixACKS2ReaxFFKokkos<DeviceType>* c_ptr):c(*c_ptr) {
+    c.cleanup_copy();
+  };
+
+  FixACKS2ReaxFFKokkosComputeHFunctor(FixACKS2ReaxFFKokkos<DeviceType> *c_ptr,
+                                  int _atoms_per_team, int _vector_length)
+      : c(*c_ptr), atoms_per_team(_atoms_per_team),
+        vector_length(_vector_length) {
+    c.cleanup_copy();
+  };
+
+  KOKKOS_INLINE_FUNCTION
+  void operator()(const int ii, int &m_fill, const bool &final) const {
+    c.template compute_h_item<NEIGHFLAG>(ii,m_fill,final);
+  }
+
+  KOKKOS_INLINE_FUNCTION
+  void operator()(
+      const typename Kokkos::TeamPolicy<DeviceType>::member_type &team) const {
+    c.template compute_h_team<NEIGHFLAG>(team, atoms_per_team, vector_length);
+  }
+
+  size_t team_shmem_size(int team_size) const {
+    size_t shmem_size =
+        Kokkos::View<int *, scratch_space, Kokkos::MemoryUnmanaged>::shmem_size(
+            atoms_per_team) + // s_ilist
+        Kokkos::View<int *, scratch_space, Kokkos::MemoryUnmanaged>::shmem_size(
+            atoms_per_team) + // s_numnbrs
+        Kokkos::View<int *, scratch_space, Kokkos::MemoryUnmanaged>::shmem_size(
+            atoms_per_team) + // s_firstnbr
+        Kokkos::View<int **, scratch_space, Kokkos::MemoryUnmanaged>::
+            shmem_size(atoms_per_team, vector_length) + // s_jtype
+        Kokkos::View<int **, scratch_space, Kokkos::MemoryUnmanaged>::
+            shmem_size(atoms_per_team, vector_length) + // s_j
+        Kokkos::View<F_FLOAT **, scratch_space,
+                     Kokkos::MemoryUnmanaged>::shmem_size(atoms_per_team,
+                                                          vector_length); // s_r
+    return shmem_size;
+  }
+};
+
+template <class DeviceType, int NEIGHFLAG>
+struct FixACKS2ReaxFFKokkosComputeXFunctor {
+  int atoms_per_team, vector_length;
+  typedef int value_type;
+  typedef Kokkos::ScratchMemorySpace<DeviceType> scratch_space;
+  FixACKS2ReaxFFKokkos<DeviceType> c;
+
+  FixACKS2ReaxFFKokkosComputeXFunctor(FixACKS2ReaxFFKokkos<DeviceType>* c_ptr):c(*c_ptr) {
+    c.cleanup_copy();
+  };
+
+  FixACKS2ReaxFFKokkosComputeXFunctor(FixACKS2ReaxFFKokkos<DeviceType> *c_ptr,
+                                  int _atoms_per_team, int _vector_length)
+      : c(*c_ptr), atoms_per_team(_atoms_per_team),
+        vector_length(_vector_length) {
+    c.cleanup_copy();
+  };
+
+  KOKKOS_INLINE_FUNCTION
+  void operator()(const int ii, int &m_fill, const bool &final) const {
+    c.template compute_x_item<NEIGHFLAG>(ii,m_fill,final);
+  }
+
+  KOKKOS_INLINE_FUNCTION
+  void operator()(
+      const typename Kokkos::TeamPolicy<DeviceType>::member_type &team) const {
+    c.template compute_x_team<NEIGHFLAG>(team, atoms_per_team, vector_length);
+  }
+
+  size_t team_shmem_size(int team_size) const {
+    size_t shmem_size =
+        Kokkos::View<int *, scratch_space, Kokkos::MemoryUnmanaged>::shmem_size(
+            atoms_per_team) + // s_ilist
+        Kokkos::View<int *, scratch_space, Kokkos::MemoryUnmanaged>::shmem_size(
+            atoms_per_team) + // s_numnbrs
+        Kokkos::View<int *, scratch_space, Kokkos::MemoryUnmanaged>::shmem_size(
+            atoms_per_team) + // s_firstnbr
+        Kokkos::View<int **, scratch_space, Kokkos::MemoryUnmanaged>::
+            shmem_size(atoms_per_team, vector_length) + // s_jtype
+        Kokkos::View<int **, scratch_space, Kokkos::MemoryUnmanaged>::
+            shmem_size(atoms_per_team, vector_length) + // s_j
+        Kokkos::View<F_FLOAT **, scratch_space,
+                     Kokkos::MemoryUnmanaged>::shmem_size(atoms_per_team,
+                                                          vector_length); // s_r
+    return shmem_size;
+  }
+};
+
+}
+
+#endif
+#endif

src/KOKKOS/fix_qeq_reaxff_kokkos.cpp

@@ -75,6 +75,7 @@ FixQEqReaxFFKokkos<DeviceType>::~FixQEqReaxFFKokkos()
 
   memoryKK->destroy_kokkos(k_s_hist,s_hist);
   memoryKK->destroy_kokkos(k_t_hist,t_hist);
+  memoryKK->destroy_kokkos(k_chi_field,chi_field);
 }
 
 /* ---------------------------------------------------------------------- */
@@ -256,6 +257,9 @@ void FixQEqReaxFFKokkos<DeviceType>::pre_force(int /*vflag*/)
 
   // init_matvec
 
+  if (field_flag)
+    get_chi_field();
+
   k_s_hist.template sync<DeviceType>();
   k_t_hist.template sync<DeviceType>();
   FixQEqReaxFFKokkosMatVecFunctor<DeviceType> matvec_functor(this);
@@ -373,10 +377,16 @@ void FixQEqReaxFFKokkos<DeviceType>::allocate_array()
     k_d = DAT::tdual_ffloat_1d("qeq/kk:d",nmax);
     d_d = k_d.template view<DeviceType>();
     h_d = k_d.h_view;
+
+    memoryKK->create_kokkos(k_chi_field,chi_field,nmax,"acks2/kk:chi_field");
+    d_chi_field = k_chi_field.template view<DeviceType>();
   }
 
   // init_storage
 
+  if (field_flag)
+    get_chi_field();
+
   FixQEqReaxFFKokkosZeroFunctor<DeviceType> zero_functor(this);
   Kokkos::parallel_for(ignum,zero_functor);
 
@@ -392,7 +402,7 @@ void FixQEqReaxFFKokkos<DeviceType>::zero_item(int ii) const
 
   if (mask[i] & groupbit) {
     d_Hdia_inv[i] = 1.0 / params(itype).eta;
-    d_b_s[i] = -params(itype).chi;
+    d_b_s[i] = -params(itype).chi - d_chi_field[i];
     d_b_t[i] = -1.0;
     d_s[i] = 0.0;
     d_t[i] = 0.0;
@@ -711,7 +721,7 @@ void FixQEqReaxFFKokkos<DeviceType>::matvec_item(int ii) const
 
   if (mask[i] & groupbit) {
     d_Hdia_inv[i] = 1.0 / params(itype).eta;
-    d_b_s[i] = -params(itype).chi;
+    d_b_s[i] = -params(itype).chi - d_chi_field[i];
     d_b_t[i] = -1.0;
     d_t[i] = d_t_hist(i,2) + 3*(d_t_hist(i,0) - d_t_hist(i,1));
     d_s[i] = 4*(d_s_hist(i,0)+d_s_hist(i,2))-(6*d_s_hist(i,1)+d_s_hist(i,3));
@@ -1575,6 +1585,16 @@ int FixQEqReaxFFKokkos<DeviceType>::unpack_exchange(int nlocal, double *buf)
 
 /* ---------------------------------------------------------------------- */
 
+template<class DeviceType>
+void FixQEqReaxFFKokkos<DeviceType>::get_chi_field()
+{
+  FixQEqReaxFF::get_chi_field();
+  k_chi_field.modify_host();
+  k_chi_field.sync_device();
+}
+
+/* ---------------------------------------------------------------------- */
+
 namespace LAMMPS_NS {
 template class FixQEqReaxFFKokkos<LMPDeviceType>;
 #ifdef LMP_KOKKOS_GPU

src/KOKKOS/fix_qeq_reaxff_kokkos.h

@@ -201,8 +201,8 @@ class FixQEqReaxFFKokkos : public FixQEqReaxFF, public KokkosBase {
   typename AT::t_int_1d d_jlist;
   typename AT::t_ffloat_1d d_val;
 
-  DAT::tdual_ffloat_1d k_t, k_s;
+  DAT::tdual_ffloat_1d k_t, k_s, k_chi_field;
-  typename AT::t_ffloat_1d d_Hdia_inv, d_b_s, d_b_t, d_t, d_s;
+  typename AT::t_ffloat_1d d_Hdia_inv, d_b_s, d_b_t, d_t, d_s, d_chi_field;
   HAT::t_ffloat_1d h_t, h_s;
   typename AT::t_ffloat_1d_randomread r_b_s, r_b_t, r_t, r_s;
 
@@ -241,6 +241,7 @@ class FixQEqReaxFFKokkos : public FixQEqReaxFF, public KokkosBase {
   void copy_arrays(int, int, int);
   int pack_exchange(int, double *);
   int unpack_exchange(int, double *);
+  void get_chi_field();
 };
 
 template <class DeviceType>

src/KOKKOS/kokkos.h

@@ -61,9 +61,11 @@ class KokkosLMP : protected Pointers {
   int neigh_count(int);
 
   template<class DeviceType>
-  int need_dup()
+  int need_dup(int qeq_flag = 0)
   {
     int value = 0;
+    int neighflag = this->neighflag;
+    if (qeq_flag) neighflag = this->neighflag_qeq;
 
     if (neighflag == HALFTHREAD)
       value = std::is_same<typename NeedDup<HALFTHREAD,DeviceType>::value,Kokkos::Experimental::ScatterDuplicated>::value;

src/KOKKOS/pair_reaxff_kokkos.cpp

@@ -23,6 +23,7 @@
 #include "comm.h"
 #include "error.h"
 #include "force.h"
+#include "fix_acks2_reaxff_kokkos.h"
 #include "kokkos.h"
 #include "math_const.h"
 #include "math_special.h"
@@ -141,6 +142,25 @@ void PairReaxFFKokkos<DeviceType>::init_style()
   if (fix_reaxff) modify->delete_fix(fix_id); // not needed in the Kokkos version
   fix_reaxff = nullptr;
 
+  acks2_flag = api->system->acks2_flag;
+  if (acks2_flag) {
+    int ifix = modify->find_fix_by_style("^acks2/reax");
+    Fix* fix = modify->fix[ifix];
+    if (!fix->kokkosable)
+      error->all(FLERR,"Must use Kokkos version of acks2/reax with pair reaxc/kk");
+    if (fix->execution_space == Host) {
+      FixACKS2ReaxFFKokkos<LMPHostType>* acks2_fix = (FixACKS2ReaxFFKokkos<LMPHostType>*) modify->fix[ifix];
+      auto k_s = acks2_fix->get_s();
+      k_s.sync<DeviceType>();
+      d_s = k_s.view<DeviceType>();
+    } else {
+      FixACKS2ReaxFFKokkos<LMPDeviceType>* acks2_fix = (FixACKS2ReaxFFKokkos<LMPDeviceType>*) modify->fix[ifix];
+      auto k_s = acks2_fix->get_s();
+      k_s.sync<DeviceType>();
+      d_s = k_s.view<DeviceType>();
+    }
+  }
+
   // irequest = neigh request made by parent class
 
   neighflag = lmp->kokkos->neighflag;
@@ -230,6 +250,9 @@ void PairReaxFFKokkos<DeviceType>::setup()
     // hydrogen bond
     k_params_sing.h_view(i).p_hbond = api->system->reax_param.sbp[map[i]].p_hbond;
 
+    // acks2
+    k_params_sing.h_view(i).bcut_acks2 = api->system->reax_param.sbp[map[i]].bcut_acks2;
+
     for (j = 1; j <= n; j++) {
       if (map[j] == -1) continue;
 
@@ -690,9 +713,11 @@ void PairReaxFFKokkos<DeviceType>::compute(int eflag_in, int vflag_in)
   tag = atomKK->k_tag.view<DeviceType>();
   type = atomKK->k_type.view<DeviceType>();
   nlocal = atomKK->nlocal;
-  nall = atom->nlocal + atom->nghost;
   newton_pair = force->newton_pair;
 
+  nn = list->inum;
+  NN = list->inum + list->gnum;
+
   const int inum = list->inum;
   const int ignum = inum + list->gnum;
   NeighListKokkos<DeviceType>* k_list = static_cast<NeighListKokkos<DeviceType>*>(list);
@@ -700,6 +725,22 @@ void PairReaxFFKokkos<DeviceType>::compute(int eflag_in, int vflag_in)
   d_neighbors = k_list->d_neighbors;
   d_ilist = k_list->d_ilist;
 
+  if (acks2_flag) {
+    int ifix = modify->find_fix_by_style("^acks2/reax");
+    Fix* fix = modify->fix[ifix];
+    if (fix->execution_space == Host) {
+      FixACKS2ReaxFFKokkos<LMPHostType>* acks2_fix = (FixACKS2ReaxFFKokkos<LMPHostType>*) modify->fix[ifix];
+      auto k_s = acks2_fix->get_s();
+      k_s.sync<DeviceType>();
+      d_s = k_s.view<DeviceType>();
+    } else {
+      FixACKS2ReaxFFKokkos<LMPDeviceType>* acks2_fix = (FixACKS2ReaxFFKokkos<LMPDeviceType>*) modify->fix[ifix];
+      auto k_s = acks2_fix->get_s();
+      k_s.sync<DeviceType>();
+      d_s = k_s.view<DeviceType>();
+    }
+  }
+
   need_dup = lmp->kokkos->need_dup<DeviceType>();
 
   // allocate duplicated memory
@@ -1057,7 +1098,12 @@ void PairReaxFFKokkos<DeviceType>::operator()(PairReaxFFComputePolar<NEIGHFLAG,E
   const F_FLOAT chi = paramssing(itype).chi;
   const F_FLOAT eta = paramssing(itype).eta;
 
-  const F_FLOAT epol = KCALpMOL_to_EV*(chi*qi+(eta/2.0)*qi*qi);
+  F_FLOAT epol = KCALpMOL_to_EV*(chi*qi+(eta/2.0)*qi*qi);
+
+  /* energy due to coupling with kinetic energy potential */
+  if (acks2_flag)
+    epol += KCALpMOL_to_EV*qi*d_s[NN + i];
+
   if (eflag) ev.ecoul += epol;
   //if (eflag_atom) this->template ev_tally<NEIGHFLAG>(ev,i,i,epol,0.0,0.0,0.0,0.0);
   if (eflag_atom) this->template e_tally_single<NEIGHFLAG>(ev,i,epol);
@@ -1195,8 +1241,47 @@ void PairReaxFFKokkos<DeviceType>::operator()(PairReaxFFComputeLJCoulomb<NEIGHFL
     const F_FLOAT shld = paramstwbp(itype,jtype).gamma;
     const F_FLOAT denom1 = rij * rij * rij + shld;
     const F_FLOAT denom3 = pow(denom1,0.3333333333333);
-    const F_FLOAT ecoul = C_ele * qi*qj*Tap/denom3;
+    F_FLOAT ecoul = C_ele * qi*qj*Tap/denom3;
-    const F_FLOAT fcoul = C_ele * qi*qj*(dTap-Tap*rij/denom1)/denom3;
+    F_FLOAT fcoul = C_ele * qi*qj*(dTap-Tap*rij/denom1)/denom3;
+
+    /* contribution to energy and gradients (atoms and cell)
+     * due to geometry-dependent terms in the ACKS2
+     * kinetic energy */
+    if (acks2_flag) {
+
+      /* kinetic energy terms */
+      double xcut = 0.5 * (paramssing(itype).bcut_acks2
+                          + paramssing(jtype).bcut_acks2);
+
+      if (rij <= xcut) {
+        const F_FLOAT d = rij / xcut;
+        const F_FLOAT bond_softness = gp[34] * pow( d, 3.0 )
+                                    * pow( 1.0 - d, 6.0 );
+
+        if (bond_softness > 0.0) {
+          /* Coulombic energy contribution */
+          const F_FLOAT effpot_diff = d_s[NN + i]
+                                    - d_s[NN + j];
+          const F_FLOAT e_ele = -0.5 * KCALpMOL_to_EV * bond_softness
+                                     * SQR( effpot_diff );
+
+          ecoul += e_ele;
+
+          /* forces contribution */
+          F_FLOAT d_bond_softness;
+          d_bond_softness = gp[34]
+                          * 3.0 / xcut * pow( d, 2.0 )
+                          * pow( 1.0 - d, 5.0 ) * (1.0 - 3.0 * d);
+          d_bond_softness = -0.5 * d_bond_softness
+                          * SQR( effpot_diff );
+          d_bond_softness = KCALpMOL_to_EV * d_bond_softness
+                          / rij;
+
+          fcoul += d_bond_softness;
+        }
+      }
+    }
+
 
     const F_FLOAT ftotal = fvdwl + fcoul;
     fxtmp += delx*ftotal;
@@ -1296,15 +1381,53 @@ void PairReaxFFKokkos<DeviceType>::operator()(PairReaxFFComputeTabulatedLJCoulom
     const F_FLOAT evdwl = ((vdW.d*dif + vdW.c)*dif + vdW.b)*dif +
       vdW.a;
 
-    const F_FLOAT ecoul = (((ele.d*dif + ele.c)*dif + ele.b)*dif +
+    F_FLOAT ecoul = (((ele.d*dif + ele.c)*dif + ele.b)*dif +
       ele.a)*qi*qj;
 
     const F_FLOAT fvdwl = ((CEvd.d*dif + CEvd.c)*dif + CEvd.b)*dif +
       CEvd.a;
 
-    const F_FLOAT fcoul = (((CEclmb.d*dif+CEclmb.c)*dif+CEclmb.b)*dif +
+    F_FLOAT fcoul = (((CEclmb.d*dif+CEclmb.c)*dif+CEclmb.b)*dif +
       CEclmb.a)*qi*qj;
 
+    /* contribution to energy and gradients (atoms and cell)
+     * due to geometry-dependent terms in the ACKS2
+     * kinetic energy */
+    if (acks2_flag) {
+
+      /* kinetic energy terms */
+      double xcut = 0.5 * (paramssing(itype).bcut_acks2
+                          + paramssing(jtype).bcut_acks2);
+
+      if (rij <= xcut) {
+        const F_FLOAT d = rij / xcut;
+        const F_FLOAT bond_softness = gp[34] * pow( d, 3.0 )
+                                    * pow( 1.0 - d, 6.0 );
+
+        if (bond_softness > 0.0) {
+          /* Coulombic energy contribution */
+          const F_FLOAT effpot_diff = d_s[NN + i]
+                                    - d_s[NN + j];
+          const F_FLOAT e_ele = -0.5 * KCALpMOL_to_EV * bond_softness
+                                     * SQR( effpot_diff );
+
+          ecoul += e_ele;
+
+          /* forces contribution */
+          F_FLOAT d_bond_softness;
+          d_bond_softness = gp[34]
+                          * 3.0 / xcut * pow( d, 2.0 )
+                          * pow( 1.0 - d, 5.0 ) * (1.0 - 3.0 * d);
+          d_bond_softness = -0.5 * d_bond_softness
+                          * SQR( effpot_diff );
+          d_bond_softness = KCALpMOL_to_EV * d_bond_softness
+                          / rij;
+
+          fcoul += d_bond_softness;
+        }
+      }
+    }
+
     const F_FLOAT ftotal = fvdwl + fcoul;
     fxtmp += delx*ftotal;
     fytmp += dely*ftotal;

src/KOKKOS/pair_reaxff_kokkos.h

@@ -252,12 +252,12 @@ class PairReaxFFKokkos : public PairReaxFF {
   struct params_sing{
     KOKKOS_INLINE_FUNCTION
     params_sing() {mass=0;chi=0;eta=0;r_s=0;r_pi=0;r_pi2=0;valency=0;valency_val=0;valency_e=0;valency_boc=0;nlp_opt=0;
-      p_lp2=0;p_ovun2=0;p_ovun5=0;p_val3=0;p_val5=0;p_hbond=0;};
+      p_lp2=0;p_ovun2=0;p_ovun5=0;p_val3=0;p_val5=0;p_hbond=0;bcut_acks2=0;};
     KOKKOS_INLINE_FUNCTION
     params_sing(int /*i*/) {mass=0;chi=0;eta=0;r_s=0;r_pi=0;r_pi2=0;valency=0;valency_val=0;valency_e=0;valency_boc=0;nlp_opt=0;
-      p_lp2=0;p_ovun2=0;p_ovun5=0;p_val3=0;p_val5=0;p_hbond=0;};
+      p_lp2=0;p_ovun2=0;p_ovun5=0;p_val3=0;p_val5=0;p_hbond=0;bcut_acks2=0;};
     F_FLOAT mass,chi,eta,r_s,r_pi,r_pi2,valency,valency_val,valency_e,valency_boc,nlp_opt,
-      p_lp2,p_ovun2,p_ovun5, p_val3, p_val5, p_hbond;
+      p_lp2,p_ovun2,p_ovun5, p_val3, p_val5, p_hbond, bcut_acks2;
   };
 
   struct params_twbp{
@@ -377,7 +377,7 @@ class PairReaxFFKokkos : public PairReaxFF {
   typename AT::t_float_1d d_tap;
   HAT::t_float_1d h_tap;
 
-  typename AT::t_float_1d d_bo_rij, d_hb_rsq, d_Deltap, d_Deltap_boc, d_total_bo;
+  typename AT::t_float_1d d_bo_rij, d_hb_rsq, d_Deltap, d_Deltap_boc, d_total_bo, d_s;
   typename AT::t_float_1d d_Delta, d_Delta_boc, d_Delta_lp, d_dDelta_lp, d_Delta_lp_temp, d_CdDelta;
   typename AT::t_ffloat_2d_dl d_BO, d_BO_s, d_BO_pi, d_BO_pi2, d_dBOp;
   typename AT::t_ffloat_2d_dl d_dln_BOp_pix, d_dln_BOp_piy, d_dln_BOp_piz;
@@ -426,7 +426,7 @@ class PairReaxFFKokkos : public PairReaxFF {
   typename AT::t_ffloat_2d_dl d_dBOpx, d_dBOpy, d_dBOpz;
 
   int neighflag, newton_pair, maxnumneigh, maxhb, maxbo;
-  int nlocal,nall,eflag,vflag;
+  int nlocal,nn,NN,eflag,vflag,acks2_flag;
   F_FLOAT cut_nbsq, cut_hbsq, cut_bosq, bo_cut, thb_cut, thb_cutsq;
   F_FLOAT bo_cut_bond;
 

src/OPENMP/pair_reaxff_omp.cpp

@@ -101,13 +101,17 @@ PairReaxFFOMP::~PairReaxFFOMP()
 
 void PairReaxFFOMP::init_style()
 {
-  if (!atom->q_flag)
-    error->all(FLERR,"Pair style reaxff/omp requires atom attribute q");
-
   bool have_qeq = ((modify->find_fix_by_style("^qeq/reax") != -1)
-                   || (modify->find_fix_by_style("^qeq/shielded") != -1));
+                   || (modify->find_fix_by_style("^qeq/shielded") != -1)
+                   || (modify->find_fix_by_style("^acks2/reax") != -1));
   if (!have_qeq && qeqflag == 1)
-    error->all(FLERR,"Pair reaxff/omp requires use of fix qeq/reaxff or qeq/shielded");
+    error->all(FLERR,"Pair reaxff/omp requires use of fix qeq/reaxff or qeq/shielded"
+                       " or fix acks2/reaxff");
+
+  int have_acks2 = (modify->find_fix_by_style("^acks2/reax") != -1);
+  api->system->acks2_flag = have_acks2;
+  if (api->system->acks2_flag)
+    error->all(FLERR,"Cannot (yet) use ACKS2 with OPENMP ReaxFF");
 
   api->system->n = atom->nlocal; // my atoms
   api->system->N = atom->nlocal + atom->nghost; // mine + ghosts

src/OPENMP/reaxff_nonbonded_omp.cpp

@@ -221,7 +221,7 @@ namespace ReaxFF {
     data->my_en.e_vdW = total_EvdW;
     data->my_en.e_ele = total_Eele;
 
-    Compute_Polarization_Energy(system, data);
+    Compute_Polarization_Energy(system, data, workspace);
   }
 
 /* ---------------------------------------------------------------------- */
@@ -343,6 +343,6 @@ namespace ReaxFF {
     data->my_en.e_vdW = total_EvdW;
     data->my_en.e_ele = total_Eele;
 
-    Compute_Polarization_Energy(system, data);
+    Compute_Polarization_Energy(system, data, workspace);
   }
 }

src/REAXFF/fix_acks2_reaxff.h

@@ -0,0 +1,90 @@
+/* -*- c++ -*- ----------------------------------------------------------
+   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+   http://lammps.sandia.gov, Sandia National Laboratories
+   Steve Plimpton, sjplimp@sandia.gov
+
+   Copyright (2003) Sandia Corporation.  Under the terms of Contract
+   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+   certain rights in this software.  This software is distributed under
+   the GNU General Public License.
+
+   See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+#ifdef FIX_CLASS
+
+FixStyle(acks2/reax,FixACKS2ReaxFF)
+FixStyle(acks2/reaxff,FixACKS2ReaxFF)
+
+#else
+
+#ifndef LMP_FIX_ACKS2_REAXFF_H
+#define LMP_FIX_ACKS2_REAXFF_H
+
+#include "fix_qeq_reaxff.h"
+
+namespace LAMMPS_NS {
+
+class FixACKS2ReaxFF : public FixQEqReaxFF {
+ public:
+  FixACKS2ReaxFF(class LAMMPS *, int, char **);
+  virtual ~FixACKS2ReaxFF();
+  void post_constructor();
+  virtual void init();
+  void init_storage();
+  virtual void pre_force(int);
+
+  double* get_s() {return s;}
+
+ protected:
+  double **s_hist_X,**s_hist_last;
+
+  double *bcut_acks2,bond_softness,**bcut; // acks2 parameters
+
+  sparse_matrix X;
+  double *Xdia_inv;
+  double *X_diag;
+
+  //BiCGStab storage
+  double *g, *q_hat, *r_hat, *y, *z;
+
+  void pertype_parameters(char*);
+  void init_bondcut();
+  void allocate_storage();
+  void deallocate_storage();
+  void allocate_matrix();
+  void deallocate_matrix();
+
+  void init_matvec();
+  void compute_X();
+  double calculate_X(double,double);
+  void calculate_Q();
+
+  int BiCGStab(double*,double*);
+  void sparse_matvec_acks2(sparse_matrix*,sparse_matrix*,double*,double*);
+
+  int pack_forward_comm(int, int *, double *, int, int *);
+  void unpack_forward_comm(int, int, double *);
+  int pack_reverse_comm(int, int, double *);
+  void unpack_reverse_comm(int, int *, double *);
+  void more_forward_comm(double *);
+  void more_reverse_comm(double *);
+  double memory_usage();
+  virtual void grow_arrays(int);
+  virtual void copy_arrays(int, int, int);
+  virtual int pack_exchange(int, double *);
+  virtual int unpack_exchange(int, double *);
+
+  double parallel_norm( double*, int );
+  double parallel_dot( double*, double*, int );
+  double parallel_vector_acc( double*, int );
+
+  void vector_sum(double*,double,double*,double,double*,int);
+  void vector_add(double*, double, double*,int);
+  void vector_copy(double*, double*,int);
+};
+
+}
+
+#endif
+#endif

src/REAXFF/fix_qeq_reaxff.cpp

@@ -25,6 +25,7 @@
 #include "citeme.h"
 #include "comm.h"
 #include "error.h"
+#include "fix_efield.h"
 #include "force.h"
 #include "group.h"
 #include "memory.h"
@@ -97,10 +98,10 @@ FixQEqReaxFF::FixQEqReaxFF(LAMMPS *lmp, int narg, char **arg) :
     else if (strcmp(arg[iarg],"nowarn") == 0) maxwarn = 0;
     else if (strcmp(arg[iarg],"maxiter") == 0) {
       if (iarg+1 > narg-1)
-        error->all(FLERR,"Illegal fix qeq/reaxff command");
+        error->all(FLERR,fmt::format("Illegal fix {} command", style));
       imax = utils::numeric(FLERR,arg[iarg+1],false,lmp);
       iarg++;
-    } else error->all(FLERR,"Illegal fix qeq/reaxff command");
+    } else error->all(FLERR,fmt::format("Illegal fix {} command", style));
     iarg++;
   }
   shld = nullptr;
@@ -115,6 +116,7 @@ FixQEqReaxFF::FixQEqReaxFF(LAMMPS *lmp, int narg, char **arg) :
 
   Hdia_inv = nullptr;
   b_s = nullptr;
+  chi_field = nullptr;
   b_t = nullptr;
   b_prc = nullptr;
   b_prm = nullptr;
@@ -271,6 +273,7 @@ void FixQEqReaxFF::allocate_storage()
 
   memory->create(Hdia_inv,nmax,"qeq:Hdia_inv");
   memory->create(b_s,nmax,"qeq:b_s");
+  memory->create(chi_field,nmax,"qeq:chi_field");
   memory->create(b_t,nmax,"qeq:b_t");
   memory->create(b_prc,nmax,"qeq:b_prc");
   memory->create(b_prm,nmax,"qeq:b_prm");
@@ -297,6 +300,7 @@ void FixQEqReaxFF::deallocate_storage()
   memory->destroy(b_t);
   memory->destroy(b_prc);
   memory->destroy(b_prm);
+  memory->destroy(chi_field);
 
   memory->destroy(p);
   memory->destroy(q);
@@ -378,6 +382,11 @@ void FixQEqReaxFF::init()
   if (group->count(igroup) == 0)
     error->all(FLERR,"Fix qeq/reaxff group has no atoms");
 
+  field_flag = 0;
+  for (int n = 0; n < modify->nfix; n++)
+    if (utils::strmatch(modify->fix[n]->style,"^efield"))
+      field_flag = 1;
+
   // need a half neighbor list w/ Newton off and ghost neighbors
   // built whenever re-neighboring occurs
 
@@ -502,22 +511,16 @@ void FixQEqReaxFF::min_setup_pre_force(int vflag)
 
 void FixQEqReaxFF::init_storage()
 {
-  int NN;
+  if (field_flag)
-  int *ilist;
+    get_chi_field();
-
-  if (reaxff) {
-    NN = reaxff->list->inum + reaxff->list->gnum;
-    ilist = reaxff->list->ilist;
-  } else {
-    NN = list->inum + list->gnum;
-    ilist = list->ilist;
-  }
 
   for (int ii = 0; ii < NN; ii++) {
     int i = ilist[ii];
     if (atom->mask[i] & groupbit) {
       Hdia_inv[i] = 1. / eta[atom->type[i]];
       b_s[i] = -chi[atom->type[i]];
+      if (field_flag)
+        b_s[i] -= chi_field[i];
       b_t[i] = -1.0;
       b_prc[i] = 0;
       b_prm[i] = 0;
@@ -555,6 +558,9 @@ void FixQEqReaxFF::pre_force(int /*vflag*/)
   if (n > n_cap*DANGER_ZONE || m_fill > m_cap*DANGER_ZONE)
     reallocate_matrix();
 
+  if (field_flag)
+    get_chi_field();
+
   init_matvec();
 
   matvecs_s = CG(b_s, s);       // CG on s - parallel
@@ -594,6 +600,8 @@ void FixQEqReaxFF::init_matvec()
       /* init pre-conditioner for H and init solution vectors */
       Hdia_inv[i] = 1. / eta[atom->type[i]];
       b_s[i]      = -chi[atom->type[i]];
+      if (field_flag)
+        b_s[i] -= chi_field[i];
       b_t[i]      = -1.0;
 
       /* quadratic extrapolation for s & t from previous solutions */
@@ -1060,3 +1068,29 @@ void FixQEqReaxFF::vector_add(double* dest, double c, double* v, int k)
   }
 }
 
+/* ---------------------------------------------------------------------- */
+
+void FixQEqReaxFF::get_chi_field()
+{
+  int nlocal = atom->nlocal;
+
+  memset(&chi_field[0],0.0,atom->nmax*sizeof(double));
+
+  if (!(strcmp(update->unit_style,"real") == 0))
+    error->all(FLERR,"Must use unit_style real with fix qeq/reax and external fields");
+
+  double factor = 1.0/force->qe2f;
+
+  // loop over all fixes, find fix efield
+
+  for (int n = 0; n < modify->nfix; n++) {
+    if (utils::strmatch(modify->fix[n]->style,"^efield")) {
+
+      FixEfield* fix_efield = (FixEfield*) modify->fix[n];
+      double* field_energy = fix_efield->get_energy(); // Real units of kcal/mol/angstrom, need to convert to eV
+
+      for (int i = 0; i < nlocal; i++)
+        chi_field[i] += field_energy[i]*factor;
+    }
+  }
+}

src/REAXFF/fix_qeq_reaxff.h

@@ -45,7 +45,7 @@ class FixQEqReaxFF : public Fix {
   virtual void init();
   void init_list(int, class NeighList *);
   virtual void init_storage();
-  void setup_pre_force(int);
+  virtual void setup_pre_force(int);
   virtual void pre_force(int);
 
   void setup_pre_force_respa(int, int);
@@ -92,6 +92,7 @@ class FixQEqReaxFF : public Fix {
   double *Hdia_inv;
   double *b_s, *b_t;
   double *b_prc, *b_prm;
+  double *chi_field;
 
   //CG storage
   double *p, *q, *r, *d;
@@ -105,7 +106,7 @@ class FixQEqReaxFF : public Fix {
   virtual void deallocate_storage();
   void reallocate_storage();
   virtual void allocate_matrix();
-  void deallocate_matrix();
+  virtual void deallocate_matrix();
   void reallocate_matrix();
 
   virtual void init_matvec();
@@ -134,9 +135,13 @@ class FixQEqReaxFF : public Fix {
   virtual void vector_sum(double *, double, double *, double, double *, int);
   virtual void vector_add(double *, double, double *, int);
 
+  virtual void get_chi_field();
+
   // dual CG support
   int dual_enabled;            // 0: Original, separate s & t optimization; 1: dual optimization
   int matvecs_s, matvecs_t;    // Iteration count for each system
+
+  int field_flag;              // 1: field enabled
 };
 
 }    // namespace LAMMPS_NS

src/REAXFF/pair_reaxff.cpp

@@ -36,6 +36,7 @@
 #include "neigh_request.h"
 #include "neighbor.h"
 #include "update.h"
+#include "fix_acks2_reaxff.h"
 
 #include <cmath>
 #include <cstring>
@@ -155,6 +156,7 @@ PairReaxFF::~PairReaxFF()
     delete [] chi;
     delete [] eta;
     delete [] gamma;
+    delete [] bcut_acks2;
   }
 
   memory->destroy(tmpid);
@@ -179,6 +181,7 @@ void PairReaxFF::allocate()
   chi = new double[n+1];
   eta = new double[n+1];
   gamma = new double[n+1];
+  bcut_acks2 = new double[n+1];
 }
 
 /* ---------------------------------------------------------------------- */
@@ -337,9 +340,19 @@ void PairReaxFF::init_style()
     error->all(FLERR,"Pair style reaxff requires atom attribute q");
 
   bool have_qeq = ((modify->find_fix_by_style("^qeq/reax") != -1)
-                   || (modify->find_fix_by_style("^qeq/shielded") != -1));
+                   || (modify->find_fix_by_style("^qeq/shielded") != -1)
+                   || (modify->find_fix_by_style("^acks2/reax") != -1));
   if (!have_qeq && qeqflag == 1)
-    error->all(FLERR,"Pair reaxff requires use of fix qeq/reaxff or qeq/shielded");
+    error->all(FLERR,"Pair reax/c requires use of fix qeq/reax or qeq/shielded"
+                       " or fix acks2/reax");
+
+  int have_acks2 = (modify->find_fix_by_style("^acks2/reax") != -1);
+  api->system->acks2_flag = have_acks2;
+  if (api->system->acks2_flag) {
+    int ifix = modify->find_fix_by_style("^acks2/reax");
+    FixACKS2ReaxFF* acks2_fix = (FixACKS2ReaxFF*) modify->fix[ifix];
+    api->workspace->s = acks2_fix->get_s();
+  }
 
   api->system->n = atom->nlocal; // my atoms
   api->system->N = atom->nlocal + atom->nghost; // mine + ghosts
@@ -466,6 +479,12 @@ void PairReaxFF::compute(int eflag, int vflag)
   api->system->N = atom->nlocal + atom->nghost; // mine + ghosts
   api->system->bigN = static_cast<int> (atom->natoms);  // all atoms in the system
 
+  if (api->system->acks2_flag) {
+    int ifix = modify->find_fix_by_style("^acks2/reax");
+    FixACKS2ReaxFF* acks2_fix = (FixACKS2ReaxFF*) modify->fix[ifix];
+    api->workspace->s = acks2_fix->get_s();
+  }
+
   // setup data structures
 
   setup();
@@ -731,6 +750,16 @@ void *PairReaxFF::extract(const char *str, int &dim)
       if (map[i] >= 0) gamma[i] = api->system->reax_param.sbp[map[i]].gamma;
       else gamma[i] = 0.0;
     return (void *) gamma;
+   }
+   if (strcmp(str,"bcut_acks2") == 0 && bcut_acks2) {
+    for (int i = 1; i <= atom->ntypes; i++)
+      if (map[i] >= 0) bcut_acks2[i] = api->system->reax_param.sbp[map[i]].bcut_acks2;
+      else bcut_acks2[i] = 0.0;
+    return (void *) bcut_acks2;
+  }
+  if (strcmp(str,"bond_softness") == 0) {
+      double* bond_softness = &api->system->reax_param.gp.l[34];
+    return (void *) bond_softness;
   }
   return nullptr;
 }

src/REAXFF/pair_reaxff.h

@@ -63,7 +63,7 @@ class PairReaxFF : public Pair {
   double cutmax;
   class FixReaxFF *fix_reaxff;
 
-  double *chi, *eta, *gamma;
+  double *chi, *eta, *gamma, *bcut_acks2;
   int qeqflag;
   int setup_flag;
   int firstwarn;

src/REAXFF/reaxff_api.h

@@ -125,7 +125,7 @@ extern void Atom_Energy(reax_system *, control_params *, simulation_data *, stor
 
 // nonbonded
 
-extern void Compute_Polarization_Energy(reax_system *, simulation_data *);
+extern void Compute_Polarization_Energy(reax_system *, simulation_data *, storage *);
 extern void vdW_Coulomb_Energy(reax_system *, control_params *, simulation_data *, storage *,
                                reax_list **);
 extern void Tabulated_vdW_Coulomb_Energy(reax_system *, control_params *, simulation_data *,

src/REAXFF/reaxff_ffield.cpp

@@ -204,6 +204,7 @@ namespace ReaxFF {
           sbp[i].b_o_131    = values.next_double();
           sbp[i].b_o_132    = values.next_double();
           sbp[i].b_o_133    = values.next_double();
+          sbp[i].bcut_acks2  = values.next_double();
 
           // line four
 

src/REAXFF/reaxff_nonbonded.cpp

@@ -32,7 +32,7 @@
 #include <cmath>
 
 namespace ReaxFF {
-  void Compute_Polarization_Energy(reax_system *system, simulation_data *data)
+  void Compute_Polarization_Energy(reax_system *system, simulation_data *data, storage *workspace)
   {
     int  i, type_i;
     double q, en_tmp;
@@ -45,6 +45,12 @@ namespace ReaxFF {
 
       en_tmp = KCALpMOL_to_EV * (system->reax_param.sbp[type_i].chi * q +
                                  (system->reax_param.sbp[type_i].eta / 2.) * SQR(q));
+
+      if (system->acks2_flag) {
+        /* energy due to coupling with kinetic energy potential */
+        en_tmp += KCALpMOL_to_EV * q * workspace->s[ system->N + i ];
+      }
+
       data->my_en.e_pol += en_tmp;
 
       /* tally energy into global or per-atom energy accumulators */
@@ -67,6 +73,7 @@ namespace ReaxFF {
     double dr3gamij_1, dr3gamij_3;
     double e_ele, e_vdW, e_core, SMALL = 0.0001;
     double e_lg, de_lg, r_ij5, r_ij6, re6;
+    double xcut, bond_softness, d_bond_softness, d, effpot_diff;
     two_body_parameters *twbp;
     far_neighbor_data *nbr_pj;
     reax_list *far_nbrs;
@@ -207,7 +214,83 @@ namespace ReaxFF {
       }
     }
 
-    Compute_Polarization_Energy(system, data);
+    /* contribution to energy and gradients (atoms and cell)
+     * due to geometry-dependent terms in the ACKS2
+     * kinetic energy */
+    if (system->acks2_flag)
+    for( i = 0; i < natoms; ++i ) {
+      if (system->my_atoms[i].type < 0) continue;
+      start_i = Start_Index(i, far_nbrs);
+      end_i   = End_Index(i, far_nbrs);
+      orig_i  = system->my_atoms[i].orig_id;
+
+      for( pj = start_i; pj < end_i; ++pj ) {
+	nbr_pj = &(far_nbrs->select.far_nbr_list[pj]);
+	j = nbr_pj->nbr;
+	if (system->my_atoms[j].type < 0) continue;
+	orig_j  = system->my_atoms[j].orig_id;
+
+	flag = 0;
+
+	/* kinetic energy terms */
+	double xcut = 0.5 * ( system->reax_param.sbp[ system->my_atoms[i].type ].bcut_acks2
+			    + system->reax_param.sbp[ system->my_atoms[j].type ].bcut_acks2 );
+
+	if(nbr_pj->d <= xcut) {
+	  if (j < natoms) flag = 1;
+	  else if (orig_i < orig_j) flag = 1;
+	  else if (orig_i == orig_j) {
+	    if (nbr_pj->dvec[2] > SMALL) flag = 1;
+	    else if (fabs(nbr_pj->dvec[2]) < SMALL) {
+	      if (nbr_pj->dvec[1] > SMALL) flag = 1;
+	      else if (fabs(nbr_pj->dvec[1]) < SMALL && nbr_pj->dvec[0] > SMALL)
+		flag = 1;
+	    }
+	  }
+	}
+
+	if (flag) {
+
+	  d = nbr_pj->d / xcut;
+	  bond_softness = system->reax_param.gp.l[34] * pow( d, 3.0 )
+			* pow( 1.0 - d, 6.0 );
+
+	  if ( bond_softness > 0.0 )
+	  {
+	    /* Coulombic energy contribution */
+	    effpot_diff = workspace->s[system->N + i]
+			- workspace->s[system->N + j];
+	    e_ele = -0.5 * KCALpMOL_to_EV * bond_softness
+			 * SQR( effpot_diff );
+
+	    data->my_en.e_ele += e_ele;
+
+	    /* forces contribution */
+	    d_bond_softness = system->reax_param.gp.l[34]
+			    * 3.0 / xcut * pow( d, 2.0 )
+			    * pow( 1.0 - d, 5.0 ) * (1.0 - 3.0 * d);
+	    d_bond_softness = -0.5 * d_bond_softness
+			    * SQR( effpot_diff );
+	    d_bond_softness = KCALpMOL_to_EV * d_bond_softness
+			    / nbr_pj->d;
+
+	    /* tally into per-atom energy */
+	    if (system->pair_ptr->evflag || system->pair_ptr->vflag_atom) {
+	      rvec_ScaledSum( delij, 1., system->my_atoms[i].x,
+				    -1., system->my_atoms[j].x );
+	      f_tmp = -d_bond_softness;
+	      system->pair_ptr->ev_tally(i,j,natoms,1,0.0,e_ele,
+				f_tmp,delij[0],delij[1],delij[2]);
+	    }
+
+	    rvec_ScaledAdd( workspace->f[i], -d_bond_softness, nbr_pj->dvec );
+	    rvec_ScaledAdd( workspace->f[j], d_bond_softness, nbr_pj->dvec );
+	  }
+	}
+      }
+    }
+
+    Compute_Polarization_Energy( system, data, workspace );
   }
 
   void Tabulated_vdW_Coulomb_Energy(reax_system *system, control_params *control,
@@ -222,6 +305,7 @@ namespace ReaxFF {
     double e_vdW, e_ele;
     double CEvd, CEclmb, SMALL = 0.0001;
     double f_tmp, delij[3];
+    double xcut, bond_softness, d_bond_softness, d, effpot_diff;
 
     far_neighbor_data *nbr_pj;
     reax_list *far_nbrs;
@@ -306,7 +390,83 @@ namespace ReaxFF {
       }
     }
 
-    Compute_Polarization_Energy(system, data);
+    /* contribution to energy and gradients (atoms and cell)
+     * due to geometry-dependent terms in the ACKS2
+     * kinetic energy */
+    if (system->acks2_flag)
+    for( i = 0; i < natoms; ++i ) {
+      if (system->my_atoms[i].type < 0) continue;
+      start_i = Start_Index(i, far_nbrs);
+      end_i   = End_Index(i, far_nbrs);
+      orig_i  = system->my_atoms[i].orig_id;
+
+      for( pj = start_i; pj < end_i; ++pj ) {
+	nbr_pj = &(far_nbrs->select.far_nbr_list[pj]);
+	j = nbr_pj->nbr;
+	if (system->my_atoms[j].type < 0) continue;
+	orig_j  = system->my_atoms[j].orig_id;
+
+	flag = 0;
+
+	/* kinetic energy terms */
+	xcut = 0.5 * ( system->reax_param.sbp[ system->my_atoms[i].type ].bcut_acks2
+		     + system->reax_param.sbp[ system->my_atoms[j].type ].bcut_acks2 );
+
+	if(nbr_pj->d <= xcut) {
+	  if (j < natoms) flag = 1;
+	  else if (orig_i < orig_j) flag = 1;
+	  else if (orig_i == orig_j) {
+	    if (nbr_pj->dvec[2] > SMALL) flag = 1;
+	    else if (fabs(nbr_pj->dvec[2]) < SMALL) {
+	      if (nbr_pj->dvec[1] > SMALL) flag = 1;
+	      else if (fabs(nbr_pj->dvec[1]) < SMALL && nbr_pj->dvec[0] > SMALL)
+		flag = 1;
+	    }
+	  }
+	}
+
+	if (flag) {
+
+	  d = nbr_pj->d / xcut;
+	  bond_softness = system->reax_param.gp.l[34] * pow( d, 3.0 )
+			* pow( 1.0 - d, 6.0 );
+
+	  if ( bond_softness > 0.0 )
+	  {
+	    /* Coulombic energy contribution */
+	    effpot_diff = workspace->s[system->N + i]
+			- workspace->s[system->N + j];
+	    e_ele = -0.5 * KCALpMOL_to_EV * bond_softness
+			 * SQR( effpot_diff );
+
+	    data->my_en.e_ele += e_ele;
+
+	    /* forces contribution */
+	    d_bond_softness = system->reax_param.gp.l[34]
+			    * 3.0 / xcut * pow( d, 2.0 )
+			    * pow( 1.0 - d, 5.0 ) * (1.0 - 3.0 * d);
+	    d_bond_softness = -0.5 * d_bond_softness
+			    * SQR( effpot_diff );
+	    d_bond_softness = KCALpMOL_to_EV * d_bond_softness
+			    / nbr_pj->d;
+
+	    /* tally into per-atom energy */
+	    if (system->pair_ptr->evflag || system->pair_ptr->vflag_atom) {
+	      rvec_ScaledSum( delij, 1., system->my_atoms[i].x,
+				    -1., system->my_atoms[j].x );
+	      f_tmp = -d_bond_softness;
+	      system->pair_ptr->ev_tally(i,j,natoms,1,0.0,e_ele,
+				f_tmp,delij[0],delij[1],delij[2]);
+	    }
+
+	    rvec_ScaledAdd( workspace->f[i], -d_bond_softness, nbr_pj->dvec );
+	    rvec_ScaledAdd( workspace->f[j], d_bond_softness, nbr_pj->dvec );
+	  }
+	}
+      }
+    }
+
+    Compute_Polarization_Energy(system, data, workspace);
   }
 
   void LR_vdW_Coulomb(reax_system *system, storage *workspace,

src/REAXFF/reaxff_types.h

@@ -77,6 +77,7 @@ struct single_body_parameters {
   double b_o_131;
   double b_o_132;
   double b_o_133;
+  double bcut_acks2; // ACKS2 bond cutoff 
 
   /* Line four in the field file */
   double p_ovun2;
@@ -212,6 +213,7 @@ struct reax_system {
   LR_lookup_table **LR;
 
   int omp_active;
+  int acks2_flag; 
 };
 
 /* system control parameters */
@@ -340,6 +342,9 @@ struct storage {
   double *CdDeltaReduction;
   int *valence_angle_atom_myoffset;
 
+  /* acks2 */
+  double *s;
+
   reallocate_data realloc;
 };
 

src/fix_efield.cpp

@@ -42,7 +42,7 @@ enum{NONE,CONSTANT,EQUAL,ATOM};
 
 FixEfield::FixEfield(LAMMPS *lmp, int narg, char **arg) :
   Fix(lmp, narg, arg), xstr(nullptr), ystr(nullptr), zstr(nullptr),
-  estr(nullptr), idregion(nullptr), efield(nullptr)
+  estr(nullptr), idregion(nullptr), efield(nullptr), energy(nullptr)
 {
   if (narg < 6) error->all(FLERR,"Illegal fix efield command");
 
@@ -111,6 +111,8 @@ FixEfield::FixEfield(LAMMPS *lmp, int narg, char **arg) :
 
   maxatom = atom->nmax;
   memory->create(efield,maxatom,4,"efield:efield");
+
+  maxatom_energy = 0;
 }
 
 /* ---------------------------------------------------------------------- */
@@ -123,6 +125,7 @@ FixEfield::~FixEfield()
   delete [] estr;
   delete [] idregion;
   memory->destroy(efield);
+  memory->destroy(energy);
 }
 
 /* ---------------------------------------------------------------------- */
@@ -449,3 +452,71 @@ double FixEfield::compute_vector(int n)
   }
   return fsum_all[n+1];
 }
+
+/* ----------------------------------------------------------------------
+   get E
+------------------------------------------------------------------------- */
+
+double* FixEfield::get_energy()
+{
+  double **x = atom->x;
+  double *q = atom->q;
+  int *mask = atom->mask;
+  imageint *image = atom->image;
+
+  int nlocal = atom->nlocal;
+
+  // reallocate energy array if necessary
+
+  if (atom->nmax > maxatom_energy) {
+    maxatom_energy = atom->nmax;
+    memory->destroy(energy);
+    memory->create(energy,maxatom_energy,"efield:energy");
+  }
+  memset(&energy[0],0.0,maxatom_energy*sizeof(double));
+
+  // update region if necessary
+
+  Region *region = NULL;
+  if (iregion >= 0) {
+    region = domain->regions[iregion];
+    region->prematch();
+  }
+
+  int warn_flag_local = 0;
+
+  // constant efield
+
+  if (varflag == CONSTANT) {
+    double unwrap[3];
+
+    // charge interactions
+    // force = qE, potential energy = F dot x in unwrapped coords
+
+    if (qflag) {
+      for (int i = 0; i < nlocal; i++) {
+        if (mask[i] & groupbit) {
+          if (region && !region->match(x[i][0],x[i][1],x[i][2])) continue;
+          const double fx = ex;
+          const double fy = ey;
+          const double fz = ez;
+
+          domain->unmap(x[i],image[i],unwrap);
+          energy[i] -= fx*unwrap[0] + fy*unwrap[1] + fz*unwrap[2];
+          if (fabs(fx*(x[i][0]-unwrap[0])) + fabs(fy*(x[i][1]-unwrap[1])) + 
+              fabs(fz*(x[i][2]-unwrap[2])) > 0.0)
+            warn_flag_local = 1;
+        }
+      }
+    }
+  } else {
+    error->all(FLERR,"Cannot yet use fix qeq/reaxff with variable efield");
+  }
+
+  int warn_flag;
+  MPI_Allreduce(&warn_flag_local,&warn_flag,1,MPI_INT,MPI_SUM,world);
+  if (warn_flag && comm->me == 0)
+    error->warning(FLERR,"Using non-zero image flags in field direction with fix qeq/reaxff");
+
+  return energy;
+}

src/fix_efield.h

@@ -38,6 +38,7 @@ class FixEfield : public Fix {
   double memory_usage();
   double compute_scalar();
   double compute_vector(int);
+  double* get_energy();
 
  private:
   double ex, ey, ez;
@@ -49,8 +50,8 @@ class FixEfield : public Fix {
   double qe2f;
   int qflag, muflag;
 
-  int maxatom;
+  int maxatom,maxatom_energy;
-  double **efield;
+  double **efield,*energy;
 
   int force_flag;
   double fsum[4], fsum_all[4];