diff --git a/src/EXTRA-MOLECULE/bond_fene_nm_split.cpp b/src/EXTRA-MOLECULE/bond_fene_nm_split.cpp index 5d9218e554..538d5ef34f 100644 --- a/src/EXTRA-MOLECULE/bond_fene_nm_split.cpp +++ b/src/EXTRA-MOLECULE/bond_fene_nm_split.cpp @@ -1,280 +1,280 @@ -/* ---------------------------------------------------------------------- - LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator - https://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. -------------------------------------------------------------------------- */ - -#include "bond_fene_nm_split.h" - -#include "atom.h" -#include "comm.h" -#include "error.h" -#include "force.h" -#include "math_const.h" -#include "memory.h" -#include "neighbor.h" -#include "update.h" - -#include -#include - -using namespace LAMMPS_NS; - -/* ---------------------------------------------------------------------- */ - -BondFENEnmSplit::BondFENEnmSplit(LAMMPS *lmp) : BondFENE(lmp) {} - -/* ---------------------------------------------------------------------- */ - -BondFENEnmSplit::~BondFENEnmSplit() -{ - if (allocated && !copymode) { - memory->destroy(nn); - memory->destroy(mm); - } -} - -/* ---------------------------------------------------------------------- */ - -void BondFENEnmSplit::compute(int eflag, int vflag) -{ - int i1, i2, n, type; - double delx, dely, delz, ebond, fbond; - double rsq, r0sq, rlogarg, sr2, sr6; - double r; - - ebond = sr6 = 0.0; - ev_init(eflag, vflag); - - double **x = atom->x; - double **f = atom->f; - int **bondlist = neighbor->bondlist; - int nbondlist = neighbor->nbondlist; - int nlocal = atom->nlocal; - int newton_bond = force->newton_bond; - - for (n = 0; n < nbondlist; n++) { - i1 = bondlist[n][0]; - i2 = bondlist[n][1]; - type = bondlist[n][2]; - - delx = x[i1][0] - x[i2][0]; - dely = x[i1][1] - x[i2][1]; - delz = x[i1][2] - x[i2][2]; - - // force from log term - - rsq = delx * delx + dely * dely + delz * delz; - r0sq = r0[type] * r0[type]; - rlogarg = 1.0 - rsq / r0sq; - - // if r -> r0, then rlogarg < 0.0 which is an error - // issue a warning and reset rlogarg = epsilon - // if r > 2*r0 something serious is wrong, abort - // change cutuff from .1 to .02 so only bond lengths > 1.485 give the warning - // and crash the run if rlogarg < -.21 rather than < -3 - // Don't print out warnings, only errors - if (rlogarg < .02) { - error->warning(FLERR, "fene/nm/split bond too long: {} {} {} {}", update->ntimestep, - atom->tag[i1], atom->tag[i2], sqrt(rsq)); - if (rlogarg <= -.21) error->one(FLERR, "Bad FENE bond"); - rlogarg = 0.02; - } - - fbond = -k[type] / rlogarg; - // force from n-m term - if (rsq < sigma[type]*sigma[type]) { - r = sqrt(rsq); - fbond += epsilon[type] * (nn[type] * mm[type] / (nn[type] - mm[type])) * - (pow(sigma[type] / r, nn[type]) - pow(sigma[type] / r, mm[type])) / rsq; - } - - // energy - - if (eflag) { - ebond = -0.5 * k[type] * r0sq * log(rlogarg); - if (rsq < sigma[type]*sigma[type]) - ebond += (epsilon[type] / (nn[type] - mm[type])) * - (mm[type] * pow(sigma[type] / r, nn[type]) - nn[type] * pow(sigma[type] / r, mm[type])); - } - // apply force to each of 2 atoms - if (newton_bond || i1 < nlocal) { - f[i1][0] += delx * fbond; - f[i1][1] += dely * fbond; - f[i1][2] += delz * fbond; - } - - if (newton_bond || i2 < nlocal) { - f[i2][0] -= delx * fbond; - f[i2][1] -= dely * fbond; - f[i2][2] -= delz * fbond; - } - - if (evflag) ev_tally(i1, i2, nlocal, newton_bond, ebond, fbond, delx, dely, delz); - } -} - -/* ---------------------------------------------------------------------- */ - -void BondFENEnmSplit::allocate() -{ - BondFENE::allocate(); - int n = atom->nbondtypes + 1; - memory->create(nn, n, "bond:nn"); - memory->create(mm, n, "bond:mm"); -} - -/* ---------------------------------------------------------------------- - set coeffs for one type -------------------------------------------------------------------------- */ - -void BondFENEnmSplit::coeff(int narg, char **arg) -{ - if (narg != 7) error->all(FLERR, "Incorrect args for bond coefficients"); - if (!allocated) allocate(); - - int ilo, ihi; - utils::bounds(FLERR, arg[0], 1, atom->nbondtypes, ilo, ihi, error); - - double k_one = utils::numeric(FLERR, arg[1], false, lmp); - double r0_one = utils::numeric(FLERR, arg[2], false, lmp); - double epsilon_one = utils::numeric(FLERR, arg[3], false, lmp); - double sigma_one = utils::numeric(FLERR, arg[4], false, lmp); - double nn_one = utils::numeric(FLERR, arg[5], false, lmp); - double mm_one = utils::numeric(FLERR, arg[6], false, lmp); - - int count = 0; - for (int i = ilo; i <= ihi; i++) { - k[i] = k_one; - r0[i] = r0_one; - epsilon[i] = epsilon_one; - sigma[i] = sigma_one; - nn[i] = nn_one; - mm[i] = mm_one; - setflag[i] = 1; - count++; - } - - if (count == 0) error->all(FLERR, "Incorrect args for bond coefficients"); -} - -/* ---------------------------------------------------------------------- - check if special_bond settings are valid -------------------------------------------------------------------------- */ - -void BondFENEnmSplit::init_style() -{ - // special bonds should be 0 1 1 - - if (force->special_lj[1] != 0.0 || force->special_lj[2] != 1.0 || force->special_lj[3] != 1.0) { - if (comm->me == 0) error->warning(FLERR, "Use special bonds = 0,1,1 with bond style fene"); - } -} - -/* ---------------------------------------------------------------------- */ - -double BondFENEnmSplit::equilibrium_distance(int i) -{ - return 0.97 * sigma[i]; -} - -/* ---------------------------------------------------------------------- - proc 0 writes to restart file -------------------------------------------------------------------------- */ - -void BondFENEnmSplit::write_restart(FILE *fp) -{ - fwrite(&k[1], sizeof(double), atom->nbondtypes, fp); - fwrite(&r0[1], sizeof(double), atom->nbondtypes, fp); - fwrite(&epsilon[1], sizeof(double), atom->nbondtypes, fp); - fwrite(&sigma[1], sizeof(double), atom->nbondtypes, fp); - fwrite(&nn[1], sizeof(double), atom->nbondtypes, fp); - fwrite(&mm[1], sizeof(double), atom->nbondtypes, fp); -} - -/* ---------------------------------------------------------------------- - proc 0 reads from restart file, bcasts -------------------------------------------------------------------------- */ - -void BondFENEnmSplit::read_restart(FILE *fp) -{ - allocate(); - - if (comm->me == 0) { - utils::sfread(FLERR, &k[1], sizeof(double), atom->nbondtypes, fp, nullptr, error); - utils::sfread(FLERR, &r0[1], sizeof(double), atom->nbondtypes, fp, nullptr, error); - utils::sfread(FLERR, &epsilon[1], sizeof(double), atom->nbondtypes, fp, nullptr, error); - utils::sfread(FLERR, &sigma[1], sizeof(double), atom->nbondtypes, fp, nullptr, error); - utils::sfread(FLERR, &nn[1], sizeof(double), atom->nbondtypes, fp, nullptr, error); - utils::sfread(FLERR, &mm[1], sizeof(double), atom->nbondtypes, fp, nullptr, error); - } - MPI_Bcast(&k[1], atom->nbondtypes, MPI_DOUBLE, 0, world); - MPI_Bcast(&r0[1], atom->nbondtypes, MPI_DOUBLE, 0, world); - MPI_Bcast(&epsilon[1], atom->nbondtypes, MPI_DOUBLE, 0, world); - MPI_Bcast(&sigma[1], atom->nbondtypes, MPI_DOUBLE, 0, world); - MPI_Bcast(&nn[1], atom->nbondtypes, MPI_DOUBLE, 0, world); - MPI_Bcast(&mm[1], atom->nbondtypes, MPI_DOUBLE, 0, world); - - for (int i = 1; i <= atom->nbondtypes; i++) setflag[i] = 1; -} - -/* ---------------------------------------------------------------------- - proc 0 writes to data file -------------------------------------------------------------------------- */ - -void BondFENEnmSplit::write_data(FILE *fp) -{ - for (int i = 1; i <= atom->nbondtypes; i++) - fprintf(fp, "%d %g %g %g %g %g %g\n", i, k[i], r0[i], epsilon[i], sigma[i], nn[i], mm[i]); -} - -/* ---------------------------------------------------------------------- */ - -double BondFENEnmSplit::single(int type, double rsq, int /*i*/, int /*j*/, double &fforce) -{ - double r0sq = r0[type] * r0[type]; - double rlogarg = 1.0 - rsq / r0sq; - double r; - // if r -> r0, then rlogarg < 0.0 which is an error - // issue a warning and reset rlogarg = epsilon - // if r > 2*r0 something serious is wrong, abort - - // change cutuff from .1 to .02 so only bond lengths > 1.485 give the warning - // and crash the run if rlogarg < -.21 rather than < -3 - // Don't print out warnings, only errors - if (rlogarg < 0.02) { - error->warning(FLERR, "FENE bond too long: {} {:.8}", update->ntimestep, sqrt(rsq)); - if (rlogarg <= -.21) error->one(FLERR, "Bad FENE bond"); - rlogarg = 0.02; - } - - double eng = -0.5 * k[type] * r0sq * log(rlogarg); - fforce = -k[type] / rlogarg; - - if (rsq < sigma[type]*sigma[type]) { - r = sqrt(rsq); - fforce += epsilon[type] * (nn[type] * mm[type] / (nn[type] - mm[type])) * - (pow(sigma[type] / r, nn[type]) - pow(sigma[type] / r, mm[type])) / rsq; - eng += (epsilon[type] / (nn[type] - mm[type])) * - (mm[type] * pow(sigma[type] / r, nn[type]) - nn[type] * pow(sigma[type] / r, mm[type])); - } - - return eng; -} - -/* ---------------------------------------------------------------------- */ - -void *BondFENEnmSplit::extract(const char *str, int &dim) -{ - dim = 1; - if (strcmp(str, "kappa") == 0) return (void *) k; - if (strcmp(str, "r0") == 0) return (void *) r0; - return nullptr; -} +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + https://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. +------------------------------------------------------------------------- */ + +#include "bond_fene_nm_split.h" + +#include "atom.h" +#include "comm.h" +#include "error.h" +#include "force.h" +#include "math_const.h" +#include "memory.h" +#include "neighbor.h" +#include "update.h" + +#include +#include + +using namespace LAMMPS_NS; + +/* ---------------------------------------------------------------------- */ + +BondFENEnmSplit::BondFENEnmSplit(LAMMPS *lmp) : BondFENE(lmp) {} + +/* ---------------------------------------------------------------------- */ + +BondFENEnmSplit::~BondFENEnmSplit() +{ + if (allocated && !copymode) { + memory->destroy(nn); + memory->destroy(mm); + } +} + +/* ---------------------------------------------------------------------- */ + +void BondFENEnmSplit::compute(int eflag, int vflag) +{ + int i1, i2, n, type; + double delx, dely, delz, ebond, fbond; + double rsq, r0sq, rlogarg, sr2, sr6; + double r; + + ebond = sr6 = 0.0; + ev_init(eflag, vflag); + + double **x = atom->x; + double **f = atom->f; + int **bondlist = neighbor->bondlist; + int nbondlist = neighbor->nbondlist; + int nlocal = atom->nlocal; + int newton_bond = force->newton_bond; + + for (n = 0; n < nbondlist; n++) { + i1 = bondlist[n][0]; + i2 = bondlist[n][1]; + type = bondlist[n][2]; + + delx = x[i1][0] - x[i2][0]; + dely = x[i1][1] - x[i2][1]; + delz = x[i1][2] - x[i2][2]; + + // force from log term + + rsq = delx * delx + dely * dely + delz * delz; + r0sq = r0[type] * r0[type]; + rlogarg = 1.0 - rsq / r0sq; + + // if r -> r0, then rlogarg < 0.0 which is an error + // issue a warning and reset rlogarg = epsilon + // if r > 2*r0 something serious is wrong, abort + // change cutuff from .1 to .02 so only bond lengths > 1.485 give the warning + // and crash the run if rlogarg < -.21 rather than < -3 + // Don't print out warnings, only errors + if (rlogarg < .02) { + error->warning(FLERR, "fene/nm/split bond too long: {} {} {} {}", update->ntimestep, + atom->tag[i1], atom->tag[i2], sqrt(rsq)); + if (rlogarg <= -.21) error->one(FLERR, "Bad FENE bond"); + rlogarg = 0.02; + } + + fbond = -k[type] / rlogarg; + // force from n-m term + if (rsq < sigma[type]*sigma[type]) { + r = sqrt(rsq); + fbond += epsilon[type] * (nn[type] * mm[type] / (nn[type] - mm[type])) * + (pow(sigma[type] / r, nn[type]) - pow(sigma[type] / r, mm[type])) / rsq; + } + + // energy + + if (eflag) { + ebond = -0.5 * k[type] * r0sq * log(rlogarg); + if (rsq < sigma[type]*sigma[type]) + ebond += (epsilon[type] / (nn[type] - mm[type])) * + (mm[type] * pow(sigma[type] / r, nn[type]) - nn[type] * pow(sigma[type] / r, mm[type])); + } + // apply force to each of 2 atoms + if (newton_bond || i1 < nlocal) { + f[i1][0] += delx * fbond; + f[i1][1] += dely * fbond; + f[i1][2] += delz * fbond; + } + + if (newton_bond || i2 < nlocal) { + f[i2][0] -= delx * fbond; + f[i2][1] -= dely * fbond; + f[i2][2] -= delz * fbond; + } + + if (evflag) ev_tally(i1, i2, nlocal, newton_bond, ebond, fbond, delx, dely, delz); + } +} + +/* ---------------------------------------------------------------------- */ + +void BondFENEnmSplit::allocate() +{ + BondFENE::allocate(); + int n = atom->nbondtypes + 1; + memory->create(nn, n, "bond:nn"); + memory->create(mm, n, "bond:mm"); +} + +/* ---------------------------------------------------------------------- + set coeffs for one type +------------------------------------------------------------------------- */ + +void BondFENEnmSplit::coeff(int narg, char **arg) +{ + if (narg != 7) error->all(FLERR, "Incorrect args for bond coefficients"); + if (!allocated) allocate(); + + int ilo, ihi; + utils::bounds(FLERR, arg[0], 1, atom->nbondtypes, ilo, ihi, error); + + double k_one = utils::numeric(FLERR, arg[1], false, lmp); + double r0_one = utils::numeric(FLERR, arg[2], false, lmp); + double epsilon_one = utils::numeric(FLERR, arg[3], false, lmp); + double sigma_one = utils::numeric(FLERR, arg[4], false, lmp); + double nn_one = utils::numeric(FLERR, arg[5], false, lmp); + double mm_one = utils::numeric(FLERR, arg[6], false, lmp); + + int count = 0; + for (int i = ilo; i <= ihi; i++) { + k[i] = k_one; + r0[i] = r0_one; + epsilon[i] = epsilon_one; + sigma[i] = sigma_one; + nn[i] = nn_one; + mm[i] = mm_one; + setflag[i] = 1; + count++; + } + + if (count == 0) error->all(FLERR, "Incorrect args for bond coefficients"); +} + +/* ---------------------------------------------------------------------- + check if special_bond settings are valid +------------------------------------------------------------------------- */ + +void BondFENEnmSplit::init_style() +{ + // special bonds should be 0 1 1 + + if (force->special_lj[1] != 0.0 || force->special_lj[2] != 1.0 || force->special_lj[3] != 1.0) { + if (comm->me == 0) error->warning(FLERR, "Use special bonds = 0,1,1 with bond style fene"); + } +} + +/* ---------------------------------------------------------------------- */ + +double BondFENEnmSplit::equilibrium_distance(int i) +{ + return 0.97 * sigma[i]; +} + +/* ---------------------------------------------------------------------- + proc 0 writes to restart file +------------------------------------------------------------------------- */ + +void BondFENEnmSplit::write_restart(FILE *fp) +{ + fwrite(&k[1], sizeof(double), atom->nbondtypes, fp); + fwrite(&r0[1], sizeof(double), atom->nbondtypes, fp); + fwrite(&epsilon[1], sizeof(double), atom->nbondtypes, fp); + fwrite(&sigma[1], sizeof(double), atom->nbondtypes, fp); + fwrite(&nn[1], sizeof(double), atom->nbondtypes, fp); + fwrite(&mm[1], sizeof(double), atom->nbondtypes, fp); +} + +/* ---------------------------------------------------------------------- + proc 0 reads from restart file, bcasts +------------------------------------------------------------------------- */ + +void BondFENEnmSplit::read_restart(FILE *fp) +{ + allocate(); + + if (comm->me == 0) { + utils::sfread(FLERR, &k[1], sizeof(double), atom->nbondtypes, fp, nullptr, error); + utils::sfread(FLERR, &r0[1], sizeof(double), atom->nbondtypes, fp, nullptr, error); + utils::sfread(FLERR, &epsilon[1], sizeof(double), atom->nbondtypes, fp, nullptr, error); + utils::sfread(FLERR, &sigma[1], sizeof(double), atom->nbondtypes, fp, nullptr, error); + utils::sfread(FLERR, &nn[1], sizeof(double), atom->nbondtypes, fp, nullptr, error); + utils::sfread(FLERR, &mm[1], sizeof(double), atom->nbondtypes, fp, nullptr, error); + } + MPI_Bcast(&k[1], atom->nbondtypes, MPI_DOUBLE, 0, world); + MPI_Bcast(&r0[1], atom->nbondtypes, MPI_DOUBLE, 0, world); + MPI_Bcast(&epsilon[1], atom->nbondtypes, MPI_DOUBLE, 0, world); + MPI_Bcast(&sigma[1], atom->nbondtypes, MPI_DOUBLE, 0, world); + MPI_Bcast(&nn[1], atom->nbondtypes, MPI_DOUBLE, 0, world); + MPI_Bcast(&mm[1], atom->nbondtypes, MPI_DOUBLE, 0, world); + + for (int i = 1; i <= atom->nbondtypes; i++) setflag[i] = 1; +} + +/* ---------------------------------------------------------------------- + proc 0 writes to data file +------------------------------------------------------------------------- */ + +void BondFENEnmSplit::write_data(FILE *fp) +{ + for (int i = 1; i <= atom->nbondtypes; i++) + fprintf(fp, "%d %g %g %g %g %g %g\n", i, k[i], r0[i], epsilon[i], sigma[i], nn[i], mm[i]); +} + +/* ---------------------------------------------------------------------- */ + +double BondFENEnmSplit::single(int type, double rsq, int /*i*/, int /*j*/, double &fforce) +{ + double r0sq = r0[type] * r0[type]; + double rlogarg = 1.0 - rsq / r0sq; + double r; + // if r -> r0, then rlogarg < 0.0 which is an error + // issue a warning and reset rlogarg = epsilon + // if r > 2*r0 something serious is wrong, abort + + // change cutuff from .1 to .02 so only bond lengths > 1.485 give the warning + // and crash the run if rlogarg < -.21 rather than < -3 + // Don't print out warnings, only errors + if (rlogarg < 0.02) { + error->warning(FLERR, "FENE bond too long: {} {:.8}", update->ntimestep, sqrt(rsq)); + if (rlogarg <= -.21) error->one(FLERR, "Bad FENE bond"); + rlogarg = 0.02; + } + + double eng = -0.5 * k[type] * r0sq * log(rlogarg); + fforce = -k[type] / rlogarg; + + if (rsq < sigma[type]*sigma[type]) { + r = sqrt(rsq); + fforce += epsilon[type] * (nn[type] * mm[type] / (nn[type] - mm[type])) * + (pow(sigma[type] / r, nn[type]) - pow(sigma[type] / r, mm[type])) / rsq; + eng += (epsilon[type] / (nn[type] - mm[type])) * + (mm[type] * pow(sigma[type] / r, nn[type]) - nn[type] * pow(sigma[type] / r, mm[type])); + } + + return eng; +} + +/* ---------------------------------------------------------------------- */ + +void *BondFENEnmSplit::extract(const char *str, int &dim) +{ + dim = 1; + if (strcmp(str, "kappa") == 0) return (void *) k; + if (strcmp(str, "r0") == 0) return (void *) r0; + return nullptr; +} diff --git a/src/EXTRA-PAIR/pair_nm_cut_split.cpp b/src/EXTRA-PAIR/pair_nm_cut_split.cpp index e8719b8479..97a4d38079 100644 --- a/src/EXTRA-PAIR/pair_nm_cut_split.cpp +++ b/src/EXTRA-PAIR/pair_nm_cut_split.cpp @@ -1,156 +1,156 @@ -// clang-format off -/* ---------------------------------------------------------------------- - LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator - https://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 Author: Julien Devemy (ICCF), Robert S. Hoy (USF), Joseph D. Dietz (USF) -------------------------------------------------------------------------- */ - -#include "pair_nm_cut_split.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 -#include - -using namespace LAMMPS_NS; -using MathSpecial::powint; - -/* ---------------------------------------------------------------------- */ -PairNMCutSplit::PairNMCutSplit(LAMMPS *lmp) : PairNMCut(lmp) -{ - writedata = 1; -} - -void PairNMCutSplit::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 r,forcenm,rminv,rninv; - int *ilist,*jlist,*numneigh,**firstneigh; - - evdwl = 0.0; - ev_init(eflag,vflag); - - 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; - r = sqrt(rsq); - -// r < r0 --> use generalized LJ - if (rsq < r0[itype][jtype]*r0[itype][jtype]) { - forcenm = e0nm[itype][jtype]*nm[itype][jtype]* - (r0n[itype][jtype]/pow(r,nn[itype][jtype]) - -r0m[itype][jtype]/pow(r,mm[itype][jtype])); - } -// r > r0 --> use standard LJ (m = 6 n = 12) - else forcenm =(e0[itype][jtype]/6.0)*72.0*(4.0/powint(r,12)-2.0/powint(r,6)); - - fpair = factor_lj*forcenm*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) { - if (rsq < r0[itype][jtype]*r0[itype][jtype]) { - rminv = pow(r2inv,mm[itype][jtype]/2.0); - rninv = pow(r2inv,nn[itype][jtype]/2.0); - - evdwl = e0nm[itype][jtype]*(mm[itype][jtype]*r0n[itype][jtype]*rninv - - nn[itype][jtype]*r0m[itype][jtype]*rminv) - - offset[itype][jtype]; - } else evdwl = (e0[itype][jtype]/6.0)*(24.0*powint(r2inv,6) - 24.0*pow(r2inv,3)); - 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(); -} - -/* ---------------------------------------------------------------------- */ - -double PairNMCutSplit::single(int /*i*/, int /*j*/, int itype, int jtype, double rsq, double /*factor_coul*/, double factor_lj, double &fforce) -{ - double r2inv,rminv,rninv,r,forcenm,phinm; - - r2inv = 1.0/rsq; - r = sqrt(rsq); - rminv = pow(r2inv,mm[itype][jtype]/2.0); - rninv = pow(r2inv,nn[itype][jtype]/2.0); - // r < 2^1/6, use generalized LJ - if (rsq < r0[itype][jtype]*r0[itype][jtype]) { // note the addition of the r0 factor - forcenm = e0nm[itype][jtype]*nm[itype][jtype]* - (r0n[itype][jtype]/pow(r,nn[itype][jtype])-r0m[itype][jtype]/pow(r,mm[itype][jtype])); - phinm = e0nm[itype][jtype]*(mm[itype][jtype]*r0n[itype][jtype]*rninv - -nn[itype][jtype]*r0m[itype][jtype]*rminv)-offset[itype][jtype]; - - } - // r > 2^1/6 --> use standard LJ (m = 6 n = 12) - else{forcenm = (e0[itype][jtype]/6.0)*72.0*(4.0/powint(r,12)-2.0/powint(r,6)); - phinm = (e0[itype][jtype]/6.0)*(24.0*powint(r2inv,6)-24.0*powint(r2inv,3)); - } - - fforce = factor_lj*forcenm*r2inv; - return factor_lj*phinm; -} - +// clang-format off +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + https://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 Author: Julien Devemy (ICCF), Robert S. Hoy (USF), Joseph D. Dietz (USF) +------------------------------------------------------------------------- */ + +#include "pair_nm_cut_split.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 +#include + +using namespace LAMMPS_NS; +using MathSpecial::powint; + +/* ---------------------------------------------------------------------- */ +PairNMCutSplit::PairNMCutSplit(LAMMPS *lmp) : PairNMCut(lmp) +{ + writedata = 1; +} + +void PairNMCutSplit::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 r,forcenm,rminv,rninv; + int *ilist,*jlist,*numneigh,**firstneigh; + + evdwl = 0.0; + ev_init(eflag,vflag); + + 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; + r = sqrt(rsq); + +// r < r0 --> use generalized LJ + if (rsq < r0[itype][jtype]*r0[itype][jtype]) { + forcenm = e0nm[itype][jtype]*nm[itype][jtype]* + (r0n[itype][jtype]/pow(r,nn[itype][jtype]) + -r0m[itype][jtype]/pow(r,mm[itype][jtype])); + } +// r > r0 --> use standard LJ (m = 6 n = 12) + else forcenm =(e0[itype][jtype]/6.0)*72.0*(4.0/powint(r,12)-2.0/powint(r,6)); + + fpair = factor_lj*forcenm*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) { + if (rsq < r0[itype][jtype]*r0[itype][jtype]) { + rminv = pow(r2inv,mm[itype][jtype]/2.0); + rninv = pow(r2inv,nn[itype][jtype]/2.0); + + evdwl = e0nm[itype][jtype]*(mm[itype][jtype]*r0n[itype][jtype]*rninv - + nn[itype][jtype]*r0m[itype][jtype]*rminv) - + offset[itype][jtype]; + } else evdwl = (e0[itype][jtype]/6.0)*(24.0*powint(r2inv,6) - 24.0*pow(r2inv,3)); + 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(); +} + +/* ---------------------------------------------------------------------- */ + +double PairNMCutSplit::single(int /*i*/, int /*j*/, int itype, int jtype, double rsq, double /*factor_coul*/, double factor_lj, double &fforce) +{ + double r2inv,rminv,rninv,r,forcenm,phinm; + + r2inv = 1.0/rsq; + r = sqrt(rsq); + rminv = pow(r2inv,mm[itype][jtype]/2.0); + rninv = pow(r2inv,nn[itype][jtype]/2.0); + // r < 2^1/6, use generalized LJ + if (rsq < r0[itype][jtype]*r0[itype][jtype]) { // note the addition of the r0 factor + forcenm = e0nm[itype][jtype]*nm[itype][jtype]* + (r0n[itype][jtype]/pow(r,nn[itype][jtype])-r0m[itype][jtype]/pow(r,mm[itype][jtype])); + phinm = e0nm[itype][jtype]*(mm[itype][jtype]*r0n[itype][jtype]*rninv + -nn[itype][jtype]*r0m[itype][jtype]*rminv)-offset[itype][jtype]; + + } + // r > 2^1/6 --> use standard LJ (m = 6 n = 12) + else{forcenm = (e0[itype][jtype]/6.0)*72.0*(4.0/powint(r,12)-2.0/powint(r,6)); + phinm = (e0[itype][jtype]/6.0)*(24.0*powint(r2inv,6)-24.0*powint(r2inv,3)); + } + + fforce = factor_lj*forcenm*r2inv; + return factor_lj*phinm; +} +