diff --git a/doc/src/Commands_pair.rst b/doc/src/Commands_pair.rst
index 2ca0e88729..f5b1ef9b38 100644
--- a/doc/src/Commands_pair.rst
+++ b/doc/src/Commands_pair.rst
@@ -263,6 +263,7 @@ OPT.
    * :doc:`ufm (got) <pair_ufm>`
    * :doc:`vashishta (gko) <pair_vashishta>`
    * :doc:`vashishta/table (o) <pair_vashishta>`
+   * :doc:`wf/cut <pair_wf_cut>`
    * :doc:`yukawa (gko) <pair_yukawa>`
    * :doc:`yukawa/colloid (go) <pair_yukawa_colloid>`
    * :doc:`zbl (gko) <pair_zbl>`
diff --git a/doc/src/JPG/WF_LJ.jpg b/doc/src/JPG/WF_LJ.jpg
new file mode 100644
index 0000000000..fdacb8e80a
Binary files /dev/null and b/doc/src/JPG/WF_LJ.jpg differ
diff --git a/doc/src/pair_style.rst b/doc/src/pair_style.rst
index a5fbb824f5..524581d2c4 100644
--- a/doc/src/pair_style.rst
+++ b/doc/src/pair_style.rst
@@ -327,6 +327,7 @@ accelerated styles exist.
 * :doc:`ufm <pair_ufm>` -
 * :doc:`vashishta <pair_vashishta>` - Vashishta 2-body and 3-body potential
 * :doc:`vashishta/table <pair_vashishta>` -
+* :doc:`wf/cut <pair_wf_cut>` - Wang-Frenkel Potential for short-ranged interactions
 * :doc:`yukawa <pair_yukawa>` - Yukawa potential
 * :doc:`yukawa/colloid <pair_yukawa_colloid>` - screened Yukawa potential for finite-size particles
 * :doc:`zbl <pair_zbl>` - Ziegler-Biersack-Littmark potential
diff --git a/doc/src/pair_wf_cut.rst b/doc/src/pair_wf_cut.rst
new file mode 100644
index 0000000000..e69b982d8a
--- /dev/null
+++ b/doc/src/pair_wf_cut.rst
@@ -0,0 +1,117 @@
+.. index:: pair_style wf/cut
+
+pair_style wf/cut command
+===========================
+
+Syntax
+""""""
+
+
+.. code-block:: LAMMPS
+
+   pair_style wf/cut cutoff
+
+* cutoff = cutoff for wf interactions (distance units)
+
+Examples
+""""""""
+
+
+.. code-block:: LAMMPS
+
+   pair_style        wf/cut 2.0
+   pair_coeff        1 1 1.0 1.0 1 1  2.0
+
+Description
+"""""""""""
+
+The *wf/cut* (Wang-Frenkel) style computes LJ-like potentials as
+described in :ref:`Wang2020 <Wang2020>`. This potential is by
+construction finite ranged and it vanishes quadratically at the cutoff
+distance, avoiding truncation, shifting, interpolation and other typical
+procedures with the LJ potential. The *wf/cut* can be used when a
+typical short-ranged potential with attraction is required. The
+potential is given by which is given by:
+
+.. math::
+  \phi(r)= \epsilon \alpha \left(\left[{\sigma\over r}\right]^{2\mu} -1 \right)\left(\left[{r_c\over r}\right]^{2\mu}-1\right)^{2\nu}
+
+with
+
+.. math::
+  \alpha=2\nu\left(\frac{r_c}{\sigma}\right)^{2\mu}\left[\frac{1+2\nu}{2\nu\left[(r_c/\sigma)^{2\mu}-1\right]}\right]^{2\nu+1}
+
+and
+
+.. math::
+  r_{min}=r_c\left[\frac{1+2\nu}{1+2\nu(r_c/\sigma)^{2\nu}}\right]^{1/{2\nu}}
+
+:math:`r_c` is the cutoff.
+
+Comparison of the non-truncated Lennard-Jones 12-6 potential (red curve),
+and the WF potentials with :math:`\mu=1` and :math:`\nu=1` are shown in
+the figure below. The blue curve has :math:`r_c =2.0` and the green
+curve has :math:`r_c =1.2` and can be used to describe colloidal
+interactions.
+
+.. image:: JPG/WF_LJ.jpg
+   :align: center
+   :scale: 33%
+
+
+The following coefficients must be defined for each pair of atoms
+types via the :doc:`pair_coeff <pair_coeff>` command as in the example
+above, or in the data file or restart files read by the
+:doc:`read_data <read_data>` or :doc:`read_restart <read_restart>`
+commands:
+
+* :math:`\epsilon` (energy units)
+* :math:`\sigma` (distance units)
+* :math:`\nu`
+* :math:`\mu`
+* :math:`r_c` (distance units)
+
+The last coefficient is optional. If not specified, the global cutoff
+given in the pair_style command is used.  The exponents :math:`\nu` and
+:math:`\mu` are positive integers, usually set to 1. There is usually
+little to be gained by choosing other values of :math:`\nu` and
+:math:`\mu` (See discussion in :ref:`Wang2020 <Wang2020>`)
+
+----------
+
+**Mixing, shift, table, tail correction, restart, rRESPA info**\ :
+
+This pair style does not support the :doc:`pair_modify <pair_modify>`
+mixing and table options.
+
+The :doc:`pair_modify <pair_modify>` tail option is not relevant
+for this pair style as it goes to zero at the cut-off radius.
+
+This pair style writes its information to :doc:`binary restart files
+<restart>`, so pair_style and pair_coeff commands do not need to be
+specified in an input script that reads a restart file.
+
+This pair style does not support the use of the *inner*\ , *middle*\ ,
+and *outer* keywords of the :doc:`run_style respa <run_style>` command.
+
+----------
+
+Restrictions
+""""""""""""
+This pair style can only be used if LAMMPS was built with the
+USER-MISC package.  See the :doc:`Build package <Build_package>` doc
+page for more info.
+
+Related commands
+""""""""""""""""
+
+:doc:`pair_coeff <pair_coeff>`
+
+**Default:** none
+
+
+----------
+
+.. _Wang2020:
+
+**(Wang2020)** X. Wang, S. Ramirez-Hinestrosa, J. Dobnikar, and D. Frenkel, Phys. Chem. Chem. Phys. 22, 10624 (2020).
diff --git a/doc/utils/sphinx-config/false_positives.txt b/doc/utils/sphinx-config/false_positives.txt
index 7e55fe3638..21582f11f1 100644
--- a/doc/utils/sphinx-config/false_positives.txt
+++ b/doc/utils/sphinx-config/false_positives.txt
@@ -692,6 +692,7 @@ dmi
 dnf
 DNi
 Dobson
