tx258/lammps-gran-kokkos
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Finally added pair_nm_cut_split fene_nm_cut_split

Browse Source
This commit is contained in:
jddietz
2021-09-23 13:23:29 -04:00
parent 7792a4db6b
commit 5b02f704cc
6 changed files with 626 additions and 9 deletions
Download Patch File Download Diff File Expand all files Collapse all files

25
doc/src/bond_fene.rst
View File

@ -1,4 +1,5 @@
.. index:: bond_style fene
.. index:: bond_style fene/nm/cut
.. index:: bond_style fene/intel
.. index:: bond_style fene/kk
.. index:: bond_style fene/omp
@ -8,13 +9,18 @@ bond_style fene command
Accelerator Variants: *fene/intel*, *fene/kk*, *fene/omp*
bond_style fene/nm/cut command
==============================
Syntax
""""""
.. code-block:: LAMMPS
bond_style fene
bond_style fene/nm/cut
Examples
""""""""
@ -23,6 +29,9 @@ Examples
bond_style fene
bond_coeff 1 30.0 1.5 1.0 1.0
bond_style fene/nm/cut
bond_coeff 1 2.25344 1.5 1.0 1.12246 2 6
Description
"""""""""""
@ -34,10 +43,15 @@ The *fene* bond style uses the potential
to define a finite extensible nonlinear elastic (FENE) potential
:ref:`(Kremer) <fene-Kremer>`, used for bead-spring polymer models. The first
term is attractive, the second Lennard-Jones term is repulsive. The
term is attractive, the second Lennard-Jones (LJ) term is repulsive. The
first term extends to :math:`R_0`, the maximum extent of the bond. The second
term is cutoff at :math:`2^\frac{1}{6} \sigma`, the minimum of the LJ potential.
The *fene/nm/cut* bond style substitutes the standard LJ potential with the generalized LJ potential. The bond energy is then given by
.. math::
E = -0.5 K R_0^2 \ln \left[ 1 - \left(\frac{r}{R_0}\right)^2\right] + \frac{\epsilon}{(n-m)} \left[ m \left(\frac{r_0}{r}\right)^n - n \left(\frac{r_0}{r}\right)^m \right]
The following coefficients must be defined for each bond type via the
:doc:`bond_coeff <bond_coeff>` command as in the example above, or in
the data file or restart files read by the :doc:`read_data <read_data>`
@ -48,6 +62,11 @@ or :doc:`read_restart <read_restart>` commands:
* :math:`\epsilon` (energy)
* :math:`\sigma` (distance)
For the *fene/nm/cut* style, the following additional coefficients are needed. Note, the standard LJ potential is recovered for (n=12 m=6).
* :math:`\n` (unitless)
* :math:`\m` (unitless)
----------
.. include:: accel_styles.rst
@ -58,7 +77,7 @@ Restrictions
""""""""""""
This bond style can only be used if LAMMPS was built with the MOLECULE
package. See the :doc:`Build package <Build_package>` page for more
package. See the :doc:`Build package <Build_package>` doc page for more
info.
You typically should specify :doc:`special_bonds fene <special_bonds>`

25
doc/src/pair_nm.rst
View File

@ -1,4 +1,5 @@
.. index:: pair_style nm/cut
.. index:: pair_style nm/cut/split
.. index:: pair_style nm/cut/coul/cut
.. index:: pair_style nm/cut/coul/long
.. index:: pair_style nm/cut/omp
@ -10,6 +11,11 @@ pair_style nm/cut command
Accelerator Variants: *nm/cut/omp*
pair_style nm/cut/split command
===============================
Accelerator Variants:*nm/cut/split/omp*
pair_style nm/cut/coul/cut command
==================================
@ -27,13 +33,15 @@ Syntax
pair_style style args
* style = *nm/cut* or *nm/cut/coul/cut* or *nm/cut/coul/long*
* style = *nm/cut* or *nm/cut/split* or *nm/cut/coul/cut* or *nm/cut/coul/long*
* args = list of arguments for a particular style
.. parsed-literal::
*nm/cut* args = cutoff
cutoff = global cutoff for Pair interactions (distance units)
*nm/cut/split* args = cutoff
cutoff = global cutoff for Pair interactions (distance units)
*nm/cut/coul/cut* args = cutoff (cutoff2)
cutoff = global cutoff for Pair (and Coulombic if only 1 arg) (distance units)
cutoff2 = global cutoff for Coulombic (optional) (distance units)
@ -50,6 +58,10 @@ Examples
pair_coeff * * 0.01 5.4 8.0 7.0
pair_coeff 1 1 0.01 4.4 7.0 6.0
pair_style nm/cut/split 1.12246
pair_coeff 1 1 1.0 1.1246 12 6
pair_coeff * * 1.0 1.1246 11 6
pair_style nm/cut/coul/cut 12.0 15.0
pair_coeff * * 0.01 5.4 8.0 7.0
pair_coeff 1 1 0.01 4.4 7.0 6.0
@ -73,6 +85,8 @@ interaction has the following form:
where :math:`r_c` is the cutoff.
Style *nm/cut/split* applies the standard LJ (6-12) potential abouve :math:`2^\frac{1}{6} \sigma`.
Style *nm/cut/coul/cut* adds a Coulombic pairwise interaction given by
.. math::
@ -114,7 +128,7 @@ N-M and Coulombic cutoffs specified in the pair_style command are used.
If only one cutoff is specified, it is used as the cutoff for both N-M
and Coulombic interactions for this type pair. If both coefficients
are specified, they are used as the N-M and Coulombic cutoffs for this
type pair. You cannot specify 2 cutoffs for style *nm*, since it
type pair. You cannot specify 2 cutoffs for style *nm*\ , since it
has no Coulombic terms.
For *nm/cut/coul/long* only the N-M cutoff can be specified since a
@ -147,7 +161,7 @@ to be specified in an input script that reads a restart file.
All of the *nm* pair styles can only be used via the *pair* keyword of
the :doc:`run_style respa <run_style>` command. They do not support the
*inner*, *middle*, *outer* keywords.
*inner*\ , *middle*\ , *outer* keywords.
----------
@ -156,9 +170,8 @@ the :doc:`run_style respa <run_style>` command. They do not support the
Restrictions
""""""""""""
These pair styles are part of the EXTRA-PAIR package. They are only enabled if
LAMMPS was built with that package. See the
:doc:`Build package <Build_package>` page for more info.
These pair styles are part of the MISC package. It is only enabled if
LAMMPS was built with that package. See the :doc:`Build package <Build_package>` doc page for more info.
Related commands
""""""""""""""""

