diff --git a/doc/src/pair_gauss.rst b/doc/src/pair_gauss.rst index 091bd52d45..f6d8134467 100644 --- a/doc/src/pair_gauss.rst +++ b/doc/src/pair_gauss.rst @@ -137,16 +137,19 @@ interaction. The :doc:`pair_modify ` table and tail options are not relevant for these pair styles. -These pair styles write their 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. +These pair styles write their 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. These pair styles can only be used via the *pair* keyword of the :doc:`run_style respa ` 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 ` 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 ` 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 ` page for more info. -The *gauss* style does not apply :doc:`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 `) that do -apply those factors. +.. versionchanged:: TBD + +Prior to this version, the *gauss* pair style did not apply +:doc:`special_bonds ` factors. Related commands """""""""""""""" diff --git a/lib/gpu/lal_gauss.cpp b/lib/gpu/lal_gauss.cpp index 6d8f0f02aa..5a7bd24b85 100644 --- a/lib/gpu/lal_gauss.cpp +++ b/lib/gpu/lal_gauss.cpp @@ -133,13 +133,13 @@ int GaussT::loop(const int eflag, const int vflag) { this->time_pair.start(); if (shared_types) { this->k_pair_sel->set_size(GX,BX); - this->k_pair_sel->run(&this->atom->x, &gauss1, + this->k_pair_sel->run(&this->atom->x, &gauss1, &sp_lj, &this->nbor->dev_nbor, &this->_nbor_data->begin(), &this->ans->force, &this->ans->engv, &eflag, &vflag, &ainum, &nbor_pitch, &this->_threads_per_atom); } else { this->k_pair.set_size(GX,BX); - this->k_pair.run(&this->atom->x, &gauss1, &_lj_types, + this->k_pair.run(&this->atom->x, &gauss1, &_lj_types, &sp_lj, &this->nbor->dev_nbor, &this->_nbor_data->begin(), &this->ans->force, &this->ans->engv, &eflag, &vflag, &ainum, &nbor_pitch, &this->_threads_per_atom); diff --git a/lib/gpu/lal_gauss.cu b/lib/gpu/lal_gauss.cu index 2540b8492f..cb6b72db30 100644 --- a/lib/gpu/lal_gauss.cu +++ b/lib/gpu/lal_gauss.cu @@ -27,6 +27,7 @@ _texture_2d( pos_tex,int4); __kernel void k_gauss(const __global numtyp4 *restrict x_, const __global numtyp4 *restrict gauss1, const int lj_types, + const __global numtyp *restrict sp_lj, const __global int *dev_nbor, const __global int *dev_packed, __global acctyp4 *restrict ans, @@ -56,9 +57,11 @@ __kernel void k_gauss(const __global numtyp4 *restrict x_, numtyp4 ix; fetch4(ix,i,pos_tex); //x_[i]; int itype=ix.w; + numtyp factor_lj; for ( ; nbor -#include #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 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; } diff --git a/src/OPENMP/pair_gauss_omp.cpp b/src/OPENMP/pair_gauss_omp.cpp index 12c197e768..18263d4d15 100644 --- a/src/OPENMP/pair_gauss_omp.cpp +++ b/src/OPENMP/pair_gauss_omp.cpp @@ -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 -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; diff --git a/unittest/force-styles/tests/mol-pair-gauss.yaml b/unittest/force-styles/tests/mol-pair-gauss.yaml index 5df85e0c3f..88ea376d5f 100644 --- a/unittest/force-styles/tests/mol-pair-gauss.yaml +++ b/unittest/force-styles/tests/mol-pair-gauss.yaml @@ -1,6 +1,6 @@ --- -lammps_version: 17 Feb 2022 -date_generated: Fri Mar 18 22:17:30 2022 +lammps_version: 8 Feb 2023 +date_generated: Fri Feb 10 04:35:42 2023 epsilon: 5e-14 skip_tests: prerequisites: ! | @@ -21,72 +21,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 ...