Merge branch 'develop' into region-lookup-refactor

src/.gitignore

@@ -453,6 +453,8 @@
 /compute_basal_atom.h
 /compute_body_local.cpp
 /compute_body_local.h
+/compute_born_matrix.cpp
+/compute_born_matrix.h
 /compute_cnp_atom.cpp
 /compute_cnp_atom.h
 /compute_damage_atom.cpp

src/EXTRA-COMPUTE/compute_born_matrix.h

@@ -0,0 +1,97 @@
+/* ----------------------------------------------------------------------
+   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+   http://www.lammps.org/, Sandia National Laboratories
+   Steve Plimpton, sjplimp@sandia.gov
+
+   Copyright (2003) Sandia Corporation.  Under the terms of Contract
+   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+   certain rights in this software.  This software is distributed under
+   the GNU General Public License.
+
+   See the README file in the top-level LAMMPS directory.
+   ------------------------------------------------------------------------- */
+
+/*------------------------------------------------------------------------
+  Contributing Authors : Germain Clavier (TUe), Aidan Thompson (Sandia)
+  --------------------------------------------------------------------------*/
+
+#ifdef COMPUTE_CLASS
+// clang-format off
+ComputeStyle(born/matrix,ComputeBornMatrix);
+// clang-format on
+#else
+
+#ifndef LMP_COMPUTE_BORN_MATRIX_H
+#define LMP_COMPUTE_BORN_MATRIX_H
+
+#include "compute.h"
+
+namespace LAMMPS_NS {
+
+class ComputeBornMatrix : public Compute {
+ public:
+  ComputeBornMatrix(class LAMMPS *, int, char **);
+  ~ComputeBornMatrix() override;
+  void init() override;
+  void init_list(int, class NeighList *) override;
+  void compute_vector() override;
+  double memory_usage() override;
+
+ private:
+  // Born matrix contributions
+
+  void compute_pairs();                     // pair and manybody
+  void compute_bonds();                     // bonds
+  void compute_angles();                    // angles
+  void compute_dihedrals();                 // dihedrals
+  void compute_numdiff();                   // stress virial finite differences
+  void displace_atoms(int, int, double);    // displace atoms
+  void force_clear(int);                    // zero out force array
+  void update_virial();                     // recalculate the virial
+  void restore_atoms(int, int);             // restore atom positions
+  void virial_addon();                      // restore atom positions
+  void reallocate();                        // grow the atom arrays
+
+  int me;             // process rank
+  int nvalues;        // length of elastic tensor
+  int numflag;        // 1 if using finite differences
+  double numdelta;    // size of finite strain
+  int maxatom;        // allocated size of atom arrays
+
+  int pairflag, bondflag, angleflag;
+  int dihedflag, impflag, kspaceflag;
+
+  double *values_local, *values_global;
+  double pos, pos1, dt, nktv2p, ftm2v;
+  class NeighList *list;
+
+  char *id_virial;                  // name of virial compute
+  class Compute *compute_virial;    // pointer to virial compute
+
+  static constexpr int NDIR_VIRIAL = 6;    // dimension of virial and strain vectors
+  static constexpr int NXYZ_VIRIAL = 3;    // number of Cartesian coordinates
+  int revalbe[NDIR_VIRIAL][NDIR_VIRIAL];
+  int virialVtoV[NDIR_VIRIAL];
+  double **temp_x;                   // original coords
+  double **temp_f;                   // original forces
+  double fixedpoint[NXYZ_VIRIAL];    // displacement field origin
+  int dirlist[NDIR_VIRIAL][2];       // strain cartesian indices
+};
+}    // namespace LAMMPS_NS
+
+#endif
+#endif
+
+/* ERROR/WARNING messages:
+
+ E: Illegal ... command
+
+ Self-explanatory.  Check the input script syntax and compare to the
+ documentation for the command.  You can use -echo screen as a
+ command-line option when running LAMMPS to see the offending line.
+
+ E: ... style does not support compute born/matrix
+
+ Some component of the force field (pair, bond, angle...) does not provide
+ a function to return the Born term contribution.
+ */

src/EXTRA-MOLECULE/dihedral_nharmonic.cpp

@@ -25,6 +25,7 @@
 #include "force.h"
 #include "memory.h"
 #include "neighbor.h"
+#include "domain.h"
 
 #include <cmath>
 