+Dobnikar
 Dodds
 docenv
 dodgerblue
@@ -1232,6 +1233,7 @@ Hibbs
 Higdon
 Hijazi
 Hilger
+Hinestrosa
 histo
 histogrammed
 histogramming
@@ -3406,6 +3408,7 @@ WeinanE
 Wennberg
 Westmere
 Westview
+wf
 wget
 Whelan
 whitesmoke
diff --git a/src/USER-MISC/README b/src/USER-MISC/README
index 538fdb7952..314fe6146e 100644
--- a/src/USER-MISC/README
+++ b/src/USER-MISC/README
@@ -107,5 +107,6 @@ pair_style morse/smooth/linear, Stefan Paquay (TU Eindhoven), stefanpaquay at gm
 pair_style srp, Tim Sirk, tim.sirk at us.army.mil, 21 Nov 14
 pair_style tersoff/table, Luca Ferraro, luca.ferraro@caspur.it, 1 Dec 11
 pair_style momb, Kristen Fichthorn, Tonnam Balankura, Ya Zhou, fichthorn@psu.edu, 18 Mar 17
+pair_style wf_cut, Simon Ramirez-Hinestrosa, Xipeng Wang, sr802 at cam.ac.uk, 4 Nov 20
 temper/grem, David Stelter, dstelter@bu.edu, 22 Nov 16
 temper/npt, Amulya K. Pervaje, Cody K. Addington, amulyapervaje@gmail.com , 31 Aug 17