310
src/EXTRA-MOLECULE/bond_fene_nm_split.cpp Normal file
View File

@ -0,0 +1,310 @@
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://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.
------------------------------------------------------------------------- */
#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 <cmath>
#include <cstring>
using namespace LAMMPS_NS;
using MathConst::MY_CUBEROOT2;
/* ---------------------------------------------------------------------- */
BondFENEnm::BondFENEnm(LAMMPS *lmp) : BondFENE(lmp)
{
// MY_CUBEROOT2 = pow(2.0,(1.0/3.0));
}
/* ---------------------------------------------------------------------- */
BondFENEnm::~BondFENEnm()
{
if (allocated && !copymode) {
memory->destroy(setflag);
memory->destroy(k);
memory->destroy(r0);
memory->destroy(epsilon);
memory->destroy(sigma);
memory->destroy(nn);
memory->destroy(mm);
}
}
/* ---------------------------------------------------------------------- */
void BondFENEnm::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) {
char str[128];
sprintf(str,"FENE bond too long: " BIGINT_FORMAT " "
TAGINT_FORMAT " " TAGINT_FORMAT " %g",
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
// MY_CUBEROOT2 cutoff assumes sigma = 2^{1/6}
if (rsq < MY_CUBEROOT2) {
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 < MY_CUBEROOT2) 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 BondFENEnm::allocate()
{
allocated = 1;
int n = atom->nbondtypes;
memory->create(k,n+1,"bond:k");
memory->create(r0,n+1,"bond:r0");
memory->create(epsilon,n+1,"bond:epsilon");
memory->create(sigma,n+1,"bond:sigma");
memory->create(nn,n+1,"bond:nn");
memory->create(mm,n+1,"bond:mm");
memory->create(setflag,n+1,"bond:setflag");
for (int i = 1; i <= n; i++) setflag[i] = 0;
}
/* ----------------------------------------------------------------------
set coeffs for one type
------------------------------------------------------------------------- */
void BondFENEnm::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 BondFENEnm::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 BondFENEnm::equilibrium_distance(int i)
{
return 0.97*sigma[i];
}
/* ----------------------------------------------------------------------
proc 0 writes to restart file
------------------------------------------------------------------------- */
void BondFENEnm::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 BondFENEnm::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 BondFENEnm::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 BondFENEnm::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) {
char str[128];
sprintf(str,"FENE bond too long: " BIGINT_FORMAT " %g",
update->ntimestep,sqrt(rsq));
// error->warning(FLERR,str,0);
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;
// MY_CUBEROOT2 cutoff assumes sigma = 2^1/6
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 *BondFENEnm::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;
}

76
src/EXTRA-MOLECULE/bond_fene_nm_split.h Executable file
View File

@ -0,0 +1,76 @@
/* -*- 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 BOND_CLASS
BondStyle(fene/nm/split,BondFENEnm)
#else
#ifndef LMP_BOND_FENE_NM_H
#define LMP_BOND_FENE_NM_H
#include "bond_fene.h"
namespace LAMMPS_NS {
class BondFENEnm : public BondFENE {
public :
BondFENEnm(class LAMMPS *);
virtual ~BondFENEnm();
virtual void compute(int, int);
virtual void coeff(int, char **);
void init_style();
double equilibrium_distance(int);
virtual void write_restart(FILE *);
void read_restart(FILE *);
void write_data(FILE *);
double single(int, double, int, int, double &);
virtual void *extract(const char *, int &);
protected:
double TWO_1_3;
double *k,*r0,*epsilon,*sigma,*nn,*mm;
virtual void allocate();
};
}
#endif
#endif
/* ERROR/WARNING messages:
W: FENE bond too long: %ld %d %d %g
A FENE bond has stretched dangerously far. It's interaction strength
will be truncated to attempt to prevent the bond from blowing up.
E: Bad FENE bond
Two atoms in a FENE bond have become so far apart that the bond cannot
be computed.
E: Incorrect args for bond coefficients
Self-explanatory. Check the input script or data file.
W: Use special bonds = 0,1,1 with bond style fene
Most FENE models need this setting for the special_bonds command.
W: FENE bond too long: %ld %g
A FENE bond has stretched dangerously far. It's interaction strength
will be truncated to attempt to prevent the bond from blowing up.
*/

148
src/EXTRA-PAIR/pair_nm_cut_split.cpp Normal file
View File

@ -0,0 +1,148 @@
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://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: Julien Devemy (ICCF), Robert S. Hoy (USF), Joseph D. Dietz (USF)
------------------------------------------------------------------------- */
#include "pair_nm_cut_split.h"
#include <cmath>
#include <cstring>
#include "atom.h"
#include "comm.h"
#include "force.h"
#include "neigh_list.h"
#include "math_const.h"
#include "memory.h"
#include "error.h"
using namespace LAMMPS_NS;
using namespace MathConst;
/* ---------------------------------------------------------------------- */
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);
rminv = pow(r2inv,mm[itype][jtype]/2.0);
rninv = pow(r2inv,nn[itype][jtype]/2.0);
// r < 2^1/6, use generlaized LJ
if (rsq < 1.25992105) {
forcenm = e0nm[itype][jtype]*nm[itype][jtype]*(r0n[itype][jtype]/pow(r,nn[itype][jtype])-r0m[itype][jtype]/pow(r,mm[itype][jtype]));
}
// r > 2^1/6 --> use standard LJ (m = 6 n = 12)
else{forcenm =(e0[itype][jtype]/6)*72*(4/pow(r,12)-2/pow(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 < 1.25992105) {
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)*(6*4*pow(r2inv,6) - 12*2*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 generlaized LJ
if (rsq < 1.25992105) {
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)*72*(4/pow(r,12)-2/pow(r,6));
phinm = (e0[itype][jtype]/6)*(6*4*pow(r2inv,6)-12*2*pow(r2inv,3));
}
fforce = factor_lj*forcenm*r2inv;
return factor_lj*phinm;
}

51
src/EXTRA-PAIR/pair_nm_cut_split.h Normal file
View File

@ -0,0 +1,51 @@
/* -*- c++ -*- ----------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://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(nm/cut/split,PairNMCutSplit)
#else
#ifndef LMP_PAIR_NM_CUT_SPLIT_H
#define LMP_PAIR_NM_CUT_SPLIT_H
#include "pair_nm_cut.h"
namespace LAMMPS_NS {
class PairNMCutSplit : public PairNMCut {
public :
PairNMCutSplit(class LAMMPS *);
double single(int, int, int, int, double, double, double, double &);
virtual void compute(int, int);
};
}
#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.
*/