@@ -39,6 +40,7 @@ DihedralNHarmonic::DihedralNHarmonic(LAMMPS *lmp) : Dihedral(lmp)
 {
   writedata = 1;
   a = nullptr;
+  born_matrix_enable = 1;
 }
 
 /* ---------------------------------------------------------------------- */
@@ -336,3 +338,64 @@ void DihedralNHarmonic::write_data(FILE *fp)
   }
 
 }
+
+/* ----------------------------------------------------------------------*/
+
+void DihedralNHarmonic::born_matrix(int nd, int i1, int i2, int i3, int i4,
+                             double &dudih, double &du2dih) {
+  int i,type;
+  double vb1x,vb1y,vb1z,vb2x,vb2y,vb2z,vb3x,vb3y,vb3z,vb2xm,vb2ym,vb2zm;
+  double ax,ay,az,bx,by,bz,rasq,rbsq,rgsq,rg,rginv,ra2inv,rb2inv,rabinv;
+  double c,s,kf;
+
+  int **dihedrallist = neighbor->dihedrallist;
+  double **x = atom->x;
+
+  int ndihedrallist = neighbor->ndihedrallist;
+  type = dihedrallist[nd][4];
+
+  vb1x = x[i1][0] - x[i2][0];
+  vb1y = x[i1][1] - x[i2][1];
+  vb1z = x[i1][2] - x[i2][2];
+
+  vb2x = x[i3][0] - x[i2][0];
+  vb2y = x[i3][1] - x[i2][1];
+  vb2z = x[i3][2] - x[i2][2];
+
+  vb2xm = -vb2x;
+  vb2ym = -vb2y;
+  vb2zm = -vb2z;
+
+  vb3x = x[i4][0] - x[i3][0];
+  vb3y = x[i4][1] - x[i3][1];
+  vb3z = x[i4][2] - x[i3][2];
+
+  // c,s calculation
+
+  ax = vb1y*vb2zm - vb1z*vb2ym;
+  ay = vb1z*vb2xm - vb1x*vb2zm;
+  az = vb1x*vb2ym - vb1y*vb2xm;
+  bx = vb3y*vb2zm - vb3z*vb2ym;
+  by = vb3z*vb2xm - vb3x*vb2zm;
+  bz = vb3x*vb2ym - vb3y*vb2xm;
+
+  rasq = ax*ax + ay*ay + az*az;
+  rbsq = bx*bx + by*by + bz*bz;
+
+  ra2inv = rb2inv = 0.0;
+  if (rasq > 0) ra2inv = 1.0/rasq;
+  if (rbsq > 0) rb2inv = 1.0/rbsq;
+  rabinv = sqrt(ra2inv*rb2inv);
+
+  c = (ax*bx + ay*by + az*bz)*rabinv;
+
+  dudih = 0.0;
+  du2dih = 0.0;
+  for (i = 1; i < nterms[type]; i++) {
+    dudih += this->a[type][i]*i*pow(c,i-1);
+  }
+
+  for (i = 2; i < nterms[type]; i++) {
+    du2dih += this->a[type][i]*i*(i-1)*pow(c, i-2);
+  }
+}

src/EXTRA-MOLECULE/dihedral_nharmonic.h

@@ -33,6 +33,7 @@ class DihedralNHarmonic : public Dihedral {
   void write_restart(FILE *) override;
   void read_restart(FILE *) override;
   void write_data(FILE *) override;
+  void born_matrix(int /*dtype*/, int, int, int, int, double &, double &) override;
 
  protected:
   int *nterms;

src/MOLECULE/angle_cosine.cpp

@@ -29,7 +29,10 @@ using MathConst::MY_PI;
 
 /* ---------------------------------------------------------------------- */
 
-AngleCosine::AngleCosine(LAMMPS *_lmp) : Angle(_lmp) {}
+AngleCosine::AngleCosine(LAMMPS *_lmp) : Angle(_lmp)
+{
+  born_matrix_enable = 1;
+}
 
 /* ---------------------------------------------------------------------- */
 
@@ -235,6 +238,14 @@ double AngleCosine::single(int type, int i1, int i2, int i3)
   return k[type] * (1.0 + c);
 }
 
+/* ---------------------------------------------------------------------- */
+
+void AngleCosine::born_matrix(int type, int i1, int i2, int i3, double &du, double &du2)
+{
+  du2 = 0;
+  du = k[type];
+}
+
 /* ----------------------------------------------------------------------
    return ptr to internal members upon request
 ------------------------------------------------------------------------ */

src/MOLECULE/angle_cosine.h

