consistently support special_bonds settings in pair style gauss
This commit is contained in:
@ -137,16 +137,19 @@ interaction.
|
||||
The :doc:`pair_modify <pair_modify>` table and tail options are not
|
||||
relevant for these pair styles.
|
||||
|
||||
These pair styles write their 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.
|
||||
These pair styles write their 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.
|
||||
|
||||
These 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.
|
||||
|
||||
The *gauss* pair style tallies an "occupancy" count of how many Gaussian-well
|
||||
sites have an atom within the distance at which the force is a maximum
|
||||
= sqrt(0.5/b). This quantity can be accessed via the :doc:`compute pair <compute_pair>` command as a vector of values of length 1.
|
||||
The *gauss* pair style tallies an "occupancy" count of how many
|
||||
Gaussian-well sites have an atom within the distance at which the force
|
||||
is a maximum = sqrt(0.5/b). This quantity can be accessed via the
|
||||
:doc:`compute pair <compute_pair>` command as a vector of values of
|
||||
length 1.
|
||||
|
||||
To print this quantity to the log file (with a descriptive column
|
||||
heading) the following commands could be included in an input script:
|
||||
@ -166,14 +169,10 @@ The *gauss* and *gauss/cut* styles are part of the EXTRA-PAIR package.
|
||||
They are only enabled if LAMMPS is build with that package. See the
|
||||
:doc:`Build package <Build_package>` page for more info.
|
||||
|
||||
The *gauss* style does not apply :doc:`special_bonds <special_bonds>`
|
||||
factors. When using this pair style on a system that has bonds, the
|
||||
special_bonds factors, if using the default setting of 0.0, may need to
|
||||
be adjusted to some very small number (e.g. 1.0e-100), so that those
|
||||
special pairs are not completely excluded from the neighbor lists, but
|
||||
won't contribute forces or energies from styles (e.g. when used in
|
||||
combination with a :doc:`hybrid pair style <pair_hybrid>`) that do
|
||||
apply those factors.
|
||||
.. versionchanged:: TBD
|
||||
|
||||
Prior to this version, the *gauss* pair style did not apply
|
||||
:doc:`special_bonds <special_bonds>` factors.
|
||||
|
||||
Related commands
|
||||
""""""""""""""""
|
||||
|
||||
@ -1,4 +1,3 @@
|
||||
// clang-format off
|
||||
/* ----------------------------------------------------------------------
|
||||
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
|
||||
https://www.lammps.org/, Sandia National Laboratories
|
||||
@ -18,20 +17,17 @@
|
||||
|
||||
#include "pair_gauss.h"
|
||||
|
||||
#include <cmath>
|
||||
#include <cstring>
|
||||
#include "atom.h"
|
||||
#include "comm.h"
|
||||
#include "force.h"
|
||||
#include "neigh_list.h"
|
||||
#include "memory.h"
|
||||
#include "error.h"
|
||||
#include "force.h"
|
||||
#include "memory.h"
|
||||
#include "neigh_list.h"
|
||||
|
||||
#include <cmath>
|
||||
|
||||
using namespace LAMMPS_NS;
|
||||
|
||||
#define EPSILON 1.0e-10
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
PairGauss::PairGauss(LAMMPS *lmp) : Pair(lmp)
|
||||
@ -45,7 +41,7 @@ PairGauss::PairGauss(LAMMPS *lmp) : Pair(lmp)
|
||||
|
||||
PairGauss::~PairGauss()
|
||||
{
|
||||
delete [] pvector;
|
||||
delete[] pvector;
|
||||
|
||||
if (allocated) {
|
||||
memory->destroy(setflag);
|
||||
@ -62,19 +58,19 @@ PairGauss::~PairGauss()
|
||||
|
||||
void PairGauss::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;
|
||||
int *ilist,*jlist,*numneigh,**firstneigh;
|
||||
int i, j, ii, jj, inum, jnum, itype, jtype;
|
||||
double xtmp, ytmp, ztmp, rsq, delx, dely, delz, evdwl, fpair, factor_lj;
|
||||
int *ilist, *jlist, *numneigh, **firstneigh;
|
||||
|
||||
evdwl = 0.0;
|
||||
ev_init(eflag,vflag);
|
||||
ev_init(eflag, vflag);
|
||||
int occ = 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;
|
||||
@ -95,36 +91,39 @@ void PairGauss::compute(int eflag, int vflag)
|
||||
|
||||
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;
|
||||
rsq = delx * delx + dely * dely + delz * delz;
|
||||
jtype = type[j];
|
||||
|
||||
// define a Gaussian well to be occupied if
|
||||
// the site it interacts with is within the force maximum
|
||||
|
||||
if (eflag_global && rsq < 0.5/b[itype][jtype]) occ++;
|
||||
if (eflag_global && rsq < 0.5 / b[itype][jtype]) occ++;
|
||||
|
||||
if (rsq < cutsq[itype][jtype]) {
|
||||
fpair = -2.0*a[itype][jtype]*b[itype][jtype] * exp(-b[itype][jtype]*rsq);
|
||||
fpair = -2.0 * a[itype][jtype] * b[itype][jtype] * exp(-b[itype][jtype] * rsq);
|
||||
fpair *= factor_lj;
|
||||
|
||||
f[i][0] += delx*fpair;
|
||||
f[i][1] += dely*fpair;
|
||||
f[i][2] += delz*fpair;
|
||||
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;
|
||||
f[j][0] -= delx * fpair;
|
||||
f[j][1] -= dely * fpair;
|
||||
f[j][2] -= delz * fpair;
|
||||
}
|
||||
|
||||
if (eflag)
|
||||
evdwl = -(a[itype][jtype]*exp(-b[itype][jtype]*rsq) - offset[itype][jtype]);
|
||||
if (eflag) {
|
||||
evdwl = -(a[itype][jtype] * exp(-b[itype][jtype] * rsq) - offset[itype][jtype]);
|
||||
evdwl *= factor_lj;
|
||||
}
|
||||
|
||||
if (evflag) ev_tally(i,j,nlocal,newton_pair,
|
||||
evdwl,0.0,fpair,delx,dely,delz);
|
||||
if (evflag) ev_tally(i, j, nlocal, newton_pair, evdwl, 0.0, fpair, delx, dely, delz);
|
||||
}
|
||||
}
|
||||
}
|
||||
@ -140,19 +139,18 @@ void PairGauss::compute(int eflag, int vflag)
|
||||
void PairGauss::allocate()
|
||||
{
|
||||
allocated = 1;
|
||||
int n = atom->ntypes;
|
||||
int np1 = atom->ntypes + 1;
|
||||
|
||||
memory->create(setflag,n+1,n+1,"pair:setflag");
|
||||
for (int i = 1; i <= n; i++)
|
||||
for (int j = 1; j <= n; j++)
|
||||
setflag[i][j] = 0;
|
||||
memory->create(setflag, np1, np1, "pair:setflag");
|
||||
for (int i = 1; i < np1; i++)
|
||||
for (int j = 1; j < np1; j++) setflag[i][j] = 0;
|
||||
|
||||
memory->create(cutsq,n+1,n+1,"pair:cutsq");
|
||||
memory->create(cutsq, np1, np1, "pair:cutsq");
|
||||
|
||||
memory->create(cut,n+1,n+1,"pair:cut_gauss");
|
||||
memory->create(a,n+1,n+1,"pair:a");
|
||||
memory->create(b,n+1,n+1,"pair:b");
|
||||
memory->create(offset,n+1,n+1,"pair:offset");
|
||||
memory->create(cut, np1, np1, "pair:cut_gauss");
|
||||
memory->create(a, np1, np1, "pair:a");
|
||||
memory->create(b, np1, np1, "pair:b");
|
||||
memory->create(offset, np1, np1, "pair:offset");
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
@ -161,14 +159,14 @@ void PairGauss::allocate()
|
||||
|
||||
void PairGauss::settings(int narg, char **arg)
|
||||
{
|
||||
if (narg != 1) error->all(FLERR,"Illegal pair_style command");
|
||||
if (narg != 1) error->all(FLERR, "Illegal pair_style command");
|
||||
|
||||
cut_global = utils::numeric(FLERR,arg[0],false,lmp);
|
||||
cut_global = utils::numeric(FLERR, arg[0], false, lmp);
|
||||
|
||||
// reset cutoffs that have been explicitly set
|
||||
|
||||
if (allocated) {
|
||||
int i,j;
|
||||
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;
|
||||
@ -181,23 +179,22 @@ void PairGauss::settings(int narg, char **arg)
|
||||
|
||||
void PairGauss::coeff(int narg, char **arg)
|
||||
{
|
||||
if (narg < 4 || narg > 5)
|
||||
error->all(FLERR,"Incorrect args for pair coefficients");
|
||||
if (narg < 4 || narg > 5) 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);
|
||||
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 a_one = utils::numeric(FLERR,arg[2],false,lmp);
|
||||
double b_one = utils::numeric(FLERR,arg[3],false,lmp);
|
||||
double a_one = utils::numeric(FLERR, arg[2], false, lmp);
|
||||
double b_one = utils::numeric(FLERR, arg[3], false, lmp);
|
||||
|
||||
double cut_one = cut_global;
|
||||
if (narg == 5) cut_one = utils::numeric(FLERR,arg[4],false,lmp);
|
||||
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++) {
|
||||
for (int j = MAX(jlo, i); j <= jhi; j++) {
|
||||
a[i][j] = a_one;
|
||||
b[i][j] = b_one;
|
||||
cut[i][j] = cut_one;
|
||||
@ -206,7 +203,7 @@ void PairGauss::coeff(int narg, char **arg)
|
||||
}
|
||||
}
|
||||
|
||||
if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients");
|
||||
if (count == 0) error->all(FLERR, "Incorrect args for pair coefficients");
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
@ -218,10 +215,10 @@ double PairGauss::init_one(int i, int j)
|
||||
if (setflag[i][j] == 0) {
|
||||
double sign_bi = (b[i][i] >= 0.0) ? 1.0 : -1.0;
|
||||
double sign_bj = (b[j][j] >= 0.0) ? 1.0 : -1.0;
|
||||
double si = sqrt(0.5/fabs(b[i][i]));
|
||||
double sj = sqrt(0.5/fabs(b[j][j]));
|
||||
double si = sqrt(0.5 / fabs(b[i][i]));
|
||||
double sj = sqrt(0.5 / fabs(b[j][j]));
|
||||
double sij = mix_distance(si, sj);
|
||||
b[i][j] = 0.5 / (sij*sij);
|
||||
b[i][j] = 0.5 / (sij * sij);
|
||||
b[i][j] *= MAX(sign_bi, sign_bj);
|
||||
|
||||
// Negative "a" values are useful for simulating repulsive particles.
|
||||
@ -232,12 +229,14 @@ double PairGauss::init_one(int i, int j)
|
||||
a[i][j] = mix_energy(fabs(a[i][i]), fabs(a[j][j]), si, sj);
|
||||
a[i][j] *= MIN(sign_ai, sign_aj);
|
||||
|
||||
cut[i][j] = mix_distance(cut[i][i],cut[j][j]);
|
||||
cut[i][j] = mix_distance(cut[i][i], cut[j][j]);
|
||||
}
|
||||
|
||||
// cutoff correction to energy
|
||||
if (offset_flag) offset[i][j] = a[i][j]*exp(-b[i][j]*cut[i][j]*cut[i][j]);
|
||||
else offset[i][j] = 0.0;
|
||||
if (offset_flag)
|
||||
offset[i][j] = a[i][j] * exp(-b[i][j] * cut[i][j] * cut[i][j]);
|
||||
else
|
||||
offset[i][j] = 0.0;
|
||||
|
||||
a[j][i] = a[i][j];
|
||||
b[j][i] = b[i][j];
|
||||
@ -254,14 +253,14 @@ void PairGauss::write_restart(FILE *fp)
|
||||
{
|
||||
write_restart_settings(fp);
|
||||
|
||||
int i,j;
|
||||
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);
|
||||
fwrite(&setflag[i][j], sizeof(int), 1, fp);
|
||||
if (setflag[i][j]) {
|
||||
fwrite(&a[i][j],sizeof(double),1,fp);
|
||||
fwrite(&b[i][j],sizeof(double),1,fp);
|
||||
fwrite(&cut[i][j],sizeof(double),1,fp);
|
||||
fwrite(&a[i][j], sizeof(double), 1, fp);
|
||||
fwrite(&b[i][j], sizeof(double), 1, fp);
|
||||
fwrite(&cut[i][j], sizeof(double), 1, fp);
|
||||
}
|
||||
}
|
||||
}
|
||||
@ -275,21 +274,21 @@ void PairGauss::read_restart(FILE *fp)
|
||||
read_restart_settings(fp);
|
||||
allocate();
|
||||
|
||||
int i,j;
|
||||
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 (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,&a[i][j],sizeof(double),1,fp,nullptr,error);
|
||||
utils::sfread(FLERR,&b[i][j],sizeof(double),1,fp,nullptr,error);
|
||||
utils::sfread(FLERR,&cut[i][j],sizeof(double),1,fp,nullptr,error);
|
||||
utils::sfread(FLERR, &a[i][j], sizeof(double), 1, fp, nullptr, error);
|
||||
utils::sfread(FLERR, &b[i][j], sizeof(double), 1, fp, nullptr, error);
|
||||
utils::sfread(FLERR, &cut[i][j], sizeof(double), 1, fp, nullptr, error);
|
||||
}
|
||||
MPI_Bcast(&a[i][j],1,MPI_DOUBLE,0,world);
|
||||
MPI_Bcast(&b[i][j],1,MPI_DOUBLE,0,world);
|
||||
MPI_Bcast(&cut[i][j],1,MPI_DOUBLE,0,world);
|
||||
MPI_Bcast(&a[i][j], 1, MPI_DOUBLE, 0, world);
|
||||
MPI_Bcast(&b[i][j], 1, MPI_DOUBLE, 0, world);
|
||||
MPI_Bcast(&cut[i][j], 1, MPI_DOUBLE, 0, world);
|
||||
}
|
||||
}
|
||||
}
|
||||
@ -300,9 +299,9 @@ void PairGauss::read_restart(FILE *fp)
|
||||
|
||||
void PairGauss::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);
|
||||
fwrite(&cut_global, sizeof(double), 1, fp);
|
||||
fwrite(&offset_flag, sizeof(int), 1, fp);
|
||||
fwrite(&mix_flag, sizeof(int), 1, fp);
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
@ -312,13 +311,13 @@ void PairGauss::write_restart_settings(FILE *fp)
|
||||
void PairGauss::read_restart_settings(FILE *fp)
|
||||
{
|
||||
if (comm->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);
|
||||
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);
|
||||
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);
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
@ -327,8 +326,7 @@ void PairGauss::read_restart_settings(FILE *fp)
|
||||
|
||||
void PairGauss::write_data(FILE *fp)
|
||||
{
|
||||
for (int i = 1; i <= atom->ntypes; i++)
|
||||
fprintf(fp,"%d %g %g\n",i,a[i][i],b[i][i]);
|
||||
for (int i = 1; i <= atom->ntypes; i++) fprintf(fp, "%d %g %g\n", i, a[i][i], b[i][i]);
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
@ -339,18 +337,18 @@ void PairGauss::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 %g\n",i,j,a[i][j],b[i][j],cut[i][j]);
|
||||
fprintf(fp, "%d %d %g %g %g\n", i, j, a[i][j], b[i][j], cut[i][j]);
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
double PairGauss::single(int /*i*/, int /*j*/, int itype, int jtype, double rsq,
|
||||
double /*factor_coul*/, double /*factor_lj*/,
|
||||
double &fforce)
|
||||
double /*factor_coul*/, double factor_lj, double &fforce)
|
||||
{
|
||||
double philj =
|
||||
-(a[itype][jtype]*exp(-b[itype][jtype]*rsq) - offset[itype][jtype]);
|
||||
fforce = -2.0*a[itype][jtype]*b[itype][jtype] * exp(-b[itype][jtype]*rsq);
|
||||
double philj = -(a[itype][jtype] * exp(-b[itype][jtype] * rsq) - offset[itype][jtype]);
|
||||
philj *= factor_lj;
|
||||
fforce = -2.0 * a[itype][jtype] * b[itype][jtype] * exp(-b[itype][jtype] * rsq);
|
||||
fforce *= factor_lj;
|
||||
return philj;
|
||||
}
|
||||
|
||||
@ -359,6 +357,6 @@ double PairGauss::single(int /*i*/, int /*j*/, int itype, int jtype, double rsq,
|
||||
void *PairGauss::extract(const char *str, int &dim)
|
||||
{
|
||||
dim = 2;
|
||||
if (strcmp(str,"a") == 0) return (void *) a;
|
||||
if (strcmp(str, "a") == 0) return (void *) a;
|
||||
return nullptr;
|
||||
}
|
||||
|
||||
@ -1,4 +1,3 @@
|
||||
// clang-format off
|
||||
/* ----------------------------------------------------------------------
|
||||
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
|
||||
https://www.lammps.org/, Sandia National Laboratories
|
||||
@ -26,11 +25,9 @@
|
||||
#include "omp_compat.h"
|
||||
using namespace LAMMPS_NS;
|
||||
|
||||
#define EPSILON 1.0e-10
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
PairGaussOMP::PairGaussOMP(LAMMPS *lmp) :
|
||||
PairGauss(lmp), ThrOMP(lmp, THR_PAIR)
|
||||
PairGaussOMP::PairGaussOMP(LAMMPS *lmp) : PairGauss(lmp), ThrOMP(lmp, THR_PAIR)
|
||||
{
|
||||
suffix_flag |= Suffix::OMP;
|
||||
respa_enable = 0;
|
||||
@ -40,7 +37,7 @@ PairGaussOMP::PairGaussOMP(LAMMPS *lmp) :
|
||||
|
||||
void PairGaussOMP::compute(int eflag, int vflag)
|
||||
{
|
||||
ev_init(eflag,vflag);
|
||||
ev_init(eflag, vflag);
|
||||
|
||||
const int nall = atom->nlocal + atom->nghost;
|
||||
const int nthreads = comm->nthreads;
|
||||
@ -48,7 +45,7 @@ void PairGaussOMP::compute(int eflag, int vflag)
|
||||
double occ = 0.0;
|
||||
|
||||
#if defined(_OPENMP)
|
||||
#pragma omp parallel LMP_DEFAULT_NONE LMP_SHARED(eflag,vflag) reduction(+:occ)
|
||||
#pragma omp parallel LMP_DEFAULT_NONE LMP_SHARED(eflag, vflag) reduction(+ : occ)
|
||||
#endif
|
||||
{
|
||||
int ifrom, ito, tid;
|
||||
@ -60,39 +57,46 @@ void PairGaussOMP::compute(int eflag, int vflag)
|
||||
|
||||
if (evflag) {
|
||||
if (eflag) {
|
||||
if (force->newton_pair) occ = eval<1,1,1>(ifrom, ito, thr);
|
||||
else occ = eval<1,1,0>(ifrom, ito, thr);
|
||||
if (force->newton_pair)
|
||||
occ = eval<1, 1, 1>(ifrom, ito, thr);
|
||||
else
|
||||
occ = eval<1, 1, 0>(ifrom, ito, thr);
|
||||
} else {
|
||||
if (force->newton_pair) occ = eval<1,0,1>(ifrom, ito, thr);
|
||||
else occ = eval<1,0,0>(ifrom, ito, thr);
|
||||
if (force->newton_pair)
|
||||
occ = eval<1, 0, 1>(ifrom, ito, thr);
|
||||
else
|
||||
occ = eval<1, 0, 0>(ifrom, ito, thr);
|
||||
}
|
||||
} else {
|
||||
if (force->newton_pair) occ = eval<0,0,1>(ifrom, ito, thr);
|
||||
else occ = eval<0,0,0>(ifrom, ito, thr);
|
||||
if (force->newton_pair)
|
||||
occ = eval<0, 0, 1>(ifrom, ito, thr);
|
||||
else
|
||||
occ = eval<0, 0, 0>(ifrom, ito, thr);
|
||||
}
|
||||
|
||||
thr->timer(Timer::PAIR);
|
||||
reduce_thr(this, eflag, vflag, thr);
|
||||
} // end of omp parallel region
|
||||
} // end of omp parallel region
|
||||
|
||||
if (eflag_global) pvector[0] = occ;
|
||||
}
|
||||
|
||||
template <int EVFLAG, int EFLAG, int NEWTON_PAIR>
|
||||
double PairGaussOMP::eval(int iifrom, int iito, ThrData * const thr)
|
||||
double PairGaussOMP::eval(int iifrom, int iito, ThrData *const thr)
|
||||
{
|
||||
int i,j,ii,jj,jnum,itype,jtype;
|
||||
double xtmp,ytmp,ztmp,delx,dely,delz,evdwl,fpair,rsq;
|
||||
int *ilist,*jlist,*numneigh,**firstneigh;
|
||||
int i, j, ii, jj, jnum, itype, jtype;
|
||||
double xtmp, ytmp, ztmp, delx, dely, delz, evdwl, fpair, rsq;
|
||||
int *ilist, *jlist, *numneigh, **firstneigh;
|
||||
int occ = 0;
|
||||
|
||||
evdwl = 0.0;
|
||||
|
||||
const auto * _noalias const x = (dbl3_t *) atom->x[0];
|
||||
auto * _noalias const f = (dbl3_t *) thr->get_f()[0];
|
||||
const int * _noalias const type = atom->type;
|
||||
const auto *_noalias const x = (dbl3_t *) atom->x[0];
|
||||
auto *_noalias const f = (dbl3_t *) thr->get_f()[0];
|
||||
const int *_noalias const type = atom->type;
|
||||
const double *_noalias const special_lj = force->special_lj;
|
||||
const int nlocal = atom->nlocal;
|
||||
double fxtmp,fytmp,fztmp;
|
||||
double fxtmp, fytmp, fztmp;
|
||||
|
||||
ilist = list->ilist;
|
||||
numneigh = list->numneigh;
|
||||
@ -109,42 +113,45 @@ double PairGaussOMP::eval(int iifrom, int iito, ThrData * const thr)
|
||||
itype = type[i];
|
||||
jlist = firstneigh[i];
|
||||
jnum = numneigh[i];
|
||||
fxtmp=fytmp=fztmp=0.0;
|
||||
fxtmp = fytmp = fztmp = 0.0;
|
||||
|
||||
for (jj = 0; jj < jnum; jj++) {
|
||||
j = jlist[jj];
|
||||
const double factor_lj = special_lj[sbmask(j)];
|
||||
j &= NEIGHMASK;
|
||||
|
||||
delx = xtmp - x[j].x;
|
||||
dely = ytmp - x[j].y;
|
||||
delz = ztmp - x[j].z;
|
||||
rsq = delx*delx + dely*dely + delz*delz;
|
||||
rsq = delx * delx + dely * dely + delz * delz;
|
||||
jtype = type[j];
|
||||
|
||||
// define a Gaussian well to be occupied if
|
||||
// the site it interacts with is within the force maximum
|
||||
|
||||
if (EFLAG)
|
||||
if (eflag_global && rsq < 0.5/b[itype][jtype]) occ++;
|
||||
if (eflag_global && rsq < 0.5 / b[itype][jtype]) occ++;
|
||||
|
||||
if (rsq < cutsq[itype][jtype]) {
|
||||
fpair = -2.0*a[itype][jtype]*b[itype][jtype]*exp(-b[itype][jtype]*rsq);
|
||||
fpair = -2.0 * a[itype][jtype] * b[itype][jtype] * exp(-b[itype][jtype] * rsq);
|
||||
fpair *= factor_lj;
|
||||
|
||||
fxtmp += delx*fpair;
|
||||
fytmp += dely*fpair;
|
||||
fztmp += delz*fpair;
|
||||
fxtmp += delx * fpair;
|
||||
fytmp += dely * fpair;
|
||||
fztmp += delz * fpair;
|
||||
if (NEWTON_PAIR || j < nlocal) {
|
||||
f[j].x -= delx*fpair;
|
||||
f[j].y -= dely*fpair;
|
||||
f[j].z -= delz*fpair;
|
||||
f[j].x -= delx * fpair;
|
||||
f[j].y -= dely * fpair;
|
||||
f[j].z -= delz * fpair;
|
||||
}
|
||||
|
||||
if (EFLAG)
|
||||
evdwl = -(a[itype][jtype]*exp(-b[itype][jtype]*rsq) -
|
||||
offset[itype][jtype]);
|
||||
if (EFLAG) {
|
||||
evdwl = -(a[itype][jtype] * exp(-b[itype][jtype] * rsq) - offset[itype][jtype]);
|
||||
evdwl *= factor_lj;
|
||||
}
|
||||
|
||||
if (EVFLAG) ev_tally_thr(this, i,j,nlocal,NEWTON_PAIR,
|
||||
evdwl,0.0,fpair,delx,dely,delz,thr);
|
||||
if (EVFLAG)
|
||||
ev_tally_thr(this, i, j, nlocal, NEWTON_PAIR, evdwl, 0.0, fpair, delx, dely, delz, thr);
|
||||
}
|
||||
}
|
||||
f[i].x += fxtmp;
|
||||
|
||||
@ -1,6 +1,7 @@
|
||||
---
|
||||
lammps_version: 17 Feb 2022
|
||||
date_generated: Fri Mar 18 22:17:30 2022
|
||||
lammps_version: 8 Feb 2023
|
||||
tags: generated
|
||||
date_generated: Fri Feb 10 04:35:42 2023
|
||||
epsilon: 5e-14
|
||||
skip_tests:
|
||||
prerequisites: ! |
|
||||
@ -21,72 +22,72 @@ pair_coeff: ! |
|
||||
extract: ! |
|
||||
a 2
|
||||
natoms: 29
|
||||
init_vdwl: -0.3864209232155955
|
||||
init_vdwl: -0.10744797594740513
|
||||
init_coul: 0
|
||||
init_stress: ! |-
|
||||
-2.2163698035273618e-01 -2.4625537112111151e-01 -2.1305188258578470e-01 6.1866130132872003e-02 6.6571531323592207e-03 2.4462489698253043e-02
|
||||
-1.1008516714247368e-01 -1.2186754451867259e-01 -9.8542625674496995e-02 3.1449705846190668e-02 7.6089266254867642e-03 1.8309130908749447e-02
|
||||
init_forces: ! |2
|
||||
1 6.6420790881819135e-03 -2.2656352315504419e-02 -1.1151371492916272e-02
|
||||
2 -4.2811279269361156e-03 -6.1319457093740913e-03 3.5655172207420469e-03
|
||||
3 1.4122450625899436e-02 -6.4669361374405238e-03 -3.4547856578885819e-03
|
||||
4 9.2022857012287029e-03 -9.6636513433784451e-04 5.6358203169252627e-03
|
||||
5 4.8130601322131554e-03 3.8754733051119977e-03 -9.5295063547545852e-03
|
||||
6 1.1694684866979214e-02 -1.1545676479204194e-04 8.9573160065516636e-03
|
||||
7 9.2903486979032928e-03 2.0412909476071818e-03 2.0895988031881762e-02
|
||||
8 8.2241496648224888e-03 3.2282787725485093e-03 -5.5984923599097887e-03
|
||||
9 -9.1025802092398352e-04 -2.9695228458060277e-03 -1.1310732844631449e-02
|
||||
10 -3.3216915418810807e-03 1.5165769941583293e-02 -1.2109299141022656e-03
|
||||
11 2.8660819158744734e-03 8.4275110632676224e-03 9.5049639732687780e-03
|
||||
12 -2.2598101552834836e-02 -4.9601916039375389e-03 8.4003031140324341e-03
|
||||
13 -9.4516726881566388e-03 1.4305285400726016e-03 1.5058333996407797e-03
|
||||
14 -3.4864885985541323e-04 -1.8689531090578431e-03 1.0873765307790662e-02
|
||||
15 -4.3507528925587919e-03 -9.8039134533324957e-03 -1.2974821049077257e-03
|
||||
16 -1.1020969725836332e-02 1.9784791088106480e-02 -1.1275809663597147e-02
|
||||
17 -3.3419871860029715e-03 8.2289097839403427e-03 -1.8738857919511739e-02
|
||||
18 -4.3521806872305014e-03 -7.8656826681131678e-03 6.5687503046006051e-03
|
||||
19 9.2492265214090606e-04 7.1711662533072390e-04 5.6458533331792311e-04
|
||||
20 -1.0085043451939246e-03 -9.3005491274762285e-04 5.0547579022850204e-05
|
||||
21 -1.1667863590836469e-03 4.2902548011202451e-03 -2.0499762074924696e-03
|
||||
22 3.1557118123179237e-03 9.5184863190457439e-04 3.9782856889301653e-03
|
||||
23 -4.4281270696133268e-03 -2.5346995221498434e-03 7.5785922268409691e-04
|
||||
24 -3.0622009464357529e-03 1.9851728054262412e-03 -6.5546361368900791e-03
|
||||
25 3.7063223755402488e-03 -2.7523325480129519e-04 2.6231087574901247e-03
|
||||
26 -2.9538296443943910e-03 -4.1309366766066275e-03 -6.6793595531402768e-04
|
||||
27 1.3703266730383692e-03 6.2247400367154799e-03 -2.6349027873549895e-03
|
||||
28 4.1766990025931270e-03 -1.2043760854333907e-03 2.6647558305012602e-03
|
||||
29 -3.5922837617955406e-03 -3.4710661493005325e-03 -1.0719806881093012e-03
|
||||
run_vdwl: -0.3863508651617164
|
||||
1 6.3967388935615935e-03 -8.9571493541161559e-03 -2.6243495367571021e-03
|
||||
2 -5.8925997783387395e-04 -1.6467452362593446e-03 3.9474993583852062e-04
|
||||
3 6.4541909782685350e-03 -2.7518790414616788e-03 -1.1988207163740401e-03
|
||||
4 1.7173451287781348e-03 -1.0661558546877452e-04 9.4766143245350446e-04
|
||||
5 1.3640604274609371e-03 5.2189464017124855e-04 -1.7289542395764162e-03
|
||||
6 7.5053558116135836e-03 -7.8008442916635047e-04 2.2571851729790136e-03
|
||||
7 8.3795968749890409e-03 8.4698140531401750e-05 7.9566896360105874e-03
|
||||
8 2.9576378106382336e-03 2.7674589237502013e-03 -1.3800458960768092e-03
|
||||
9 -6.8101137697709914e-05 -3.4479763208106613e-04 -2.5941356000719553e-03
|
||||
10 -2.5606662761393502e-03 5.4931707692542301e-03 -9.3119684414857094e-04
|
||||
11 4.8953489579958687e-04 1.9131103445091924e-03 1.6644396272410174e-03
|
||||
12 -1.1675327637667321e-02 9.5973244924826004e-04 1.8760227478274104e-03
|
||||
13 -1.9229408858018125e-03 -1.5273825100896492e-04 5.5179749722991656e-04
|
||||
14 -1.2828240307888158e-03 -3.8245527716947414e-04 2.7609608498598013e-03
|
||||
15 -2.6060442482268012e-03 -2.1186266400159987e-03 2.3778251062230583e-04
|
||||
16 -3.3980337566683412e-03 5.9576900525289264e-03 -4.5705974537114087e-03
|
||||
17 -3.9313325721691108e-03 5.7862524954095653e-03 -7.8476500647319345e-03
|
||||
18 -4.4024833272868361e-03 -8.0199476070002811e-03 7.0506526853925692e-03
|
||||
19 1.7435858376498004e-04 1.2618943785668791e-04 1.4777623368941996e-04
|
||||
20 -2.0763763676166395e-04 -1.8486278638647300e-04 -1.4545702140610665e-05
|
||||
21 -2.3094929681373769e-03 2.8627074594353183e-03 2.2096366480569823e-03
|
||||
22 4.4288112377976487e-04 1.5533647577931545e-04 4.5374867545560895e-04
|
||||
23 -5.7258977202143764e-04 -3.1064002433965675e-04 2.2783380609201144e-05
|
||||
24 -2.0888998040544597e-03 -1.8207426839036761e-03 -4.2932649963137327e-03
|
||||
25 1.9475456066254769e-04 -1.1833520320337950e-04 -1.6748792244510374e-04
|
||||
26 -4.1556297189798295e-04 -4.8191923887462638e-04 -1.3871041595514520e-04
|
||||
27 1.8904879208478888e-03 2.0114948956210790e-03 -1.1958555043652932e-03
|
||||
28 5.4016949423872510e-04 -8.4138794056923025e-05 3.2483661978743096e-04
|
||||
29 -4.7591550125065788e-04 -3.7805829958259914e-04 -1.7110876038516800e-04
|
||||
run_vdwl: -0.10742358116477577
|
||||
run_coul: 0
|
||||
run_stress: ! |-
|
||||
-2.2160110236171965e-01 -2.4617131296379000e-01 -2.1305840663314177e-01 6.1838681126107659e-02 6.6381976968359312e-03 2.4402003044659750e-02
|
||||
-1.1006082820490694e-01 -1.2183895033722210e-01 -9.8544056061287286e-02 3.1434066108895035e-02 7.6054332065046241e-03 1.8296316078185181e-02
|
||||
run_forces: ! |2
|
||||
1 6.6413062980785250e-03 -2.2643146328241739e-02 -1.1142942591208689e-02
|
||||
2 -4.2827175516235612e-03 -6.1293692790411374e-03 3.5544412327760560e-03
|
||||
3 1.4119948940250827e-02 -6.4770130378595647e-03 -3.4543300745505934e-03
|
||||
4 9.1950674649061044e-03 -9.6181880037213058e-04 5.6407327955521839e-03
|
||||
5 4.8139605807970181e-03 3.8756785373474725e-03 -9.5265345859732731e-03
|
||||
6 1.1686087085798024e-02 -1.2058408819619352e-04 8.9440737448735885e-03
|
||||
7 9.2903951061620062e-03 2.0470436506746757e-03 2.0893971575447097e-02
|
||||
8 8.2293078546638102e-03 3.2320594356858930e-03 -5.5960424447553192e-03
|
||||
9 -9.0531229139671548e-04 -2.9682892647793285e-03 -1.1310091180000625e-02
|
||||
10 -3.3259540649788755e-03 1.5162058066139865e-02 -1.2051475979010121e-03
|
||||
11 2.8653459452275495e-03 8.4232873380028633e-03 9.4992751088281042e-03
|
||||
12 -2.2597913046721457e-02 -4.9554168226113669e-03 8.4125389603604290e-03
|
||||
13 -9.4443544519883221e-03 1.4213463221588211e-03 1.5001514052453023e-03
|
||||
14 -3.4584663474122380e-04 -1.8589408798251432e-03 1.0867482270335016e-02
|
||||
15 -4.3501928045490360e-03 -9.8084491667171560e-03 -1.3028078079482534e-03
|
||||
16 -1.1016280904427918e-02 1.9779721304587411e-02 -1.1270569619019189e-02
|
||||
17 -3.3423858079814819e-03 8.2235538802195984e-03 -1.8736744422015096e-02
|
||||
18 -4.3608515422088765e-03 -7.8746354806103706e-03 6.5742308189911851e-03
|
||||
19 9.1906994708977803e-04 7.1287404601029134e-04 5.6400074161977465e-04
|
||||
20 -9.9447892889066537e-04 -9.1770728964396954e-04 4.8034421345478542e-05
|
||||
21 -1.1680763991239612e-03 4.2792686042234344e-03 -2.0399769668755026e-03
|
||||
22 3.1591713175596393e-03 9.5574540137949749e-04 3.9720638190945401e-03
|
||||
23 -4.4307388063199804e-03 -2.5263233108993562e-03 7.5473099454995423e-04
|
||||
24 -3.0620441606671513e-03 1.9754767097317143e-03 -6.5470110878329164e-03
|
||||
25 3.7077221581000219e-03 -2.6401622505421555e-04 2.6260279440605853e-03
|
||||
26 -2.9551493935843225e-03 -4.1316756613083479e-03 -6.7756590474090739e-04
|
||||
27 1.3660170615493567e-03 6.2236935799399491e-03 -2.6252513993441232e-03
|
||||
28 4.1796626204497374e-03 -1.2023509872189564e-03 2.6623174423427176e-03
|
||||
29 -3.5907655914288470e-03 -3.4720702537225047e-03 -1.0790575932565119e-03
|
||||
1 6.3951182607409262e-03 -8.9545730680437412e-03 -2.6242747559400842e-03
|
||||
2 -5.8940132985142637e-04 -1.6432865472356978e-03 3.9314816550261645e-04
|
||||
3 6.4493680135179910e-03 -2.7508136736323781e-03 -1.1985914252700854e-03
|
||||
4 1.7157340254504952e-03 -1.0645373608815925e-04 9.4797119947664693e-04
|
||||
5 1.3638703830031701e-03 5.2199495314776776e-04 -1.7280157917857904e-03
|
||||
6 7.5027803029913730e-03 -7.7992545467626452e-04 2.2542580363439260e-03
|
||||
7 8.3790973418566787e-03 8.6964711023362301e-05 7.9549476717222789e-03
|
||||
8 2.9600170939576595e-03 2.7677876464161800e-03 -1.3796592742572695e-03
|
||||
9 -6.4851944444991254e-05 -3.4672203705875542e-04 -2.5935699450003449e-03
|
||||
10 -2.5611753535538561e-03 5.4925092261102065e-03 -9.3045690231003775e-04
|
||||
11 4.8983736198662245e-04 1.9115748764361731e-03 1.6629308583155874e-03
|
||||
12 -1.1670568560250798e-02 9.5936201106161740e-04 1.8766858864372907e-03
|
||||
13 -1.9194289686350334e-03 -1.5306898539863079e-04 5.4961829174280503e-04
|
||||
14 -1.2795546260879820e-03 -3.8095864067789689e-04 2.7581590326053668e-03
|
||||
15 -2.6123992882847165e-03 -2.1212782856001471e-03 2.3834468823455228e-04
|
||||
16 -3.3974822936998559e-03 5.9551171086075615e-03 -4.5676348427171106e-03
|
||||
17 -3.9304986491614574e-03 5.7834906955578757e-03 -7.8464040391146503e-03
|
||||
18 -4.4036759081377110e-03 -8.0214055777730296e-03 7.0529490020863929e-03
|
||||
19 1.7246101738774175e-04 1.2467198317308374e-04 1.4771603949694447e-04
|
||||
20 -2.0504560588348406e-04 -1.8273510588466313e-04 -1.4399117733621655e-05
|
||||
21 -2.3100182076283716e-03 2.8627603983808606e-03 2.2112169537780555e-03
|
||||
22 4.4298531465963761e-04 1.5549858717522525e-04 4.5294060847149182e-04
|
||||
23 -5.7261101837172428e-04 -3.0956825334589567e-04 2.2660319228822179e-05
|
||||
24 -2.0891130362891692e-03 -1.8217407569914542e-03 -4.2926590911060896e-03
|
||||
25 1.9572941661502480e-04 -1.1608017944239833e-04 -1.6583785661338941e-04
|
||||
26 -4.1608783560330648e-04 -4.8239423005540976e-04 -1.4005216918917679e-04
|
||||
27 1.8902324250959931e-03 2.0113877029906300e-03 -1.1947874546904412e-03
|
||||
28 5.4065339019595659e-04 -8.3834857585872457e-05 3.2474503896207270e-04
|
||||
29 -4.7597172157538961e-04 -3.7828051059014876e-04 -1.7194912667675645e-04
|
||||
...
|
||||
|
||||
Reference in New Issue
Block a user