diff --git a/src/USER-MISC/pair_wf_cut.cpp b/src/USER-MISC/pair_wf_cut.cpp
new file mode 100644
index 0000000000..d29ae2a4fe
--- /dev/null
+++ b/src/USER-MISC/pair_wf_cut.cpp
@@ -0,0 +1,397 @@
+/* ----------------------------------------------------------------------
+   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+   http://lammps.sandia.gov, Sandia National Laboratories
+   Steve Plimpton, sjplimp@sandia.gov
+
+   Copyright (2003) Sandia Corporation.  Under the terms of Contract
+   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+   certain rights in this software.  This software is distributed under
+   the GNU General Public License.
+
+   See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+/* ----------------------------------------------------------------------
+   Contributing Author: Xipeng Wang, Simon Ramirez-Hinestrosa
+------------------------------------------------------------------------- */
+#include "pair_wf_cut.h"
+
+#include "atom.h"
+#include "comm.h"
+#include "error.h"
+#include "force.h"
+#include "math_const.h"
+#include "math_special.h"
+#include "memory.h"
+#include "neigh_list.h"
+
+#include <cmath>
+#include <cstdlib>
+#include <cstring>
+
+using namespace LAMMPS_NS;
+using namespace MathConst;
+using namespace MathSpecial;
+
+/* ---------------------------------------------------------------------- */
+
+PairWFCut::PairWFCut(LAMMPS *lmp) : Pair(lmp)
+{
+  writedata = 1;
+}
+
+/* ---------------------------------------------------------------------- */
+
+PairWFCut::~PairWFCut()
+{
+  if (allocated) {
+    memory->destroy(setflag);
+    memory->destroy(cutsq);
+    memory->destroy(cut);
+    memory->destroy(epsilon);
+    memory->destroy(sigma);
+    memory->destroy(nu);
+    memory->destroy(mu);
+    memory->destroy(nm);
+    memory->destroy(e0nm);  //Alpha * epsilon
+    memory->destroy(rcmu);
+    memory->destroy(sigma_mu);
+    memory->destroy(offset);
+  }
+}
+
+/* ---------------------------------------------------------------------- */
+
+void PairWFCut::compute(int eflag, int vflag)
+{
+  int i,j,ii,jj,inum,jnum,itype,jtype;
+  double xtmp,ytmp,ztmp,delx,dely,delz,evdwl,fpair;
+  double rsq,r2inv,factor_lj;
+  double forcenm,rminv, rm, rn;
+  int *ilist,*jlist,*numneigh,**firstneigh;
+
+  evdwl = 0.0;
+  if (eflag || vflag) ev_setup(eflag,vflag);
+  else evflag = vflag_fdotr = 0;
+
+  double **x = atom->x;
+  double **f = atom->f;
+  int *type = atom->type;
+  int nlocal = atom->nlocal;
+  double *special_lj = force->special_lj;
+  int newton_pair = force->newton_pair;
+
+  inum = list->inum;
+  ilist = list->ilist;
+  numneigh = list->numneigh;
+  firstneigh = list->firstneigh;
+
+  // loop over neighbors of my atoms
+
+  for (ii = 0; ii < inum; ii++) {
+    i = ilist[ii];
+    xtmp = x[i][0];
+    ytmp = x[i][1];
+    ztmp = x[i][2];
+    itype = type[i];
+    jlist = firstneigh[i];
+    jnum = numneigh[i];
+
+    for (jj = 0; jj < jnum; jj++) {
+      j = jlist[jj];
+      factor_lj = special_lj[sbmask(j)];
+      j &= NEIGHMASK;
+
+      delx = xtmp - x[j][0];
+      dely = ytmp - x[j][1];
+      delz = ztmp - x[j][2];
+      rsq = delx*delx + dely*dely + delz*delz;
+      jtype = type[j];
+
+      if (rsq < cutsq[itype][jtype]) {
+        r2inv = 1.0/rsq;
+        rminv = powint(r2inv,mu[itype][jtype]);
+        rm = sigma_mu[itype][jtype]*rminv - 1.0;
+        rn = rcmu[itype][jtype]*rminv - 1.0;
+
+        forcenm = 2.0*mu[itype][jtype] *sigma_mu[itype][jtype]*powint(rn,2*nu[itype][jtype])
+                + 4.0*nm[itype][jtype] *rcmu[itype][jtype]*rm*powint(rn,2*nu[itype][jtype]-1);
+        fpair = factor_lj*e0nm[itype][jtype]*forcenm*powint(r2inv,mu[itype][jtype]+1);
+
+        f[i][0] += delx*fpair;
+        f[i][1] += dely*fpair;
+        f[i][2] += delz*fpair;
+        if (newton_pair || j < nlocal) {
+          f[j][0] -= delx*fpair;
+          f[j][1] -= dely*fpair;
+          f[j][2] -= delz*fpair;
+        }
+
+        if (eflag) {
+          evdwl = e0nm[itype][jtype] *
+            (rm*powint(rn,2*nu[itype][jtype])) - offset[itype][jtype];
+          evdwl *= factor_lj;
+        }
+
+        if (evflag) ev_tally(i,j,nlocal,newton_pair,
+                             evdwl,0.0,fpair,delx,dely,delz);
+      }
+    }
+  }
+
+  if (vflag_fdotr) virial_fdotr_compute();
+}
+
+/* ----------------------------------------------------------------------
+   allocate all arrays
+------------------------------------------------------------------------- */
+
+void PairWFCut::allocate()
+{
+  allocated = 1;
+  int n = atom->ntypes;
+
+  memory->create(setflag,n+1,n+1,"pair:setflag");
+  for (int i = 1; i <= n; i++)
+    for (int j = i; j <= n; j++)
+      setflag[i][j] = 0;
+
+  memory->create(cutsq,n+1,n+1,"pair:cutsq");
+
+  memory->create(cut,n+1,n+1,"pair:cut");
+  memory->create(epsilon,n+1,n+1,"pair:epsilon");
+  memory->create(sigma,n+1,n+1,"pair:sigma");
+  memory->create(nu,n+1,n+1,"pair:nu");
+  memory->create(mu,n+1,n+1,"pair:mu");
+  memory->create(nm,n+1,n+1,"pair:nm");
+  memory->create(e0nm,n+1,n+1,"pair:e0nm");
+  memory->create(rcmu,n+1,n+1,"pair:rcmu");
+  memory->create(sigma_mu,n+1,n+1,"pair:sigma_mu");
+  memory->create(offset,n+1,n+1,"pair:offset");
+}
+
+/* ----------------------------------------------------------------------
+   global settings
+------------------------------------------------------------------------- */
+
+void PairWFCut::settings(int narg, char **arg)
+{
+  if (narg != 1) error->all(FLERR,"Illegal pair_style command");
+
+  cut_global = utils::numeric(FLERR,arg[0],false,lmp);
+
+  // reset cutoffs that have been explicitly set
+
+  if (allocated) {
+    int i,j;
+    for (i = 1; i <= atom->ntypes; i++)
+      for (j = i; j <= atom->ntypes; j++)
+        if (setflag[i][j]) cut[i][j] = cut_global;
+  }
+}
+
+/* ----------------------------------------------------------------------
+   set coeffs for one or more type pairs
+------------------------------------------------------------------------- */
+
+void PairWFCut::coeff(int narg, char **arg)
+{
+  if (narg < 6 || narg > 7)
+    error->all(FLERR,"Incorrect args for pair coefficients");
+  if (!allocated) allocate();
+
+  int ilo,ihi,jlo,jhi;
+  utils::bounds(FLERR,arg[0],1,atom->ntypes,ilo,ihi,error);
+  utils::bounds(FLERR,arg[1],1,atom->ntypes,jlo,jhi,error);
+
+  double epsilon_one = utils::numeric(FLERR,arg[2],false,lmp);
+  double sigma_one = utils::numeric(FLERR,arg[3],false,lmp);
+  int nu_one = utils::inumeric(FLERR,arg[4],false,lmp);
+  int mu_one = utils::inumeric(FLERR,arg[5],false,lmp);
+
+  double cut_one = cut_global;
+  if (narg == 7) cut_one = utils::numeric(FLERR,arg[6],false,lmp);
+
+  int count = 0;
+  for (int i = ilo; i <= ihi; i++) {
+    for (int j = MAX(jlo,i); j <= jhi; j++) {
+      epsilon[i][j] = epsilon_one;
+      sigma[i][j] = sigma_one;
+      nu[i][j] = nu_one;
+      mu[i][j] = mu_one;
+      cut[i][j] = cut_one;
+      setflag[i][j] = 1;
+      count++;
+    }
+  }
+
+  if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients");
+}
+
+/* ----------------------------------------------------------------------
+   init for one type pair i,j and corresponding j,i
+------------------------------------------------------------------------- */
+
+double PairWFCut::init_one(int i, int j)
+{
+  if (setflag[i][j] == 0) error->all(FLERR,"All pair coeffs are not set");
+
+  nm[i][j] = nu[i][j]*mu[i][j];
+  e0nm[i][j] = epsilon[i][j]*2.0*nu[i][j]*powint(cut[i][j]/sigma[i][j],2*mu[i][j])
+                       *powint((1+2.0*nu[i][j])/(2.0*nu[i][j])/(MathSpecial::powint(cut[i][j]/sigma[i][j],2*mu[i][j])-1.0),
+                              2*nu[i][j]+1);
+  rcmu[i][j] = powint(cut[i][j],2*mu[i][j]);
+  sigma_mu[i][j] = powint(sigma[i][j], 2*mu[i][j]);
+
+  if (offset_flag && (cut[i][j] > 0.0)) {
+    offset[i][j] = 0.0;
+  } else offset[i][j] = 0.0;
+
+  epsilon[j][i] = epsilon[i][j];
+  nu[j][i] = nu[i][j];
+  mu[j][i] = mu[i][j];
+  nm[j][i] = nm[i][j];
+  sigma[j][i] = sigma[i][j];
+  e0nm[j][i] = e0nm[i][j];
+  rcmu[j][i] = rcmu[i][j];
+  sigma_mu[j][i] = sigma_mu[i][j];
+  offset[j][i] = offset[i][j];
+
+  return cut[i][j];
+}
+
+/* ----------------------------------------------------------------------
+  proc 0 writes to restart file
+------------------------------------------------------------------------- */
+
+void PairWFCut::write_restart(FILE *fp)
+{
+  write_restart_settings(fp);
+
+  int i,j;
+  for (i = 1; i <= atom->ntypes; i++)
+    for (j = i; j <= atom->ntypes; j++) {
+      fwrite(&setflag[i][j],sizeof(int),1,fp);
+      if (setflag[i][j]) {
+        fwrite(&epsilon[i][j],sizeof(double),1,fp);
+        fwrite(&sigma[i][j],sizeof(double),1,fp);
+        fwrite(&nu[i][j],sizeof(int),1,fp);
+        fwrite(&mu[i][j],sizeof(int),1,fp);
+        fwrite(&cut[i][j],sizeof(double),1,fp);
+      }
+    }
+}
+
+/* ----------------------------------------------------------------------
+  proc 0 reads from restart file, bcasts
+------------------------------------------------------------------------- */
+
+void PairWFCut::read_restart(FILE *fp)
+{
+  read_restart_settings(fp);
+  allocate();
+
+  int i,j;
+  int me = comm->me;
+  for (i = 1; i <= atom->ntypes; i++)
+    for (j = i; j <= atom->ntypes; j++) {
+      if (me == 0) fread(&setflag[i][j],sizeof(int),1,fp);
+      MPI_Bcast(&setflag[i][j],1,MPI_INT,0,world);
+      if (setflag[i][j]) {
+        if (me == 0) {
+          fread(&epsilon[i][j],sizeof(double),1,fp);
+          fread(&sigma[i][j],sizeof(double),1,fp);
+          fread(&nu[i][j],sizeof(int),1,fp);
+          fread(&mu[i][j],sizeof(int),1,fp);
+          fread(&cut[i][j],sizeof(double),1,fp);
+        }
+        MPI_Bcast(&epsilon[i][j],1,MPI_DOUBLE,0,world);
+        MPI_Bcast(&sigma[i][j],1,MPI_DOUBLE,0,world);
+        MPI_Bcast(&nu[i][j],1,MPI_INT,0,world);
+        MPI_Bcast(&mu[i][j],1,MPI_INT,0,world);
+        MPI_Bcast(&cut[i][j],1,MPI_DOUBLE,0,world);
+      }
+    }
+}
+
+/* ----------------------------------------------------------------------
+  proc 0 writes to restart file
+------------------------------------------------------------------------- */
+
+void PairWFCut::write_restart_settings(FILE *fp)
+{
+  fwrite(&cut_global,sizeof(double),1,fp);
+  fwrite(&offset_flag,sizeof(int),1,fp);
+  fwrite(&mix_flag,sizeof(int),1,fp);
+}
+
+/* ----------------------------------------------------------------------
+  proc 0 reads from restart file, bcasts
+------------------------------------------------------------------------- */
+
+void PairWFCut::read_restart_settings(FILE *fp)
+{
+  if (comm->me == 0) {
+    fread(&cut_global,sizeof(double),1,fp);
+    fread(&offset_flag,sizeof(int),1,fp);
+    fread(&mix_flag,sizeof(int),1,fp);
+  }
+  MPI_Bcast(&cut_global,1,MPI_DOUBLE,0,world);
+  MPI_Bcast(&offset_flag,1,MPI_INT,0,world);
+  MPI_Bcast(&mix_flag,1,MPI_INT,0,world);
+}
+
+/* ----------------------------------------------------------------------
+   proc 0 writes to data file
+------------------------------------------------------------------------- */
+
+void PairWFCut::write_data(FILE *fp)
+{
+  for (int i = 1; i <= atom->ntypes; i++)
+    fprintf(fp,"%d %g %g %d %d\n",i,epsilon[i][i],sigma[i][i],nu[i][i],mu[i][i]);
+}
+
+/* ----------------------------------------------------------------------
+   proc 0 writes all pairs to data file
+------------------------------------------------------------------------- */
+
+void PairWFCut::write_data_all(FILE *fp)
+{
+  for (int i = 1; i <= atom->ntypes; i++)
+    for (int j = i; j <= atom->ntypes; j++)
+      fprintf(fp,"%d %d %g %g %d %d %g\n",i,j,
+              epsilon[i][j],sigma[i][j],nu[i][j],mu[i][j],cut[i][j]);
+}
+
+/* ---------------------------------------------------------------------- */
+
+double PairWFCut::single(int /*i*/, int /*j*/, int itype, int jtype,
+                      double rsq, double /*factor_coul*/, double factor_lj,
+                      double &fforce)
+{
+  double r2inv,rminv,rm,rn,forcenm,phinm;
+
+  r2inv = 1.0/rsq;
+  rminv =powint(r2inv,mu[itype][jtype]);
+  rm = sigma_mu[itype][jtype]*rminv - 1.0;
+  rn = rcmu[itype][jtype]*rminv - 1.0;
+  forcenm = 2.0*mu[itype][jtype] *sigma_mu[itype][jtype]*powint(rn,2*nu[itype][jtype])
+                + 4.0*nm[itype][jtype] *rcmu[itype][jtype]*rm*powint(rn,2*nu[itype][jtype]-1);
+  fforce = factor_lj*e0nm[itype][jtype]*forcenm*powint(r2inv,mu[itype][jtype]+1);
+
+  phinm = e0nm[itype][jtype] * rm*powint(rn,2*nu[itype][jtype]) -
+    offset[itype][jtype];
+  return factor_lj*phinm;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void *PairWFCut::extract(const char *str, int &dim)
+{
+  dim = 2;
+  if (strcmp(str,"epsilon") == 0) return (void *) epsilon;
+  if (strcmp(str,"sigma") == 0) return (void *) sigma;
+  if (strcmp(str,"nu") == 0) return (void *) nu;
+  if (strcmp(str,"mu") == 0) return (void *) mu;
+  return NULL;
+}
diff --git a/src/USER-MISC/pair_wf_cut.h b/src/USER-MISC/pair_wf_cut.h
new file mode 100644
index 0000000000..d69efdf006
--- /dev/null
+++ b/src/USER-MISC/pair_wf_cut.h
@@ -0,0 +1,77 @@
+/* -*- 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 PAIR_CLASS
+
+PairStyle(wf/cut,PairWFCut)
+
+#else
+
+#ifndef LMP_PAIR_WF_CUT_H
+#define LMP_PAIR_WF_CUT_H
+
+#include "pair.h"
+
+namespace LAMMPS_NS {
+
+class PairWFCut : public Pair {
+ public:
+  PairWFCut(class LAMMPS *);
+  virtual ~PairWFCut();
+
+  virtual void compute(int, int);
+  void settings(int, char **);
+  void coeff(int, char **);
+  double init_one(int, int);
+  void write_restart(FILE *);
+  void read_restart(FILE *);
+  void write_restart_settings(FILE *);
+  void read_restart_settings(FILE *);
+  void write_data(FILE *);
+  void write_data_all(FILE *);
+  double single(int, int, int, int, double, double, double, double &);
+  void *extract(const char *, int &);
+
+ protected:
+  int    **nu,**mu;
+  double cut_global;
+  double **cut;
+  double **epsilon,**sigma;
+  double **nm,**e0nm,**rcmu,**sigma_mu,**offset;
+
+  void allocate();
+};
+
+}
+
+#endif
+#endif
+
+/* ERROR/WARNING messages:
+
+E: Illegal ... command
+
+Self-explanatory.  Check the input script syntax and compare to the
+documentation for the command.  You can use -echo screen as a
+command-line option when running LAMMPS to see the offending line.
+
+E: Incorrect args for pair coefficients
+
+Self-explanatory.  Check the input script or data file.
+
+E: All pair coeffs are not set
+
+All pair coefficients must be set in the data file or by the
+pair_coeff command before running a simulation.
+
+*/
diff --git a/unittest/force-styles/tests/mol-pair-wf_cut.yaml b/unittest/force-styles/tests/mol-pair-wf_cut.yaml
new file mode 100644
index 0000000000..359c6d7006
--- /dev/null
+++ b/unittest/force-styles/tests/mol-pair-wf_cut.yaml
@@ -0,0 +1,103 @@
+---
+lammps_version: 24 Dec 2020
+date_generated: Thu Jan 28 03:47:32 202
+epsilon: 5e-14
+prerequisites: ! |
+  atom full
+  pair wf/cut
+pre_commands: ! |
+  variable write_data_pair index ij
+post_commands: ! ""
+input_file: in.fourmol
+pair_style: wf/cut 8.0
+pair_coeff: ! |
+  1 1 0.02 2.5 1 1
+  1 2 0.01 1.75 1 1
+  1 3 0.02 2.85 1 1 8.0
+  1 4 0.0173205 2.8 1 1
+  1 5 0.0173205 2.8 1 1
+  2 2 0.005 1.0 1 1
+  2 3 0.01 2.1 1 1
+  2 4 0.005 0.5 1 1
+  2 5 0.00866025 2.05 1 1
+  3 3 0.02 3.2 1 1
+  3 4 0.0173205 3.15 1 1
+  3 5 0.0173205 3.15 1 1
+  4 4 0.015 3.1 1 1
+  4 5 0.015 3.1 1 1
+  5 5 0.015 3.1 1 1
+extract: ! |
+  epsilon 2
+  sigma 2
+  mu 2
+  nu 2
+natoms: 29
+init_vdwl: 13.331701868543
+init_coul: 0
+init_stress: ! |2-
+   2.5178354370947517e+01  2.5642248474815997e+01  4.2120178706774510e+01 -6.3752199820262074e+00  2.4423821790103029e+00  4.0948131051558949e+00
+init_forces: ! |2
+    1  2.0952781234870946e-02  3.8845446670745600e+00  4.0266101185891161e+00
+    2  1.7887164462329612e+00  1.4609613965099468e+00 -2.0999998459119902e+00
+    3 -1.8691796751627372e+00 -4.4796677668687161e+00 -1.7170104063474256e+00
+    4 -3.8589820127510022e-01  8.0372789351205343e-02 -2.7929310936892798e-01
+    5 -1.2955829232072852e-01 -2.4706552738630230e-01  5.3701075594699055e-01
+    6 -6.6148738809131027e+00  7.1811520725339770e+00  6.3379187807855732e+00
+    7 -1.7981150332175785e-01 -1.9777540538053207e+00 -1.2725481634272681e+01
+    8  3.7391579476987608e-01 -4.5422703098451167e-01  4.3004390815756706e+00
+    9  7.8558519258875970e-01  8.6126695536264386e-01  3.4956293730696517e+00
+   10  4.4874125886822629e+00 -5.8119138927398062e+00 -1.7734559280979514e+00
+   11 -9.9403778094523679e-02 -2.7520368302139697e-01 -4.3804356637156183e-01
+   12  1.0868312455520435e+00  4.0910511978376651e-01 -4.5424665758568461e-01
+   13  3.8952669951265123e-01 -1.6043031845021499e-01 -5.1542590980086338e-03
+   14 -1.6611954727058934e-01  3.3384477515331905e-02 -4.1278689987033085e-01
+   15 -4.5653517060236324e-03  4.1957529431250451e-01  1.4273175604750324e-01
+   16  3.8705064540182841e+00 -3.0790107392412147e+00 -8.4106346959871434e+00
+   17 -3.3157104185347768e+00  2.1771887621745951e+00  9.4592735301519415e+00
+   18 -1.5492810751018058e-02 -3.2610168906092168e-02  3.2838321862976375e-02
+   19 -1.0221075635570953e-04 -2.4819963720238802e-03  6.3588088449151167e-03
+   20 -2.5170248769922246e-03 -3.1219689267982533e-03  2.0971562342533825e-03
+   21 -8.2441875163539302e-01 -9.3159913536522632e-01  2.7408941662865209e+00
+   22 -1.3111629321445617e+00 -3.1620494060282150e-01 -2.0559535894812671e+00
+   23  2.1151339870601102e+00  1.2713008492411884e+00 -6.6308294784684352e-01
+   24  5.0599831218561930e-01 -2.5174724379982512e+00  1.3636099540694366e+00
+   25 -1.8335394725098328e+00  2.9153335061855473e-01 -1.5462990744267575e+00
+   26  1.3115091747032070e+00  2.2039747754639123e+00  1.4641348730469819e-01
+   27  5.1192148870411813e-01 -2.7035848557139470e+00  1.0102720504380645e+00
+   28 -2.0671223896069999e+00  8.9019054155521471e-01 -1.3998051355791850e+00
+   29  1.5714660756357337e+00  1.8277974648852422e+00  3.7915040903844482e-01
+run_vdwl: 13.3325799070269
+run_coul: 0
+run_stress: ! |2-
+   2.5226505134166011e+01  2.5665098713276272e+01  4.2052612927488731e+01 -6.3480044192300973e+00  2.4548005748737984e+00  4.1205804855222423e+00
+run_forces: ! |2
+    1  2.5680403549306934e-02  3.8827166426800459e+00  4.0194469167946600e+00
+    2  1.7880086746222570e+00  1.4633410618480136e+00 -2.0917783930712170e+00
+    3 -1.8746243827554767e+00 -4.4803529797494974e+00 -1.7183386746157114e+00
+    4 -3.8419882068706390e-01  7.9679467451782574e-02 -2.7873843176033536e-01
+    5 -1.2933215119796876e-01 -2.4678741398042314e-01  5.3604662977634610e-01
+    6 -6.5911667151219158e+00  7.1628788517392143e+00  6.2986239762979812e+00
+    7 -1.8049056720639833e-01 -1.9732770892371352e+00 -1.2664364406289327e+01
+    8  3.4117839418858381e-01 -4.2686258225758750e-01  4.2905525151899671e+00
+    9  7.8492784707757846e-01  8.5653579701851090e-01  3.4886297945287605e+00
+   10  4.4990483538894663e+00 -5.8243793398035351e+00 -1.7801144827973929e+00
+   11 -9.8434450383760988e-02 -2.7265166825530213e-01 -4.3384298520494036e-01
+   12  1.0862014078554925e+00  4.0749467141919393e-01 -4.5810984874699312e-01
+   13  3.8840851136656979e-01 -1.5918511941530258e-01 -5.0435188450808476e-03
+   14 -1.6498465402555887e-01  3.2555700014858351e-02 -4.0969046354316158e-01
+   15 -5.1746183244565733e-03  4.2030050245425560e-01  1.4369692057790906e-01
+   16  3.8716092530809463e+00 -3.0834565917976824e+00 -8.4237297910920166e+00
+   17 -3.3183293718136500e+00  2.1837078960805343e+00  9.4702404317304740e+00
+   18 -1.5490237139094618e-02 -3.2605923516534234e-02  3.2842553892020329e-02
+   19 -1.1022993717618882e-04 -2.4919998809809947e-03  6.3698702514743471e-03
+   20 -2.5070784591324888e-03 -3.1118152231195575e-03  2.0924415709628569e-03
+   21 -8.3146680786849725e-01 -9.3100084977990383e-01  2.7490008132217891e+00
+   22 -1.3200074013367948e+00 -3.2075876797175989e-01 -2.0624746416557977e+00
+   23  2.1310227732736915e+00  1.2752677600846360e+00 -6.6466450635212104e-01
+   24  5.2174019224112633e-01 -2.5432192692104154e+00  1.3836830641468214e+00
+   25 -1.8640901009223221e+00  2.9226034502451703e-01 -1.5724314889166486e+00
+   26  1.3263176769235392e+00  2.2290027561589696e+00  1.5247509279036245e-01
+   27  5.2006996521896620e-01 -2.7148614233414463e+00  1.0110231119667439e+00
+   28 -2.0796779375288135e+00  8.9431296914727643e-01 -1.4054509798002364e+00
+   29  1.5758720714205563e+00  1.8349484122988167e+00  3.8404847995470470e-01
+...