diff --git a/doc/src/Commands_pair.rst b/doc/src/Commands_pair.rst index e659a429f0..5328d1163f 100644 --- a/doc/src/Commands_pair.rst +++ b/doc/src/Commands_pair.rst @@ -185,6 +185,7 @@ OPT. * :doc:`lj/long/tip4p/long (o) ` * :doc:`lj/mdf ` * :doc:`lj/relres (o) ` + * :doc:`lj/sphere ` * :doc:`lj/spica (gko) ` * :doc:`lj/spica/coul/long (go) ` * :doc:`lj/spica/coul/msm (o) ` diff --git a/doc/src/pair_lj_sphere.rst b/doc/src/pair_lj_sphere.rst new file mode 100644 index 0000000000..6f692400da --- /dev/null +++ b/doc/src/pair_lj_sphere.rst @@ -0,0 +1,117 @@ +.. index:: pair_style lj/sphere + +pair_style lj/sphere command +============================ + +Syntax +"""""" + +.. code-block:: LAMMPS + + pair_style style args + +* style = *lj/sphere* +* args = list of arguments for a particular style + +.. parsed-literal:: + + *lj/sphere* args = cutoff + cutoff = global cutoff for Lennard Jones interactions (distance units) + +Examples +"""""""" + +.. code-block:: LAMMPS + + pair_style lj/sphere 2.5 + pair_coeff * * 1.0 + pair_coeff 1 1 1.1 2.8 + +Description +""""""""""" + +The *lj/sphere* styles compute the standard 12/6 Lennard-Jones potential, +given by + +.. math:: + + E = 4 \epsilon \left[ \left(\frac{\sigma_{ij}}{r}\right)^{12} - + \left(\frac{\sigma_{ij}}{r}\right)^6 \right] + \qquad r < r_c + +:math:`r_c` is the cutoff. + +This is the same potential function as used by the :doc:`lj/cut +` pair style, but the :math:`\sigma_{ij}` parameter is not set +as a per-type parameter via the :doc:`pair_coeff command `, +but taken from the per-atom radius attribute of :doc:`atom_style sphere +`. The individual value of :math:`\sigma_{ij}` is computed +using the mixing rule for pair coefficients as set by the +:doc:`pair_modify mix ` command. + +Note that :math:`\sigma_{ij}` is defined in the LJ formula above as the +zero-crossing distance for the potential, *not* as the energy minimum which +is at :math:`2^{\frac{1}{6}} \sigma_{ij}`. + +Coefficients +"""""""""""" + +The following coefficients must be defined for each pair of atoms types via the +:doc:`pair_coeff ` command as in the examples above, or in the data +file or restart files read by the :doc:`read_data ` or +:doc:`read_restart ` commands, or by mixing as described below: + +* :math:`\epsilon` (energy units) +* LJ cutoff (distance units) (optional) + +The last coefficient is optional. If not specified, the global +LJ cutoff specified in the pair_style command is used. + +---------- + +Mixing, shift, table, tail correction, restart, rRESPA info +""""""""""""""""""""""""""""""""""""""""""""""""""""""""""" + +For atom type pairs I,J and I != J, the epsilon coefficients and cutoff +distance for the *lj/sphere* pair style can be mixed. The default mix value +is *geometric*. See the "pair_modify" command for details. + +The *lj/sphere* pair style supports the :doc:`pair_modify ` +shift option for the energy of the Lennard-Jones portion of the pair +interaction. + +The *lj/sphere* pair style does *not* support the :doc:`pair_modify +` tail option for adding a long-range tail corrections to +the energy and pressure. + +The *lj/sphere* pair style writes its information to :doc:`binary +restart files `, 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 can only be used via the *pair* keyword of the +:doc:`run_style respa ` command. It does *not* support the +*inner*, *middle*, *outer* keywords. + +---------- + +Restrictions +"""""""""""" + +The *lj/sphere* pair style is only enabled if LAMMPS was built with the +EXTRA-PAIR package. See the :doc:`Build package ` page +for more info. + +The *lj/sphere* pair style does not support the *sixthpower* mixing rule. + +---------- + +Related commands +"""""""""""""""" + +* :doc:`pair_coeff ` +* :doc:`pair_style lj/cut ` + +Default +""""""" + +none diff --git a/doc/src/pair_style.rst b/doc/src/pair_style.rst index a1c1dee33e..a0fb6eeefd 100644 --- a/doc/src/pair_style.rst +++ b/doc/src/pair_style.rst @@ -263,6 +263,7 @@ accelerated styles exist. * :doc:`lj/long/tip4p/long ` - long-range LJ and long-range Coulomb for TIP4P water * :doc:`lj/mdf ` - LJ potential with a taper function * :doc:`lj/relres ` - LJ using multiscale Relative Resolution (RelRes) methodology :ref:`(Chaimovich) `. +* :doc:`lj/sphere ` - LJ where per-atom radius is used as LJ sigma * :doc:`lj/spica ` - LJ for SPICA coarse-graining * :doc:`lj/spica/coul/long ` - LJ for SPICA coarse-graining with long-range Coulomb * :doc:`lj/spica/coul/msm ` - LJ for SPICA coarse-graining with long-range Coulomb via MSM diff --git a/src/EXTRA-PAIR/pair_lj_sphere.cpp b/src/EXTRA-PAIR/pair_lj_sphere.cpp new file mode 100644 index 0000000000..c8c0211d50 --- /dev/null +++ b/src/EXTRA-PAIR/pair_lj_sphere.cpp @@ -0,0 +1,354 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + https://www.lammps.org/, Sandia National Laboratories + LAMMPS development team: developers@lammps.org + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#include "pair_lj_sphere.h" + +#include "atom.h" +#include "comm.h" +#include "error.h" +#include "force.h" +#include "math_special.h" +#include "memory.h" +#include "neigh_list.h" +#include "neighbor.h" +#include "update.h" + +#include +#include + +using namespace LAMMPS_NS; +using MathSpecial::powint; + +/* ---------------------------------------------------------------------- */ + +PairLJSphere::PairLJSphere(LAMMPS *lmp) : Pair(lmp) +{ + writedata = 1; +} + +/* ---------------------------------------------------------------------- */ + +PairLJSphere::~PairLJSphere() +{ + if (allocated) { + memory->destroy(setflag); + memory->destroy(cutsq); + + memory->destroy(cut); + memory->destroy(epsilon); + } +} + +/* ---------------------------------------------------------------------- */ + +void PairLJSphere::compute(int eflag, int vflag) +{ + int i, j, ii, jj, inum, jnum, itype, jtype; + double xtmp, ytmp, ztmp, rtmp, delx, dely, delz, evdwl, sigma6, fpair; + double rsq, r2inv, r6inv, forcelj, factor_lj; + int *ilist, *jlist, *numneigh, **firstneigh; + + evdwl = 0.0; + ev_init(eflag, vflag); + + double **x = atom->x; + double **f = atom->f; + double *radius = atom->radius; + 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]; + rtmp = radius[i]; + 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; + r6inv = r2inv * r2inv * r2inv; + + sigma6 = powint(mix_distance(rtmp, radius[j]), 6); + forcelj = r6inv * 24.0 * epsilon[itype][jtype] * (2.0 * sigma6 * sigma6 * r6inv - sigma6); + fpair = factor_lj * forcelj * r2inv; + + 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 = r6inv * 4.0 * epsilon[itype][jtype]; + evdwl *= sigma6 * sigma6 * r6inv - sigma6; + if (offset_flag && (cutsq[itype][jtype] > 0.0)) { + double ratio6 = sigma6 / powint(cutsq[itype][jtype], 3); + evdwl -= 4.0 * epsilon[itype][jtype] * (ratio6 * ratio6 - ratio6); + } + 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 PairLJSphere::allocate() +{ + allocated = 1; + int n = atom->ntypes + 1; + + memory->create(setflag, n, n, "pair:setflag"); + for (int i = 1; i < n; i++) + for (int j = i; j < n; j++) setflag[i][j] = 0; + + memory->create(cutsq, n, n, "pair:cutsq"); + + memory->create(cut, n, n, "pair:cut"); + memory->create(epsilon, n, n, "pair:epsilon"); +} + +/* ---------------------------------------------------------------------- + global settings +------------------------------------------------------------------------- */ + +void PairLJSphere::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 PairLJSphere::coeff(int narg, char **arg) +{ + if (narg < 3 || narg > 4) 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 cut_one = cut_global; + if (narg == 5) cut_one = utils::numeric(FLERR, arg[4], 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; + cut[i][j] = cut_one; + setflag[i][j] = 1; + count++; + } + } + + if (count == 0) error->all(FLERR, "Incorrect args for pair coefficients"); +} + +/* ---------------------------------------------------------------------- + init specific to this pair style +------------------------------------------------------------------------- */ + +void PairLJSphere::init_style() +{ + Pair::init_style(); + + if (!atom->radius_flag) error->all(FLERR, "Pair style lj/sphere requires atom attribute radius"); + if (mix_flag == SIXTHPOWER) + error->all(FLERR, "Pair_modify mix sixthpower is not compatible with pair style lj/sphere"); +} + +/* ---------------------------------------------------------------------- + init for one type pair i,j and corresponding j,i +------------------------------------------------------------------------- */ + +double PairLJSphere::init_one(int i, int j) +{ + if (setflag[i][j] == 0) { + epsilon[i][j] = mix_energy(epsilon[i][i], epsilon[j][j], 0.0, 0.0); + cut[i][j] = mix_distance(cut[i][i], cut[j][j]); + } + + epsilon[j][i] = epsilon[i][j]; + cut[j][i] = cut[i][j]; + return cut[i][j]; +} + +/* ---------------------------------------------------------------------- + proc 0 writes to restart file +------------------------------------------------------------------------- */ + +void PairLJSphere::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(&cut[i][j], sizeof(double), 1, fp); + } + } +} + +/* ---------------------------------------------------------------------- + proc 0 reads from restart file, bcasts +------------------------------------------------------------------------- */ + +void PairLJSphere::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) utils::sfread(FLERR, &setflag[i][j], sizeof(int), 1, fp, nullptr, error); + MPI_Bcast(&setflag[i][j], 1, MPI_INT, 0, world); + if (setflag[i][j]) { + if (me == 0) { + utils::sfread(FLERR, &epsilon[i][j], sizeof(double), 1, fp, nullptr, error); + utils::sfread(FLERR, &cut[i][j], sizeof(double), 1, fp, nullptr, error); + } + MPI_Bcast(&epsilon[i][j], 1, MPI_DOUBLE, 0, world); + MPI_Bcast(&cut[i][j], 1, MPI_DOUBLE, 0, world); + } + } +} + +/* ---------------------------------------------------------------------- + proc 0 writes to restart file +------------------------------------------------------------------------- */ + +void PairLJSphere::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 PairLJSphere::read_restart_settings(FILE *fp) +{ + int me = comm->me; + if (me == 0) { + utils::sfread(FLERR, &cut_global, sizeof(double), 1, fp, nullptr, error); + utils::sfread(FLERR, &offset_flag, sizeof(int), 1, fp, nullptr, error); + utils::sfread(FLERR, &mix_flag, sizeof(int), 1, fp, nullptr, error); + } + 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 PairLJSphere::write_data(FILE *fp) +{ + for (int i = 1; i <= atom->ntypes; i++) fprintf(fp, "%d %g\n", i, epsilon[i][i]); +} + +/* ---------------------------------------------------------------------- + proc 0 writes all pairs to data file +------------------------------------------------------------------------- */ + +void PairLJSphere::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\n", i, j, epsilon[i][j], cut[i][j]); +} + +/* ---------------------------------------------------------------------- */ + +double PairLJSphere::single(int i, int j, int itype, int jtype, double rsq, double /*factor_coul*/, + double factor_lj, double &fforce) +{ + double r2inv, r6inv, sigma6, forcelj, philj; + + sigma6 = powint(mix_distance(atom->radius[i], atom->radius[j]), 6); + r2inv = 1.0 / rsq; + r6inv = r2inv * r2inv * r2inv; + forcelj = r6inv * 24.0 * epsilon[itype][jtype] * (sigma6 * sigma6 * r6inv - sigma6); + fforce = factor_lj * forcelj * r2inv; + + philj = r6inv * 4.0 * epsilon[itype][jtype] * (sigma6 * sigma6 * r6inv - sigma6); + if (offset_flag && (cut[itype][jtype] > 0.0)) { + double ratio6 = sigma6 / powint(cutsq[itype][jtype], 3); + philj -= 4.0 * epsilon[itype][jtype] * (ratio6 * ratio6 - ratio6); + } + return factor_lj * philj; +} + +/* ---------------------------------------------------------------------- */ + +void *PairLJSphere::extract(const char *str, int &dim) +{ + dim = 2; + if (strcmp(str, "epsilon") == 0) return (void *) epsilon; + return nullptr; +} diff --git a/src/EXTRA-PAIR/pair_lj_sphere.h b/src/EXTRA-PAIR/pair_lj_sphere.h new file mode 100644 index 0000000000..d5125b690c --- /dev/null +++ b/src/EXTRA-PAIR/pair_lj_sphere.h @@ -0,0 +1,56 @@ +/* -*- c++ -*- ---------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + https://www.lammps.org/, Sandia National Laboratories + LAMMPS development team: developers@lammps.org + + 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 +// clang-format off +PairStyle(lj/sphere,PairLJSphere); +// clang-format on +#else + +#ifndef LMP_PAIR_LJ_SPHERE_H +#define LMP_PAIR_LJ_SPHERE_H + +#include "pair.h" + +namespace LAMMPS_NS { + +class PairLJSphere : public Pair { + public: + PairLJSphere(class LAMMPS *); + ~PairLJSphere() override; + void compute(int, int) override; + void settings(int, char **) override; + void coeff(int, char **) override; + void init_style() override; + double init_one(int, int) override; + void write_restart(FILE *) override; + void read_restart(FILE *) override; + void write_restart_settings(FILE *) override; + void read_restart_settings(FILE *) override; + void write_data(FILE *) override; + void write_data_all(FILE *) override; + double single(int, int, int, int, double, double, double, double &) override; + void *extract(const char *, int &) override; + + protected: + double cut_global; + double **cut; + double **epsilon; + + virtual void allocate(); +}; + +} // namespace LAMMPS_NS + +#endif +#endif