@@ -35,6 +35,7 @@ class AngleCosine : public Angle {
   void read_restart(FILE *) override;
   void write_data(FILE *) override;
   double single(int, int, int, int) override;
+  void born_matrix(int type, int i1, int i2, int i3, double& du, double& du2) override;
   void *extract(const char *, int &) override;
 
  protected:

src/MOLECULE/angle_cosine_squared.cpp

@@ -38,6 +38,7 @@ AngleCosineSquared::AngleCosineSquared(LAMMPS *_lmp) : Angle(_lmp)
 {
   k = nullptr;
   theta0 = nullptr;
+  born_matrix_enable = 1;
 }
 
 /* ---------------------------------------------------------------------- */
@@ -262,3 +263,32 @@ double AngleCosineSquared::single(int type, int i1, int i2, int i3)
   double tk = k[type] * dcostheta;
   return tk * dcostheta;
 }
+
+/* ---------------------------------------------------------------------- */
+
+void AngleCosineSquared::born_matrix(int type, int i1, int i2, int i3, double& du, double& du2)
+{
+  double **x = atom->x;
+
+  double delx1 = x[i1][0] - x[i2][0];
+  double dely1 = x[i1][1] - x[i2][1];
+  double delz1 = x[i1][2] - x[i2][2];
+  domain->minimum_image(delx1,dely1,delz1);
+  double r1 = sqrt(delx1*delx1 + dely1*dely1 + delz1*delz1);
+
+  double delx2 = x[i3][0] - x[i2][0];
+  double dely2 = x[i3][1] - x[i2][1];
+  double delz2 = x[i3][2] - x[i2][2];
+  domain->minimum_image(delx2,dely2,delz2);
+  double r2 = sqrt(delx2*delx2 + dely2*dely2 + delz2*delz2);
+
+  double c = delx1*delx2 + dely1*dely2 + delz1*delz2;
+  c /= r1*r2;
+  if (c > 1.0) c = 1.0;
+  if (c < -1.0) c = -1.0;
+
+  double dcostheta = c - cos(theta0[type]);
+  double tk = k[type] * dcostheta;
+  du2 = 2*k[type];
+  du = du2*dcostheta;
+}

src/MOLECULE/angle_cosine_squared.h

@@ -35,6 +35,7 @@ class AngleCosineSquared : public Angle {
   void read_restart(FILE *) override;
   void write_data(FILE *) override;
   double single(int, int, int, int) override;
+  void born_matrix(int type, int i1, int i2, int i3, double& du, double& du2) override;
 
  protected:
   double *k, *theta0;

src/MOLECULE/bond_harmonic.cpp

@@ -27,7 +27,10 @@ using namespace LAMMPS_NS;
 
 /* ---------------------------------------------------------------------- */
 
-BondHarmonic::BondHarmonic(LAMMPS *_lmp) : Bond(_lmp) {}
+BondHarmonic::BondHarmonic(LAMMPS *_lmp) : Bond(_lmp)
+{
+  born_matrix_enable = 1;
+}
 
 /* ---------------------------------------------------------------------- */
 
@@ -197,6 +200,20 @@ double BondHarmonic::single(int type, double rsq, int /*i*/, int /*j*/, double &
   return rk * dr;
 }
 
+
+/* ---------------------------------------------------------------------- */
+
+void BondHarmonic::born_matrix(int type, double rsq, int /*i*/, int /*j*/,
+                            double &du, double& du2)
+{
+  double r = sqrt(rsq);
+  double dr = r - r0[type];
+  du2 = 0.0;
+  du = 0.0;
+  du2 = 2*k[type];
+  if (r > 0.0) du = du2*dr;
+}
+
 /* ----------------------------------------------------------------------
    return ptr to internal members upon request
 ------------------------------------------------------------------------ */

src/MOLECULE/bond_harmonic.h

@@ -35,6 +35,7 @@ class BondHarmonic : public Bond {
   void read_restart(FILE *) override;
   void write_data(FILE *) override;
   double single(int, double, int, int, double &) override;
+  void born_matrix(int, double, int, int, double &, double &) override;
   void *extract(const char *, int &) override;
 
  protected:

src/angle.cpp

@@ -34,6 +34,7 @@ Angle::Angle(LAMMPS *_lmp) : Pointers(_lmp)
 
   allocated = 0;
   suffix_flag = Suffix::NONE;
+  born_matrix_enable = 0;
 
   maxeatom = maxvatom = maxcvatom = 0;
   eatom = nullptr;

src/angle.h

@@ -25,7 +25,8 @@ class Angle : protected Pointers {
  public:
   int allocated;
   int *setflag;
-  int writedata;             // 1 if writes coeffs to data file
+  int writedata;    // 1 if writes coeffs to data file
+  int born_matrix_enable;
   double energy;             // accumulated energies
   double virial[6];          // accumulated virial: xx,yy,zz,xy,xz,yz
   double *eatom, **vatom;    // accumulated per-atom energy/virial
@@ -59,6 +60,12 @@ class Angle : protected Pointers {
   virtual void read_restart_settings(FILE *){};
   virtual void write_data(FILE *) {}
   virtual double single(int, int, int, int) = 0;
+  virtual void born_matrix(int /*atype*/, int /*at1*/, int /*at2*/, int /*at3*/, double &du,
+                           double &du2)
+  {
+    du = 0.0;
+    du2 = 0.0;
+  }
   virtual double memory_usage();
   virtual void *extract(const char *, int &) { return nullptr; }
   void reinit();

src/bond.cpp

@@ -49,6 +49,7 @@ Bond::Bond(LAMMPS *_lmp) : Pointers(_lmp)
 
   allocated = 0;
   suffix_flag = Suffix::NONE;
+  born_matrix_enable = 0;
   partial_flag = 0;
 
   maxeatom = maxvatom = 0;

src/bond.h

@@ -33,6 +33,8 @@ class Bond : protected Pointers {
   double virial[6];          // accumulated virial: xx,yy,zz,xy,xz,yz
   double *eatom, **vatom;    // accumulated per-atom energy/virial
 
+  int born_matrix_enable;
+
   int comm_forward;        // size of forward communication (0 if none)
   int comm_reverse;        // size of reverse communication (0 if none)
   int comm_reverse_off;    // size of reverse comm even if newton off
@@ -69,6 +71,13 @@ class Bond : protected Pointers {
   virtual int pack_reverse_comm(int, int, double *) { return 0; }
   virtual void unpack_reverse_comm(int, int *, double *) {}
 
+  virtual void born_matrix(int /*btype*/, double /*rsq*/, int /*at1*/, int /*at2*/, double &du,
+                           double &du2)
+  {
+    du = 0.0;
+    du2 = 0.0;
+  }
+
   void write_file(int, char **);
 
  protected:

src/compute.h

@@ -108,7 +108,7 @@ class Compute : protected Pointers {
   Compute(class LAMMPS *, int, char **);
   ~Compute() override;
   void modify_params(int, char **);
-  void reset_extra_dof();
+  virtual void reset_extra_dof();
 
   virtual void init() = 0;
   virtual void init_list(int, class NeighList *) {}

src/compute_temp_profile.cpp

@@ -119,6 +119,9 @@ ComputeTempProfile::ComputeTempProfile(LAMMPS *lmp, int narg, char **arg) :
   nbins = nbinx*nbiny*nbinz;
   if (nbins <= 0) error->all(FLERR,"Illegal compute temp/profile command");
 
+  nstreaming = (xflag==0 ? 0 : 1) + (yflag==0 ? 0 : 1) + (zflag==0 ? 0 : 1);
+  reset_extra_dof();
+
   memory->create(vbin,nbins,ncount,"temp/profile:vbin");
   memory->create(binave,nbins,ncount,"temp/profile:binave");
 
@@ -197,9 +200,10 @@ void ComputeTempProfile::dof_compute()
   natoms_temp = group->count(igroup);
   dof = domain->dimension * natoms_temp;
 
-  // subtract additional d*Nbins DOF, as in Evans and Morriss paper
+  // subtract additional Nbins DOF for each adjusted direction,
+  // as in Evans and Morriss paper
 
-  dof -= extra_dof + fix_dof + domain->dimension*nbins;
+  dof -= extra_dof + fix_dof + nstreaming*nbins;
   if (dof > 0) tfactor = force->mvv2e / (dof * force->boltz);
   else tfactor = 0.0;
 }
@@ -334,14 +338,19 @@ void ComputeTempProfile::compute_array()
 
   MPI_Allreduce(tbin,tbinall,nbins,MPI_DOUBLE,MPI_SUM,world);
 
-  int nper = domain->dimension;
+  double totcount = 0.0;
   for (i = 0; i < nbins; i++) {
     array[i][0] = binave[i][ncount-1];
+    totcount += array[i][0];
+  }
+  double nper = domain->dimension - (extra_dof + fix_dof)/totcount;
+  double dofbin, tfactorbin;
+  for (i = 0; i < nbins; i++) {
     if (array[i][0] > 0.0) {
-      dof = nper*array[i][0] - extra_dof;
-      if (dof > 0) tfactor = force->mvv2e / (dof * force->boltz);
-      else tfactor = 0.0;
-      array[i][1] = tfactor*tbinall[i];
+      dofbin = nper*array[i][0] - nstreaming;
+      if (dofbin > 0) tfactorbin = force->mvv2e / (dofbin * force->boltz);
+      else tfactorbin = 0.0;
+      array[i][1] = tfactorbin*tbinall[i];
     } else array[i][1] = 0.0;
   }
 }
@@ -576,6 +585,12 @@ void ComputeTempProfile::bin_assign()
 
 /* ---------------------------------------------------------------------- */
 
+void ComputeTempProfile::reset_extra_dof() {
+  extra_dof = domain->dimension - nstreaming;
+}
+
+/* ---------------------------------------------------------------------- */
+
 double ComputeTempProfile::memory_usage()
 {
   double bytes = (double)maxatom * sizeof(int);

src/compute_temp_profile.h

@@ -34,6 +34,7 @@ class ComputeTempProfile : public Compute {
   void compute_vector() override;
   void compute_array() override;
 
+  void reset_extra_dof() override;
   void remove_bias(int, double *) override;
   void remove_bias_thr(int, double *, double *) override;
   void remove_bias_all() override;
@@ -47,6 +48,7 @@ class ComputeTempProfile : public Compute {
   int nbinx, nbiny, nbinz, nbins;
   int ivx, ivy, ivz;
   double tfactor;
+  double nstreaming;
 
   int box_change, triclinic;
   int *periodicity;

src/dihedral.cpp

@@ -35,6 +35,7 @@ Dihedral::Dihedral(LAMMPS *_lmp) : Pointers(_lmp)
 
   allocated = 0;
   suffix_flag = Suffix::NONE;
+  born_matrix_enable = 0;
 
   maxeatom = maxvatom = maxcvatom = 0;
   eatom = nullptr;

src/dihedral.h

@@ -25,7 +25,8 @@ class Dihedral : protected Pointers {
  public:
   int allocated;
   int *setflag;
-  int writedata;             // 1 if writes coeffs to data file
+  int writedata;    // 1 if writes coeffs to data file
+  int born_matrix_enable;
   double energy;             // accumulated energy
   double virial[6];          // accumulated virial: xx,yy,zz,xy,xz,yz
   double *eatom, **vatom;    // accumulated per-atom energy/virial
@@ -55,6 +56,12 @@ class Dihedral : protected Pointers {
   virtual void read_restart_settings(FILE *){};
   virtual void write_data(FILE *) {}
   virtual double memory_usage();
+  virtual void born_matrix(int /*dtype*/, int /*at1*/, int /*at2*/, int /*at3*/, int /*at4*/,
+                           double &du, double &du2)
+  {
+    du = 0.0;
+    du2 = 0.0;
+  }
 
  protected:
   int suffix_flag;    // suffix compatibility flag

src/fix_bond_history.cpp

@@ -251,7 +251,7 @@ void FixBondHistory::post_neighbor()
 
 double FixBondHistory::memory_usage()
 {
-  return maxbond * ndata * sizeof(double);
+  return (double) maxbond * ndata * sizeof(double);
 }
 
 /* ---------------------------------------------------------------------- */

src/improper.cpp

@@ -33,6 +33,7 @@ Improper::Improper(LAMMPS *_lmp) : Pointers(_lmp)
 
   allocated = 0;
   suffix_flag = Suffix::NONE;
+  born_matrix_enable = 0;
 
   maxeatom = maxvatom = maxcvatom = 0;
   eatom = nullptr;

src/improper.h

@@ -25,7 +25,8 @@ class Improper : protected Pointers {
  public:
   int allocated;
   int *setflag;
-  int writedata;             // 1 if writes coeffs to data file
+  int writedata;    // 1 if writes coeffs to data file
+  int born_matrix_enable;
   double energy;             // accumulated energies
   double virial[6];          // accumulated virial: xx,yy,zz,xy,xz,yz
   double *eatom, **vatom;    // accumulated per-atom energy/virial
@@ -55,6 +56,12 @@ class Improper : protected Pointers {
   virtual void read_restart_settings(FILE *){};
   virtual void write_data(FILE *) {}
   virtual double memory_usage();
+  virtual void born_matrix(int /*dtype*/, int /*at1*/, int /*at2*/, int /*at3*/, int /*at4*/,
+                           double &du, double &du2)
+  {
+    du = 0.0;
+    du2 = 0.0;
+  }
 
  protected:
   int suffix_flag;    // suffix compatibility flag

src/pair.cpp

@@ -58,6 +58,7 @@ Pair::Pair(LAMMPS *lmp) : Pointers(lmp)
   comm_forward = comm_reverse = comm_reverse_off = 0;
 
   single_enable = 1;
+  born_matrix_enable = 0;
   single_hessian_enable = 0;
   restartinfo = 1;
   respa_enable = 0;

src/pair.h

@@ -52,6 +52,7 @@ class Pair : protected Pointers {
   int comm_reverse_off;    // size of reverse comm even if newton off
 
   int single_enable;              // 1 if single() routine exists
+  int born_matrix_enable;         // 1 if born_matrix() routine exists
   int single_hessian_enable;      // 1 if single_hessian() routine exists
   int restartinfo;                // 1 if pair style writes restart info
   int respa_enable;               // 1 if inner/middle/outer rRESPA routines
@@ -169,6 +170,12 @@ class Pair : protected Pointers {
     return 0.0;
   }
 
+  virtual void born_matrix(int /*i*/, int /*j*/, int /*itype*/, int /*jtype*/, double /*rsq*/,
+                           double /*factor_coul*/, double /*factor_lj*/, double &du, double &du2)
+  {
+    du = du2 = 0.0;
+  }
+
   virtual void settings(int, char **) = 0;
   virtual void coeff(int, char **) = 0;
 

src/pair_hybrid.cpp

@@ -391,6 +391,7 @@ void PairHybrid::flags()
   // single_enable = 1 if all sub-styles are set
   // respa_enable = 1 if all sub-styles are set
   // manybody_flag = 1 if any sub-style is set
+  // born_matrix_enable = 1 if all sub-styles are set
   // no_virial_fdotr_compute = 1 if any sub-style is set
   // ghostneigh = 1 if any sub-style is set
   // ewaldflag, pppmflag, msmflag, dipoleflag, dispersionflag, tip4pflag = 1
@@ -401,11 +402,13 @@ void PairHybrid::flags()
   compute_flag = 0;
   respa_enable = 0;
   restartinfo = 0;
+  born_matrix_enable = 0;
 
   for (m = 0; m < nstyles; m++) {
     if (styles[m]->single_enable) ++single_enable;
     if (styles[m]->respa_enable) ++respa_enable;
     if (styles[m]->restartinfo) ++restartinfo;
+    if (styles[m]->born_matrix_enable) ++born_matrix_enable;
     if (styles[m]->manybody_flag) manybody_flag = 1;
     if (styles[m]->no_virial_fdotr_compute) no_virial_fdotr_compute = 1;
     if (styles[m]->ghostneigh) ghostneigh = 1;
@@ -422,6 +425,7 @@ void PairHybrid::flags()
   single_enable = (single_enable == nstyles) ? 1 : 0;
   respa_enable = (respa_enable == nstyles) ? 1 : 0;
   restartinfo = (restartinfo == nstyles) ? 1 : 0;
+  born_matrix_enable = (born_matrix_enable == nstyles) ? 1 : 0;
   init_svector();
 
   // set centroidstressflag for pair hybrid
@@ -867,6 +871,40 @@ double PairHybrid::single(int i, int j, int itype, int jtype,
   return esum;
 }
 
+/* ----------------------------------------------------------------------
+   call sub-style to compute born matrix interaction
+   error if sub-style does not support born_matrix call
+   since overlay could have multiple sub-styles, sum results explicitly
+------------------------------------------------------------------------- */
+
+void PairHybrid::born_matrix(int i, int j, int itype, int jtype, double rsq,
+                             double factor_coul, double factor_lj,
+                             double &dupair, double &du2pair)
+{
+  if (nmap[itype][jtype] == 0)
+    error->one(FLERR,"Invoked pair born_matrix on pair style none");
+
+  double du, du2;
+  dupair = du2pair = 0.0;
+
+  for (int m = 0; m < nmap[itype][jtype]; m++) {
+    if (rsq < styles[map[itype][jtype][m]]->cutsq[itype][jtype]) {
+      if (styles[map[itype][jtype][m]]->born_matrix_enable == 0)
+        error->one(FLERR,"Pair hybrid sub-style does not support born_matrix call");
+
+      if ((special_lj[map[itype][jtype][m]] != nullptr) ||
+          (special_coul[map[itype][jtype][m]] != nullptr))
+        error->one(FLERR,"Pair hybrid born_matrix calls do not support"
+                   " per sub-style special bond values");
+
+      du = du2 = 0.0;
+      styles[map[itype][jtype][m]]->born_matrix(i,j,itype,jtype,rsq,factor_coul,factor_lj,du,du2);
+      dupair += du;
+      du2pair += du2;
+    }
+  }
+}
+
 /* ----------------------------------------------------------------------
    copy Pair::svector data
 ------------------------------------------------------------------------- */

src/pair_hybrid.h

@@ -49,6 +49,8 @@ class PairHybrid : public Pair {
   void write_restart(FILE *) override;
   void read_restart(FILE *) override;
   double single(int, int, int, int, double, double, double, double &) override;
+  void born_matrix(int, int, int, int, double, double, double, double &, double &) override;
+
   void modify_params(int narg, char **arg) override;
   double memory_usage() override;
 

src/pair_hybrid_scaled.cpp

@@ -413,7 +413,7 @@ double PairHybridScaled::single(int i, int j, int itype, int jtype, double rsq,
           (special_coul[map[itype][jtype][m]] != nullptr))
         error->one(FLERR, "Pair hybrid single() does not support per sub-style special_bond");
 
-      scale = scaleval[map[itype][jtype][m]];
+      double scale = scaleval[map[itype][jtype][m]];
       esum += scale * pstyle->single(i, j, itype, jtype, rsq, factor_coul, factor_lj, fone);
       fforce += scale * fone;
     }
@@ -423,6 +423,57 @@ double PairHybridScaled::single(int i, int j, int itype, int jtype, double rsq,
   return esum;
 }
 
+/* ----------------------------------------------------------------------
+   call sub-style to compute born matrix interaction
+   error if sub-style does not support born_matrix call
+   since overlay could have multiple sub-styles, sum results explicitly
+------------------------------------------------------------------------- */
+
+void PairHybridScaled::born_matrix(int i, int j, int itype, int jtype, double rsq,
+                                   double factor_coul, double factor_lj, double &dupair,
+                                   double &du2pair)
+{
+  if (nmap[itype][jtype] == 0) error->one(FLERR, "Invoked pair born_matrix on pair style none");
+
+  // update scale values from variables where needed
+
+  const int nvars = scalevars.size();
+  if (nvars > 0) {
+    double *vals = new double[nvars];
+    for (int k = 0; k < nvars; ++k) {
+      int m = input->variable->find(scalevars[k].c_str());
+      if (m < 0)
+        error->all(FLERR, "Variable '{}' not found when updating scale factors", scalevars[k]);
+      vals[k] = input->variable->compute_equal(m);
+    }
+    for (int k = 0; k < nstyles; ++k) {
+      if (scaleidx[k] >= 0) scaleval[k] = vals[scaleidx[k]];
+    }
+    delete[] vals;
+  }
+
+  double du, du2, scale;
+  dupair = du2pair = scale = 0.0;
+
+  for (int m = 0; m < nmap[itype][jtype]; m++) {
+    auto pstyle = styles[map[itype][jtype][m]];
+    if (rsq < pstyle->cutsq[itype][jtype]) {
+      if (pstyle->born_matrix_enable == 0)
+        error->one(FLERR, "Pair hybrid sub-style does not support born_matrix call");
+
+      if ((special_lj[map[itype][jtype][m]] != nullptr) ||
+          (special_coul[map[itype][jtype][m]] != nullptr))
+        error->one(FLERR, "Pair hybrid born_matrix() does not support per sub-style special_bond");
+
+      du = du2 = 0.0;
+      scale = scaleval[map[itype][jtype][m]];
+      pstyle->born_matrix(i, j, itype, jtype, rsq, factor_coul, factor_lj, du, du2);
+      dupair += scale * du;
+      du2pair += scale * du2;
+    }
+  }
+}
+
 /* ----------------------------------------------------------------------
    set coeffs for one or more type pairs
 ------------------------------------------------------------------------- */

src/pair_hybrid_scaled.h

@@ -38,6 +38,7 @@ class PairHybridScaled : public PairHybrid {
   void write_restart(FILE *) override;
   void read_restart(FILE *) override;
   double single(int, int, int, int, double, double, double, double &) override;
+  void born_matrix(int, int, int, int, double, double, double, double &, double &) override;
 
   void init_svector() override;
   void copy_svector(int, int) override;

src/pair_lj_cut.cpp

@@ -39,6 +39,7 @@ using namespace MathConst;
 PairLJCut::PairLJCut(LAMMPS *lmp) : Pair(lmp)
 {
   respa_enable = 1;
+  born_matrix_enable = 1;
   writedata = 1;
 }
 
@@ -672,6 +673,28 @@ double PairLJCut::single(int /*i*/, int /*j*/, int itype, int jtype, double rsq,
 
 /* ---------------------------------------------------------------------- */
 
+void PairLJCut::born_matrix(int /*i*/, int /*j*/, int itype, int jtype, double rsq,
+                            double /*factor_coul*/, double factor_lj, double &dupair,
+                            double &du2pair)
+{
+  double rinv, r2inv, r6inv, du, du2;
+
+  r2inv = 1.0 / rsq;
+  rinv = sqrt(r2inv);
+  r6inv = r2inv * r2inv * r2inv;
+
+  // Reminder: lj1 = 48*e*s^12, lj2 = 24*e*s^6
+  // so dupair = -forcelj/r = -fforce*r (forcelj from single method)
+
+  du = r6inv * rinv * (lj2[itype][jtype] - lj1[itype][jtype] * r6inv);
+  du2 = r6inv * r2inv * (13 * lj1[itype][jtype] * r6inv - 7 * lj2[itype][jtype]);
+
+  dupair = factor_lj * du;
+  du2pair = factor_lj * du2;
+}
+
+/* ---------------------------------------------------------------------- */
+
 void *PairLJCut::extract(const char *str, int &dim)
 {
   dim = 2;

src/pair_lj_cut.h

@@ -40,6 +40,7 @@ class PairLJCut : public Pair {
   void write_data(FILE *) override;
   void write_data_all(FILE *) override;
   double single(int, int, int, int, double, double, double, double &) override;
+  void born_matrix(int, int, int, int, double, double, double, double &, double &) override;
   void *extract(const char *, int &) override;
 
   void compute_inner() override;

src/utils.cpp

@@ -120,6 +120,12 @@ std::string utils::strfind(const std::string &text, const std::string &pattern)
     return "";
 }
 
+void utils::missing_cmd_args(const std::string &file, int line, const std::string &cmd,
+                             Error *error)
+{
+  if (error) error->all(file, line, "Illegal {} command: missing argument(s)", cmd);
+}
+
 /* specialization for the case of just a single string argument */
 
 void utils::logmesg(LAMMPS *lmp, const std::string &mesg)
@@ -571,14 +577,14 @@ void utils::bounds(const char *file, int line, const std::string &str,
       error->all(file, line, fmt::format("Invalid range string: {}", str));
 
     if (nlo < nmin)
-      error->all(file, line, fmt::format("Numeric index {} is out of bounds "
-                             "({}-{})", nlo, nmin, nmax));
+      error->all(file, line, fmt::format("Numeric index {} is out of bounds ({}-{})",
+                                         nlo, nmin, nmax));
     else if (nhi > nmax)
-      error->all(file, line, fmt::format("Numeric index {} is out of bounds "
-                             "({}-{})", nhi, nmin, nmax));
+      error->all(file, line, fmt::format("Numeric index {} is out of bounds ({}-{})",
+                                         nhi, nmin, nmax));
     else if (nlo > nhi)
-      error->all(file, line, fmt::format("Numeric index {} is out of bounds "
-                             "({}-{})", nlo, nmin, nhi));
+      error->all(file, line, fmt::format("Numeric index {} is out of bounds ({}-{})",
+                                         nlo, nmin, nhi));
   }
 }
 

src/utils.h

@@ -47,6 +47,18 @@ namespace utils {
 
   std::string strfind(const std::string &text, const std::string &pattern);
 
+  /*! Print error message about missing arguments for command
+   *
+   * This function simplifies the repetitive reporting missing arguments to a command.
+   *
+   *  \param file     name of source file for error message
+   *  \param line     line number in source file for error message
+   *  \param cmd      name of the failing command
+   *  \param error    pointer to Error class instance (for abort) or nullptr */
+
+  [[noreturn]] void missing_cmd_args(const std::string &file, int line, const std::string &cmd,
+                                     Error *error);
+
   /* Internal function handling the argument list for logmesg(). */
 
   void fmtargs_logmesg(LAMMPS *lmp, fmt::string_view format, fmt::format_args args);