From d2b5f557379fdd8f7abb4ed6886aeb4028968084 Mon Sep 17 00:00:00 2001 From: EiPi Fun Date: Sun, 1 Sep 2024 16:07:52 +0800 Subject: [PATCH 01/29] Copy and rename base files --- .../pair_hbond_dreiding_lj_angleoffset.cpp | 572 ++++++++++++++++++ .../pair_hbond_dreiding_lj_angleoffset.h | 66 ++ .../pair_hbond_dreiding_morse_angleoffset.cpp | 474 +++++++++++++++ .../pair_hbond_dreiding_morse_angleoffset.h | 40 ++ ...pair_hbond_dreiding_lj_angleoffset_omp.cpp | 317 ++++++++++ .../pair_hbond_dreiding_lj_angleoffset_omp.h | 52 ++ ...r_hbond_dreiding_morse_angleoffset_omp.cpp | 316 ++++++++++ ...air_hbond_dreiding_morse_angleoffset_omp.h | 52 ++ 8 files changed, 1889 insertions(+) create mode 100644 src/EXTRA-MOLECULE/pair_hbond_dreiding_lj_angleoffset.cpp create mode 100644 src/EXTRA-MOLECULE/pair_hbond_dreiding_lj_angleoffset.h create mode 100644 src/EXTRA-MOLECULE/pair_hbond_dreiding_morse_angleoffset.cpp create mode 100644 src/EXTRA-MOLECULE/pair_hbond_dreiding_morse_angleoffset.h create mode 100644 src/OPENMP/pair_hbond_dreiding_lj_angleoffset_omp.cpp create mode 100644 src/OPENMP/pair_hbond_dreiding_lj_angleoffset_omp.h create mode 100644 src/OPENMP/pair_hbond_dreiding_morse_angleoffset_omp.cpp create mode 100644 src/OPENMP/pair_hbond_dreiding_morse_angleoffset_omp.h diff --git a/src/EXTRA-MOLECULE/pair_hbond_dreiding_lj_angleoffset.cpp b/src/EXTRA-MOLECULE/pair_hbond_dreiding_lj_angleoffset.cpp new file mode 100644 index 0000000000..274f8bc2a3 --- /dev/null +++ b/src/EXTRA-MOLECULE/pair_hbond_dreiding_lj_angleoffset.cpp @@ -0,0 +1,572 @@ +// clang-format off +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + https://www.lammps.org/, Sandia National Laboratories + LAMMPS development team: developers@lammps.org + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +/* ---------------------------------------------------------------------- + Contributing author: Tod A Pascal (Caltech) +------------------------------------------------------------------------- */ + +#include "pair_hbond_dreiding_lj.h" + +#include "atom.h" +#include "atom_vec.h" +#include "domain.h" +#include "error.h" +#include "force.h" +#include "math_const.h" +#include "math_special.h" +#include "memory.h" +#include "molecule.h" +#include "neigh_list.h" +#include "neighbor.h" + +#include +#include + +using namespace LAMMPS_NS; +using namespace MathConst; +using namespace MathSpecial; + +static constexpr double SMALL = 0.001; +static constexpr int CHUNK = 8; + +/* ---------------------------------------------------------------------- */ + +PairHbondDreidingLJ::PairHbondDreidingLJ(LAMMPS *lmp) : Pair(lmp) +{ + // hbond cannot compute virial as F dot r + // due to using map() to find bonded H atoms which are not near donor atom + + no_virial_fdotr_compute = 1; + restartinfo = 0; + + nparams = maxparam = 0; + params = nullptr; + + nextra = 2; + pvector = new double[2]; +} + +/* ---------------------------------------------------------------------- */ + +PairHbondDreidingLJ::~PairHbondDreidingLJ() +{ + memory->sfree(params); + delete [] pvector; + + if (allocated) { + memory->destroy(setflag); + memory->destroy(cutsq); + + delete [] donor; + delete [] acceptor; + memory->destroy(type2param); + } +} + +/* ---------------------------------------------------------------------- */ + +void PairHbondDreidingLJ::compute(int eflag, int vflag) +{ + int i,j,k,m,ii,jj,kk,inum,jnum,knum,itype,jtype,ktype,iatom,imol; + tagint tagprev; + double delx,dely,delz,rsq,rsq1,rsq2,r1,r2; + double factor_hb,force_angle,force_kernel,evdwl,eng_lj,ehbond,force_switch; + double c,s,a,b,ac,a11,a12,a22,vx1,vx2,vy1,vy2,vz1,vz2,d; + double fi[3],fj[3],delr1[3],delr2[3]; + double r2inv,r10inv; + double switch1,switch2; + int *ilist,*jlist,*numneigh,**firstneigh; + tagint *klist; + + evdwl = ehbond = 0.0; + ev_init(eflag,vflag); + + double **x = atom->x; + double **f = atom->f; + tagint *tag = atom->tag; + int *molindex = atom->molindex; + int *molatom = atom->molatom; + tagint **special = atom->special; + int **nspecial = atom->nspecial; + int *type = atom->type; + double *special_lj = force->special_lj; + int molecular = atom->molecular; + Molecule **onemols = atom->avec->onemols; + + inum = list->inum; + ilist = list->ilist; + numneigh = list->numneigh; + firstneigh = list->firstneigh; + + // ii = loop over donors + // jj = loop over acceptors + // kk = loop over hydrogens bonded to donor + + int hbcount = 0; + + for (ii = 0; ii < inum; ii++) { + i = ilist[ii]; + itype = type[i]; + if (!donor[itype]) continue; + if (molecular == Atom::MOLECULAR) { + klist = special[i]; + knum = nspecial[i][0]; + } else { + if (molindex[i] < 0) continue; + imol = molindex[i]; + iatom = molatom[i]; + klist = onemols[imol]->special[iatom]; + knum = onemols[imol]->nspecial[iatom][0]; + tagprev = tag[i] - iatom - 1; + } + jlist = firstneigh[i]; + jnum = numneigh[i]; + + for (jj = 0; jj < jnum; jj++) { + j = jlist[jj]; + factor_hb = special_lj[sbmask(j)]; + j &= NEIGHMASK; + + jtype = type[j]; + if (!acceptor[jtype]) continue; + + delx = x[i][0] - x[j][0]; + dely = x[i][1] - x[j][1]; + delz = x[i][2] - x[j][2]; + rsq = delx*delx + dely*dely + delz*delz; + + for (kk = 0; kk < knum; kk++) { + if (molecular == Atom::MOLECULAR) k = atom->map(klist[kk]); + else k = atom->map(klist[kk]+tagprev); + if (k < 0) continue; + ktype = type[k]; + m = type2param[itype][jtype][ktype]; + if (m < 0) continue; + const Param &pm = params[m]; + + if (rsq < pm.cut_outersq) { + delr1[0] = x[i][0] - x[k][0]; + delr1[1] = x[i][1] - x[k][1]; + delr1[2] = x[i][2] - x[k][2]; + domain->minimum_image(delr1); + rsq1 = delr1[0]*delr1[0] + delr1[1]*delr1[1] + delr1[2]*delr1[2]; + r1 = sqrt(rsq1); + + delr2[0] = x[j][0] - x[k][0]; + delr2[1] = x[j][1] - x[k][1]; + delr2[2] = x[j][2] - x[k][2]; + domain->minimum_image(delr2); + rsq2 = delr2[0]*delr2[0] + delr2[1]*delr2[1] + delr2[2]*delr2[2]; + r2 = sqrt(rsq2); + + // angle (cos and sin) + + c = delr1[0]*delr2[0] + delr1[1]*delr2[1] + delr1[2]*delr2[2]; + c /= r1*r2; + if (c > 1.0) c = 1.0; + if (c < -1.0) c = -1.0; + ac = acos(c); + + if (ac > pm.cut_angle && ac < (2.0*MY_PI - pm.cut_angle)) { + s = sqrt(1.0 - c*c); + if (s < SMALL) s = SMALL; + + // LJ-specific kernel + + r2inv = 1.0/rsq; + r10inv = r2inv*r2inv*r2inv*r2inv*r2inv; + force_kernel = r10inv*(pm.lj1*r2inv - pm.lj2)*r2inv * + powint(c,pm.ap); + force_angle = pm.ap * r10inv*(pm.lj3*r2inv - pm.lj4) * + powint(c,pm.ap-1)*s; + + eng_lj = r10inv*(pm.lj3*r2inv - pm.lj4); + + force_switch=0.0; + + if (rsq > pm.cut_innersq) { + switch1 = (pm.cut_outersq-rsq) * (pm.cut_outersq-rsq) * + (pm.cut_outersq + 2.0*rsq - 3.0*pm.cut_innersq) / + pm.denom_vdw; + switch2 = 12.0*rsq * (pm.cut_outersq-rsq) * + (rsq-pm.cut_innersq) / pm.denom_vdw; + + force_kernel *= switch1; + force_angle *= switch1; + force_switch = eng_lj*switch2/rsq; + eng_lj *= switch1; + } + + if (eflag) { + evdwl = eng_lj * powint(c,pm.ap); + evdwl *= factor_hb; + ehbond += evdwl; + } + + a = factor_hb*force_angle/s; + b = factor_hb*force_kernel; + d = factor_hb*force_switch; + + a11 = a*c / rsq1; + a12 = -a / (r1*r2); + a22 = a*c / rsq2; + + vx1 = a11*delr1[0] + a12*delr2[0]; + vx2 = a22*delr2[0] + a12*delr1[0]; + vy1 = a11*delr1[1] + a12*delr2[1]; + vy2 = a22*delr2[1] + a12*delr1[1]; + vz1 = a11*delr1[2] + a12*delr2[2]; + vz2 = a22*delr2[2] + a12*delr1[2]; + + fi[0] = vx1 + b*delx + d*delx; + fi[1] = vy1 + b*dely + d*dely; + fi[2] = vz1 + b*delz + d*delz; + fj[0] = vx2 - b*delx - d*delx; + fj[1] = vy2 - b*dely - d*dely; + fj[2] = vz2 - b*delz - d*delz; + + f[i][0] += fi[0]; + f[i][1] += fi[1]; + f[i][2] += fi[2]; + + f[j][0] += fj[0]; + f[j][1] += fj[1]; + f[j][2] += fj[2]; + + f[k][0] -= vx1 + vx2; + f[k][1] -= vy1 + vy2; + f[k][2] -= vz1 + vz2; + + // KIJ instead of IJK b/c delr1/delr2 are both with respect to k + + if (evflag) ev_tally3(k,i,j,evdwl,0.0,fi,fj,delr1,delr2); + + hbcount++; + } + } + } + } + } + + if (eflag_global) { + pvector[0] = hbcount; + pvector[1] = ehbond; + } +} + +/* ---------------------------------------------------------------------- + allocate all arrays +------------------------------------------------------------------------- */ + +void PairHbondDreidingLJ::allocate() +{ + allocated = 1; + int n = atom->ntypes; + + // mark all setflag as set, since don't require pair_coeff of all I,J + + memory->create(setflag,n+1,n+1,"pair:setflag"); + for (int i = 1; i <= n; i++) + for (int j = i; j <= n; j++) + setflag[i][j] = 1; + + memory->create(cutsq,n+1,n+1,"pair:cutsq"); + + donor = new int[n+1]; + acceptor = new int[n+1]; + memory->create(type2param,n+1,n+1,n+1,"pair:type2param"); + + int i,j,k; + for (i = 1; i <= n; i++) + for (j = 1; j <= n; j++) + for (k = 1; k <= n; k++) + type2param[i][j][k] = -1; +} + +/* ---------------------------------------------------------------------- + global settings +------------------------------------------------------------------------- */ + +void PairHbondDreidingLJ::settings(int narg, char **arg) +{ + if (narg != 4) error->all(FLERR,"Illegal pair_style command"); + + ap_global = utils::inumeric(FLERR,arg[0],false,lmp); + cut_inner_global = utils::numeric(FLERR,arg[1],false,lmp); + cut_outer_global = utils::numeric(FLERR,arg[2],false,lmp); + cut_angle_global = utils::numeric(FLERR,arg[3],false,lmp) * MY_PI/180.0; +} + +/* ---------------------------------------------------------------------- + set coeffs for one or more type pairs +------------------------------------------------------------------------- */ + +void PairHbondDreidingLJ::coeff(int narg, char **arg) +{ + if (narg < 6 || narg > 10) + error->all(FLERR,"Incorrect args for pair coefficients"); + if (!allocated) allocate(); + + int ilo,ihi,jlo,jhi,klo,khi; + utils::bounds(FLERR, arg[0], 1, atom->ntypes, ilo, ihi, error); + utils::bounds(FLERR, arg[1], 1, atom->ntypes, jlo, jhi, error); + utils::bounds_typelabel(FLERR, arg[2], 1, atom->ntypes, klo, khi, lmp, Atom::ATOM); + + int donor_flag; + if (strcmp(arg[3],"i") == 0) donor_flag = 0; + else if (strcmp(arg[3],"j") == 0) donor_flag = 1; + else error->all(FLERR,"Incorrect args for pair coefficients"); + + double epsilon_one = utils::numeric(FLERR, arg[4], false, lmp); + double sigma_one = utils::numeric(FLERR, arg[5], false, lmp); + + int ap_one = ap_global; + if (narg > 6) ap_one = utils::inumeric(FLERR, arg[6], false, lmp); + double cut_inner_one = cut_inner_global; + double cut_outer_one = cut_outer_global; + if (narg > 8) { + cut_inner_one = utils::numeric(FLERR, arg[7], false, lmp); + cut_outer_one = utils::numeric(FLERR, arg[8], false, lmp); + } + if (cut_inner_one>cut_outer_one) + error->all(FLERR,"Pair inner cutoff >= Pair outer cutoff"); + double cut_angle_one = cut_angle_global; + if (narg == 10) cut_angle_one = utils::numeric(FLERR, arg[9], false, lmp) * MY_PI/180.0; + // grow params array if necessary + + if (nparams == maxparam) { + maxparam += CHUNK; + params = (Param *) memory->srealloc(params, maxparam*sizeof(Param), + "pair:params"); + + // make certain all addional allocated storage is initialized + // to avoid false positives when checking with valgrind + + memset(params + nparams, 0, CHUNK*sizeof(Param)); + } + + params[nparams].epsilon = epsilon_one; + params[nparams].sigma = sigma_one; + params[nparams].ap = ap_one; + params[nparams].cut_inner = cut_inner_one; + params[nparams].cut_outer = cut_outer_one; + params[nparams].cut_innersq = cut_inner_one*cut_inner_one; + params[nparams].cut_outersq = cut_outer_one*cut_outer_one; + params[nparams].cut_angle = cut_angle_one; + params[nparams].denom_vdw = + (params[nparams].cut_outersq-params[nparams].cut_innersq) * + (params[nparams].cut_outersq-params[nparams].cut_innersq) * + (params[nparams].cut_outersq-params[nparams].cut_innersq); + + // flag type2param with either i,j = D,A or j,i = D,A + + int count = 0; + for (int i = ilo; i <= ihi; i++) + for (int j = MAX(jlo,i); j <= jhi; j++) + for (int k = klo; k <= khi; k++) { + if (donor_flag == 0) type2param[i][j][k] = nparams; + else type2param[j][i][k] = nparams; + count++; + } + nparams++; + + if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients"); +} + +/* ---------------------------------------------------------------------- + init specific to this pair style +------------------------------------------------------------------------- */ + +void PairHbondDreidingLJ::init_style() +{ + // molecular system required to use special list to find H atoms + // tags required to use special list + // pair newton on required since are looping over D atoms + // and computing forces on A,H which may be on different procs + + if (atom->molecular == Atom::ATOMIC) + error->all(FLERR,"Pair style hbond/dreiding requires molecular system"); + if (atom->tag_enable == 0) + error->all(FLERR,"Pair style hbond/dreiding requires atom IDs"); + if (atom->map_style == Atom::MAP_NONE) + error->all(FLERR,"Pair style hbond/dreiding requires an atom map, " + "see atom_modify"); + if (force->newton_pair == 0) + error->all(FLERR,"Pair style hbond/dreiding requires newton pair on"); + + // set donor[M]/acceptor[M] if any atom of type M is a donor/acceptor + + int anyflag = 0; + int n = atom->ntypes; + for (int m = 1; m <= n; m++) donor[m] = acceptor[m] = 0; + for (int i = 1; i <= n; i++) + for (int j = 1; j <= n; j++) + for (int k = 1; k <= n; k++) + if (type2param[i][j][k] >= 0) { + anyflag = 1; + donor[i] = 1; + acceptor[j] = 1; + } + + if (!anyflag) error->all(FLERR,"No pair hbond/dreiding coefficients set"); + + // set additional param values + // offset is for LJ only, angle term is not included + + for (int m = 0; m < nparams; m++) { + params[m].lj1 = 60.0*params[m].epsilon*pow(params[m].sigma,12.0); + params[m].lj2 = 60.0*params[m].epsilon*pow(params[m].sigma,10.0); + params[m].lj3 = 5.0*params[m].epsilon*pow(params[m].sigma,12.0); + params[m].lj4 = 6.0*params[m].epsilon*pow(params[m].sigma,10.0); + + /* + if (offset_flag) { + double ratio = params[m].sigma / params[m].cut_outer; + params[m].offset = params[m].epsilon * + ((2.0*pow(ratio,9.0)) - (3.0*pow(ratio,6.0))); + } else params[m].offset = 0.0; + */ + } + + // full neighbor list request + + neighbor->add_request(this, NeighConst::REQ_FULL); +} + +/* ---------------------------------------------------------------------- + init for one type pair i,j and corresponding j,i +------------------------------------------------------------------------- */ + +double PairHbondDreidingLJ::init_one(int i, int j) +{ + int m; + + // return maximum cutoff for any K with I,J = D,A or J,I = D,A + // donor/acceptor is not symmetric, IJ interaction != JI interaction + + double cut = 0.0; + for (int k = 1; k <= atom->ntypes; k++) { + m = type2param[i][j][k]; + if (m >= 0) cut = MAX(cut,params[m].cut_outer); + m = type2param[j][i][k]; + if (m >= 0) cut = MAX(cut,params[m].cut_outer); + } + return cut; +} + +/* ---------------------------------------------------------------------- */ + +double PairHbondDreidingLJ::single(int i, int j, int itype, int jtype, + double rsq, + double /*factor_coul*/, double /*factor_lj*/, + double &fforce) +{ + int k,kk,ktype,knum,m; + tagint tagprev; + double eng,eng_lj,force_kernel,force_angle; + double rsq1,rsq2,r1,r2,c,s,ac,r2inv,r10inv,factor_hb; + double switch1,switch2; + double delr1[3],delr2[3]; + tagint *klist; + + double **x = atom->x; + int *type = atom->type; + double *special_lj = force->special_lj; + + eng = 0.0; + fforce = 0; + + // sanity check + + if (!donor[itype]) return 0.0; + if (!acceptor[jtype]) return 0.0; + + int molecular = atom->molecular; + if (molecular == Atom::MOLECULAR) { + klist = atom->special[i]; + knum = atom->nspecial[i][0]; + } else { + if (atom->molindex[i] < 0) return 0.0; + int imol = atom->molindex[i]; + int iatom = atom->molatom[i]; + Molecule **onemols = atom->avec->onemols; + klist = onemols[imol]->special[iatom]; + knum = onemols[imol]->nspecial[iatom][0]; + tagprev = atom->tag[i] - iatom - 1; + } + + factor_hb = special_lj[sbmask(j)]; + + for (kk = 0; kk < knum; kk++) { + if (molecular == Atom::MOLECULAR) k = atom->map(klist[kk]); + else k = atom->map(klist[kk]+tagprev); + + if (k < 0) continue; + ktype = type[k]; + m = type2param[itype][jtype][ktype]; + if (m < 0) continue; + const Param &pm = params[m]; + + delr1[0] = x[i][0] - x[k][0]; + delr1[1] = x[i][1] - x[k][1]; + delr1[2] = x[i][2] - x[k][2]; + domain->minimum_image(delr1); + rsq1 = delr1[0]*delr1[0] + delr1[1]*delr1[1] + delr1[2]*delr1[2]; + r1 = sqrt(rsq1); + + delr2[0] = x[j][0] - x[k][0]; + delr2[1] = x[j][1] - x[k][1]; + delr2[2] = x[j][2] - x[k][2]; + domain->minimum_image(delr2); + rsq2 = delr2[0]*delr2[0] + delr2[1]*delr2[1] + delr2[2]*delr2[2]; + r2 = sqrt(rsq2); + + // angle (cos and sin) + + c = delr1[0]*delr2[0] + delr1[1]*delr2[1] + delr1[2]*delr2[2]; + c /= r1*r2; + if (c > 1.0) c = 1.0; + if (c < -1.0) c = -1.0; + ac = acos(c); + + if (ac < pm.cut_angle || ac > (2.0*MY_PI - pm.cut_angle)) return 0.0; + s = sqrt(1.0 - c*c); + if (s < SMALL) s = SMALL; + + // LJ-specific kernel + + r2inv = 1.0/rsq; + r10inv = r2inv*r2inv*r2inv*r2inv*r2inv; + force_kernel = r10inv*(pm.lj1*r2inv - pm.lj2)*r2inv * powint(c,pm.ap); + force_angle = pm.ap * r10inv*(pm.lj3*r2inv - pm.lj4) * + powint(c,pm.ap-1)*s; + + // only lj part for now + + eng_lj = r10inv*(pm.lj3*r2inv - pm.lj4); + if (rsq > pm.cut_innersq) { + switch1 = (pm.cut_outersq-rsq) * (pm.cut_outersq-rsq) * + (pm.cut_outersq + 2.0*rsq - 3.0*pm.cut_innersq) / pm.denom_vdw; + switch2 = 12.0*rsq * (pm.cut_outersq-rsq) * + (rsq-pm.cut_innersq) / pm.denom_vdw; + force_kernel = force_kernel*switch1 + eng_lj*switch2; + eng_lj *= switch1; + } + + fforce += force_kernel*powint(c,pm.ap) + eng_lj*force_angle; + eng += eng_lj * powint(c,pm.ap) * factor_hb; + } + + return eng; +} diff --git a/src/EXTRA-MOLECULE/pair_hbond_dreiding_lj_angleoffset.h b/src/EXTRA-MOLECULE/pair_hbond_dreiding_lj_angleoffset.h new file mode 100644 index 0000000000..91a906c5cb --- /dev/null +++ b/src/EXTRA-MOLECULE/pair_hbond_dreiding_lj_angleoffset.h @@ -0,0 +1,66 @@ +/* -*- c++ -*- ---------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + https://www.lammps.org/, Sandia National Laboratories + LAMMPS development team: developers@lammps.org + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#ifdef PAIR_CLASS +// clang-format off +PairStyle(hbond/dreiding/lj,PairHbondDreidingLJ); +// clang-format on +#else + +#ifndef LMP_PAIR_HBOND_DREIDING_LJ_H +#define LMP_PAIR_HBOND_DREIDING_LJ_H + +#include "pair.h" + +namespace LAMMPS_NS { + +class PairHbondDreidingLJ : public Pair { + public: + PairHbondDreidingLJ(class LAMMPS *); + ~PairHbondDreidingLJ() override; + void compute(int, int) override; + void settings(int, char **) override; + void coeff(int, char **) override; + void init_style() override; + double init_one(int, int) override; + double single(int, int, int, int, double, double, double, double &) override; + + protected: + double cut_inner_global, cut_outer_global, cut_angle_global; + int ap_global; + + struct Param { + double epsilon, sigma; + double lj1, lj2, lj3, lj4; + double d0, alpha, r0; + double morse1; + double denom_vdw; + double cut_inner, cut_outer, cut_innersq, cut_outersq, cut_angle, offset; + int ap; + }; + + Param *params; // parameter set for an I-J-K interaction + int nparams; // number of parameters read + int maxparam; + + int *donor; // 1 if this type is ever a donor, else 0 + int *acceptor; // 1 if this type is ever an acceptor, else 0 + int ***type2param; // mapping from D,A,H to params, -1 if no map + + void allocate(); +}; + +} // namespace LAMMPS_NS + +#endif +#endif diff --git a/src/EXTRA-MOLECULE/pair_hbond_dreiding_morse_angleoffset.cpp b/src/EXTRA-MOLECULE/pair_hbond_dreiding_morse_angleoffset.cpp new file mode 100644 index 0000000000..c8bc0a627d --- /dev/null +++ b/src/EXTRA-MOLECULE/pair_hbond_dreiding_morse_angleoffset.cpp @@ -0,0 +1,474 @@ +// clang-format off +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + https://www.lammps.org/, Sandia National Laboratories + LAMMPS development team: developers@lammps.org + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. + ------------------------------------------------------------------------- */ + +/* ---------------------------------------------------------------------- + Contributing author: Tod A Pascal (Caltech) +------------------------------------------------------------------------- */ + +#include "pair_hbond_dreiding_morse.h" + +#include "atom.h" +#include "atom_vec.h" +#include "domain.h" +#include "error.h" +#include "force.h" +#include "math_const.h" +#include "math_special.h" +#include "memory.h" +#include "molecule.h" +#include "neigh_list.h" +#include "neighbor.h" + +#include +#include + +using namespace LAMMPS_NS; +using namespace MathConst; +using namespace MathSpecial; + +static constexpr double SMALL = 0.001; +static constexpr int CHUNK = 8; + +/* ---------------------------------------------------------------------- */ + +PairHbondDreidingMorse::PairHbondDreidingMorse(LAMMPS *lmp) : + PairHbondDreidingLJ(lmp) {} + +/* ---------------------------------------------------------------------- */ + +void PairHbondDreidingMorse::compute(int eflag, int vflag) +{ + int i,j,k,m,ii,jj,kk,inum,jnum,knum,itype,jtype,ktype,imol,iatom; + tagint tagprev; + double delx,dely,delz,rsq,rsq1,rsq2,r1,r2; + double factor_hb,force_angle,force_kernel,force_switch,evdwl,ehbond; + double c,s,a,b,d,ac,a11,a12,a22,vx1,vx2,vy1,vy2,vz1,vz2; + double fi[3],fj[3],delr1[3],delr2[3]; + double r,dr,dexp,eng_morse,switch1,switch2; + int *ilist,*jlist,*numneigh,**firstneigh; + tagint *klist; + + evdwl = ehbond = 0.0; + ev_init(eflag,vflag); + + double **x = atom->x; + double **f = atom->f; + tagint *tag = atom->tag; + int *molindex = atom->molindex; + int *molatom = atom->molatom; + tagint **special = atom->special; + int **nspecial = atom->nspecial; + int *type = atom->type; + double *special_lj = force->special_lj; + int molecular = atom->molecular; + Molecule **onemols = atom->avec->onemols; + + inum = list->inum; + ilist = list->ilist; + numneigh = list->numneigh; + firstneigh = list->firstneigh; + + // ii = loop over donors + // jj = loop over acceptors + // kk = loop over hydrogens bonded to donor + + int hbcount = 0; + + for (ii = 0; ii < inum; ii++) { + i = ilist[ii]; + itype = type[i]; + if (!donor[itype]) continue; + if (molecular == Atom::MOLECULAR) { + klist = special[i]; + knum = nspecial[i][0]; + } else { + if (molindex[i] < 0) continue; + imol = molindex[i]; + iatom = molatom[i]; + klist = onemols[imol]->special[iatom]; + knum = onemols[imol]->nspecial[iatom][0]; + tagprev = tag[i] - iatom - 1; + } + jlist = firstneigh[i]; + jnum = numneigh[i]; + + for (jj = 0; jj < jnum; jj++) { + j = jlist[jj]; + factor_hb = special_lj[sbmask(j)]; + j &= NEIGHMASK; + + jtype = type[j]; + if (!acceptor[jtype]) continue; + + delx = x[i][0] - x[j][0]; + dely = x[i][1] - x[j][1]; + delz = x[i][2] - x[j][2]; + rsq = delx*delx + dely*dely + delz*delz; + + for (kk = 0; kk < knum; kk++) { + if (molecular == Atom::MOLECULAR) k = atom->map(klist[kk]); + else k = atom->map(klist[kk]+tagprev); + if (k < 0) continue; + ktype = type[k]; + m = type2param[itype][jtype][ktype]; + if (m < 0) continue; + const Param &pm = params[m]; + + if (rsq < pm.cut_outersq) { + delr1[0] = x[i][0] - x[k][0]; + delr1[1] = x[i][1] - x[k][1]; + delr1[2] = x[i][2] - x[k][2]; + domain->minimum_image(delr1); + rsq1 = delr1[0]*delr1[0] + delr1[1]*delr1[1] + delr1[2]*delr1[2]; + r1 = sqrt(rsq1); + + delr2[0] = x[j][0] - x[k][0]; + delr2[1] = x[j][1] - x[k][1]; + delr2[2] = x[j][2] - x[k][2]; + domain->minimum_image(delr2); + rsq2 = delr2[0]*delr2[0] + delr2[1]*delr2[1] + delr2[2]*delr2[2]; + r2 = sqrt(rsq2); + + // angle (cos and sin) + + c = delr1[0]*delr2[0] + delr1[1]*delr2[1] + delr1[2]*delr2[2]; + c /= r1*r2; + if (c > 1.0) c = 1.0; + if (c < -1.0) c = -1.0; + ac = acos(c); + + if (ac > pm.cut_angle && ac < (2.0*MY_PI - pm.cut_angle)) { + s = sqrt(1.0 - c*c); + if (s < SMALL) s = SMALL; + + // Morse-specific kernel + + r = sqrt(rsq); + dr = r - pm.r0; + dexp = exp(-pm.alpha * dr); + eng_morse = pm.d0 * (dexp*dexp - 2.0*dexp); + force_kernel = pm.morse1*(dexp*dexp - dexp)/r * powint(c,pm.ap); + force_angle = pm.ap * eng_morse * powint(c,pm.ap-1)*s; + force_switch = 0.0; + + if (rsq > pm.cut_innersq) { + switch1 = (pm.cut_outersq-rsq) * (pm.cut_outersq-rsq) * + (pm.cut_outersq + 2.0*rsq - 3.0*pm.cut_innersq) / + pm.denom_vdw; + switch2 = 12.0*rsq * (pm.cut_outersq-rsq) * + (rsq-pm.cut_innersq) / pm.denom_vdw; + + force_kernel *= switch1; + force_angle *= switch1; + force_switch = eng_morse*switch2/rsq; + eng_morse *= switch1; + } + + if (eflag) { + evdwl = eng_morse * powint(c,pm.ap); + evdwl *= factor_hb; + ehbond += evdwl; + } + + a = factor_hb*force_angle/s; + b = factor_hb*force_kernel; + d = factor_hb*force_switch; + + a11 = a*c / rsq1; + a12 = -a / (r1*r2); + a22 = a*c / rsq2; + + vx1 = a11*delr1[0] + a12*delr2[0]; + vx2 = a22*delr2[0] + a12*delr1[0]; + vy1 = a11*delr1[1] + a12*delr2[1]; + vy2 = a22*delr2[1] + a12*delr1[1]; + vz1 = a11*delr1[2] + a12*delr2[2]; + vz2 = a22*delr2[2] + a12*delr1[2]; + + fi[0] = vx1 + (b+d)*delx; + fi[1] = vy1 + (b+d)*dely; + fi[2] = vz1 + (b+d)*delz; + fj[0] = vx2 - (b+d)*delx; + fj[1] = vy2 - (b+d)*dely; + fj[2] = vz2 - (b+d)*delz; + + f[i][0] += fi[0]; + f[i][1] += fi[1]; + f[i][2] += fi[2]; + + f[j][0] += fj[0]; + f[j][1] += fj[1]; + f[j][2] += fj[2]; + + f[k][0] -= vx1 + vx2; + f[k][1] -= vy1 + vy2; + f[k][2] -= vz1 + vz2; + + // KIJ instead of IJK b/c delr1/delr2 are both with respect to k + + if (evflag) ev_tally3(k,i,j,evdwl,0.0,fi,fj,delr1,delr2); + + hbcount++; + } + } + } + } + } + + if (eflag_global) { + pvector[0] = hbcount; + pvector[1] = ehbond; + } +} + +/* ---------------------------------------------------------------------- + set coeffs for one or more type pairs +------------------------------------------------------------------------- */ + +void PairHbondDreidingMorse::coeff(int narg, char **arg) +{ + if (narg < 7 || narg > 11) + error->all(FLERR,"Incorrect args for pair coefficients"); + if (!allocated) allocate(); + + int ilo,ihi,jlo,jhi,klo,khi; + utils::bounds(FLERR, arg[0], 1, atom->ntypes, ilo, ihi, error); + utils::bounds(FLERR, arg[1], 1, atom->ntypes, jlo, jhi, error); + utils::bounds_typelabel(FLERR, arg[2], 1, atom->ntypes, klo, khi, lmp, Atom::ATOM); + + int donor_flag; + if (strcmp(arg[3],"i") == 0) donor_flag = 0; + else if (strcmp(arg[3],"j") == 0) donor_flag = 1; + else error->all(FLERR,"Incorrect args for pair coefficients"); + + double d0_one = utils::numeric(FLERR, arg[4], false, lmp); + double alpha_one = utils::numeric(FLERR, arg[5], false, lmp); + double r0_one = utils::numeric(FLERR, arg[6], false, lmp); + + int ap_one = ap_global; + if (narg > 7) ap_one = utils::inumeric(FLERR, arg[7], false, lmp); + double cut_inner_one = cut_inner_global; + double cut_outer_one = cut_outer_global; + if (narg > 9) { + cut_inner_one = utils::numeric(FLERR, arg[8], false, lmp); + cut_outer_one = utils::numeric(FLERR, arg[9], false, lmp); + } + if (cut_inner_one>cut_outer_one) + error->all(FLERR,"Pair inner cutoff >= Pair outer cutoff"); + double cut_angle_one = cut_angle_global; + if (narg > 10) cut_angle_one = utils::numeric(FLERR, arg[10], false, lmp) * MY_PI/180.0; + + // grow params array if necessary + + if (nparams == maxparam) { + maxparam += CHUNK; + params = (Param *) memory->srealloc(params, maxparam*sizeof(Param), + "pair:params"); + + // make certain all addional allocated storage is initialized + // to avoid false positives when checking with valgrind + + memset(params + nparams, 0, CHUNK*sizeof(Param)); + } + + params[nparams].d0 = d0_one; + params[nparams].alpha = alpha_one; + params[nparams].r0 = r0_one; + params[nparams].ap = ap_one; + params[nparams].cut_inner = cut_inner_one; + params[nparams].cut_outer = cut_outer_one; + params[nparams].cut_innersq = cut_inner_one*cut_inner_one; + params[nparams].cut_outersq = cut_outer_one*cut_outer_one; + params[nparams].cut_angle = cut_angle_one; + params[nparams].denom_vdw = + (params[nparams].cut_outersq-params[nparams].cut_innersq) * + (params[nparams].cut_outersq-params[nparams].cut_innersq) * + (params[nparams].cut_outersq-params[nparams].cut_innersq); + + // flag type2param with either i,j = D,A or j,i = D,A + + int count = 0; + for (int i = ilo; i <= ihi; i++) + for (int j = MAX(jlo,i); j <= jhi; j++) + for (int k = klo; k <= khi; k++) { + if (donor_flag == 0) type2param[i][j][k] = nparams; + else type2param[j][i][k] = nparams; + count++; + } + nparams++; + + if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients"); +} + +/* ---------------------------------------------------------------------- + init specific to this pair style +------------------------------------------------------------------------- */ + +void PairHbondDreidingMorse::init_style() +{ + // molecular system required to use special list to find H atoms + // tags required to use special list + // pair newton on required since are looping over D atoms + // and computing forces on A,H which may be on different procs + + if (atom->molecular == Atom::ATOMIC) + error->all(FLERR,"Pair style hbond/dreiding requires molecular system"); + if (atom->tag_enable == 0) + error->all(FLERR,"Pair style hbond/dreiding requires atom IDs"); + if (atom->map_style == Atom::MAP_NONE) + error->all(FLERR,"Pair style hbond/dreiding requires an atom map, " + "see atom_modify"); + if (force->newton_pair == 0) + error->all(FLERR,"Pair style hbond/dreiding requires newton pair on"); + + // set donor[M]/acceptor[M] if any atom of type M is a donor/acceptor + + int anyflag = 0; + int n = atom->ntypes; + for (int m = 1; m <= n; m++) donor[m] = acceptor[m] = 0; + for (int i = 1; i <= n; i++) + for (int j = 1; j <= n; j++) + for (int k = 1; k <= n; k++) + if (type2param[i][j][k] >= 0) { + anyflag = 1; + donor[i] = 1; + acceptor[j] = 1; + } + + if (!anyflag) error->all(FLERR,"No pair hbond/dreiding coefficients set"); + + // set additional param values + // offset is for Morse only, angle term is not included + + for (int m = 0; m < nparams; m++) { + params[m].morse1 = 2.0*params[m].d0*params[m].alpha; + + /* + if (offset_flag) { + double alpha_dr = -params[m].alpha * (params[m].cut - params[m].r0); + params[m].offset = params[m].d0 * + ((exp(2.0*alpha_dr)) - (2.0*exp(alpha_dr))); + } else params[m].offset = 0.0; + */ + } + + // full neighbor list request + + neighbor->add_request(this, NeighConst::REQ_FULL); +} + +/* ---------------------------------------------------------------------- */ + +double PairHbondDreidingMorse::single(int i, int j, int itype, int jtype, + double rsq, + double /*factor_coul*/, double /*factor_lj*/, + double &fforce) +{ + int k,kk,ktype,knum,m; + tagint tagprev; + double eng,eng_morse,force_kernel,force_angle; + double rsq1,rsq2,r1,r2,c,s,ac,r,dr,dexp,factor_hb; + double switch1,switch2; + double delr1[3],delr2[3]; + tagint *klist; + + double **x = atom->x; + int *type = atom->type; + double *special_lj = force->special_lj; + + eng = 0.0; + fforce = 0; + + //sanity check + + if (!donor[itype]) return 0.0; + if (!acceptor[jtype]) return 0.0; + + int molecular = atom->molecular; + if (molecular == Atom::MOLECULAR) { + klist = atom->special[i]; + knum = atom->nspecial[i][0]; + } else { + if (atom->molindex[i] < 0) return 0.0; + int imol = atom->molindex[i]; + int iatom = atom->molatom[i]; + Molecule **onemols = atom->avec->onemols; + klist = onemols[imol]->special[iatom]; + knum = onemols[imol]->nspecial[iatom][0]; + tagprev = atom->tag[i] - iatom - 1; + } + + factor_hb = special_lj[sbmask(j)]; + + for (kk = 0; kk < knum; kk++) { + if (molecular == Atom::MOLECULAR) k = atom->map(klist[kk]); + else k = atom->map(klist[kk]+tagprev); + + if (k < 0) continue; + ktype = type[k]; + m = type2param[itype][jtype][ktype]; + if (m < 0) continue; + const Param &pm = params[m]; + + delr1[0] = x[i][0] - x[k][0]; + delr1[1] = x[i][1] - x[k][1]; + delr1[2] = x[i][2] - x[k][2]; + domain->minimum_image(delr1); + rsq1 = delr1[0]*delr1[0] + delr1[1]*delr1[1] + delr1[2]*delr1[2]; + r1 = sqrt(rsq1); + + delr2[0] = x[j][0] - x[k][0]; + delr2[1] = x[j][1] - x[k][1]; + delr2[2] = x[j][2] - x[k][2]; + domain->minimum_image(delr2); + rsq2 = delr2[0]*delr2[0] + delr2[1]*delr2[1] + delr2[2]*delr2[2]; + r2 = sqrt(rsq2); + + // angle (cos and sin) + + c = delr1[0]*delr2[0] + delr1[1]*delr2[1] + delr1[2]*delr2[2]; + c /= r1*r2; + if (c > 1.0) c = 1.0; + if (c < -1.0) c = -1.0; + ac = acos(c); + + if (ac < pm.cut_angle || ac > (2.0*MY_PI - pm.cut_angle)) return 0.0; + s = sqrt(1.0 - c*c); + if (s < SMALL) s = SMALL; + + // Morse-specific kernel + + r = sqrt(rsq); + dr = r - pm.r0; + dexp = exp(-pm.alpha * dr); + eng_morse = pm.d0 * (dexp*dexp - 2.0*dexp); //<-- BUGFIX 2012-11-14 + force_kernel = pm.morse1*(dexp*dexp - dexp)/r * powint(c,pm.ap); + force_angle = pm.ap * eng_morse * powint(c,pm.ap-1)*s; + + if (rsq > pm.cut_innersq) { + switch1 = (pm.cut_outersq-rsq) * (pm.cut_outersq-rsq) * + (pm.cut_outersq + 2.0*rsq - 3.0*pm.cut_innersq) / + pm.denom_vdw; + switch2 = 12.0*rsq * (pm.cut_outersq-rsq) * + (rsq-pm.cut_innersq) / pm.denom_vdw; + force_kernel = force_kernel*switch1 + eng_morse*switch2; + eng_morse *= switch1; + } + + eng += eng_morse * powint(c,pm.ap)* factor_hb; + fforce += force_kernel*powint(c,pm.ap) + eng_morse*force_angle; + } + + return eng; +} diff --git a/src/EXTRA-MOLECULE/pair_hbond_dreiding_morse_angleoffset.h b/src/EXTRA-MOLECULE/pair_hbond_dreiding_morse_angleoffset.h new file mode 100644 index 0000000000..8d4646451e --- /dev/null +++ b/src/EXTRA-MOLECULE/pair_hbond_dreiding_morse_angleoffset.h @@ -0,0 +1,40 @@ +/* -*- c++ -*- ---------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + https://www.lammps.org/, Sandia National Laboratories + LAMMPS development team: developers@lammps.org + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#ifdef PAIR_CLASS +// clang-format off +PairStyle(hbond/dreiding/morse,PairHbondDreidingMorse); +// clang-format on +#else + +#ifndef LMP_PAIR_HBOND_DREIDING_MORSE_H +#define LMP_PAIR_HBOND_DREIDING_MORSE_H + +#include "pair_hbond_dreiding_lj.h" + +namespace LAMMPS_NS { + +class PairHbondDreidingMorse : public PairHbondDreidingLJ { + public: + PairHbondDreidingMorse(class LAMMPS *); + + void compute(int, int) override; + void coeff(int, char **) override; + void init_style() override; + double single(int, int, int, int, double, double, double, double &) override; +}; + +} // namespace LAMMPS_NS + +#endif +#endif diff --git a/src/OPENMP/pair_hbond_dreiding_lj_angleoffset_omp.cpp b/src/OPENMP/pair_hbond_dreiding_lj_angleoffset_omp.cpp new file mode 100644 index 0000000000..b0f6dcfb5b --- /dev/null +++ b/src/OPENMP/pair_hbond_dreiding_lj_angleoffset_omp.cpp @@ -0,0 +1,317 @@ +// clang-format off +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + https://www.lammps.org/, Sandia National Laboratories + LAMMPS development team: developers@lammps.org + + This software is distributed under the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +/* ---------------------------------------------------------------------- + Contributing author: Axel Kohlmeyer (Temple U) +------------------------------------------------------------------------- */ + +#include "pair_hbond_dreiding_lj_omp.h" + +#include "atom.h" +#include "atom_vec.h" +#include "comm.h" +#include "domain.h" +#include "force.h" +#include "math_const.h" +#include "math_special.h" +#include "molecule.h" +#include "neigh_list.h" +#include "suffix.h" + +#include + +#include "omp_compat.h" +using namespace LAMMPS_NS; +using namespace MathConst; +using namespace MathSpecial; + +static constexpr double SMALL = 0.001; + +/* ---------------------------------------------------------------------- */ + +PairHbondDreidingLJOMP::PairHbondDreidingLJOMP(LAMMPS *lmp) : + PairHbondDreidingLJ(lmp), ThrOMP(lmp, THR_PAIR) +{ + suffix_flag |= Suffix::OMP; + respa_enable = 0; + hbcount_thr = hbeng_thr = nullptr; +} + +/* ---------------------------------------------------------------------- */ + +PairHbondDreidingLJOMP::~PairHbondDreidingLJOMP() +{ + if (hbcount_thr) { + delete[] hbcount_thr; + delete[] hbeng_thr; + } +} + +/* ---------------------------------------------------------------------- */ + +void PairHbondDreidingLJOMP::compute(int eflag, int vflag) +{ + ev_init(eflag,vflag); + + const int nall = atom->nlocal + atom->nghost; + const int nthreads = comm->nthreads; + const int inum = list->inum; + + if (!hbcount_thr) { + hbcount_thr = new double[nthreads]; + hbeng_thr = new double[nthreads]; + } + + for (int i=0; i < nthreads; ++i) { + hbcount_thr[i] = 0.0; + hbeng_thr[i] = 0.0; + } + +#if defined(_OPENMP) +#pragma omp parallel LMP_DEFAULT_NONE LMP_SHARED(eflag,vflag) +#endif + { + int ifrom, ito, tid; + + loop_setup_thr(ifrom, ito, tid, inum, nthreads); + ThrData *thr = fix->get_thr(tid); + thr->timer(Timer::START); + ev_setup_thr(eflag, vflag, nall, eatom, vatom, nullptr, thr); + + if (evflag) { + if (eflag) { + if (force->newton_pair) eval<1,1,1>(ifrom, ito, thr); + else eval<1,1,0>(ifrom, ito, thr); + } else { + if (force->newton_pair) eval<1,0,1>(ifrom, ito, thr); + else eval<1,0,0>(ifrom, ito, thr); + } + } else { + if (force->newton_pair) eval<0,0,1>(ifrom, ito, thr); + else eval<0,0,0>(ifrom, ito, thr); + } + + thr->timer(Timer::PAIR); + reduce_thr(this, eflag, vflag, thr); + } // end of omp parallel region + + // reduce per thread hbond data + if (eflag_global) { + pvector[0] = 0.0; + pvector[1] = 0.0; + for (int i=0; i < nthreads; ++i) { + pvector[0] += hbcount_thr[i]; + pvector[1] += hbeng_thr[i]; + } + } +} + +template +void PairHbondDreidingLJOMP::eval(int iifrom, int iito, ThrData * const thr) +{ + int i,j,k,m,ii,jj,kk,jnum,knum,itype,jtype,ktype,iatom,imol; + tagint tagprev; + double xtmp,ytmp,ztmp,delx,dely,delz,rsq,rsq1,rsq2,r1,r2; + double factor_hb,force_angle,force_kernel,evdwl,eng_lj; + double c,s,a,b,ac,a11,a12,a22,vx1,vx2,vy1,vy2,vz1,vz2; + double fi[3],fj[3],delr1[3],delr2[3]; + double r2inv,r10inv; + double switch1,switch2; + int *ilist,*jlist,*numneigh,**firstneigh; + const tagint *klist; + + evdwl = 0.0; + + const auto * _noalias const x = (dbl3_t *) atom->x[0]; + auto * _noalias const f = (dbl3_t *) thr->get_f()[0]; + const tagint * _noalias const tag = atom->tag; + const int * _noalias const molindex = atom->molindex; + const int * _noalias const molatom = atom->molatom; + const int * _noalias const type = atom->type; + const double * _noalias const special_lj = force->special_lj; + const int * const * const nspecial = atom->nspecial; + const tagint * const * const special = atom->special; + const int molecular = atom->molecular; + Molecule * const * const onemols = atom->avec->onemols; + double fxtmp,fytmp,fztmp; + + ilist = list->ilist; + numneigh = list->numneigh; + firstneigh = list->firstneigh; + + // ii = loop over donors + // jj = loop over acceptors + // kk = loop over hydrogens bonded to donor + + int hbcount = 0; + double hbeng = 0.0; + + for (ii = iifrom; ii < iito; ++ii) { + i = ilist[ii]; + itype = type[i]; + if (!donor[itype]) continue; + if (molecular == Atom::MOLECULAR) { + klist = special[i]; + knum = nspecial[i][0]; + } else { + if (molindex[i] < 0) continue; + imol = molindex[i]; + iatom = molatom[i]; + klist = onemols[imol]->special[iatom]; + knum = onemols[imol]->nspecial[iatom][0]; + tagprev = tag[i] - iatom - 1; + } + jlist = firstneigh[i]; + jnum = numneigh[i]; + fxtmp=fytmp=fztmp=0.0; + + xtmp = x[i].x; + ytmp = x[i].y; + ztmp = x[i].z; + + for (jj = 0; jj < jnum; jj++) { + j = jlist[jj]; + factor_hb = special_lj[sbmask(j)]; + j &= NEIGHMASK; + + jtype = type[j]; + if (!acceptor[jtype]) continue; + + delx = xtmp - x[j].x; + dely = ytmp - x[j].y; + delz = ztmp - x[j].z; + rsq = delx*delx + dely*dely + delz*delz; + + for (kk = 0; kk < knum; kk++) { + if (molecular == Atom::MOLECULAR) k = atom->map(klist[kk]); + else k = atom->map(klist[kk]+tagprev); + if (k < 0) continue; + ktype = type[k]; + m = type2param[itype][jtype][ktype]; + if (m < 0) continue; + const Param &pm = params[m]; + + if (rsq < pm.cut_outersq) { + delr1[0] = xtmp - x[k].x; + delr1[1] = ytmp - x[k].y; + delr1[2] = ztmp - x[k].z; + domain->minimum_image(delr1); + rsq1 = delr1[0]*delr1[0] + delr1[1]*delr1[1] + delr1[2]*delr1[2]; + r1 = sqrt(rsq1); + + delr2[0] = x[j].x - x[k].x; + delr2[1] = x[j].y - x[k].y; + delr2[2] = x[j].z - x[k].z; + domain->minimum_image(delr2); + rsq2 = delr2[0]*delr2[0] + delr2[1]*delr2[1] + delr2[2]*delr2[2]; + r2 = sqrt(rsq2); + + // angle (cos and sin) + + c = delr1[0]*delr2[0] + delr1[1]*delr2[1] + delr1[2]*delr2[2]; + c /= r1*r2; + if (c > 1.0) c = 1.0; + if (c < -1.0) c = -1.0; + ac = acos(c); + + if (ac > pm.cut_angle && ac < (2.0*MY_PI - pm.cut_angle)) { + s = sqrt(1.0 - c*c); + if (s < SMALL) s = SMALL; + + // LJ-specific kernel + + r2inv = 1.0/rsq; + r10inv = r2inv*r2inv*r2inv*r2inv*r2inv; + force_kernel = r10inv*(pm.lj1*r2inv - pm.lj2)*r2inv * + powint(c,pm.ap); + force_angle = pm.ap * r10inv*(pm.lj3*r2inv - pm.lj4) * + powint(c,pm.ap-1)*s; + + eng_lj = r10inv*(pm.lj3*r2inv - pm.lj4); + if (rsq > pm.cut_innersq) { + switch1 = (pm.cut_outersq-rsq) * (pm.cut_outersq-rsq) * + (pm.cut_outersq + 2.0*rsq - 3.0*pm.cut_innersq) / + pm.denom_vdw; + switch2 = 12.0*rsq * (pm.cut_outersq-rsq) * + (rsq-pm.cut_innersq) / pm.denom_vdw; + force_kernel = force_kernel*switch1 + eng_lj*switch2/rsq; + force_angle *= switch1; + eng_lj *= switch1; + } + + if (EFLAG) { + evdwl = eng_lj * powint(c,pm.ap); + evdwl *= factor_hb; + } + + a = factor_hb*force_angle/s; + b = factor_hb*force_kernel; + + a11 = a*c / rsq1; + a12 = -a / (r1*r2); + a22 = a*c / rsq2; + + vx1 = a11*delr1[0] + a12*delr2[0]; + vx2 = a22*delr2[0] + a12*delr1[0]; + vy1 = a11*delr1[1] + a12*delr2[1]; + vy2 = a22*delr2[1] + a12*delr1[1]; + vz1 = a11*delr1[2] + a12*delr2[2]; + vz2 = a22*delr2[2] + a12*delr1[2]; + + fi[0] = vx1 + b*delx; + fi[1] = vy1 + b*dely; + fi[2] = vz1 + b*delz; + fj[0] = vx2 - b*delx; + fj[1] = vy2 - b*dely; + fj[2] = vz2 - b*delz; + + fxtmp += fi[0]; + fytmp += fi[1]; + fztmp += fi[2]; + + f[j].x += fj[0]; + f[j].y += fj[1]; + f[j].z += fj[2]; + + f[k].x -= vx1 + vx2; + f[k].y -= vy1 + vy2; + f[k].z -= vz1 + vz2; + + // KIJ instead of IJK b/c delr1/delr2 are both with respect to k + + if (EVFLAG) ev_tally3_thr(this,k,i,j,evdwl,0.0,fi,fj,delr1,delr2,thr); + if (EFLAG) { + hbcount++; + hbeng += evdwl; + } + } + } + } + } + f[i].x += fxtmp; + f[i].y += fytmp; + f[i].z += fztmp; + } + const int tid = thr->get_tid(); + hbcount_thr[tid] = static_cast(hbcount); + hbeng_thr[tid] = hbeng; +} + +/* ---------------------------------------------------------------------- */ + +double PairHbondDreidingLJOMP::memory_usage() +{ + double bytes = memory_usage_thr(); + bytes += (double)comm->nthreads * 2 * sizeof(double); + bytes += PairHbondDreidingLJ::memory_usage(); + + return bytes; +} diff --git a/src/OPENMP/pair_hbond_dreiding_lj_angleoffset_omp.h b/src/OPENMP/pair_hbond_dreiding_lj_angleoffset_omp.h new file mode 100644 index 0000000000..792a9a6f8c --- /dev/null +++ b/src/OPENMP/pair_hbond_dreiding_lj_angleoffset_omp.h @@ -0,0 +1,52 @@ +/* -*- c++ -*- ---------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + https://www.lammps.org/, Sandia National Laboratories + LAMMPS development team: developers@lammps.org + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +/* ---------------------------------------------------------------------- + Contributing author: Axel Kohlmeyer (Temple U) +------------------------------------------------------------------------- */ + +#ifdef PAIR_CLASS +// clang-format off +PairStyle(hbond/dreiding/lj/omp,PairHbondDreidingLJOMP); +// clang-format on +#else + +#ifndef LMP_PAIR_HBOND_DREIDING_LJ_OMP_H +#define LMP_PAIR_HBOND_DREIDING_LJ_OMP_H + +#include "pair_hbond_dreiding_lj.h" +#include "thr_omp.h" + +namespace LAMMPS_NS { + +class PairHbondDreidingLJOMP : public PairHbondDreidingLJ, public ThrOMP { + + public: + PairHbondDreidingLJOMP(class LAMMPS *); + ~PairHbondDreidingLJOMP() override; + + void compute(int, int) override; + double memory_usage() override; + + protected: + double *hbcount_thr, *hbeng_thr; + + private: + template + void eval(int ifrom, int ito, ThrData *const thr); +}; + +} // namespace LAMMPS_NS + +#endif +#endif diff --git a/src/OPENMP/pair_hbond_dreiding_morse_angleoffset_omp.cpp b/src/OPENMP/pair_hbond_dreiding_morse_angleoffset_omp.cpp new file mode 100644 index 0000000000..0e43e2a037 --- /dev/null +++ b/src/OPENMP/pair_hbond_dreiding_morse_angleoffset_omp.cpp @@ -0,0 +1,316 @@ +// clang-format off +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + https://www.lammps.org/, Sandia National Laboratories + LAMMPS development team: developers@lammps.org + + This software is distributed under the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +/* ---------------------------------------------------------------------- + Contributing author: Axel Kohlmeyer (Temple U) +------------------------------------------------------------------------- */ + +#include "pair_hbond_dreiding_morse_omp.h" + +#include "atom.h" +#include "atom_vec.h" +#include "comm.h" +#include "domain.h" +#include "force.h" +#include "math_const.h" +#include "math_special.h" +#include "molecule.h" +#include "neigh_list.h" +#include "suffix.h" + +#include + +#include "omp_compat.h" +using namespace LAMMPS_NS; +using namespace MathConst; +using namespace MathSpecial; + +static constexpr double SMALL = 0.001; + +/* ---------------------------------------------------------------------- */ + +PairHbondDreidingMorseOMP::PairHbondDreidingMorseOMP(LAMMPS *lmp) : + PairHbondDreidingMorse(lmp), ThrOMP(lmp, THR_PAIR) +{ + suffix_flag |= Suffix::OMP; + respa_enable = 0; + hbcount_thr = hbeng_thr = nullptr; +} + +/* ---------------------------------------------------------------------- */ + +PairHbondDreidingMorseOMP::~PairHbondDreidingMorseOMP() +{ + if (hbcount_thr) { + delete[] hbcount_thr; + delete[] hbeng_thr; + } +} + +/* ---------------------------------------------------------------------- */ + +void PairHbondDreidingMorseOMP::compute(int eflag, int vflag) +{ + ev_init(eflag,vflag); + + const int nall = atom->nlocal + atom->nghost; + const int nthreads = comm->nthreads; + const int inum = list->inum; + + if (!hbcount_thr) { + hbcount_thr = new double[nthreads]; + hbeng_thr = new double[nthreads]; + } + + for (int i=0; i < nthreads; ++i) { + hbcount_thr[i] = 0.0; + hbeng_thr[i] = 0.0; + } + +#if defined(_OPENMP) +#pragma omp parallel LMP_DEFAULT_NONE LMP_SHARED(eflag,vflag) +#endif + { + int ifrom, ito, tid; + + loop_setup_thr(ifrom, ito, tid, inum, nthreads); + ThrData *thr = fix->get_thr(tid); + thr->timer(Timer::START); + ev_setup_thr(eflag, vflag, nall, eatom, vatom, nullptr, thr); + + if (evflag) { + if (eflag) { + if (force->newton_pair) eval<1,1,1>(ifrom, ito, thr); + else eval<1,1,0>(ifrom, ito, thr); + } else { + if (force->newton_pair) eval<1,0,1>(ifrom, ito, thr); + else eval<1,0,0>(ifrom, ito, thr); + } + } else { + if (force->newton_pair) eval<0,0,1>(ifrom, ito, thr); + else eval<0,0,0>(ifrom, ito, thr); + } + + thr->timer(Timer::PAIR); + reduce_thr(this, eflag, vflag, thr); + } // end of omp parallel region + + // reduce per thread hbond data + if (eflag_global) { + pvector[0] = 0.0; + pvector[1] = 0.0; + for (int i=0; i < nthreads; ++i) { + pvector[0] += hbcount_thr[i]; + pvector[1] += hbeng_thr[i]; + } + } +} + +template +void PairHbondDreidingMorseOMP::eval(int iifrom, int iito, ThrData * const thr) +{ + int i,j,k,m,ii,jj,kk,jnum,knum,itype,jtype,ktype,imol,iatom; + tagint tagprev; + double xtmp,ytmp,ztmp,delx,dely,delz,rsq,rsq1,rsq2,r1,r2; + double factor_hb,force_angle,force_kernel,evdwl; + double c,s,a,b,ac,a11,a12,a22,vx1,vx2,vy1,vy2,vz1,vz2; + double fi[3],fj[3],delr1[3],delr2[3]; + double r,dr,dexp,eng_morse,switch1,switch2; + int *ilist,*jlist,*numneigh,**firstneigh; + const tagint *klist; + + evdwl = 0.0; + + const auto * _noalias const x = (dbl3_t *) atom->x[0]; + auto * _noalias const f = (dbl3_t *) thr->get_f()[0]; + const tagint * _noalias const tag = atom->tag; + const int * _noalias const type = atom->type; + const int * _noalias const molindex = atom->molindex; + const int * _noalias const molatom = atom->molatom; + const double * _noalias const special_lj = force->special_lj; + const int * const * const nspecial = atom->nspecial; + const tagint * const * const special = atom->special; + const int molecular = atom->molecular; + Molecule * const * const onemols = atom->avec->onemols; + double fxtmp,fytmp,fztmp; + + ilist = list->ilist; + numneigh = list->numneigh; + firstneigh = list->firstneigh; + + // ii = loop over donors + // jj = loop over acceptors + // kk = loop over hydrogens bonded to donor + + int hbcount = 0; + double hbeng = 0.0; + + for (ii = iifrom; ii < iito; ++ii) { + + i = ilist[ii]; + itype = type[i]; + if (!donor[itype]) continue; + if (molecular == Atom::MOLECULAR) { + klist = special[i]; + knum = nspecial[i][0]; + } else { + if (molindex[i] < 0) continue; + imol = molindex[i]; + iatom = molatom[i]; + klist = onemols[imol]->special[iatom]; + knum = onemols[imol]->nspecial[iatom][0]; + tagprev = tag[i] - iatom - 1; + } + jlist = firstneigh[i]; + jnum = numneigh[i]; + fxtmp=fytmp=fztmp=0.0; + + xtmp = x[i].x; + ytmp = x[i].y; + ztmp = x[i].z; + + for (jj = 0; jj < jnum; jj++) { + j = jlist[jj]; + factor_hb = special_lj[sbmask(j)]; + j &= NEIGHMASK; + + jtype = type[j]; + if (!acceptor[jtype]) continue; + + delx = xtmp - x[j].x; + dely = ytmp - x[j].y; + delz = ztmp - x[j].z; + rsq = delx*delx + dely*dely + delz*delz; + + for (kk = 0; kk < knum; kk++) { + if (molecular == Atom::MOLECULAR) k = atom->map(klist[kk]); + else k = atom->map(klist[kk]+tagprev); + if (k < 0) continue; + ktype = type[k]; + m = type2param[itype][jtype][ktype]; + if (m < 0) continue; + const Param &pm = params[m]; + + if (rsq < pm.cut_outersq) { + delr1[0] = xtmp - x[k].x; + delr1[1] = ytmp - x[k].y; + delr1[2] = ztmp - x[k].z; + domain->minimum_image(delr1); + rsq1 = delr1[0]*delr1[0] + delr1[1]*delr1[1] + delr1[2]*delr1[2]; + r1 = sqrt(rsq1); + + delr2[0] = x[j].x - x[k].x; + delr2[1] = x[j].y - x[k].y; + delr2[2] = x[j].z - x[k].z; + domain->minimum_image(delr2); + rsq2 = delr2[0]*delr2[0] + delr2[1]*delr2[1] + delr2[2]*delr2[2]; + r2 = sqrt(rsq2); + + // angle (cos and sin) + + c = delr1[0]*delr2[0] + delr1[1]*delr2[1] + delr1[2]*delr2[2]; + c /= r1*r2; + if (c > 1.0) c = 1.0; + if (c < -1.0) c = -1.0; + ac = acos(c); + + if (ac > pm.cut_angle && ac < (2.0*MY_PI - pm.cut_angle)) { + s = sqrt(1.0 - c*c); + if (s < SMALL) s = SMALL; + + // Morse-specific kernel + + r = sqrt(rsq); + dr = r - pm.r0; + dexp = exp(-pm.alpha * dr); + eng_morse = pm.d0 * (dexp*dexp - 2.0*dexp); + force_kernel = pm.morse1*(dexp*dexp - dexp)/r * powint(c,pm.ap); + force_angle = pm.ap * eng_morse * powint(c,pm.ap-1)*s; + + if (rsq > pm.cut_innersq) { + switch1 = (pm.cut_outersq-rsq) * (pm.cut_outersq-rsq) * + (pm.cut_outersq + 2.0*rsq - 3.0*pm.cut_innersq) / + pm.denom_vdw; + switch2 = 12.0*rsq * (pm.cut_outersq-rsq) * + (rsq-pm.cut_innersq) / pm.denom_vdw; + force_kernel = force_kernel*switch1 + eng_morse*switch2/rsq; + force_angle *= switch1; + eng_morse *= switch1; + } + + if (EFLAG) { + evdwl = eng_morse * powint(c,pm.ap); + evdwl *= factor_hb; + } + + a = factor_hb*force_angle/s; + b = factor_hb*force_kernel; + + a11 = a*c / rsq1; + a12 = -a / (r1*r2); + a22 = a*c / rsq2; + + vx1 = a11*delr1[0] + a12*delr2[0]; + vx2 = a22*delr2[0] + a12*delr1[0]; + vy1 = a11*delr1[1] + a12*delr2[1]; + vy2 = a22*delr2[1] + a12*delr1[1]; + vz1 = a11*delr1[2] + a12*delr2[2]; + vz2 = a22*delr2[2] + a12*delr1[2]; + + fi[0] = vx1 + b*delx; + fi[1] = vy1 + b*dely; + fi[2] = vz1 + b*delz; + fj[0] = vx2 - b*delx; + fj[1] = vy2 - b*dely; + fj[2] = vz2 - b*delz; + + fxtmp += fi[0]; + fytmp += fi[1]; + fztmp += fi[2]; + + f[j].x += fj[0]; + f[j].y += fj[1]; + f[j].z += fj[2]; + + f[k].x -= vx1 + vx2; + f[k].y -= vy1 + vy2; + f[k].z -= vz1 + vz2; + + // KIJ instead of IJK b/c delr1/delr2 are both with respect to k + + if (EVFLAG) ev_tally3_thr(this,k,i,j,evdwl,0.0,fi,fj,delr1,delr2,thr); + if (EFLAG) { + hbcount++; + hbeng += evdwl; + } + } + } + } + } + f[i].x += fxtmp; + f[i].y += fytmp; + f[i].z += fztmp; + } + const int tid = thr->get_tid(); + hbcount_thr[tid] = static_cast(hbcount); + hbeng_thr[tid] = hbeng; +} + +/* ---------------------------------------------------------------------- */ + +double PairHbondDreidingMorseOMP::memory_usage() +{ + double bytes = memory_usage_thr(); + bytes += (double)comm->nthreads * 2 * sizeof(double); + bytes += PairHbondDreidingMorse::memory_usage(); + + return bytes; +} diff --git a/src/OPENMP/pair_hbond_dreiding_morse_angleoffset_omp.h b/src/OPENMP/pair_hbond_dreiding_morse_angleoffset_omp.h new file mode 100644 index 0000000000..016154b762 --- /dev/null +++ b/src/OPENMP/pair_hbond_dreiding_morse_angleoffset_omp.h @@ -0,0 +1,52 @@ +/* -*- c++ -*- ---------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + https://www.lammps.org/, Sandia National Laboratories + LAMMPS development team: developers@lammps.org + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +/* ---------------------------------------------------------------------- + Contributing author: Axel Kohlmeyer (Temple U) +------------------------------------------------------------------------- */ + +#ifdef PAIR_CLASS +// clang-format off +PairStyle(hbond/dreiding/morse/omp,PairHbondDreidingMorseOMP); +// clang-format on +#else + +#ifndef LMP_PAIR_HBOND_DREIDING_MORSE_OMP_H +#define LMP_PAIR_HBOND_DREIDING_MORSE_OMP_H + +#include "pair_hbond_dreiding_morse.h" +#include "thr_omp.h" + +namespace LAMMPS_NS { + +class PairHbondDreidingMorseOMP : public PairHbondDreidingMorse, public ThrOMP { + + public: + PairHbondDreidingMorseOMP(class LAMMPS *); + ~PairHbondDreidingMorseOMP() override; + + void compute(int, int) override; + double memory_usage() override; + + protected: + double *hbcount_thr, *hbeng_thr; + + private: + template + void eval(int ifrom, int ito, ThrData *const thr); +}; + +} // namespace LAMMPS_NS + +#endif +#endif From 898d97e603f7ddce47456e0a8e5f5979c5ba6f27 Mon Sep 17 00:00:00 2001 From: EiPi Fun Date: Sun, 1 Sep 2024 16:13:20 +0800 Subject: [PATCH 02/29] Add angle offset for pair_hbond_dreiding --- .../pair_hbond_dreiding_lj_angleoffset.cpp | 45 +++++++++++++------ .../pair_hbond_dreiding_lj_angleoffset.h | 16 +++---- .../pair_hbond_dreiding_morse_angleoffset.cpp | 33 ++++++++++---- .../pair_hbond_dreiding_morse_angleoffset.h | 12 ++--- ...pair_hbond_dreiding_lj_angleoffset_omp.cpp | 23 ++++++---- .../pair_hbond_dreiding_lj_angleoffset_omp.h | 16 +++---- ...r_hbond_dreiding_morse_angleoffset_omp.cpp | 23 ++++++---- ...air_hbond_dreiding_morse_angleoffset_omp.h | 16 +++---- 8 files changed, 113 insertions(+), 71 deletions(-) diff --git a/src/EXTRA-MOLECULE/pair_hbond_dreiding_lj_angleoffset.cpp b/src/EXTRA-MOLECULE/pair_hbond_dreiding_lj_angleoffset.cpp index 274f8bc2a3..ab15d7c0b5 100644 --- a/src/EXTRA-MOLECULE/pair_hbond_dreiding_lj_angleoffset.cpp +++ b/src/EXTRA-MOLECULE/pair_hbond_dreiding_lj_angleoffset.cpp @@ -13,10 +13,10 @@ ------------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- - Contributing author: Tod A Pascal (Caltech) + Contributing author: Tod A Pascal (Caltech), Don Xu/EiPi Fun ------------------------------------------------------------------------- */ -#include "pair_hbond_dreiding_lj.h" +#include "pair_hbond_dreiding_lj_angleoffset.h" #include "atom.h" #include "atom_vec.h" @@ -42,7 +42,7 @@ static constexpr int CHUNK = 8; /* ---------------------------------------------------------------------- */ -PairHbondDreidingLJ::PairHbondDreidingLJ(LAMMPS *lmp) : Pair(lmp) +PairHbondDreidingLJangleoffset::PairHbondDreidingLJangleoffset(LAMMPS *lmp) : Pair(lmp) { // hbond cannot compute virial as F dot r // due to using map() to find bonded H atoms which are not near donor atom @@ -59,7 +59,7 @@ PairHbondDreidingLJ::PairHbondDreidingLJ(LAMMPS *lmp) : Pair(lmp) /* ---------------------------------------------------------------------- */ -PairHbondDreidingLJ::~PairHbondDreidingLJ() +PairHbondDreidingLJangleoffset::~PairHbondDreidingLJangleoffset() { memory->sfree(params); delete [] pvector; @@ -76,7 +76,7 @@ PairHbondDreidingLJ::~PairHbondDreidingLJ() /* ---------------------------------------------------------------------- */ -void PairHbondDreidingLJ::compute(int eflag, int vflag) +void PairHbondDreidingLJangleoffset::compute(int eflag, int vflag) { int i,j,k,m,ii,jj,kk,inum,jnum,knum,itype,jtype,ktype,iatom,imol; tagint tagprev; @@ -178,6 +178,11 @@ void PairHbondDreidingLJ::compute(int eflag, int vflag) if (c < -1.0) c = -1.0; ac = acos(c); + ac = ac + pm.angle_offset; + c = cos(ac); + if (c > 1.0) c = 1.0; + if (c < -1.0) c = -1.0; + if (ac > pm.cut_angle && ac < (2.0*MY_PI - pm.cut_angle)) { s = sqrt(1.0 - c*c); if (s < SMALL) s = SMALL; @@ -269,7 +274,7 @@ void PairHbondDreidingLJ::compute(int eflag, int vflag) allocate all arrays ------------------------------------------------------------------------- */ -void PairHbondDreidingLJ::allocate() +void PairHbondDreidingLJangleoffset::allocate() { allocated = 1; int n = atom->ntypes; @@ -298,23 +303,24 @@ void PairHbondDreidingLJ::allocate() global settings ------------------------------------------------------------------------- */ -void PairHbondDreidingLJ::settings(int narg, char **arg) +void PairHbondDreidingLJangleoffset::settings(int narg, char **arg) { - if (narg != 4) error->all(FLERR,"Illegal pair_style command"); + if (narg != 5) error->all(FLERR,"Illegal pair_style command"); ap_global = utils::inumeric(FLERR,arg[0],false,lmp); cut_inner_global = utils::numeric(FLERR,arg[1],false,lmp); cut_outer_global = utils::numeric(FLERR,arg[2],false,lmp); cut_angle_global = utils::numeric(FLERR,arg[3],false,lmp) * MY_PI/180.0; + angle_offset_global = (180.0 - utils::numeric(FLERR,arg[4],false,lmp)) * MY_PI/180.0; } /* ---------------------------------------------------------------------- set coeffs for one or more type pairs ------------------------------------------------------------------------- */ -void PairHbondDreidingLJ::coeff(int narg, char **arg) +void PairHbondDreidingLJangleoffset::coeff(int narg, char **arg) { - if (narg < 6 || narg > 10) + if (narg < 6 || narg > 11) error->all(FLERR,"Incorrect args for pair coefficients"); if (!allocated) allocate(); @@ -342,7 +348,12 @@ void PairHbondDreidingLJ::coeff(int narg, char **arg) if (cut_inner_one>cut_outer_one) error->all(FLERR,"Pair inner cutoff >= Pair outer cutoff"); double cut_angle_one = cut_angle_global; - if (narg == 10) cut_angle_one = utils::numeric(FLERR, arg[9], false, lmp) * MY_PI/180.0; + if (narg > 9) cut_angle_one = utils::numeric(FLERR, arg[9], false, lmp) * MY_PI/180.0; + double angle_offset_one = angle_offset_global; + if (narg == 11) angle_offset_one = (180.0 - utils::numeric(FLERR, arg[10], false, lmp)) * MY_PI/180.0; + if (angle_offset_one < 0.0 || angle_offset_one > 90.0 * MY_PI/180.0) + error->all(FLERR,"Illegal angle offset"); + // grow params array if necessary if (nparams == maxparam) { @@ -364,6 +375,7 @@ void PairHbondDreidingLJ::coeff(int narg, char **arg) params[nparams].cut_innersq = cut_inner_one*cut_inner_one; params[nparams].cut_outersq = cut_outer_one*cut_outer_one; params[nparams].cut_angle = cut_angle_one; + params[nparams].angle_offset = angle_offset_one; params[nparams].denom_vdw = (params[nparams].cut_outersq-params[nparams].cut_innersq) * (params[nparams].cut_outersq-params[nparams].cut_innersq) * @@ -388,7 +400,7 @@ void PairHbondDreidingLJ::coeff(int narg, char **arg) init specific to this pair style ------------------------------------------------------------------------- */ -void PairHbondDreidingLJ::init_style() +void PairHbondDreidingLJangleoffset::init_style() { // molecular system required to use special list to find H atoms // tags required to use special list @@ -448,7 +460,7 @@ void PairHbondDreidingLJ::init_style() init for one type pair i,j and corresponding j,i ------------------------------------------------------------------------- */ -double PairHbondDreidingLJ::init_one(int i, int j) +double PairHbondDreidingLJangleoffset::init_one(int i, int j) { int m; @@ -467,7 +479,7 @@ double PairHbondDreidingLJ::init_one(int i, int j) /* ---------------------------------------------------------------------- */ -double PairHbondDreidingLJ::single(int i, int j, int itype, int jtype, +double PairHbondDreidingLJangleoffset::single(int i, int j, int itype, int jtype, double rsq, double /*factor_coul*/, double /*factor_lj*/, double &fforce) @@ -540,6 +552,11 @@ double PairHbondDreidingLJ::single(int i, int j, int itype, int jtype, if (c < -1.0) c = -1.0; ac = acos(c); + ac = ac + pm.angle_offset; + c = cos(ac); + if (c > 1.0) c = 1.0; + if (c < -1.0) c = -1.0; + if (ac < pm.cut_angle || ac > (2.0*MY_PI - pm.cut_angle)) return 0.0; s = sqrt(1.0 - c*c); if (s < SMALL) s = SMALL; diff --git a/src/EXTRA-MOLECULE/pair_hbond_dreiding_lj_angleoffset.h b/src/EXTRA-MOLECULE/pair_hbond_dreiding_lj_angleoffset.h index 91a906c5cb..331d8e2ac1 100644 --- a/src/EXTRA-MOLECULE/pair_hbond_dreiding_lj_angleoffset.h +++ b/src/EXTRA-MOLECULE/pair_hbond_dreiding_lj_angleoffset.h @@ -13,21 +13,21 @@ #ifdef PAIR_CLASS // clang-format off -PairStyle(hbond/dreiding/lj,PairHbondDreidingLJ); +PairStyle(hbond/dreiding/lj/angleoffset,PairHbondDreidingLJangleoffset); // clang-format on #else -#ifndef LMP_PAIR_HBOND_DREIDING_LJ_H -#define LMP_PAIR_HBOND_DREIDING_LJ_H +#ifndef LMP_PAIR_HBOND_DREIDING_LJ_ANGLEOFFSET_H +#define LMP_PAIR_HBOND_DREIDING_LJ_ANGLEOFFSET_H #include "pair.h" namespace LAMMPS_NS { -class PairHbondDreidingLJ : public Pair { +class PairHbondDreidingLJangleoffset : public Pair { public: - PairHbondDreidingLJ(class LAMMPS *); - ~PairHbondDreidingLJ() override; + PairHbondDreidingLJangleoffset(class LAMMPS *); + ~PairHbondDreidingLJangleoffset() override; void compute(int, int) override; void settings(int, char **) override; void coeff(int, char **) override; @@ -36,7 +36,7 @@ class PairHbondDreidingLJ : public Pair { double single(int, int, int, int, double, double, double, double &) override; protected: - double cut_inner_global, cut_outer_global, cut_angle_global; + double cut_inner_global, cut_outer_global, cut_angle_global, angle_offset_global; int ap_global; struct Param { @@ -45,7 +45,7 @@ class PairHbondDreidingLJ : public Pair { double d0, alpha, r0; double morse1; double denom_vdw; - double cut_inner, cut_outer, cut_innersq, cut_outersq, cut_angle, offset; + double cut_inner, cut_outer, cut_innersq, cut_outersq, cut_angle, offset, angle_offset; int ap; }; diff --git a/src/EXTRA-MOLECULE/pair_hbond_dreiding_morse_angleoffset.cpp b/src/EXTRA-MOLECULE/pair_hbond_dreiding_morse_angleoffset.cpp index c8bc0a627d..1ee4f7ba08 100644 --- a/src/EXTRA-MOLECULE/pair_hbond_dreiding_morse_angleoffset.cpp +++ b/src/EXTRA-MOLECULE/pair_hbond_dreiding_morse_angleoffset.cpp @@ -13,10 +13,10 @@ ------------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- - Contributing author: Tod A Pascal (Caltech) + Contributing author: Tod A Pascal (Caltech), Don Xu/EiPi Fun ------------------------------------------------------------------------- */ -#include "pair_hbond_dreiding_morse.h" +#include "pair_hbond_dreiding_morse_angleoffset.h" #include "atom.h" #include "atom_vec.h" @@ -42,12 +42,12 @@ static constexpr int CHUNK = 8; /* ---------------------------------------------------------------------- */ -PairHbondDreidingMorse::PairHbondDreidingMorse(LAMMPS *lmp) : - PairHbondDreidingLJ(lmp) {} +PairHbondDreidingMorseAngleoffset::PairHbondDreidingMorseAngleoffset(LAMMPS *lmp) : + PairHbondDreidingLJangleoffset(lmp) {} /* ---------------------------------------------------------------------- */ -void PairHbondDreidingMorse::compute(int eflag, int vflag) +void PairHbondDreidingMorseAngleoffset::compute(int eflag, int vflag) { int i,j,k,m,ii,jj,kk,inum,jnum,knum,itype,jtype,ktype,imol,iatom; tagint tagprev; @@ -148,6 +148,11 @@ void PairHbondDreidingMorse::compute(int eflag, int vflag) if (c < -1.0) c = -1.0; ac = acos(c); + ac = ac + pm.angle_offset; + c = cos(ac); + if (c > 1.0) c = 1.0; + if (c < -1.0) c = -1.0; + if (ac > pm.cut_angle && ac < (2.0*MY_PI - pm.cut_angle)) { s = sqrt(1.0 - c*c); if (s < SMALL) s = SMALL; @@ -236,9 +241,9 @@ void PairHbondDreidingMorse::compute(int eflag, int vflag) set coeffs for one or more type pairs ------------------------------------------------------------------------- */ -void PairHbondDreidingMorse::coeff(int narg, char **arg) +void PairHbondDreidingMorseAngleoffset::coeff(int narg, char **arg) { - if (narg < 7 || narg > 11) + if (narg < 7 || narg > 12) error->all(FLERR,"Incorrect args for pair coefficients"); if (!allocated) allocate(); @@ -268,6 +273,10 @@ void PairHbondDreidingMorse::coeff(int narg, char **arg) error->all(FLERR,"Pair inner cutoff >= Pair outer cutoff"); double cut_angle_one = cut_angle_global; if (narg > 10) cut_angle_one = utils::numeric(FLERR, arg[10], false, lmp) * MY_PI/180.0; + double angle_offset_one = angle_offset_global; + if (narg == 12) angle_offset_one = (180.0 - utils::numeric(FLERR,arg[11],false,lmp)) * MY_PI/180.0; + if (angle_offset_one < 0.0 || angle_offset_one > 90.0 * MY_PI/180.0) + error->all(FLERR,"Illegal angle offset"); // grow params array if necessary @@ -291,6 +300,7 @@ void PairHbondDreidingMorse::coeff(int narg, char **arg) params[nparams].cut_innersq = cut_inner_one*cut_inner_one; params[nparams].cut_outersq = cut_outer_one*cut_outer_one; params[nparams].cut_angle = cut_angle_one; + params[nparams].angle_offset = angle_offset_one; params[nparams].denom_vdw = (params[nparams].cut_outersq-params[nparams].cut_innersq) * (params[nparams].cut_outersq-params[nparams].cut_innersq) * @@ -315,7 +325,7 @@ void PairHbondDreidingMorse::coeff(int narg, char **arg) init specific to this pair style ------------------------------------------------------------------------- */ -void PairHbondDreidingMorse::init_style() +void PairHbondDreidingMorseAngleoffset::init_style() { // molecular system required to use special list to find H atoms // tags required to use special list @@ -370,7 +380,7 @@ void PairHbondDreidingMorse::init_style() /* ---------------------------------------------------------------------- */ -double PairHbondDreidingMorse::single(int i, int j, int itype, int jtype, +double PairHbondDreidingMorseAngleoffset::single(int i, int j, int itype, int jtype, double rsq, double /*factor_coul*/, double /*factor_lj*/, double &fforce) @@ -443,6 +453,11 @@ double PairHbondDreidingMorse::single(int i, int j, int itype, int jtype, if (c < -1.0) c = -1.0; ac = acos(c); + ac = ac + pm.angle_offset; + c = cos(ac); + if (c > 1.0) c = 1.0; + if (c < -1.0) c = -1.0; + if (ac < pm.cut_angle || ac > (2.0*MY_PI - pm.cut_angle)) return 0.0; s = sqrt(1.0 - c*c); if (s < SMALL) s = SMALL; diff --git a/src/EXTRA-MOLECULE/pair_hbond_dreiding_morse_angleoffset.h b/src/EXTRA-MOLECULE/pair_hbond_dreiding_morse_angleoffset.h index 8d4646451e..c1b0cf2f30 100644 --- a/src/EXTRA-MOLECULE/pair_hbond_dreiding_morse_angleoffset.h +++ b/src/EXTRA-MOLECULE/pair_hbond_dreiding_morse_angleoffset.h @@ -13,20 +13,20 @@ #ifdef PAIR_CLASS // clang-format off -PairStyle(hbond/dreiding/morse,PairHbondDreidingMorse); +PairStyle(hbond/dreiding/morse/angleoffset,PairHbondDreidingMorseAngleoffset); // clang-format on #else -#ifndef LMP_PAIR_HBOND_DREIDING_MORSE_H -#define LMP_PAIR_HBOND_DREIDING_MORSE_H +#ifndef LMP_PAIR_HBOND_DREIDING_MORSE_ANGLEOFFSET_H +#define LMP_PAIR_HBOND_DREIDING_MORSE_ANGLEOFFSET_H -#include "pair_hbond_dreiding_lj.h" +#include "pair_hbond_dreiding_lj_angleoffset.h" namespace LAMMPS_NS { -class PairHbondDreidingMorse : public PairHbondDreidingLJ { +class PairHbondDreidingMorseAngleoffset : public PairHbondDreidingLJangleoffset { public: - PairHbondDreidingMorse(class LAMMPS *); + PairHbondDreidingMorseAngleoffset(class LAMMPS *); void compute(int, int) override; void coeff(int, char **) override; diff --git a/src/OPENMP/pair_hbond_dreiding_lj_angleoffset_omp.cpp b/src/OPENMP/pair_hbond_dreiding_lj_angleoffset_omp.cpp index b0f6dcfb5b..094beca70a 100644 --- a/src/OPENMP/pair_hbond_dreiding_lj_angleoffset_omp.cpp +++ b/src/OPENMP/pair_hbond_dreiding_lj_angleoffset_omp.cpp @@ -10,10 +10,10 @@ ------------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- - Contributing author: Axel Kohlmeyer (Temple U) + Contributing author: Axel Kohlmeyer (Temple U), Don Xu/EiPi Fun ------------------------------------------------------------------------- */ -#include "pair_hbond_dreiding_lj_omp.h" +#include "pair_hbond_dreiding_lj_angleoffset_omp.h" #include "atom.h" #include "atom_vec.h" @@ -37,8 +37,8 @@ static constexpr double SMALL = 0.001; /* ---------------------------------------------------------------------- */ -PairHbondDreidingLJOMP::PairHbondDreidingLJOMP(LAMMPS *lmp) : - PairHbondDreidingLJ(lmp), ThrOMP(lmp, THR_PAIR) +PairHbondDreidingLJangleoffsetOMP::PairHbondDreidingLJangleoffsetOMP(LAMMPS *lmp) : + PairHbondDreidingLJangleoffset(lmp), ThrOMP(lmp, THR_PAIR) { suffix_flag |= Suffix::OMP; respa_enable = 0; @@ -47,7 +47,7 @@ PairHbondDreidingLJOMP::PairHbondDreidingLJOMP(LAMMPS *lmp) : /* ---------------------------------------------------------------------- */ -PairHbondDreidingLJOMP::~PairHbondDreidingLJOMP() +PairHbondDreidingLJangleoffsetOMP::~PairHbondDreidingLJangleoffsetOMP() { if (hbcount_thr) { delete[] hbcount_thr; @@ -57,7 +57,7 @@ PairHbondDreidingLJOMP::~PairHbondDreidingLJOMP() /* ---------------------------------------------------------------------- */ -void PairHbondDreidingLJOMP::compute(int eflag, int vflag) +void PairHbondDreidingLJangleoffsetOMP::compute(int eflag, int vflag) { ev_init(eflag,vflag); @@ -115,7 +115,7 @@ void PairHbondDreidingLJOMP::compute(int eflag, int vflag) } template -void PairHbondDreidingLJOMP::eval(int iifrom, int iito, ThrData * const thr) +void PairHbondDreidingLJangleoffsetOMP::eval(int iifrom, int iito, ThrData * const thr) { int i,j,k,m,ii,jj,kk,jnum,knum,itype,jtype,ktype,iatom,imol; tagint tagprev; @@ -222,6 +222,11 @@ void PairHbondDreidingLJOMP::eval(int iifrom, int iito, ThrData * const thr) if (c < -1.0) c = -1.0; ac = acos(c); + ac = ac + pm.angle_offset; + c = cos(ac); + if (c > 1.0) c = 1.0; + if (c < -1.0) c = -1.0; + if (ac > pm.cut_angle && ac < (2.0*MY_PI - pm.cut_angle)) { s = sqrt(1.0 - c*c); if (s < SMALL) s = SMALL; @@ -307,11 +312,11 @@ void PairHbondDreidingLJOMP::eval(int iifrom, int iito, ThrData * const thr) /* ---------------------------------------------------------------------- */ -double PairHbondDreidingLJOMP::memory_usage() +double PairHbondDreidingLJangleoffsetOMP::memory_usage() { double bytes = memory_usage_thr(); bytes += (double)comm->nthreads * 2 * sizeof(double); - bytes += PairHbondDreidingLJ::memory_usage(); + bytes += PairHbondDreidingLJangleoffset::memory_usage(); return bytes; } diff --git a/src/OPENMP/pair_hbond_dreiding_lj_angleoffset_omp.h b/src/OPENMP/pair_hbond_dreiding_lj_angleoffset_omp.h index 792a9a6f8c..0a4e14f68d 100644 --- a/src/OPENMP/pair_hbond_dreiding_lj_angleoffset_omp.h +++ b/src/OPENMP/pair_hbond_dreiding_lj_angleoffset_omp.h @@ -12,28 +12,28 @@ ------------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- - Contributing author: Axel Kohlmeyer (Temple U) + Contributing author: Axel Kohlmeyer (Temple U), Don Xu/EiPi Fun ------------------------------------------------------------------------- */ #ifdef PAIR_CLASS // clang-format off -PairStyle(hbond/dreiding/lj/omp,PairHbondDreidingLJOMP); +PairStyle(hbond/dreiding/lj/angleoffset/omp,PairHbondDreidingLJangleoffsetOMP); // clang-format on #else -#ifndef LMP_PAIR_HBOND_DREIDING_LJ_OMP_H -#define LMP_PAIR_HBOND_DREIDING_LJ_OMP_H +#ifndef LMP_PAIR_HBOND_DREIDING_LJ_ANGLEOFFSET_OMP_H +#define LMP_PAIR_HBOND_DREIDING_LJ_ANGLEOFFSET_OMP_H -#include "pair_hbond_dreiding_lj.h" +#include "pair_hbond_dreiding_lj_angleoffset.h" #include "thr_omp.h" namespace LAMMPS_NS { -class PairHbondDreidingLJOMP : public PairHbondDreidingLJ, public ThrOMP { +class PairHbondDreidingLJangleoffsetOMP : public PairHbondDreidingLJangleoffset, public ThrOMP { public: - PairHbondDreidingLJOMP(class LAMMPS *); - ~PairHbondDreidingLJOMP() override; + PairHbondDreidingLJangleoffsetOMP(class LAMMPS *); + ~PairHbondDreidingLJangleoffsetOMP() override; void compute(int, int) override; double memory_usage() override; diff --git a/src/OPENMP/pair_hbond_dreiding_morse_angleoffset_omp.cpp b/src/OPENMP/pair_hbond_dreiding_morse_angleoffset_omp.cpp index 0e43e2a037..ad5c43f7b7 100644 --- a/src/OPENMP/pair_hbond_dreiding_morse_angleoffset_omp.cpp +++ b/src/OPENMP/pair_hbond_dreiding_morse_angleoffset_omp.cpp @@ -10,10 +10,10 @@ ------------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- - Contributing author: Axel Kohlmeyer (Temple U) + Contributing author: Axel Kohlmeyer (Temple U), Don Xu/EiPi Fun ------------------------------------------------------------------------- */ -#include "pair_hbond_dreiding_morse_omp.h" +#include "pair_hbond_dreiding_morse_angleoffset_omp.h" #include "atom.h" #include "atom_vec.h" @@ -37,8 +37,8 @@ static constexpr double SMALL = 0.001; /* ---------------------------------------------------------------------- */ -PairHbondDreidingMorseOMP::PairHbondDreidingMorseOMP(LAMMPS *lmp) : - PairHbondDreidingMorse(lmp), ThrOMP(lmp, THR_PAIR) +PairHbondDreidingMorseAngleoffsetOMP::PairHbondDreidingMorseAngleoffsetOMP(LAMMPS *lmp) : + PairHbondDreidingMorseAngleoffset(lmp), ThrOMP(lmp, THR_PAIR) { suffix_flag |= Suffix::OMP; respa_enable = 0; @@ -47,7 +47,7 @@ PairHbondDreidingMorseOMP::PairHbondDreidingMorseOMP(LAMMPS *lmp) : /* ---------------------------------------------------------------------- */ -PairHbondDreidingMorseOMP::~PairHbondDreidingMorseOMP() +PairHbondDreidingMorseAngleoffsetOMP::~PairHbondDreidingMorseAngleoffsetOMP() { if (hbcount_thr) { delete[] hbcount_thr; @@ -57,7 +57,7 @@ PairHbondDreidingMorseOMP::~PairHbondDreidingMorseOMP() /* ---------------------------------------------------------------------- */ -void PairHbondDreidingMorseOMP::compute(int eflag, int vflag) +void PairHbondDreidingMorseAngleoffsetOMP::compute(int eflag, int vflag) { ev_init(eflag,vflag); @@ -115,7 +115,7 @@ void PairHbondDreidingMorseOMP::compute(int eflag, int vflag) } template -void PairHbondDreidingMorseOMP::eval(int iifrom, int iito, ThrData * const thr) +void PairHbondDreidingMorseAngleoffsetOMP::eval(int iifrom, int iito, ThrData * const thr) { int i,j,k,m,ii,jj,kk,jnum,knum,itype,jtype,ktype,imol,iatom; tagint tagprev; @@ -222,6 +222,11 @@ void PairHbondDreidingMorseOMP::eval(int iifrom, int iito, ThrData * const thr) if (c < -1.0) c = -1.0; ac = acos(c); + ac = ac + pm.angle_offset; + c = cos(ac); + if (c > 1.0) c = 1.0; + if (c < -1.0) c = -1.0; + if (ac > pm.cut_angle && ac < (2.0*MY_PI - pm.cut_angle)) { s = sqrt(1.0 - c*c); if (s < SMALL) s = SMALL; @@ -306,11 +311,11 @@ void PairHbondDreidingMorseOMP::eval(int iifrom, int iito, ThrData * const thr) /* ---------------------------------------------------------------------- */ -double PairHbondDreidingMorseOMP::memory_usage() +double PairHbondDreidingMorseAngleoffsetOMP::memory_usage() { double bytes = memory_usage_thr(); bytes += (double)comm->nthreads * 2 * sizeof(double); - bytes += PairHbondDreidingMorse::memory_usage(); + bytes += PairHbondDreidingMorseAngleoffset::memory_usage(); return bytes; } diff --git a/src/OPENMP/pair_hbond_dreiding_morse_angleoffset_omp.h b/src/OPENMP/pair_hbond_dreiding_morse_angleoffset_omp.h index 016154b762..fb95988c6f 100644 --- a/src/OPENMP/pair_hbond_dreiding_morse_angleoffset_omp.h +++ b/src/OPENMP/pair_hbond_dreiding_morse_angleoffset_omp.h @@ -12,28 +12,28 @@ ------------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- - Contributing author: Axel Kohlmeyer (Temple U) + Contributing author: Axel Kohlmeyer (Temple U), Don Xu/EiPi Fun ------------------------------------------------------------------------- */ #ifdef PAIR_CLASS // clang-format off -PairStyle(hbond/dreiding/morse/omp,PairHbondDreidingMorseOMP); +PairStyle(hbond/dreiding/morse/angleoffset/omp,PairHbondDreidingMorseAngleoffsetOMP); // clang-format on #else -#ifndef LMP_PAIR_HBOND_DREIDING_MORSE_OMP_H -#define LMP_PAIR_HBOND_DREIDING_MORSE_OMP_H +#ifndef LMP_PAIR_HBOND_DREIDING_MORSE_ANGLEOFFSET_OMP_H +#define LMP_PAIR_HBOND_DREIDING_MORSE_ANGLEOFFSET_OMP_H -#include "pair_hbond_dreiding_morse.h" +#include "pair_hbond_dreiding_morse_angleoffset.h" #include "thr_omp.h" namespace LAMMPS_NS { -class PairHbondDreidingMorseOMP : public PairHbondDreidingMorse, public ThrOMP { +class PairHbondDreidingMorseAngleoffsetOMP : public PairHbondDreidingMorseAngleoffset, public ThrOMP { public: - PairHbondDreidingMorseOMP(class LAMMPS *); - ~PairHbondDreidingMorseOMP() override; + PairHbondDreidingMorseAngleoffsetOMP(class LAMMPS *); + ~PairHbondDreidingMorseAngleoffsetOMP() override; void compute(int, int) override; double memory_usage() override; From bb1624b20d772129d93c18b9f3d0e5c4da491c91 Mon Sep 17 00:00:00 2001 From: EiPi Fun Date: Sun, 1 Sep 2024 16:16:47 +0800 Subject: [PATCH 03/29] Add documentation for modified pair_hbond_dreiding with angleoffset --- doc/src/pair_hbond_dreiding_angleoffset.rst | 130 ++++++++++++++++++++ 1 file changed, 130 insertions(+) create mode 100644 doc/src/pair_hbond_dreiding_angleoffset.rst diff --git a/doc/src/pair_hbond_dreiding_angleoffset.rst b/doc/src/pair_hbond_dreiding_angleoffset.rst new file mode 100644 index 0000000000..e827a644ac --- /dev/null +++ b/doc/src/pair_hbond_dreiding_angleoffset.rst @@ -0,0 +1,130 @@ +.. index:: pair_style hbond/dreiding/lj/angleoffset +.. index:: pair_style hbond/dreiding/lj/angleoffset/omp +.. index:: pair_style hbond/dreiding/morse/angleoffset +.. index:: pair_style hbond/dreiding/morse/angleoffest/omp + +pair_style hbond/dreiding/lj/angleoffset command +==================================== + +Accelerator Variants: *hbond/dreiding/lj/angleoffset/omp* + +pair_style hbond/dreiding/morse/angleoffset command +======================================= + +Accelerator Variants: *hbond/dreiding/morse/angleoffset/omp* + +Syntax +"""""" + +.. code-block:: LAMMPS + + pair_style style N inner_distance_cutoff outer_distance_cutoff angle_cutoff equilibrium_angle + +* style = *hbond/dreiding/lj/angleoffset* or *hbond/dreiding/morse/angleoffset* +* N = power of cosine of sum of angle theta and angle offset (integer) +* inner_distance_cutoff = global inner cutoff for Donor-Acceptor interactions (distance units) +* outer_distance_cutoff = global cutoff for Donor-Acceptor interactions (distance units) +* angle_cutoff = global angle cutoff for Acceptor-Hydrogen-Donor interactions (degrees) +* equilibrium_angle = global equilibrium angle for Acceptor-Hydrogen-Donor interactions (degrees) + +(Tips: angle offset is the supplementary angle of equilibrium angle. It means if equilibrium angle is 166.6, the angle offset is 13.4) + +Examples +"""""""" + +.. code-block:: LAMMPS + + pair_style hybrid/overlay lj/cut 10.0 hbond/dreiding/lj/angleoffset 4 9.0 11.0 90.0 170.0 + pair_coeff 1 2 hbond/dreiding/lj 3 i 9.5 2.75 4 9.0 11.0 90.0 160.0 + + pair_style hybrid/overlay lj/cut 10.0 hbond/dreiding/morse/angleoffset 2 9.0 11.0 90.0 170.0 + pair_coeff 1 2 hbond/dreiding/morse 3 i 3.88 1.7241379 2.9 2 9.0 11.0 90.0 160.0 + + labelmap atom 1 C 2 O 3 H + pair_coeff C O hbond/dreiding/morse H i 3.88 1.7241379 2.9 2 9.0 11.0 90.0 160.0 + +Description +""""""""""" + +The *hbond/dreiding/\*/angleoffset* styles are modified version of *hbond/dreiding* styles. + +In some cases, the angle of acceptor-hydrogen-donor in the equilibrium state could not achieve 180 degree especially in some coarse grained models. +In these cases, an angle offset is required to ensure that equilibrium state could be the minimum energy state. + +The *hbond/dreiding/\*/angleoffset* styles compute the Acceptor-Hydrogen-Donor (AHD) +3-body hydrogen bond interaction for the :doc:`DREIDING ` +force field with an angle offset, given by: + +.. math:: + + E = & \left[LJ(r) | Morse(r) \right] \qquad \qquad \qquad r < r_{\rm in} and \theta + \theta_{offset} < \theta_c \\ + = & S(r) * \left[LJ(r) | Morse(r) \right] \qquad \qquad r_{\rm in} < r < r_{\rm out} and \theta + \theta_{offset} < \theta_c \\ + = & 0 \qquad \qquad \qquad \qquad \qquad \qquad \qquad r > r_{\rm out} and \theta + \theta_{offset} < \theta_c \\ + LJ(r) = & AR^{-12}-BR^{-10}cos^n(\theta + \theta_{offset})= + \epsilon\left\lbrace 5\left[ \frac{\sigma}{r}\right]^{12}- + 6\left[ \frac{\sigma}{r}\right]^{10} \right\rbrace cos^n(\theta + \theta_{offset})\\ + Morse(r) = & D_0\left\lbrace \chi^2 - 2\chi\right\rbrace cos^n(\theta + \theta_{offset})= + D_{0}\left\lbrace e^{- 2 \alpha (r - r_0)} - 2 e^{- \alpha (r - r_0)} + \right\rbrace cos^n(\theta + \theta_{offset})\ + S(r) = & \frac{ \left[r_{\rm out}^2 - r^2\right]^2 + \left[r_{\rm out}^2 + 2r^2 - 3{r_{\rm in}^2}\right]} + { \left[r_{\rm out}^2 - {r_{\rm in}}^2\right]^3 } + +where :math:`r_{\rm in}` is the inner spline distance cutoff, +:math:`r_{\rm out}` is the outer distance cutoff, :math:`\theta_c` is +the angle cutoff, :math:`\theta_offset` is the angle offset, and :math:`n` is the power of the cosine of the sum of the angle :math:`\theta` and the angle offset :math:`\theta_offset` + +Here, *r* is the radial distance between the donor (D) and acceptor +(A) atoms and :math:`\theta` is the bond angle between the acceptor, the +hydrogen (H) and the donor atoms: + +.. image:: JPG/dreiding_hbond.jpg + :align: center + +For the *hbond/dreiding/lj/angleoffset* style the list of coefficients is as +follows: + +* K = hydrogen atom type = 1 to Ntypes, or type label +* donor flag = *i* or *j* +* :math:`\epsilon` (energy units) +* :math:`\sigma` (distance units) +* *n* = exponent in formula above +* distance cutoff :math:`r_{\rm in}` (distance units) +* distance cutoff :math:`r_{\rm out}` (distance units) +* angle cutoff (degrees) +* equilibrium angle (degrees) + +(Tips: angle offset is the supplementary angle of equilibrium angle) + +For the *hbond/dreiding/morse/angleoffset* style the list of coefficients is as +follows: + +* K = hydrogen atom type = 1 to Ntypes, or type label +* donor flag = *i* or *j* +* :math:`D_0` (energy units) +* :math:`\alpha` (1/distance units) +* :math:`r_0` (distance units) +* *n* = exponent in formula above +* distance cutoff :math:`r_{\rm in}` (distance units) +* distance cutoff :math:`r_{out}` (distance units) +* angle cutoff (degrees) +* equilibrium angle (degrees) + +(Tips: angle offset is the supplementary angle of equilibrium angle) + +---------- + +Additional Information +"""""""""""" + +For more information about DREIDING force field and other notes, please refer to the documentation of *hbond/dreiding* styles. + +---------- + +Restrictions +"""""""""""" + +This pair style can only be used if LAMMPS was built with the +EXTRA-MOLECULE package. See the :doc:`Build package ` doc page +for more info. + From 9f73494c91f82cd46079ac3496fdee006be26a77 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Mon, 2 Sep 2024 20:42:05 -0400 Subject: [PATCH 04/29] correct ReStructuredText and LaTeX formatting issues --- doc/src/pair_hbond_dreiding_angleoffset.rst | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/doc/src/pair_hbond_dreiding_angleoffset.rst b/doc/src/pair_hbond_dreiding_angleoffset.rst index e827a644ac..5aa8950961 100644 --- a/doc/src/pair_hbond_dreiding_angleoffset.rst +++ b/doc/src/pair_hbond_dreiding_angleoffset.rst @@ -4,12 +4,12 @@ .. index:: pair_style hbond/dreiding/morse/angleoffest/omp pair_style hbond/dreiding/lj/angleoffset command -==================================== +================================================ Accelerator Variants: *hbond/dreiding/lj/angleoffset/omp* pair_style hbond/dreiding/morse/angleoffset command -======================================= +=================================================== Accelerator Variants: *hbond/dreiding/morse/angleoffset/omp* @@ -65,7 +65,7 @@ force field with an angle offset, given by: 6\left[ \frac{\sigma}{r}\right]^{10} \right\rbrace cos^n(\theta + \theta_{offset})\\ Morse(r) = & D_0\left\lbrace \chi^2 - 2\chi\right\rbrace cos^n(\theta + \theta_{offset})= D_{0}\left\lbrace e^{- 2 \alpha (r - r_0)} - 2 e^{- \alpha (r - r_0)} - \right\rbrace cos^n(\theta + \theta_{offset})\ + \right\rbrace cos^n(\theta + \theta_{offset})\\ S(r) = & \frac{ \left[r_{\rm out}^2 - r^2\right]^2 \left[r_{\rm out}^2 + 2r^2 - 3{r_{\rm in}^2}\right]} { \left[r_{\rm out}^2 - {r_{\rm in}}^2\right]^3 } @@ -115,7 +115,7 @@ follows: ---------- Additional Information -"""""""""""" +"""""""""""""""""""""" For more information about DREIDING force field and other notes, please refer to the documentation of *hbond/dreiding* styles. @@ -125,6 +125,6 @@ Restrictions """""""""""" This pair style can only be used if LAMMPS was built with the -EXTRA-MOLECULE package. See the :doc:`Build package ` doc page -for more info. +EXTRA-MOLECULE package. See the :doc:`Build package ` +doc page for more info. From 7c13562c4101635face95b9ea34a3c9db91c4607 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Mon, 2 Sep 2024 20:42:15 -0400 Subject: [PATCH 05/29] spelling --- doc/utils/sphinx-config/false_positives.txt | 1 + 1 file changed, 1 insertion(+) diff --git a/doc/utils/sphinx-config/false_positives.txt b/doc/utils/sphinx-config/false_positives.txt index cfbddbe5f6..7b9ff540f2 100644 --- a/doc/utils/sphinx-config/false_positives.txt +++ b/doc/utils/sphinx-config/false_positives.txt @@ -107,6 +107,7 @@ Andrienko Andzelm Ang anglegrad +angleoffset angletangrad angmom angmomx From 9f2e542c807f1c0802ec172eddde5f873a534410 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Mon, 2 Sep 2024 20:42:33 -0400 Subject: [PATCH 06/29] integrate into manual build --- doc/src/Commands_pair.rst | 2 ++ doc/src/pair_style.rst | 2 ++ 2 files changed, 4 insertions(+) diff --git a/doc/src/Commands_pair.rst b/doc/src/Commands_pair.rst index dfed8f7485..a813bb771d 100644 --- a/doc/src/Commands_pair.rst +++ b/doc/src/Commands_pair.rst @@ -113,7 +113,9 @@ OPT. * :doc:`gw/zbl ` * :doc:`harmonic/cut (o) ` * :doc:`hbond/dreiding/lj (o) ` + * :doc:`hbond/dreiding/lj/angleoffset (o) ` * :doc:`hbond/dreiding/morse (o) ` + * :doc:`hbond/dreiding/morse/angleoffset (o) ` * :doc:`hdnnp ` * :doc:`hippo (g) ` * :doc:`ilp/graphene/hbn (t) ` diff --git a/doc/src/pair_style.rst b/doc/src/pair_style.rst index 51350c86a4..6e59d697a8 100644 --- a/doc/src/pair_style.rst +++ b/doc/src/pair_style.rst @@ -205,7 +205,9 @@ accelerated styles exist. * :doc:`gw/zbl ` - Gao-Weber potential with a repulsive ZBL core * :doc:`harmonic/cut ` - repulsive-only harmonic potential * :doc:`hbond/dreiding/lj ` - DREIDING hydrogen bonding LJ potential +* :doc:`hbond/dreiding/lj/angleoffset ` - DREIDING hydrogen bonding LJ potential with offset for hbond angle * :doc:`hbond/dreiding/morse ` - DREIDING hydrogen bonding Morse potential +* :doc:`hbond/dreiding/morse/angleoffset ` - DREIDING hydrogen bonding Morse potential with offset for hbond angle * :doc:`hdnnp ` - High-dimensional neural network potential * :doc:`hippo ` - * :doc:`ilp/graphene/hbn ` - registry-dependent interlayer potential (ILP) From 6aa592d28696a39b326c64c8458b9c5baddbdfc1 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Mon, 2 Sep 2024 20:43:48 -0400 Subject: [PATCH 07/29] build system integration --- src/.gitignore | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/src/.gitignore b/src/.gitignore index c26eaaba30..d706ef4a7d 100644 --- a/src/.gitignore +++ b/src/.gitignore @@ -1262,6 +1262,10 @@ /pair_hbond_dreiding_lj.h /pair_hbond_dreiding_morse.cpp /pair_hbond_dreiding_morse.h +/pair_hbond_dreiding_lj_angleoffset.cpp +/pair_hbond_dreiding_lj_angleoffset.h +/pair_hbond_dreiding_morse_angleoffset.cpp +/pair_hbond_dreiding_morse_angleoffset.h /pair_hdnnp.cpp /pair_hdnnp.h /pair_ilp_graphene_hbn.cpp From e763f9e052e893c6e79f1624132f9d0506e96cec Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Mon, 2 Sep 2024 20:49:41 -0400 Subject: [PATCH 08/29] use correct style in error messages --- .../pair_hbond_dreiding_lj_angleoffset.cpp | 18 ++++---------- .../pair_hbond_dreiding_morse_angleoffset.cpp | 24 ++++++------------- 2 files changed, 12 insertions(+), 30 deletions(-) diff --git a/src/EXTRA-MOLECULE/pair_hbond_dreiding_lj_angleoffset.cpp b/src/EXTRA-MOLECULE/pair_hbond_dreiding_lj_angleoffset.cpp index ab15d7c0b5..3b30baa2b3 100644 --- a/src/EXTRA-MOLECULE/pair_hbond_dreiding_lj_angleoffset.cpp +++ b/src/EXTRA-MOLECULE/pair_hbond_dreiding_lj_angleoffset.cpp @@ -408,14 +408,14 @@ void PairHbondDreidingLJangleoffset::init_style() // and computing forces on A,H which may be on different procs if (atom->molecular == Atom::ATOMIC) - error->all(FLERR,"Pair style hbond/dreiding requires molecular system"); + error->all(FLERR,"Pair style hbond/dreiding/lj/angleoffset requires molecular system"); if (atom->tag_enable == 0) - error->all(FLERR,"Pair style hbond/dreiding requires atom IDs"); + error->all(FLERR,"Pair style hbond/dreiding/lj/angleoffset requires atom IDs"); if (atom->map_style == Atom::MAP_NONE) - error->all(FLERR,"Pair style hbond/dreiding requires an atom map, " + error->all(FLERR,"Pair style hbond/dreiding/lj/angleoffset requires an atom map, " "see atom_modify"); if (force->newton_pair == 0) - error->all(FLERR,"Pair style hbond/dreiding requires newton pair on"); + error->all(FLERR,"Pair style hbond/dreiding/lj/angleoffset requires newton pair on"); // set donor[M]/acceptor[M] if any atom of type M is a donor/acceptor @@ -431,7 +431,7 @@ void PairHbondDreidingLJangleoffset::init_style() acceptor[j] = 1; } - if (!anyflag) error->all(FLERR,"No pair hbond/dreiding coefficients set"); + if (!anyflag) error->all(FLERR,"No pair hbond/dreiding/lj/angleoffset coefficients set"); // set additional param values // offset is for LJ only, angle term is not included @@ -441,14 +441,6 @@ void PairHbondDreidingLJangleoffset::init_style() params[m].lj2 = 60.0*params[m].epsilon*pow(params[m].sigma,10.0); params[m].lj3 = 5.0*params[m].epsilon*pow(params[m].sigma,12.0); params[m].lj4 = 6.0*params[m].epsilon*pow(params[m].sigma,10.0); - - /* - if (offset_flag) { - double ratio = params[m].sigma / params[m].cut_outer; - params[m].offset = params[m].epsilon * - ((2.0*pow(ratio,9.0)) - (3.0*pow(ratio,6.0))); - } else params[m].offset = 0.0; - */ } // full neighbor list request diff --git a/src/EXTRA-MOLECULE/pair_hbond_dreiding_morse_angleoffset.cpp b/src/EXTRA-MOLECULE/pair_hbond_dreiding_morse_angleoffset.cpp index 1ee4f7ba08..173eac605f 100644 --- a/src/EXTRA-MOLECULE/pair_hbond_dreiding_morse_angleoffset.cpp +++ b/src/EXTRA-MOLECULE/pair_hbond_dreiding_morse_angleoffset.cpp @@ -282,8 +282,7 @@ void PairHbondDreidingMorseAngleoffset::coeff(int narg, char **arg) if (nparams == maxparam) { maxparam += CHUNK; - params = (Param *) memory->srealloc(params, maxparam*sizeof(Param), - "pair:params"); + params = (Param *) memory->srealloc(params, maxparam*sizeof(Param),"pair:params"); // make certain all addional allocated storage is initialized // to avoid false positives when checking with valgrind @@ -301,8 +300,7 @@ void PairHbondDreidingMorseAngleoffset::coeff(int narg, char **arg) params[nparams].cut_outersq = cut_outer_one*cut_outer_one; params[nparams].cut_angle = cut_angle_one; params[nparams].angle_offset = angle_offset_one; - params[nparams].denom_vdw = - (params[nparams].cut_outersq-params[nparams].cut_innersq) * + params[nparams].denom_vdw = (params[nparams].cut_outersq-params[nparams].cut_innersq) * (params[nparams].cut_outersq-params[nparams].cut_innersq) * (params[nparams].cut_outersq-params[nparams].cut_innersq); @@ -333,14 +331,14 @@ void PairHbondDreidingMorseAngleoffset::init_style() // and computing forces on A,H which may be on different procs if (atom->molecular == Atom::ATOMIC) - error->all(FLERR,"Pair style hbond/dreiding requires molecular system"); + error->all(FLERR,"Pair style hbond/dreiding/morse/angleoffset requires molecular system"); if (atom->tag_enable == 0) - error->all(FLERR,"Pair style hbond/dreiding requires atom IDs"); + error->all(FLERR,"Pair style hbond/dreiding/morse/angleoffset requires atom IDs"); if (atom->map_style == Atom::MAP_NONE) - error->all(FLERR,"Pair style hbond/dreiding requires an atom map, " + error->all(FLERR,"Pair style hbond/dreiding/morse/angleoffset requires an atom map, " "see atom_modify"); if (force->newton_pair == 0) - error->all(FLERR,"Pair style hbond/dreiding requires newton pair on"); + error->all(FLERR,"Pair style hbond/dreiding/morse/angleoffset requires newton pair on"); // set donor[M]/acceptor[M] if any atom of type M is a donor/acceptor @@ -356,21 +354,13 @@ void PairHbondDreidingMorseAngleoffset::init_style() acceptor[j] = 1; } - if (!anyflag) error->all(FLERR,"No pair hbond/dreiding coefficients set"); + if (!anyflag) error->all(FLERR,"No pair hbond/dreiding/morse/angleoffset coefficients set"); // set additional param values // offset is for Morse only, angle term is not included for (int m = 0; m < nparams; m++) { params[m].morse1 = 2.0*params[m].d0*params[m].alpha; - - /* - if (offset_flag) { - double alpha_dr = -params[m].alpha * (params[m].cut - params[m].r0); - params[m].offset = params[m].d0 * - ((exp(2.0*alpha_dr)) - (2.0*exp(alpha_dr))); - } else params[m].offset = 0.0; - */ } // full neighbor list request From 9cfbf3dcdd46e91570fb24b7102ef3fb06c8e1f0 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Mon, 2 Sep 2024 20:58:40 -0400 Subject: [PATCH 09/29] fix link and grammar --- doc/src/pair_hbond_dreiding_angleoffset.rst | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/doc/src/pair_hbond_dreiding_angleoffset.rst b/doc/src/pair_hbond_dreiding_angleoffset.rst index 5aa8950961..888be1bb10 100644 --- a/doc/src/pair_hbond_dreiding_angleoffset.rst +++ b/doc/src/pair_hbond_dreiding_angleoffset.rst @@ -27,7 +27,8 @@ Syntax * angle_cutoff = global angle cutoff for Acceptor-Hydrogen-Donor interactions (degrees) * equilibrium_angle = global equilibrium angle for Acceptor-Hydrogen-Donor interactions (degrees) -(Tips: angle offset is the supplementary angle of equilibrium angle. It means if equilibrium angle is 166.6, the angle offset is 13.4) +(Tips: angle offset is the supplementary angle of the equilibrium angle. It means if the +equilibrium angle is 166.6, the angle offset is 13.4) Examples """""""" @@ -117,7 +118,8 @@ follows: Additional Information """""""""""""""""""""" -For more information about DREIDING force field and other notes, please refer to the documentation of *hbond/dreiding* styles. +For more information about DREIDING force field and other notes, please refer +to the :doc:`documentation of the *hbond/dreiding* styles `. ---------- From 61ffe1ece17e09c01dd0d2e585f801f5d794ddfe Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Mon, 2 Sep 2024 21:00:02 -0400 Subject: [PATCH 10/29] correct link --- doc/src/pair_hbond_dreiding_angleoffset.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/doc/src/pair_hbond_dreiding_angleoffset.rst b/doc/src/pair_hbond_dreiding_angleoffset.rst index 888be1bb10..8262937db2 100644 --- a/doc/src/pair_hbond_dreiding_angleoffset.rst +++ b/doc/src/pair_hbond_dreiding_angleoffset.rst @@ -119,7 +119,7 @@ Additional Information """""""""""""""""""""" For more information about DREIDING force field and other notes, please refer -to the :doc:`documentation of the *hbond/dreiding* styles `. +to the :doc:`documentation of the *hbond/dreiding* pair styles `. ---------- From 4b56e81b66e5ed603ee3876b3e3be61ed6cc3e57 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Mon, 2 Sep 2024 21:03:16 -0400 Subject: [PATCH 11/29] fix typo --- doc/src/pair_hbond_dreiding_angleoffset.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/doc/src/pair_hbond_dreiding_angleoffset.rst b/doc/src/pair_hbond_dreiding_angleoffset.rst index 8262937db2..682dd7564a 100644 --- a/doc/src/pair_hbond_dreiding_angleoffset.rst +++ b/doc/src/pair_hbond_dreiding_angleoffset.rst @@ -1,7 +1,7 @@ .. index:: pair_style hbond/dreiding/lj/angleoffset .. index:: pair_style hbond/dreiding/lj/angleoffset/omp .. index:: pair_style hbond/dreiding/morse/angleoffset -.. index:: pair_style hbond/dreiding/morse/angleoffest/omp +.. index:: pair_style hbond/dreiding/morse/angleoffset/omp pair_style hbond/dreiding/lj/angleoffset command ================================================ From a48b67baae4b037d68fd67229fcd83a3a83168ab Mon Sep 17 00:00:00 2001 From: EiPiFun <132569654+EiPiFun@users.noreply.github.com> Date: Tue, 3 Sep 2024 11:09:35 +0800 Subject: [PATCH 12/29] Fix angle cutoff logical error Angle should > angle cutoff instead of < --- doc/src/pair_hbond_dreiding_angleoffset.rst | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/doc/src/pair_hbond_dreiding_angleoffset.rst b/doc/src/pair_hbond_dreiding_angleoffset.rst index 682dd7564a..dff6040218 100644 --- a/doc/src/pair_hbond_dreiding_angleoffset.rst +++ b/doc/src/pair_hbond_dreiding_angleoffset.rst @@ -58,9 +58,9 @@ force field with an angle offset, given by: .. math:: - E = & \left[LJ(r) | Morse(r) \right] \qquad \qquad \qquad r < r_{\rm in} and \theta + \theta_{offset} < \theta_c \\ - = & S(r) * \left[LJ(r) | Morse(r) \right] \qquad \qquad r_{\rm in} < r < r_{\rm out} and \theta + \theta_{offset} < \theta_c \\ - = & 0 \qquad \qquad \qquad \qquad \qquad \qquad \qquad r > r_{\rm out} and \theta + \theta_{offset} < \theta_c \\ + E = & \left[LJ(r) | Morse(r) \right] \qquad \qquad \qquad r < r_{\rm in} and \theta + \theta_{offset} > \theta_c \\ + = & S(r) * \left[LJ(r) | Morse(r) \right] \qquad \qquad r_{\rm in} < r < r_{\rm out} and \theta + \theta_{offset} > \theta_c \\ + = & 0 \qquad \qquad \qquad \qquad \qquad \qquad \qquad r > r_{\rm out} and \theta + \theta_{offset} > \theta_c \\ LJ(r) = & AR^{-12}-BR^{-10}cos^n(\theta + \theta_{offset})= \epsilon\left\lbrace 5\left[ \frac{\sigma}{r}\right]^{12}- 6\left[ \frac{\sigma}{r}\right]^{10} \right\rbrace cos^n(\theta + \theta_{offset})\\ From 787c49d8415a15d07e7565762b25b7ef0c00ad73 Mon Sep 17 00:00:00 2001 From: EiPi Fun Date: Tue, 3 Sep 2024 21:14:10 +0800 Subject: [PATCH 13/29] Add information about mixing and shift into doc and improve code format --- doc/src/pair_hbond_dreiding_angleoffset.rst | 6 ++++++ .../pair_hbond_dreiding_lj_angleoffset.cpp | 12 ++++++------ .../pair_hbond_dreiding_morse_angleoffset.cpp | 2 +- 3 files changed, 13 insertions(+), 7 deletions(-) diff --git a/doc/src/pair_hbond_dreiding_angleoffset.rst b/doc/src/pair_hbond_dreiding_angleoffset.rst index dff6040218..97bf93ffe4 100644 --- a/doc/src/pair_hbond_dreiding_angleoffset.rst +++ b/doc/src/pair_hbond_dreiding_angleoffset.rst @@ -118,6 +118,12 @@ follows: Additional Information """""""""""""""""""""" +These pair styles do not support mixing. You must explicitly identify +each donor/acceptor type pair. + +These styles do not support the :doc:`pair_modify ` shift +option for the energy of the interactions. + For more information about DREIDING force field and other notes, please refer to the :doc:`documentation of the *hbond/dreiding* pair styles `. diff --git a/src/EXTRA-MOLECULE/pair_hbond_dreiding_lj_angleoffset.cpp b/src/EXTRA-MOLECULE/pair_hbond_dreiding_lj_angleoffset.cpp index 3b30baa2b3..68cf24c5ac 100644 --- a/src/EXTRA-MOLECULE/pair_hbond_dreiding_lj_angleoffset.cpp +++ b/src/EXTRA-MOLECULE/pair_hbond_dreiding_lj_angleoffset.cpp @@ -305,13 +305,13 @@ void PairHbondDreidingLJangleoffset::allocate() void PairHbondDreidingLJangleoffset::settings(int narg, char **arg) { - if (narg != 5) error->all(FLERR,"Illegal pair_style command"); + if (narg != 5) error->all(FLERR, "Illegal pair_style command"); - ap_global = utils::inumeric(FLERR,arg[0],false,lmp); - cut_inner_global = utils::numeric(FLERR,arg[1],false,lmp); - cut_outer_global = utils::numeric(FLERR,arg[2],false,lmp); - cut_angle_global = utils::numeric(FLERR,arg[3],false,lmp) * MY_PI/180.0; - angle_offset_global = (180.0 - utils::numeric(FLERR,arg[4],false,lmp)) * MY_PI/180.0; + ap_global = utils::inumeric(FLERR, arg[0], false, lmp); + cut_inner_global = utils::numeric(FLERR, arg[1], false, lmp); + cut_outer_global = utils::numeric(FLERR, arg[2], false, lmp); + cut_angle_global = utils::numeric(FLERR, arg[3], false, lmp) * MY_PI/180.0; + angle_offset_global = (180.0 - utils::numeric(FLERR, arg[4], false, lmp)) * MY_PI/180.0; } /* ---------------------------------------------------------------------- diff --git a/src/EXTRA-MOLECULE/pair_hbond_dreiding_morse_angleoffset.cpp b/src/EXTRA-MOLECULE/pair_hbond_dreiding_morse_angleoffset.cpp index 173eac605f..2a1b6f1142 100644 --- a/src/EXTRA-MOLECULE/pair_hbond_dreiding_morse_angleoffset.cpp +++ b/src/EXTRA-MOLECULE/pair_hbond_dreiding_morse_angleoffset.cpp @@ -274,7 +274,7 @@ void PairHbondDreidingMorseAngleoffset::coeff(int narg, char **arg) double cut_angle_one = cut_angle_global; if (narg > 10) cut_angle_one = utils::numeric(FLERR, arg[10], false, lmp) * MY_PI/180.0; double angle_offset_one = angle_offset_global; - if (narg == 12) angle_offset_one = (180.0 - utils::numeric(FLERR,arg[11],false,lmp)) * MY_PI/180.0; + if (narg == 12) angle_offset_one = (180.0 - utils::numeric(FLERR, arg[11], false, lmp)) * MY_PI/180.0; if (angle_offset_one < 0.0 || angle_offset_one > 90.0 * MY_PI/180.0) error->all(FLERR,"Illegal angle offset"); From c6e9e90e15a053c1766002250e5b30c711743873 Mon Sep 17 00:00:00 2001 From: megmcca Date: Mon, 16 Dec 2024 09:37:21 -0700 Subject: [PATCH 14/29] refactor to mirror airebo/morse class structure --- doc/src/pair_hbond_dreiding.rst | 40 +- doc/src/pair_hbond_dreiding_angleoffset.rst | 138 ------ .../pair_hbond_dreiding_lj_angleoffset.cpp | 464 +----------------- .../pair_hbond_dreiding_lj_angleoffset.h | 43 +- .../pair_hbond_dreiding_morse_angleoffset.cpp | 375 +------------- .../pair_hbond_dreiding_morse_angleoffset.h | 14 +- src/MOLECULE/pair_hbond_dreiding_lj.cpp | 43 +- src/MOLECULE/pair_hbond_dreiding_lj.h | 4 +- src/MOLECULE/pair_hbond_dreiding_morse.cpp | 16 + 9 files changed, 147 insertions(+), 990 deletions(-) delete mode 100644 doc/src/pair_hbond_dreiding_angleoffset.rst diff --git a/doc/src/pair_hbond_dreiding.rst b/doc/src/pair_hbond_dreiding.rst index 7e73f23b08..33f81f2439 100644 --- a/doc/src/pair_hbond_dreiding.rst +++ b/doc/src/pair_hbond_dreiding.rst @@ -1,31 +1,46 @@ .. index:: pair_style hbond/dreiding/lj .. index:: pair_style hbond/dreiding/lj/omp +.. index:: pair_style hbond/dreiding/lj/angleoffset +.. index:: pair_style hbond/dreiding/lj/angleoffset/omp .. index:: pair_style hbond/dreiding/morse .. index:: pair_style hbond/dreiding/morse/omp +.. index:: pair_style hbond/dreiding/morse/angleoffset +.. index:: pair_style hbond/dreiding/morse/angleoffset/omp pair_style hbond/dreiding/lj command ==================================== Accelerator Variants: *hbond/dreiding/lj/omp* +pair_style hbond/dreiding/lj/angleoffset command +================================================ + +Accelerator Variants: *hbond/dreiding/lj/angleoffset/omp* + pair_style hbond/dreiding/morse command ======================================= Accelerator Variants: *hbond/dreiding/morse/omp* +pair_style hbond/dreiding/morse/angleoffset command +=================================================== + +Accelerator Variants: *hbond/dreiding/morse/angleoffset/omp* + + Syntax """""" .. code-block:: LAMMPS - pair_style style N inner_distance_cutoff outer_distance_cutoff angle_cutof + pair_style style N inner_distance_cutoff outer_distance_cutoff angle_cutoff equilibrium_angle -* style = *hbond/dreiding/lj* or *hbond/dreiding/morse* +* style = *hbond/dreiding/lj* or *hbond/dreiding/morse* or *hbond/dreiding/lj/angleoffset* or *hbond/dreiding/morse/angleoffset* * n = cosine angle periodicity * inner_distance_cutoff = global inner cutoff for Donor-Acceptor interactions (distance units) * outer_distance_cutoff = global cutoff for Donor-Acceptor interactions (distance units) -* angle_cutoff = global angle cutoff for Acceptor-Hydrogen-Donor -* interactions (degrees) +* angle_cutoff = global angle cutoff for Acceptor-Hydrogen-Donor interactions (degrees) +* (with style angleoffset) equilibrium_angle = global equilibrium angle for Acceptor-Hydrogen-Donor interactions (degrees) Examples """""""" @@ -41,6 +56,9 @@ Examples labelmap atom 1 C 2 O 3 H pair_coeff C O hbond/dreiding/morse H i 3.88 1.7241379 2.9 2 9 11 90 + pair_style hybrid/overlay lj/cut 10.0 hbond/dreiding/lj 4 9.0 11.0 90 170.0 + pair_coeff 1 2 hbond/dreiding/lj 3 i 9.5 2.75 4 9.0 11.0 90.0 + Description """"""""""" @@ -91,6 +109,8 @@ potential for the Donor-Acceptor interactions. :ref:`(Liu) ` showed that the Morse form gives improved results for Dendrimer simulations, when n = 2. +The style variants *hbond/dreiding/lj/angleoffset* and *hbond/dreiding/lj/angleoffset* takes the equilibrium angle of the AHD as input, allowing it to reach 180 degrees. This variant option was added to account for cases (especially in some coarse-grained models) in which the equilibrium state of the bonds may equal the minimum energy state. + See the :doc:`Howto bioFF ` page for more information on the DREIDING force field. @@ -119,6 +139,10 @@ on the DREIDING force field. special_bonds command (e.g. "special_bonds lj 0.0 0.0 1.0") to turn these interactions on. +.. note:: + + For the *angle_offset* variants, the referenced angle offset is the supplementary angle of the equilibrium angle parameter. It means if the equilibrium angle is 166.6 degrees, the calculated angle offset is 13.4 degrees. + ---------- The following coefficients must be defined for pairs of eligible @@ -169,7 +193,10 @@ follows: * distance cutoff :math:`r_{out}` (distance units) * angle cutoff (degrees) -A single hydrogen atom type K can be specified, or a wild-card asterisk +For both the *hbond/dreiding/lj/angleoffset* and *hbond/dreiding/morse/angleoffset* styles an additional parameter is added: +* equilibrium angle (degrees) + +For all styles, a single hydrogen atom type K can be specified, or a wild-card asterisk can be used in place of or in conjunction with the K arguments to select multiple types as hydrogen atoms. This takes the form "\*" or "\*n" or "n\*" or "m\*n". See the :doc:`pair_coeff ` @@ -244,8 +271,7 @@ heading) the following commands could be included in an input script: Restrictions """""""""""" -This pair style can only be used if LAMMPS was built with the -MOLECULE package. See the :doc:`Build package ` doc page +The base pair styles can only be used if LAMMPS was built with the MOLECULE package. The *angleoffset* variant also requires the EXTRA-MOLECULE package. See the :doc:`Build package ` doc page for more info. Related commands diff --git a/doc/src/pair_hbond_dreiding_angleoffset.rst b/doc/src/pair_hbond_dreiding_angleoffset.rst deleted file mode 100644 index 97bf93ffe4..0000000000 --- a/doc/src/pair_hbond_dreiding_angleoffset.rst +++ /dev/null @@ -1,138 +0,0 @@ -.. index:: pair_style hbond/dreiding/lj/angleoffset -.. index:: pair_style hbond/dreiding/lj/angleoffset/omp -.. index:: pair_style hbond/dreiding/morse/angleoffset -.. index:: pair_style hbond/dreiding/morse/angleoffset/omp - -pair_style hbond/dreiding/lj/angleoffset command -================================================ - -Accelerator Variants: *hbond/dreiding/lj/angleoffset/omp* - -pair_style hbond/dreiding/morse/angleoffset command -=================================================== - -Accelerator Variants: *hbond/dreiding/morse/angleoffset/omp* - -Syntax -"""""" - -.. code-block:: LAMMPS - - pair_style style N inner_distance_cutoff outer_distance_cutoff angle_cutoff equilibrium_angle - -* style = *hbond/dreiding/lj/angleoffset* or *hbond/dreiding/morse/angleoffset* -* N = power of cosine of sum of angle theta and angle offset (integer) -* inner_distance_cutoff = global inner cutoff for Donor-Acceptor interactions (distance units) -* outer_distance_cutoff = global cutoff for Donor-Acceptor interactions (distance units) -* angle_cutoff = global angle cutoff for Acceptor-Hydrogen-Donor interactions (degrees) -* equilibrium_angle = global equilibrium angle for Acceptor-Hydrogen-Donor interactions (degrees) - -(Tips: angle offset is the supplementary angle of the equilibrium angle. It means if the -equilibrium angle is 166.6, the angle offset is 13.4) - -Examples -"""""""" - -.. code-block:: LAMMPS - - pair_style hybrid/overlay lj/cut 10.0 hbond/dreiding/lj/angleoffset 4 9.0 11.0 90.0 170.0 - pair_coeff 1 2 hbond/dreiding/lj 3 i 9.5 2.75 4 9.0 11.0 90.0 160.0 - - pair_style hybrid/overlay lj/cut 10.0 hbond/dreiding/morse/angleoffset 2 9.0 11.0 90.0 170.0 - pair_coeff 1 2 hbond/dreiding/morse 3 i 3.88 1.7241379 2.9 2 9.0 11.0 90.0 160.0 - - labelmap atom 1 C 2 O 3 H - pair_coeff C O hbond/dreiding/morse H i 3.88 1.7241379 2.9 2 9.0 11.0 90.0 160.0 - -Description -""""""""""" - -The *hbond/dreiding/\*/angleoffset* styles are modified version of *hbond/dreiding* styles. - -In some cases, the angle of acceptor-hydrogen-donor in the equilibrium state could not achieve 180 degree especially in some coarse grained models. -In these cases, an angle offset is required to ensure that equilibrium state could be the minimum energy state. - -The *hbond/dreiding/\*/angleoffset* styles compute the Acceptor-Hydrogen-Donor (AHD) -3-body hydrogen bond interaction for the :doc:`DREIDING ` -force field with an angle offset, given by: - -.. math:: - - E = & \left[LJ(r) | Morse(r) \right] \qquad \qquad \qquad r < r_{\rm in} and \theta + \theta_{offset} > \theta_c \\ - = & S(r) * \left[LJ(r) | Morse(r) \right] \qquad \qquad r_{\rm in} < r < r_{\rm out} and \theta + \theta_{offset} > \theta_c \\ - = & 0 \qquad \qquad \qquad \qquad \qquad \qquad \qquad r > r_{\rm out} and \theta + \theta_{offset} > \theta_c \\ - LJ(r) = & AR^{-12}-BR^{-10}cos^n(\theta + \theta_{offset})= - \epsilon\left\lbrace 5\left[ \frac{\sigma}{r}\right]^{12}- - 6\left[ \frac{\sigma}{r}\right]^{10} \right\rbrace cos^n(\theta + \theta_{offset})\\ - Morse(r) = & D_0\left\lbrace \chi^2 - 2\chi\right\rbrace cos^n(\theta + \theta_{offset})= - D_{0}\left\lbrace e^{- 2 \alpha (r - r_0)} - 2 e^{- \alpha (r - r_0)} - \right\rbrace cos^n(\theta + \theta_{offset})\\ - S(r) = & \frac{ \left[r_{\rm out}^2 - r^2\right]^2 - \left[r_{\rm out}^2 + 2r^2 - 3{r_{\rm in}^2}\right]} - { \left[r_{\rm out}^2 - {r_{\rm in}}^2\right]^3 } - -where :math:`r_{\rm in}` is the inner spline distance cutoff, -:math:`r_{\rm out}` is the outer distance cutoff, :math:`\theta_c` is -the angle cutoff, :math:`\theta_offset` is the angle offset, and :math:`n` is the power of the cosine of the sum of the angle :math:`\theta` and the angle offset :math:`\theta_offset` - -Here, *r* is the radial distance between the donor (D) and acceptor -(A) atoms and :math:`\theta` is the bond angle between the acceptor, the -hydrogen (H) and the donor atoms: - -.. image:: JPG/dreiding_hbond.jpg - :align: center - -For the *hbond/dreiding/lj/angleoffset* style the list of coefficients is as -follows: - -* K = hydrogen atom type = 1 to Ntypes, or type label -* donor flag = *i* or *j* -* :math:`\epsilon` (energy units) -* :math:`\sigma` (distance units) -* *n* = exponent in formula above -* distance cutoff :math:`r_{\rm in}` (distance units) -* distance cutoff :math:`r_{\rm out}` (distance units) -* angle cutoff (degrees) -* equilibrium angle (degrees) - -(Tips: angle offset is the supplementary angle of equilibrium angle) - -For the *hbond/dreiding/morse/angleoffset* style the list of coefficients is as -follows: - -* K = hydrogen atom type = 1 to Ntypes, or type label -* donor flag = *i* or *j* -* :math:`D_0` (energy units) -* :math:`\alpha` (1/distance units) -* :math:`r_0` (distance units) -* *n* = exponent in formula above -* distance cutoff :math:`r_{\rm in}` (distance units) -* distance cutoff :math:`r_{out}` (distance units) -* angle cutoff (degrees) -* equilibrium angle (degrees) - -(Tips: angle offset is the supplementary angle of equilibrium angle) - ----------- - -Additional Information -"""""""""""""""""""""" - -These pair styles do not support mixing. You must explicitly identify -each donor/acceptor type pair. - -These styles do not support the :doc:`pair_modify ` shift -option for the energy of the interactions. - -For more information about DREIDING force field and other notes, please refer -to the :doc:`documentation of the *hbond/dreiding* pair styles `. - ----------- - -Restrictions -"""""""""""" - -This pair style can only be used if LAMMPS was built with the -EXTRA-MOLECULE package. See the :doc:`Build package ` -doc page for more info. - diff --git a/src/EXTRA-MOLECULE/pair_hbond_dreiding_lj_angleoffset.cpp b/src/EXTRA-MOLECULE/pair_hbond_dreiding_lj_angleoffset.cpp index 68cf24c5ac..829c750536 100644 --- a/src/EXTRA-MOLECULE/pair_hbond_dreiding_lj_angleoffset.cpp +++ b/src/EXTRA-MOLECULE/pair_hbond_dreiding_lj_angleoffset.cpp @@ -33,292 +33,44 @@ #include #include -using namespace LAMMPS_NS; -using namespace MathConst; -using namespace MathSpecial; - static constexpr double SMALL = 0.001; static constexpr int CHUNK = 8; -/* ---------------------------------------------------------------------- */ - -PairHbondDreidingLJangleoffset::PairHbondDreidingLJangleoffset(LAMMPS *lmp) : Pair(lmp) -{ - // hbond cannot compute virial as F dot r - // due to using map() to find bonded H atoms which are not near donor atom - - no_virial_fdotr_compute = 1; - restartinfo = 0; - - nparams = maxparam = 0; - params = nullptr; - - nextra = 2; - pvector = new double[2]; -} +using namespace LAMMPS_NS; +using namespace MathConst; +// using namespace MathSpecial; /* ---------------------------------------------------------------------- */ -PairHbondDreidingLJangleoffset::~PairHbondDreidingLJangleoffset() -{ - memory->sfree(params); - delete [] pvector; +PairHbondDreidingLJAngleoffset::PairHbondDreidingLJAngleoffset(LAMMPS *lmp) + : PairHbondDreidingLJ(lmp) { - if (allocated) { - memory->destroy(setflag); - memory->destroy(cutsq); - - delete [] donor; - delete [] acceptor; - memory->destroy(type2param); - } -} - -/* ---------------------------------------------------------------------- */ - -void PairHbondDreidingLJangleoffset::compute(int eflag, int vflag) -{ - int i,j,k,m,ii,jj,kk,inum,jnum,knum,itype,jtype,ktype,iatom,imol; - tagint tagprev; - double delx,dely,delz,rsq,rsq1,rsq2,r1,r2; - double factor_hb,force_angle,force_kernel,evdwl,eng_lj,ehbond,force_switch; - double c,s,a,b,ac,a11,a12,a22,vx1,vx2,vy1,vy2,vz1,vz2,d; - double fi[3],fj[3],delr1[3],delr2[3]; - double r2inv,r10inv; - double switch1,switch2; - int *ilist,*jlist,*numneigh,**firstneigh; - tagint *klist; - - evdwl = ehbond = 0.0; - ev_init(eflag,vflag); - - double **x = atom->x; - double **f = atom->f; - tagint *tag = atom->tag; - int *molindex = atom->molindex; - int *molatom = atom->molatom; - tagint **special = atom->special; - int **nspecial = atom->nspecial; - int *type = atom->type; - double *special_lj = force->special_lj; - int molecular = atom->molecular; - Molecule **onemols = atom->avec->onemols; - - inum = list->inum; - ilist = list->ilist; - numneigh = list->numneigh; - firstneigh = list->firstneigh; - - // ii = loop over donors - // jj = loop over acceptors - // kk = loop over hydrogens bonded to donor - - int hbcount = 0; - - for (ii = 0; ii < inum; ii++) { - i = ilist[ii]; - itype = type[i]; - if (!donor[itype]) continue; - if (molecular == Atom::MOLECULAR) { - klist = special[i]; - knum = nspecial[i][0]; - } else { - if (molindex[i] < 0) continue; - imol = molindex[i]; - iatom = molatom[i]; - klist = onemols[imol]->special[iatom]; - knum = onemols[imol]->nspecial[iatom][0]; - tagprev = tag[i] - iatom - 1; - } - jlist = firstneigh[i]; - jnum = numneigh[i]; - - for (jj = 0; jj < jnum; jj++) { - j = jlist[jj]; - factor_hb = special_lj[sbmask(j)]; - j &= NEIGHMASK; - - jtype = type[j]; - if (!acceptor[jtype]) continue; - - delx = x[i][0] - x[j][0]; - dely = x[i][1] - x[j][1]; - delz = x[i][2] - x[j][2]; - rsq = delx*delx + dely*dely + delz*delz; - - for (kk = 0; kk < knum; kk++) { - if (molecular == Atom::MOLECULAR) k = atom->map(klist[kk]); - else k = atom->map(klist[kk]+tagprev); - if (k < 0) continue; - ktype = type[k]; - m = type2param[itype][jtype][ktype]; - if (m < 0) continue; - const Param &pm = params[m]; - - if (rsq < pm.cut_outersq) { - delr1[0] = x[i][0] - x[k][0]; - delr1[1] = x[i][1] - x[k][1]; - delr1[2] = x[i][2] - x[k][2]; - domain->minimum_image(delr1); - rsq1 = delr1[0]*delr1[0] + delr1[1]*delr1[1] + delr1[2]*delr1[2]; - r1 = sqrt(rsq1); - - delr2[0] = x[j][0] - x[k][0]; - delr2[1] = x[j][1] - x[k][1]; - delr2[2] = x[j][2] - x[k][2]; - domain->minimum_image(delr2); - rsq2 = delr2[0]*delr2[0] + delr2[1]*delr2[1] + delr2[2]*delr2[2]; - r2 = sqrt(rsq2); - - // angle (cos and sin) - - c = delr1[0]*delr2[0] + delr1[1]*delr2[1] + delr1[2]*delr2[2]; - c /= r1*r2; - if (c > 1.0) c = 1.0; - if (c < -1.0) c = -1.0; - ac = acos(c); - - ac = ac + pm.angle_offset; - c = cos(ac); - if (c > 1.0) c = 1.0; - if (c < -1.0) c = -1.0; - - if (ac > pm.cut_angle && ac < (2.0*MY_PI - pm.cut_angle)) { - s = sqrt(1.0 - c*c); - if (s < SMALL) s = SMALL; - - // LJ-specific kernel - - r2inv = 1.0/rsq; - r10inv = r2inv*r2inv*r2inv*r2inv*r2inv; - force_kernel = r10inv*(pm.lj1*r2inv - pm.lj2)*r2inv * - powint(c,pm.ap); - force_angle = pm.ap * r10inv*(pm.lj3*r2inv - pm.lj4) * - powint(c,pm.ap-1)*s; - - eng_lj = r10inv*(pm.lj3*r2inv - pm.lj4); - - force_switch=0.0; - - if (rsq > pm.cut_innersq) { - switch1 = (pm.cut_outersq-rsq) * (pm.cut_outersq-rsq) * - (pm.cut_outersq + 2.0*rsq - 3.0*pm.cut_innersq) / - pm.denom_vdw; - switch2 = 12.0*rsq * (pm.cut_outersq-rsq) * - (rsq-pm.cut_innersq) / pm.denom_vdw; - - force_kernel *= switch1; - force_angle *= switch1; - force_switch = eng_lj*switch2/rsq; - eng_lj *= switch1; - } - - if (eflag) { - evdwl = eng_lj * powint(c,pm.ap); - evdwl *= factor_hb; - ehbond += evdwl; - } - - a = factor_hb*force_angle/s; - b = factor_hb*force_kernel; - d = factor_hb*force_switch; - - a11 = a*c / rsq1; - a12 = -a / (r1*r2); - a22 = a*c / rsq2; - - vx1 = a11*delr1[0] + a12*delr2[0]; - vx2 = a22*delr2[0] + a12*delr1[0]; - vy1 = a11*delr1[1] + a12*delr2[1]; - vy2 = a22*delr2[1] + a12*delr1[1]; - vz1 = a11*delr1[2] + a12*delr2[2]; - vz2 = a22*delr2[2] + a12*delr1[2]; - - fi[0] = vx1 + b*delx + d*delx; - fi[1] = vy1 + b*dely + d*dely; - fi[2] = vz1 + b*delz + d*delz; - fj[0] = vx2 - b*delx - d*delx; - fj[1] = vy2 - b*dely - d*dely; - fj[2] = vz2 - b*delz - d*delz; - - f[i][0] += fi[0]; - f[i][1] += fi[1]; - f[i][2] += fi[2]; - - f[j][0] += fj[0]; - f[j][1] += fj[1]; - f[j][2] += fj[2]; - - f[k][0] -= vx1 + vx2; - f[k][1] -= vy1 + vy2; - f[k][2] -= vz1 + vz2; - - // KIJ instead of IJK b/c delr1/delr2 are both with respect to k - - if (evflag) ev_tally3(k,i,j,evdwl,0.0,fi,fj,delr1,delr2); - - hbcount++; - } - } - } - } - } - - if (eflag_global) { - pvector[0] = hbcount; - pvector[1] = ehbond; - } -} - -/* ---------------------------------------------------------------------- - allocate all arrays -------------------------------------------------------------------------- */ - -void PairHbondDreidingLJangleoffset::allocate() -{ - allocated = 1; - int n = atom->ntypes; - - // mark all setflag as set, since don't require pair_coeff of all I,J - - memory->create(setflag,n+1,n+1,"pair:setflag"); - for (int i = 1; i <= n; i++) - for (int j = i; j <= n; j++) - setflag[i][j] = 1; - - memory->create(cutsq,n+1,n+1,"pair:cutsq"); - - donor = new int[n+1]; - acceptor = new int[n+1]; - memory->create(type2param,n+1,n+1,n+1,"pair:type2param"); - - int i,j,k; - for (i = 1; i <= n; i++) - for (j = 1; j <= n; j++) - for (k = 1; k <= n; k++) - type2param[i][j][k] = -1; + angle_offset_flag = 1; + angle_offset_global = 0.0; + angle_offset_one = 0.0; } /* ---------------------------------------------------------------------- global settings ------------------------------------------------------------------------- */ -void PairHbondDreidingLJangleoffset::settings(int narg, char **arg) +void PairHbondDreidingLJAngleoffset::settings(int narg, char **arg) { - if (narg != 5) error->all(FLERR, "Illegal pair_style command"); + if (narg != 5) error->all(FLERR,"Illegal pair_style command"); + + // use first four settings args in parent method to update other variables + + PairHbondDreidingLJ::settings(4,arg); - ap_global = utils::inumeric(FLERR, arg[0], false, lmp); - cut_inner_global = utils::numeric(FLERR, arg[1], false, lmp); - cut_outer_global = utils::numeric(FLERR, arg[2], false, lmp); - cut_angle_global = utils::numeric(FLERR, arg[3], false, lmp) * MY_PI/180.0; angle_offset_global = (180.0 - utils::numeric(FLERR, arg[4], false, lmp)) * MY_PI/180.0; + } /* ---------------------------------------------------------------------- set coeffs for one or more type pairs ------------------------------------------------------------------------- */ -void PairHbondDreidingLJangleoffset::coeff(int narg, char **arg) +void PairHbondDreidingLJAngleoffset::coeff(int narg, char **arg) { if (narg < 6 || narg > 11) error->all(FLERR,"Incorrect args for pair coefficients"); @@ -395,187 +147,3 @@ void PairHbondDreidingLJangleoffset::coeff(int narg, char **arg) if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients"); } - -/* ---------------------------------------------------------------------- - init specific to this pair style -------------------------------------------------------------------------- */ - -void PairHbondDreidingLJangleoffset::init_style() -{ - // molecular system required to use special list to find H atoms - // tags required to use special list - // pair newton on required since are looping over D atoms - // and computing forces on A,H which may be on different procs - - if (atom->molecular == Atom::ATOMIC) - error->all(FLERR,"Pair style hbond/dreiding/lj/angleoffset requires molecular system"); - if (atom->tag_enable == 0) - error->all(FLERR,"Pair style hbond/dreiding/lj/angleoffset requires atom IDs"); - if (atom->map_style == Atom::MAP_NONE) - error->all(FLERR,"Pair style hbond/dreiding/lj/angleoffset requires an atom map, " - "see atom_modify"); - if (force->newton_pair == 0) - error->all(FLERR,"Pair style hbond/dreiding/lj/angleoffset requires newton pair on"); - - // set donor[M]/acceptor[M] if any atom of type M is a donor/acceptor - - int anyflag = 0; - int n = atom->ntypes; - for (int m = 1; m <= n; m++) donor[m] = acceptor[m] = 0; - for (int i = 1; i <= n; i++) - for (int j = 1; j <= n; j++) - for (int k = 1; k <= n; k++) - if (type2param[i][j][k] >= 0) { - anyflag = 1; - donor[i] = 1; - acceptor[j] = 1; - } - - if (!anyflag) error->all(FLERR,"No pair hbond/dreiding/lj/angleoffset coefficients set"); - - // set additional param values - // offset is for LJ only, angle term is not included - - for (int m = 0; m < nparams; m++) { - params[m].lj1 = 60.0*params[m].epsilon*pow(params[m].sigma,12.0); - params[m].lj2 = 60.0*params[m].epsilon*pow(params[m].sigma,10.0); - params[m].lj3 = 5.0*params[m].epsilon*pow(params[m].sigma,12.0); - params[m].lj4 = 6.0*params[m].epsilon*pow(params[m].sigma,10.0); - } - - // full neighbor list request - - neighbor->add_request(this, NeighConst::REQ_FULL); -} - -/* ---------------------------------------------------------------------- - init for one type pair i,j and corresponding j,i -------------------------------------------------------------------------- */ - -double PairHbondDreidingLJangleoffset::init_one(int i, int j) -{ - int m; - - // return maximum cutoff for any K with I,J = D,A or J,I = D,A - // donor/acceptor is not symmetric, IJ interaction != JI interaction - - double cut = 0.0; - for (int k = 1; k <= atom->ntypes; k++) { - m = type2param[i][j][k]; - if (m >= 0) cut = MAX(cut,params[m].cut_outer); - m = type2param[j][i][k]; - if (m >= 0) cut = MAX(cut,params[m].cut_outer); - } - return cut; -} - -/* ---------------------------------------------------------------------- */ - -double PairHbondDreidingLJangleoffset::single(int i, int j, int itype, int jtype, - double rsq, - double /*factor_coul*/, double /*factor_lj*/, - double &fforce) -{ - int k,kk,ktype,knum,m; - tagint tagprev; - double eng,eng_lj,force_kernel,force_angle; - double rsq1,rsq2,r1,r2,c,s,ac,r2inv,r10inv,factor_hb; - double switch1,switch2; - double delr1[3],delr2[3]; - tagint *klist; - - double **x = atom->x; - int *type = atom->type; - double *special_lj = force->special_lj; - - eng = 0.0; - fforce = 0; - - // sanity check - - if (!donor[itype]) return 0.0; - if (!acceptor[jtype]) return 0.0; - - int molecular = atom->molecular; - if (molecular == Atom::MOLECULAR) { - klist = atom->special[i]; - knum = atom->nspecial[i][0]; - } else { - if (atom->molindex[i] < 0) return 0.0; - int imol = atom->molindex[i]; - int iatom = atom->molatom[i]; - Molecule **onemols = atom->avec->onemols; - klist = onemols[imol]->special[iatom]; - knum = onemols[imol]->nspecial[iatom][0]; - tagprev = atom->tag[i] - iatom - 1; - } - - factor_hb = special_lj[sbmask(j)]; - - for (kk = 0; kk < knum; kk++) { - if (molecular == Atom::MOLECULAR) k = atom->map(klist[kk]); - else k = atom->map(klist[kk]+tagprev); - - if (k < 0) continue; - ktype = type[k]; - m = type2param[itype][jtype][ktype]; - if (m < 0) continue; - const Param &pm = params[m]; - - delr1[0] = x[i][0] - x[k][0]; - delr1[1] = x[i][1] - x[k][1]; - delr1[2] = x[i][2] - x[k][2]; - domain->minimum_image(delr1); - rsq1 = delr1[0]*delr1[0] + delr1[1]*delr1[1] + delr1[2]*delr1[2]; - r1 = sqrt(rsq1); - - delr2[0] = x[j][0] - x[k][0]; - delr2[1] = x[j][1] - x[k][1]; - delr2[2] = x[j][2] - x[k][2]; - domain->minimum_image(delr2); - rsq2 = delr2[0]*delr2[0] + delr2[1]*delr2[1] + delr2[2]*delr2[2]; - r2 = sqrt(rsq2); - - // angle (cos and sin) - - c = delr1[0]*delr2[0] + delr1[1]*delr2[1] + delr1[2]*delr2[2]; - c /= r1*r2; - if (c > 1.0) c = 1.0; - if (c < -1.0) c = -1.0; - ac = acos(c); - - ac = ac + pm.angle_offset; - c = cos(ac); - if (c > 1.0) c = 1.0; - if (c < -1.0) c = -1.0; - - if (ac < pm.cut_angle || ac > (2.0*MY_PI - pm.cut_angle)) return 0.0; - s = sqrt(1.0 - c*c); - if (s < SMALL) s = SMALL; - - // LJ-specific kernel - - r2inv = 1.0/rsq; - r10inv = r2inv*r2inv*r2inv*r2inv*r2inv; - force_kernel = r10inv*(pm.lj1*r2inv - pm.lj2)*r2inv * powint(c,pm.ap); - force_angle = pm.ap * r10inv*(pm.lj3*r2inv - pm.lj4) * - powint(c,pm.ap-1)*s; - - // only lj part for now - - eng_lj = r10inv*(pm.lj3*r2inv - pm.lj4); - if (rsq > pm.cut_innersq) { - switch1 = (pm.cut_outersq-rsq) * (pm.cut_outersq-rsq) * - (pm.cut_outersq + 2.0*rsq - 3.0*pm.cut_innersq) / pm.denom_vdw; - switch2 = 12.0*rsq * (pm.cut_outersq-rsq) * - (rsq-pm.cut_innersq) / pm.denom_vdw; - force_kernel = force_kernel*switch1 + eng_lj*switch2; - eng_lj *= switch1; - } - - fforce += force_kernel*powint(c,pm.ap) + eng_lj*force_angle; - eng += eng_lj * powint(c,pm.ap) * factor_hb; - } - - return eng; -} diff --git a/src/EXTRA-MOLECULE/pair_hbond_dreiding_lj_angleoffset.h b/src/EXTRA-MOLECULE/pair_hbond_dreiding_lj_angleoffset.h index 331d8e2ac1..590c5198b7 100644 --- a/src/EXTRA-MOLECULE/pair_hbond_dreiding_lj_angleoffset.h +++ b/src/EXTRA-MOLECULE/pair_hbond_dreiding_lj_angleoffset.h @@ -13,53 +13,28 @@ #ifdef PAIR_CLASS // clang-format off -PairStyle(hbond/dreiding/lj/angleoffset,PairHbondDreidingLJangleoffset); +PairStyle(hbond/dreiding/lj/angleoffset,PairHbondDreidingLJAngleoffset); // clang-format on #else #ifndef LMP_PAIR_HBOND_DREIDING_LJ_ANGLEOFFSET_H #define LMP_PAIR_HBOND_DREIDING_LJ_ANGLEOFFSET_H -#include "pair.h" +#include "pair_hbond_dreiding_lj.h" namespace LAMMPS_NS { -class PairHbondDreidingLJangleoffset : public Pair { - public: - PairHbondDreidingLJangleoffset(class LAMMPS *); - ~PairHbondDreidingLJangleoffset() override; - void compute(int, int) override; - void settings(int, char **) override; - void coeff(int, char **) override; - void init_style() override; - double init_one(int, int) override; - double single(int, int, int, int, double, double, double, double &) override; + class PairHbondDreidingLJAngleoffset : public PairHbondDreidingLJ { - protected: - double cut_inner_global, cut_outer_global, cut_angle_global, angle_offset_global; - int ap_global; + public: + PairHbondDreidingLJAngleoffset(class LAMMPS *); + void settings(int, char **) override; + void coeff(int, char **) override; - struct Param { - double epsilon, sigma; - double lj1, lj2, lj3, lj4; - double d0, alpha, r0; - double morse1; - double denom_vdw; - double cut_inner, cut_outer, cut_innersq, cut_outersq, cut_angle, offset, angle_offset; - int ap; + protected: + double angle_offset_global, angle_offset_one, cut_angle_one; }; - Param *params; // parameter set for an I-J-K interaction - int nparams; // number of parameters read - int maxparam; - - int *donor; // 1 if this type is ever a donor, else 0 - int *acceptor; // 1 if this type is ever an acceptor, else 0 - int ***type2param; // mapping from D,A,H to params, -1 if no map - - void allocate(); -}; - } // namespace LAMMPS_NS #endif diff --git a/src/EXTRA-MOLECULE/pair_hbond_dreiding_morse_angleoffset.cpp b/src/EXTRA-MOLECULE/pair_hbond_dreiding_morse_angleoffset.cpp index 2a1b6f1142..227b84cc4d 100644 --- a/src/EXTRA-MOLECULE/pair_hbond_dreiding_morse_angleoffset.cpp +++ b/src/EXTRA-MOLECULE/pair_hbond_dreiding_morse_angleoffset.cpp @@ -35,212 +35,40 @@ using namespace LAMMPS_NS; using namespace MathConst; -using namespace MathSpecial; -static constexpr double SMALL = 0.001; static constexpr int CHUNK = 8; /* ---------------------------------------------------------------------- */ PairHbondDreidingMorseAngleoffset::PairHbondDreidingMorseAngleoffset(LAMMPS *lmp) : - PairHbondDreidingLJangleoffset(lmp) {} + PairHbondDreidingLJ(lmp) { -/* ---------------------------------------------------------------------- */ - -void PairHbondDreidingMorseAngleoffset::compute(int eflag, int vflag) -{ - int i,j,k,m,ii,jj,kk,inum,jnum,knum,itype,jtype,ktype,imol,iatom; - tagint tagprev; - double delx,dely,delz,rsq,rsq1,rsq2,r1,r2; - double factor_hb,force_angle,force_kernel,force_switch,evdwl,ehbond; - double c,s,a,b,d,ac,a11,a12,a22,vx1,vx2,vy1,vy2,vz1,vz2; - double fi[3],fj[3],delr1[3],delr2[3]; - double r,dr,dexp,eng_morse,switch1,switch2; - int *ilist,*jlist,*numneigh,**firstneigh; - tagint *klist; - - evdwl = ehbond = 0.0; - ev_init(eflag,vflag); - - double **x = atom->x; - double **f = atom->f; - tagint *tag = atom->tag; - int *molindex = atom->molindex; - int *molatom = atom->molatom; - tagint **special = atom->special; - int **nspecial = atom->nspecial; - int *type = atom->type; - double *special_lj = force->special_lj; - int molecular = atom->molecular; - Molecule **onemols = atom->avec->onemols; - - inum = list->inum; - ilist = list->ilist; - numneigh = list->numneigh; - firstneigh = list->firstneigh; - - // ii = loop over donors - // jj = loop over acceptors - // kk = loop over hydrogens bonded to donor - - int hbcount = 0; - - for (ii = 0; ii < inum; ii++) { - i = ilist[ii]; - itype = type[i]; - if (!donor[itype]) continue; - if (molecular == Atom::MOLECULAR) { - klist = special[i]; - knum = nspecial[i][0]; - } else { - if (molindex[i] < 0) continue; - imol = molindex[i]; - iatom = molatom[i]; - klist = onemols[imol]->special[iatom]; - knum = onemols[imol]->nspecial[iatom][0]; - tagprev = tag[i] - iatom - 1; - } - jlist = firstneigh[i]; - jnum = numneigh[i]; - - for (jj = 0; jj < jnum; jj++) { - j = jlist[jj]; - factor_hb = special_lj[sbmask(j)]; - j &= NEIGHMASK; - - jtype = type[j]; - if (!acceptor[jtype]) continue; - - delx = x[i][0] - x[j][0]; - dely = x[i][1] - x[j][1]; - delz = x[i][2] - x[j][2]; - rsq = delx*delx + dely*dely + delz*delz; - - for (kk = 0; kk < knum; kk++) { - if (molecular == Atom::MOLECULAR) k = atom->map(klist[kk]); - else k = atom->map(klist[kk]+tagprev); - if (k < 0) continue; - ktype = type[k]; - m = type2param[itype][jtype][ktype]; - if (m < 0) continue; - const Param &pm = params[m]; - - if (rsq < pm.cut_outersq) { - delr1[0] = x[i][0] - x[k][0]; - delr1[1] = x[i][1] - x[k][1]; - delr1[2] = x[i][2] - x[k][2]; - domain->minimum_image(delr1); - rsq1 = delr1[0]*delr1[0] + delr1[1]*delr1[1] + delr1[2]*delr1[2]; - r1 = sqrt(rsq1); - - delr2[0] = x[j][0] - x[k][0]; - delr2[1] = x[j][1] - x[k][1]; - delr2[2] = x[j][2] - x[k][2]; - domain->minimum_image(delr2); - rsq2 = delr2[0]*delr2[0] + delr2[1]*delr2[1] + delr2[2]*delr2[2]; - r2 = sqrt(rsq2); - - // angle (cos and sin) - - c = delr1[0]*delr2[0] + delr1[1]*delr2[1] + delr1[2]*delr2[2]; - c /= r1*r2; - if (c > 1.0) c = 1.0; - if (c < -1.0) c = -1.0; - ac = acos(c); - - ac = ac + pm.angle_offset; - c = cos(ac); - if (c > 1.0) c = 1.0; - if (c < -1.0) c = -1.0; - - if (ac > pm.cut_angle && ac < (2.0*MY_PI - pm.cut_angle)) { - s = sqrt(1.0 - c*c); - if (s < SMALL) s = SMALL; - - // Morse-specific kernel - - r = sqrt(rsq); - dr = r - pm.r0; - dexp = exp(-pm.alpha * dr); - eng_morse = pm.d0 * (dexp*dexp - 2.0*dexp); - force_kernel = pm.morse1*(dexp*dexp - dexp)/r * powint(c,pm.ap); - force_angle = pm.ap * eng_morse * powint(c,pm.ap-1)*s; - force_switch = 0.0; - - if (rsq > pm.cut_innersq) { - switch1 = (pm.cut_outersq-rsq) * (pm.cut_outersq-rsq) * - (pm.cut_outersq + 2.0*rsq - 3.0*pm.cut_innersq) / - pm.denom_vdw; - switch2 = 12.0*rsq * (pm.cut_outersq-rsq) * - (rsq-pm.cut_innersq) / pm.denom_vdw; - - force_kernel *= switch1; - force_angle *= switch1; - force_switch = eng_morse*switch2/rsq; - eng_morse *= switch1; - } - - if (eflag) { - evdwl = eng_morse * powint(c,pm.ap); - evdwl *= factor_hb; - ehbond += evdwl; - } - - a = factor_hb*force_angle/s; - b = factor_hb*force_kernel; - d = factor_hb*force_switch; - - a11 = a*c / rsq1; - a12 = -a / (r1*r2); - a22 = a*c / rsq2; - - vx1 = a11*delr1[0] + a12*delr2[0]; - vx2 = a22*delr2[0] + a12*delr1[0]; - vy1 = a11*delr1[1] + a12*delr2[1]; - vy2 = a22*delr2[1] + a12*delr1[1]; - vz1 = a11*delr1[2] + a12*delr2[2]; - vz2 = a22*delr2[2] + a12*delr1[2]; - - fi[0] = vx1 + (b+d)*delx; - fi[1] = vy1 + (b+d)*dely; - fi[2] = vz1 + (b+d)*delz; - fj[0] = vx2 - (b+d)*delx; - fj[1] = vy2 - (b+d)*dely; - fj[2] = vz2 - (b+d)*delz; - - f[i][0] += fi[0]; - f[i][1] += fi[1]; - f[i][2] += fi[2]; - - f[j][0] += fj[0]; - f[j][1] += fj[1]; - f[j][2] += fj[2]; - - f[k][0] -= vx1 + vx2; - f[k][1] -= vy1 + vy2; - f[k][2] -= vz1 + vz2; - - // KIJ instead of IJK b/c delr1/delr2 are both with respect to k - - if (evflag) ev_tally3(k,i,j,evdwl,0.0,fi,fj,delr1,delr2); - - hbcount++; - } - } - } - } - } - - if (eflag_global) { - pvector[0] = hbcount; - pvector[1] = ehbond; - } + angle_offset_flag = 1; + angle_offset_global = 0.0; + angle_offset_one = 0.0; } + /* ---------------------------------------------------------------------- - set coeffs for one or more type pairs + global settings ------------------------------------------------------------------------- */ +void PairHbondDreidingMorseAngleoffset::settings(int narg, char **arg) +{ + if (narg != 5) error->all(FLERR,"Illegal pair_style command"); + + // use first four settings args in parent method to update other variables + + PairHbondDreidingLJ::settings(4,arg); + + angle_offset_global = (180.0 - utils::numeric(FLERR, arg[4], false, lmp)) * MY_PI/180.0; + +} + +// /* ---------------------------------------------------------------------- +// set coeffs for one or more type pairs +// ------------------------------------------------------------------------- */ + void PairHbondDreidingMorseAngleoffset::coeff(int narg, char **arg) { if (narg < 7 || narg > 12) @@ -318,162 +146,3 @@ void PairHbondDreidingMorseAngleoffset::coeff(int narg, char **arg) if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients"); } - -/* ---------------------------------------------------------------------- - init specific to this pair style -------------------------------------------------------------------------- */ - -void PairHbondDreidingMorseAngleoffset::init_style() -{ - // molecular system required to use special list to find H atoms - // tags required to use special list - // pair newton on required since are looping over D atoms - // and computing forces on A,H which may be on different procs - - if (atom->molecular == Atom::ATOMIC) - error->all(FLERR,"Pair style hbond/dreiding/morse/angleoffset requires molecular system"); - if (atom->tag_enable == 0) - error->all(FLERR,"Pair style hbond/dreiding/morse/angleoffset requires atom IDs"); - if (atom->map_style == Atom::MAP_NONE) - error->all(FLERR,"Pair style hbond/dreiding/morse/angleoffset requires an atom map, " - "see atom_modify"); - if (force->newton_pair == 0) - error->all(FLERR,"Pair style hbond/dreiding/morse/angleoffset requires newton pair on"); - - // set donor[M]/acceptor[M] if any atom of type M is a donor/acceptor - - int anyflag = 0; - int n = atom->ntypes; - for (int m = 1; m <= n; m++) donor[m] = acceptor[m] = 0; - for (int i = 1; i <= n; i++) - for (int j = 1; j <= n; j++) - for (int k = 1; k <= n; k++) - if (type2param[i][j][k] >= 0) { - anyflag = 1; - donor[i] = 1; - acceptor[j] = 1; - } - - if (!anyflag) error->all(FLERR,"No pair hbond/dreiding/morse/angleoffset coefficients set"); - - // set additional param values - // offset is for Morse only, angle term is not included - - for (int m = 0; m < nparams; m++) { - params[m].morse1 = 2.0*params[m].d0*params[m].alpha; - } - - // full neighbor list request - - neighbor->add_request(this, NeighConst::REQ_FULL); -} - -/* ---------------------------------------------------------------------- */ - -double PairHbondDreidingMorseAngleoffset::single(int i, int j, int itype, int jtype, - double rsq, - double /*factor_coul*/, double /*factor_lj*/, - double &fforce) -{ - int k,kk,ktype,knum,m; - tagint tagprev; - double eng,eng_morse,force_kernel,force_angle; - double rsq1,rsq2,r1,r2,c,s,ac,r,dr,dexp,factor_hb; - double switch1,switch2; - double delr1[3],delr2[3]; - tagint *klist; - - double **x = atom->x; - int *type = atom->type; - double *special_lj = force->special_lj; - - eng = 0.0; - fforce = 0; - - //sanity check - - if (!donor[itype]) return 0.0; - if (!acceptor[jtype]) return 0.0; - - int molecular = atom->molecular; - if (molecular == Atom::MOLECULAR) { - klist = atom->special[i]; - knum = atom->nspecial[i][0]; - } else { - if (atom->molindex[i] < 0) return 0.0; - int imol = atom->molindex[i]; - int iatom = atom->molatom[i]; - Molecule **onemols = atom->avec->onemols; - klist = onemols[imol]->special[iatom]; - knum = onemols[imol]->nspecial[iatom][0]; - tagprev = atom->tag[i] - iatom - 1; - } - - factor_hb = special_lj[sbmask(j)]; - - for (kk = 0; kk < knum; kk++) { - if (molecular == Atom::MOLECULAR) k = atom->map(klist[kk]); - else k = atom->map(klist[kk]+tagprev); - - if (k < 0) continue; - ktype = type[k]; - m = type2param[itype][jtype][ktype]; - if (m < 0) continue; - const Param &pm = params[m]; - - delr1[0] = x[i][0] - x[k][0]; - delr1[1] = x[i][1] - x[k][1]; - delr1[2] = x[i][2] - x[k][2]; - domain->minimum_image(delr1); - rsq1 = delr1[0]*delr1[0] + delr1[1]*delr1[1] + delr1[2]*delr1[2]; - r1 = sqrt(rsq1); - - delr2[0] = x[j][0] - x[k][0]; - delr2[1] = x[j][1] - x[k][1]; - delr2[2] = x[j][2] - x[k][2]; - domain->minimum_image(delr2); - rsq2 = delr2[0]*delr2[0] + delr2[1]*delr2[1] + delr2[2]*delr2[2]; - r2 = sqrt(rsq2); - - // angle (cos and sin) - - c = delr1[0]*delr2[0] + delr1[1]*delr2[1] + delr1[2]*delr2[2]; - c /= r1*r2; - if (c > 1.0) c = 1.0; - if (c < -1.0) c = -1.0; - ac = acos(c); - - ac = ac + pm.angle_offset; - c = cos(ac); - if (c > 1.0) c = 1.0; - if (c < -1.0) c = -1.0; - - if (ac < pm.cut_angle || ac > (2.0*MY_PI - pm.cut_angle)) return 0.0; - s = sqrt(1.0 - c*c); - if (s < SMALL) s = SMALL; - - // Morse-specific kernel - - r = sqrt(rsq); - dr = r - pm.r0; - dexp = exp(-pm.alpha * dr); - eng_morse = pm.d0 * (dexp*dexp - 2.0*dexp); //<-- BUGFIX 2012-11-14 - force_kernel = pm.morse1*(dexp*dexp - dexp)/r * powint(c,pm.ap); - force_angle = pm.ap * eng_morse * powint(c,pm.ap-1)*s; - - if (rsq > pm.cut_innersq) { - switch1 = (pm.cut_outersq-rsq) * (pm.cut_outersq-rsq) * - (pm.cut_outersq + 2.0*rsq - 3.0*pm.cut_innersq) / - pm.denom_vdw; - switch2 = 12.0*rsq * (pm.cut_outersq-rsq) * - (rsq-pm.cut_innersq) / pm.denom_vdw; - force_kernel = force_kernel*switch1 + eng_morse*switch2; - eng_morse *= switch1; - } - - eng += eng_morse * powint(c,pm.ap)* factor_hb; - fforce += force_kernel*powint(c,pm.ap) + eng_morse*force_angle; - } - - return eng; -} diff --git a/src/EXTRA-MOLECULE/pair_hbond_dreiding_morse_angleoffset.h b/src/EXTRA-MOLECULE/pair_hbond_dreiding_morse_angleoffset.h index c1b0cf2f30..2622aed852 100644 --- a/src/EXTRA-MOLECULE/pair_hbond_dreiding_morse_angleoffset.h +++ b/src/EXTRA-MOLECULE/pair_hbond_dreiding_morse_angleoffset.h @@ -20,18 +20,20 @@ PairStyle(hbond/dreiding/morse/angleoffset,PairHbondDreidingMorseAngleoffset); #ifndef LMP_PAIR_HBOND_DREIDING_MORSE_ANGLEOFFSET_H #define LMP_PAIR_HBOND_DREIDING_MORSE_ANGLEOFFSET_H -#include "pair_hbond_dreiding_lj_angleoffset.h" +#include "pair_hbond_dreiding_lj.h" +// #include "pair_hbond_dreiding_morse.h" namespace LAMMPS_NS { -class PairHbondDreidingMorseAngleoffset : public PairHbondDreidingLJangleoffset { +class PairHbondDreidingMorseAngleoffset : public PairHbondDreidingLJ { + public: PairHbondDreidingMorseAngleoffset(class LAMMPS *); - - void compute(int, int) override; + void settings(int, char **) override; void coeff(int, char **) override; - void init_style() override; - double single(int, int, int, int, double, double, double, double &) override; + + protected: + double angle_offset_global, angle_offset_one, cut_angle_one; }; } // namespace LAMMPS_NS diff --git a/src/MOLECULE/pair_hbond_dreiding_lj.cpp b/src/MOLECULE/pair_hbond_dreiding_lj.cpp index 274f8bc2a3..491e17230a 100644 --- a/src/MOLECULE/pair_hbond_dreiding_lj.cpp +++ b/src/MOLECULE/pair_hbond_dreiding_lj.cpp @@ -13,13 +13,14 @@ ------------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- - Contributing author: Tod A Pascal (Caltech) + Contributing authors: Tod A Pascal (Caltech), Don Xu/EiPi Fun ------------------------------------------------------------------------- */ #include "pair_hbond_dreiding_lj.h" #include "atom.h" #include "atom_vec.h" +#include "comm.h" #include "domain.h" #include "error.h" #include "force.h" @@ -55,6 +56,8 @@ PairHbondDreidingLJ::PairHbondDreidingLJ(LAMMPS *lmp) : Pair(lmp) nextra = 2; pvector = new double[2]; + + angle_offset_flag = 0; } /* ---------------------------------------------------------------------- */ @@ -178,6 +181,13 @@ void PairHbondDreidingLJ::compute(int eflag, int vflag) if (c < -1.0) c = -1.0; ac = acos(c); + if (angle_offset_flag){ + ac = ac + pm.angle_offset; + c = cos(ac); + if (c > 1.0) c = 1.0; + if (c < -1.0) c = -1.0; + } + if (ac > pm.cut_angle && ac < (2.0*MY_PI - pm.cut_angle)) { s = sqrt(1.0 - c*c); if (s < SMALL) s = SMALL; @@ -314,8 +324,14 @@ void PairHbondDreidingLJ::settings(int narg, char **arg) void PairHbondDreidingLJ::coeff(int narg, char **arg) { - if (narg < 6 || narg > 10) + // account for angleoffset variant in EXTRA-MOLECULE + int maxarg = 10; + if (angle_offset_flag == 1) maxarg = 11; + + // check settings + if (narg < 6 || narg > maxarg) error->all(FLERR,"Incorrect args for pair coefficients"); + if (!allocated) allocate(); int ilo,ihi,jlo,jhi,klo,khi; @@ -343,8 +359,14 @@ void PairHbondDreidingLJ::coeff(int narg, char **arg) error->all(FLERR,"Pair inner cutoff >= Pair outer cutoff"); double cut_angle_one = cut_angle_global; if (narg == 10) cut_angle_one = utils::numeric(FLERR, arg[9], false, lmp) * MY_PI/180.0; + + if (comm->me == 0) + utils::logmesg(lmp,"---> DEBUG BASECLASS cut_angle_one {} \n", cut_angle_one); + // grow params array if necessary + if (comm->me == 0) + utils::logmesg(lmp,"---> DEBUG BASECLASS coeff nparams {} maxparam {}\n", nparams, maxparam); if (nparams == maxparam) { maxparam += CHUNK; params = (Param *) memory->srealloc(params, maxparam*sizeof(Param), @@ -352,10 +374,12 @@ void PairHbondDreidingLJ::coeff(int narg, char **arg) // make certain all addional allocated storage is initialized // to avoid false positives when checking with valgrind - memset(params + nparams, 0, CHUNK*sizeof(Param)); } + if (comm->me == 0) + utils::logmesg(lmp,"---> DEBUG BASECLASS coeff POST MEM nparams {} maxparam {}\n", nparams, maxparam); + params[nparams].epsilon = epsilon_one; params[nparams].sigma = sigma_one; params[nparams].ap = ap_one; @@ -369,6 +393,9 @@ void PairHbondDreidingLJ::coeff(int narg, char **arg) (params[nparams].cut_outersq-params[nparams].cut_innersq) * (params[nparams].cut_outersq-params[nparams].cut_innersq); + if (comm->me == 0) + utils::logmesg(lmp,"BASECLASS PARAMS angle_offset_flag {}\n eps {}\n sigma {}\n ap {}\n cin {}\n cout {}\n cangle {}\n denom_vdw {}\n ",angle_offset_flag, params[nparams].epsilon,params[nparams].sigma, params[nparams].ap,params[nparams].cut_inner, params[nparams].cut_outer, params[nparams].cut_angle, params[nparams].denom_vdw); + // flag type2param with either i,j = D,A or j,i = D,A int count = 0; @@ -381,6 +408,9 @@ void PairHbondDreidingLJ::coeff(int narg, char **arg) } nparams++; + if (comm->me == 0) + utils::logmesg(lmp,"---> DEBUG BASECLASS coeff END nparams {} count {}\n", nparams, count); + if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients"); } @@ -540,6 +570,13 @@ double PairHbondDreidingLJ::single(int i, int j, int itype, int jtype, if (c < -1.0) c = -1.0; ac = acos(c); + if (angle_offset_flag){ + ac = ac + pm.angle_offset; + c = cos(ac); + if (c > 1.0) c = 1.0; + if (c < -1.0) c = -1.0; + } + if (ac < pm.cut_angle || ac > (2.0*MY_PI - pm.cut_angle)) return 0.0; s = sqrt(1.0 - c*c); if (s < SMALL) s = SMALL; diff --git a/src/MOLECULE/pair_hbond_dreiding_lj.h b/src/MOLECULE/pair_hbond_dreiding_lj.h index 91a906c5cb..e5cde8da69 100644 --- a/src/MOLECULE/pair_hbond_dreiding_lj.h +++ b/src/MOLECULE/pair_hbond_dreiding_lj.h @@ -25,6 +25,7 @@ PairStyle(hbond/dreiding/lj,PairHbondDreidingLJ); namespace LAMMPS_NS { class PairHbondDreidingLJ : public Pair { + public: PairHbondDreidingLJ(class LAMMPS *); ~PairHbondDreidingLJ() override; @@ -38,6 +39,7 @@ class PairHbondDreidingLJ : public Pair { protected: double cut_inner_global, cut_outer_global, cut_angle_global; int ap_global; + int angle_offset_flag; // 1 if angle offset variant used struct Param { double epsilon, sigma; @@ -45,7 +47,7 @@ class PairHbondDreidingLJ : public Pair { double d0, alpha, r0; double morse1; double denom_vdw; - double cut_inner, cut_outer, cut_innersq, cut_outersq, cut_angle, offset; + double cut_inner, cut_outer, cut_innersq, cut_outersq, cut_angle, offset, angle_offset; int ap; }; diff --git a/src/MOLECULE/pair_hbond_dreiding_morse.cpp b/src/MOLECULE/pair_hbond_dreiding_morse.cpp index c8bc0a627d..11cb9cc354 100644 --- a/src/MOLECULE/pair_hbond_dreiding_morse.cpp +++ b/src/MOLECULE/pair_hbond_dreiding_morse.cpp @@ -148,6 +148,14 @@ void PairHbondDreidingMorse::compute(int eflag, int vflag) if (c < -1.0) c = -1.0; ac = acos(c); + if (angle_offset_flag){ + ac = ac + pm.angle_offset; + c = cos(ac); + if (c > 1.0) c = 1.0; + if (c < -1.0) c = -1.0; + } + + if (ac > pm.cut_angle && ac < (2.0*MY_PI - pm.cut_angle)) { s = sqrt(1.0 - c*c); if (s < SMALL) s = SMALL; @@ -295,6 +303,7 @@ void PairHbondDreidingMorse::coeff(int narg, char **arg) (params[nparams].cut_outersq-params[nparams].cut_innersq) * (params[nparams].cut_outersq-params[nparams].cut_innersq) * (params[nparams].cut_outersq-params[nparams].cut_innersq); + // if (angle_offset_flag) params[nparams].angle_offset = angle_offset_one; // flag type2param with either i,j = D,A or j,i = D,A @@ -443,6 +452,13 @@ double PairHbondDreidingMorse::single(int i, int j, int itype, int jtype, if (c < -1.0) c = -1.0; ac = acos(c); + if (angle_offset_flag){ + ac = ac + pm.angle_offset; + c = cos(ac); + if (c > 1.0) c = 1.0; + if (c < -1.0) c = -1.0; + } + if (ac < pm.cut_angle || ac > (2.0*MY_PI - pm.cut_angle)) return 0.0; s = sqrt(1.0 - c*c); if (s < SMALL) s = SMALL; From 41555a66e995a802f627d17f210a4394ce32f2aa Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Wed, 15 Jan 2025 04:20:34 -0500 Subject: [PATCH 15/29] correct documentation and add versionadded tag --- doc/src/Commands_pair.rst | 4 +-- doc/src/pair_hbond_dreiding.rst | 48 ++++++++++++++++++++------------- doc/src/pair_style.rst | 4 +-- 3 files changed, 34 insertions(+), 22 deletions(-) diff --git a/doc/src/Commands_pair.rst b/doc/src/Commands_pair.rst index 8c0e09b446..048a54ed37 100644 --- a/doc/src/Commands_pair.rst +++ b/doc/src/Commands_pair.rst @@ -115,9 +115,9 @@ OPT. * :doc:`gw/zbl ` * :doc:`harmonic/cut (o) ` * :doc:`hbond/dreiding/lj (o) ` - * :doc:`hbond/dreiding/lj/angleoffset (o) ` + * :doc:`hbond/dreiding/lj/angleoffset (o) ` * :doc:`hbond/dreiding/morse (o) ` - * :doc:`hbond/dreiding/morse/angleoffset (o) ` + * :doc:`hbond/dreiding/morse/angleoffset (o) ` * :doc:`hdnnp ` * :doc:`hippo (g) ` * :doc:`ilp/graphene/hbn (t) ` diff --git a/doc/src/pair_hbond_dreiding.rst b/doc/src/pair_hbond_dreiding.rst index 363186f108..3b4c6e08a2 100644 --- a/doc/src/pair_hbond_dreiding.rst +++ b/doc/src/pair_hbond_dreiding.rst @@ -35,7 +35,7 @@ Syntax pair_style style N inner_distance_cutoff outer_distance_cutoff angle_cutoff equilibrium_angle -* style = *hbond/dreiding/lj* or *hbond/dreiding/morse* or *hbond/dreiding/lj/angleoffset* or *hbond/dreiding/morse/angleoffset +* style = *hbond/dreiding/lj* or *hbond/dreiding/morse* or *hbond/dreiding/lj/angleoffset* or *hbond/dreiding/morse/angleoffset* * N = power of cosine of angle theta (integer) * inner_distance_cutoff = global inner cutoff for Donor-Acceptor interactions (distance units) * outer_distance_cutoff = global cutoff for Donor-Acceptor interactions (distance units) @@ -93,44 +93,53 @@ hydrogen (H) and the donor atoms: .. image:: JPG/dreiding_hbond.jpg :align: center -These 3-body interactions can be defined for pairs of acceptor and -donor atoms, based on atom types. For each donor/acceptor atom pair, -the third atom in the interaction is a hydrogen permanently bonded to -the donor atom, e.g. in a bond list read in from a data file via the +These 3-body interactions can be defined for pairs of acceptor and donor +atoms, based on atom types. For each donor/acceptor atom pair, the +third atom in the interaction is a hydrogen permanently bonded to the +donor atom, e.g. in a bond list read in from a data file via the :doc:`read_data ` command. The atom types of possible hydrogen atoms for each donor/acceptor type pair are specified by the :doc:`pair_coeff ` command (see below). Style *hbond/dreiding/lj* is the original DREIDING potential of -:ref:`(Mayo) `. It uses a LJ 12/10 functional for the Donor-Acceptor -interactions. To match the results in the original paper, use n = 4. +:ref:`(Mayo) `. It uses a LJ 12/10 functional for the +Donor-Acceptor interactions. To match the results in the original paper, +use n = 4. Style *hbond/dreiding/morse* is an improved version using a Morse potential for the Donor-Acceptor interactions. :ref:`(Liu) ` showed that the Morse form gives improved results for Dendrimer simulations, when n = 2. -The style variants *hbond/dreiding/lj/angleoffset* and *hbond/dreiding/lj/angleoffset* take the equilibrium angle of the AHD as input, allowing it to reach 180 degrees. This variant option was added to account for cases (especially in some coarse-grained models) in which the equilibrium state of the bonds may equal the minimum energy state. +.. versionadded:: TBD + +The style variants *hbond/dreiding/lj/angleoffset* and +*hbond/dreiding/lj/angleoffset* take the equilibrium angle of the AHD as +input, allowing it to reach 180 degrees. This variant option was added +to account for cases (especially in some coarse-grained models) in which +the equilibrium state of the bonds may equal the minimum energy state. See the :doc:`Howto bioFF ` page for more information on the DREIDING force field. .. note:: - Because the Dreiding hydrogen bond potential is only one portion - of an overall force field which typically includes other pairwise - interactions, it is common to use it as a sub-style in a :doc:`pair_style hybrid/overlay ` command, where another pair style - provides the repulsive core interaction between pairs of atoms, e.g. a - 1/r\^12 Lennard-Jones repulsion. + Because the Dreiding hydrogen bond potential is only one portion of + an overall force field which typically includes other pairwise + interactions, it is common to use it as a sub-style in a + :doc:`pair_style hybrid/overlay ` command, where another + pair style provides the repulsive core interaction between pairs of + atoms, e.g. a 1/r\^12 Lennard-Jones repulsion. .. note:: - When using the hbond/dreiding pair styles with :doc:`pair_style hybrid/overlay `, you should explicitly define pair + When using the hbond/dreiding pair styles with :doc:`pair_style + hybrid/overlay `, you should explicitly define pair interactions between the donor atom and acceptor atoms, (as well as between these atoms and ALL other atoms in your system). Whenever - :doc:`pair_style hybrid/overlay ` is used, ordinary mixing - rules are not applied to atoms like the donor and acceptor atoms - because they are typically referenced in multiple pair styles. + :doc:`pair_style hybrid/overlay ` is used, ordinary + mixing rules are not applied to atoms like the donor and acceptor + atoms because they are typically referenced in multiple pair styles. Neglecting to do this can cause difficult-to-detect physics problems. .. note:: @@ -142,7 +151,10 @@ on the DREIDING force field. .. note:: - For the *angle_offset* variants, the referenced angle offset is the supplementary angle of the equilibrium angle parameter. It means if the equilibrium angle is 166.6 degrees, the calculated angle offset is 13.4 degrees. + For the *angleoffset* variants, the referenced angle offset is the + supplementary angle of the equilibrium angle parameter. It means if + the equilibrium angle is 166.6 degrees, the calculated angle offset + is 13.4 degrees. ---------- diff --git a/doc/src/pair_style.rst b/doc/src/pair_style.rst index 709383bb93..bdf06d6b66 100644 --- a/doc/src/pair_style.rst +++ b/doc/src/pair_style.rst @@ -207,9 +207,9 @@ accelerated styles exist. * :doc:`gw/zbl ` - Gao-Weber potential with a repulsive ZBL core * :doc:`harmonic/cut ` - repulsive-only harmonic potential * :doc:`hbond/dreiding/lj ` - DREIDING hydrogen bonding LJ potential -* :doc:`hbond/dreiding/lj/angleoffset ` - DREIDING hydrogen bonding LJ potential with offset for hbond angle +* :doc:`hbond/dreiding/lj/angleoffset ` - DREIDING hydrogen bonding LJ potential with offset for hbond angle * :doc:`hbond/dreiding/morse ` - DREIDING hydrogen bonding Morse potential -* :doc:`hbond/dreiding/morse/angleoffset ` - DREIDING hydrogen bonding Morse potential with offset for hbond angle +* :doc:`hbond/dreiding/morse/angleoffset ` - DREIDING hydrogen bonding Morse potential with offset for hbond angle * :doc:`hdnnp ` - High-dimensional neural network potential * :doc:`hippo ` - * :doc:`ilp/graphene/hbn ` - registry-dependent interlayer potential (ILP) From 11a790a04a7247c695e463aba87de75a4266b81e Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Wed, 15 Jan 2025 04:36:51 -0500 Subject: [PATCH 16/29] angle_offset_one and cut_angle_one are only local variables --- src/EXTRA-MOLECULE/pair_hbond_dreiding_lj_angleoffset.cpp | 1 - src/EXTRA-MOLECULE/pair_hbond_dreiding_lj_angleoffset.h | 2 +- src/EXTRA-MOLECULE/pair_hbond_dreiding_morse_angleoffset.cpp | 1 - src/EXTRA-MOLECULE/pair_hbond_dreiding_morse_angleoffset.h | 2 +- 4 files changed, 2 insertions(+), 4 deletions(-) diff --git a/src/EXTRA-MOLECULE/pair_hbond_dreiding_lj_angleoffset.cpp b/src/EXTRA-MOLECULE/pair_hbond_dreiding_lj_angleoffset.cpp index 829c750536..403377b3e0 100644 --- a/src/EXTRA-MOLECULE/pair_hbond_dreiding_lj_angleoffset.cpp +++ b/src/EXTRA-MOLECULE/pair_hbond_dreiding_lj_angleoffset.cpp @@ -47,7 +47,6 @@ PairHbondDreidingLJAngleoffset::PairHbondDreidingLJAngleoffset(LAMMPS *lmp) angle_offset_flag = 1; angle_offset_global = 0.0; - angle_offset_one = 0.0; } /* ---------------------------------------------------------------------- diff --git a/src/EXTRA-MOLECULE/pair_hbond_dreiding_lj_angleoffset.h b/src/EXTRA-MOLECULE/pair_hbond_dreiding_lj_angleoffset.h index 590c5198b7..6a286d9306 100644 --- a/src/EXTRA-MOLECULE/pair_hbond_dreiding_lj_angleoffset.h +++ b/src/EXTRA-MOLECULE/pair_hbond_dreiding_lj_angleoffset.h @@ -32,7 +32,7 @@ namespace LAMMPS_NS { void coeff(int, char **) override; protected: - double angle_offset_global, angle_offset_one, cut_angle_one; + double angle_offset_global; }; } // namespace LAMMPS_NS diff --git a/src/EXTRA-MOLECULE/pair_hbond_dreiding_morse_angleoffset.cpp b/src/EXTRA-MOLECULE/pair_hbond_dreiding_morse_angleoffset.cpp index 227b84cc4d..ba23b03775 100644 --- a/src/EXTRA-MOLECULE/pair_hbond_dreiding_morse_angleoffset.cpp +++ b/src/EXTRA-MOLECULE/pair_hbond_dreiding_morse_angleoffset.cpp @@ -45,7 +45,6 @@ PairHbondDreidingMorseAngleoffset::PairHbondDreidingMorseAngleoffset(LAMMPS *lmp angle_offset_flag = 1; angle_offset_global = 0.0; - angle_offset_one = 0.0; } diff --git a/src/EXTRA-MOLECULE/pair_hbond_dreiding_morse_angleoffset.h b/src/EXTRA-MOLECULE/pair_hbond_dreiding_morse_angleoffset.h index 2622aed852..3cfb68821f 100644 --- a/src/EXTRA-MOLECULE/pair_hbond_dreiding_morse_angleoffset.h +++ b/src/EXTRA-MOLECULE/pair_hbond_dreiding_morse_angleoffset.h @@ -33,7 +33,7 @@ class PairHbondDreidingMorseAngleoffset : public PairHbondDreidingLJ { void coeff(int, char **) override; protected: - double angle_offset_global, angle_offset_one, cut_angle_one; + double angle_offset_global; }; } // namespace LAMMPS_NS From dfd8631394b66976b2202a9026012866bcbe2261 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Wed, 15 Jan 2025 04:40:48 -0500 Subject: [PATCH 17/29] remove debug output --- src/MOLECULE/pair_hbond_dreiding_lj.cpp | 18 +----------------- 1 file changed, 1 insertion(+), 17 deletions(-) diff --git a/src/MOLECULE/pair_hbond_dreiding_lj.cpp b/src/MOLECULE/pair_hbond_dreiding_lj.cpp index 66055f754e..65b333ba02 100644 --- a/src/MOLECULE/pair_hbond_dreiding_lj.cpp +++ b/src/MOLECULE/pair_hbond_dreiding_lj.cpp @@ -20,7 +20,6 @@ #include "atom.h" #include "atom_vec.h" -#include "comm.h" #include "domain.h" #include "error.h" #include "force.h" @@ -360,26 +359,17 @@ void PairHbondDreidingLJ::coeff(int narg, char **arg) double cut_angle_one = cut_angle_global; if (narg == 10) cut_angle_one = utils::numeric(FLERR, arg[9], false, lmp) * MY_PI/180.0; - if (comm->me == 0) - utils::logmesg(lmp,"---> DEBUG BASECLASS cut_angle_one {} \n", cut_angle_one); - // grow params array if necessary - if (comm->me == 0) - utils::logmesg(lmp,"---> DEBUG BASECLASS coeff nparams {} maxparam {}\n", nparams, maxparam); if (nparams == maxparam) { maxparam += CHUNK; - params = (Param *) memory->srealloc(params, maxparam*sizeof(Param), - "pair:params"); + params = (Param *) memory->srealloc(params, maxparam*sizeof(Param), "pair:params"); // make certain all addional allocated storage is initialized // to avoid false positives when checking with valgrind memset(params + nparams, 0, CHUNK*sizeof(Param)); } - if (comm->me == 0) - utils::logmesg(lmp,"---> DEBUG BASECLASS coeff POST MEM nparams {} maxparam {}\n", nparams, maxparam); - params[nparams].epsilon = epsilon_one; params[nparams].sigma = sigma_one; params[nparams].ap = ap_one; @@ -393,9 +383,6 @@ void PairHbondDreidingLJ::coeff(int narg, char **arg) (params[nparams].cut_outersq-params[nparams].cut_innersq) * (params[nparams].cut_outersq-params[nparams].cut_innersq); - if (comm->me == 0) - utils::logmesg(lmp,"BASECLASS PARAMS angle_offset_flag {}\n eps {}\n sigma {}\n ap {}\n cin {}\n cout {}\n cangle {}\n denom_vdw {}\n ",angle_offset_flag, params[nparams].epsilon,params[nparams].sigma, params[nparams].ap,params[nparams].cut_inner, params[nparams].cut_outer, params[nparams].cut_angle, params[nparams].denom_vdw); - // flag type2param with either i,j = D,A or j,i = D,A int count = 0; @@ -408,9 +395,6 @@ void PairHbondDreidingLJ::coeff(int narg, char **arg) } nparams++; - if (comm->me == 0) - utils::logmesg(lmp,"---> DEBUG BASECLASS coeff END nparams {} count {}\n", nparams, count); - if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients"); } From 865ce67e83aba8f53960b431b3e02bbb15c52e21 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Wed, 15 Jan 2025 04:48:07 -0500 Subject: [PATCH 18/29] use correct base class --- .../pair_hbond_dreiding_morse_angleoffset.cpp | 6 +++--- src/EXTRA-MOLECULE/pair_hbond_dreiding_morse_angleoffset.h | 7 +++---- 2 files changed, 6 insertions(+), 7 deletions(-) diff --git a/src/EXTRA-MOLECULE/pair_hbond_dreiding_morse_angleoffset.cpp b/src/EXTRA-MOLECULE/pair_hbond_dreiding_morse_angleoffset.cpp index ba23b03775..0d2dbf4fcf 100644 --- a/src/EXTRA-MOLECULE/pair_hbond_dreiding_morse_angleoffset.cpp +++ b/src/EXTRA-MOLECULE/pair_hbond_dreiding_morse_angleoffset.cpp @@ -41,7 +41,7 @@ static constexpr int CHUNK = 8; /* ---------------------------------------------------------------------- */ PairHbondDreidingMorseAngleoffset::PairHbondDreidingMorseAngleoffset(LAMMPS *lmp) : - PairHbondDreidingLJ(lmp) { + PairHbondDreidingMorse(lmp) { angle_offset_flag = 1; angle_offset_global = 0.0; @@ -58,7 +58,7 @@ void PairHbondDreidingMorseAngleoffset::settings(int narg, char **arg) // use first four settings args in parent method to update other variables - PairHbondDreidingLJ::settings(4,arg); + PairHbondDreidingMorse::settings(4,arg); angle_offset_global = (180.0 - utils::numeric(FLERR, arg[4], false, lmp)) * MY_PI/180.0; @@ -103,7 +103,7 @@ void PairHbondDreidingMorseAngleoffset::coeff(int narg, char **arg) double angle_offset_one = angle_offset_global; if (narg == 12) angle_offset_one = (180.0 - utils::numeric(FLERR, arg[11], false, lmp)) * MY_PI/180.0; if (angle_offset_one < 0.0 || angle_offset_one > 90.0 * MY_PI/180.0) - error->all(FLERR,"Illegal angle offset"); + error->all(FLERR,"Illegal angle offset {}", angle_offset_one); // grow params array if necessary diff --git a/src/EXTRA-MOLECULE/pair_hbond_dreiding_morse_angleoffset.h b/src/EXTRA-MOLECULE/pair_hbond_dreiding_morse_angleoffset.h index 3cfb68821f..82c5b3bd50 100644 --- a/src/EXTRA-MOLECULE/pair_hbond_dreiding_morse_angleoffset.h +++ b/src/EXTRA-MOLECULE/pair_hbond_dreiding_morse_angleoffset.h @@ -20,12 +20,11 @@ PairStyle(hbond/dreiding/morse/angleoffset,PairHbondDreidingMorseAngleoffset); #ifndef LMP_PAIR_HBOND_DREIDING_MORSE_ANGLEOFFSET_H #define LMP_PAIR_HBOND_DREIDING_MORSE_ANGLEOFFSET_H -#include "pair_hbond_dreiding_lj.h" -// #include "pair_hbond_dreiding_morse.h" +#include "pair_hbond_dreiding_morse.h" namespace LAMMPS_NS { -class PairHbondDreidingMorseAngleoffset : public PairHbondDreidingLJ { +class PairHbondDreidingMorseAngleoffset : public PairHbondDreidingMorse { public: PairHbondDreidingMorseAngleoffset(class LAMMPS *); @@ -33,7 +32,7 @@ class PairHbondDreidingMorseAngleoffset : public PairHbondDreidingLJ { void coeff(int, char **) override; protected: - double angle_offset_global; + double angle_offset_global; }; } // namespace LAMMPS_NS From f1fb0906be91f33e943352a843d58ba42b821fdc Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Wed, 15 Jan 2025 04:49:21 -0500 Subject: [PATCH 19/29] correct class name --- .../pair_hbond_dreiding_lj_angleoffset_omp.cpp | 14 +++++++------- .../pair_hbond_dreiding_lj_angleoffset_omp.h | 8 ++++---- 2 files changed, 11 insertions(+), 11 deletions(-) diff --git a/src/OPENMP/pair_hbond_dreiding_lj_angleoffset_omp.cpp b/src/OPENMP/pair_hbond_dreiding_lj_angleoffset_omp.cpp index 094beca70a..9d4b06af43 100644 --- a/src/OPENMP/pair_hbond_dreiding_lj_angleoffset_omp.cpp +++ b/src/OPENMP/pair_hbond_dreiding_lj_angleoffset_omp.cpp @@ -37,8 +37,8 @@ static constexpr double SMALL = 0.001; /* ---------------------------------------------------------------------- */ -PairHbondDreidingLJangleoffsetOMP::PairHbondDreidingLJangleoffsetOMP(LAMMPS *lmp) : - PairHbondDreidingLJangleoffset(lmp), ThrOMP(lmp, THR_PAIR) +PairHbondDreidingLJAngleoffsetOMP::PairHbondDreidingLJAngleoffsetOMP(LAMMPS *lmp) : + PairHbondDreidingLJAngleoffset(lmp), ThrOMP(lmp, THR_PAIR) { suffix_flag |= Suffix::OMP; respa_enable = 0; @@ -47,7 +47,7 @@ PairHbondDreidingLJangleoffsetOMP::PairHbondDreidingLJangleoffsetOMP(LAMMPS *lmp /* ---------------------------------------------------------------------- */ -PairHbondDreidingLJangleoffsetOMP::~PairHbondDreidingLJangleoffsetOMP() +PairHbondDreidingLJAngleoffsetOMP::~PairHbondDreidingLJAngleoffsetOMP() { if (hbcount_thr) { delete[] hbcount_thr; @@ -57,7 +57,7 @@ PairHbondDreidingLJangleoffsetOMP::~PairHbondDreidingLJangleoffsetOMP() /* ---------------------------------------------------------------------- */ -void PairHbondDreidingLJangleoffsetOMP::compute(int eflag, int vflag) +void PairHbondDreidingLJAngleoffsetOMP::compute(int eflag, int vflag) { ev_init(eflag,vflag); @@ -115,7 +115,7 @@ void PairHbondDreidingLJangleoffsetOMP::compute(int eflag, int vflag) } template -void PairHbondDreidingLJangleoffsetOMP::eval(int iifrom, int iito, ThrData * const thr) +void PairHbondDreidingLJAngleoffsetOMP::eval(int iifrom, int iito, ThrData * const thr) { int i,j,k,m,ii,jj,kk,jnum,knum,itype,jtype,ktype,iatom,imol; tagint tagprev; @@ -312,11 +312,11 @@ void PairHbondDreidingLJangleoffsetOMP::eval(int iifrom, int iito, ThrData * con /* ---------------------------------------------------------------------- */ -double PairHbondDreidingLJangleoffsetOMP::memory_usage() +double PairHbondDreidingLJAngleoffsetOMP::memory_usage() { double bytes = memory_usage_thr(); bytes += (double)comm->nthreads * 2 * sizeof(double); - bytes += PairHbondDreidingLJangleoffset::memory_usage(); + bytes += PairHbondDreidingLJAngleoffset::memory_usage(); return bytes; } diff --git a/src/OPENMP/pair_hbond_dreiding_lj_angleoffset_omp.h b/src/OPENMP/pair_hbond_dreiding_lj_angleoffset_omp.h index 0a4e14f68d..ef2ff9ee2e 100644 --- a/src/OPENMP/pair_hbond_dreiding_lj_angleoffset_omp.h +++ b/src/OPENMP/pair_hbond_dreiding_lj_angleoffset_omp.h @@ -17,7 +17,7 @@ #ifdef PAIR_CLASS // clang-format off -PairStyle(hbond/dreiding/lj/angleoffset/omp,PairHbondDreidingLJangleoffsetOMP); +PairStyle(hbond/dreiding/lj/angleoffset/omp,PairHbondDreidingLJAngleoffsetOMP); // clang-format on #else @@ -29,11 +29,11 @@ PairStyle(hbond/dreiding/lj/angleoffset/omp,PairHbondDreidingLJangleoffsetOMP); namespace LAMMPS_NS { -class PairHbondDreidingLJangleoffsetOMP : public PairHbondDreidingLJangleoffset, public ThrOMP { +class PairHbondDreidingLJAngleoffsetOMP : public PairHbondDreidingLJAngleoffset, public ThrOMP { public: - PairHbondDreidingLJangleoffsetOMP(class LAMMPS *); - ~PairHbondDreidingLJangleoffsetOMP() override; + PairHbondDreidingLJAngleoffsetOMP(class LAMMPS *); + ~PairHbondDreidingLJAngleoffsetOMP() override; void compute(int, int) override; double memory_usage() override; From 29fca919b1dfc732c266551e636ac2f909d13a63 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Wed, 15 Jan 2025 04:49:34 -0500 Subject: [PATCH 20/29] apply clang-format --- .../pair_hbond_dreiding_lj_angleoffset.h | 16 ++++++++-------- src/MOLECULE/pair_hbond_dreiding_morse.cpp | 1 - .../pair_hbond_dreiding_morse_angleoffset_omp.h | 4 +++- 3 files changed, 11 insertions(+), 10 deletions(-) diff --git a/src/EXTRA-MOLECULE/pair_hbond_dreiding_lj_angleoffset.h b/src/EXTRA-MOLECULE/pair_hbond_dreiding_lj_angleoffset.h index 6a286d9306..8b1a93b1bc 100644 --- a/src/EXTRA-MOLECULE/pair_hbond_dreiding_lj_angleoffset.h +++ b/src/EXTRA-MOLECULE/pair_hbond_dreiding_lj_angleoffset.h @@ -24,16 +24,16 @@ PairStyle(hbond/dreiding/lj/angleoffset,PairHbondDreidingLJAngleoffset); namespace LAMMPS_NS { - class PairHbondDreidingLJAngleoffset : public PairHbondDreidingLJ { +class PairHbondDreidingLJAngleoffset : public PairHbondDreidingLJ { - public: - PairHbondDreidingLJAngleoffset(class LAMMPS *); - void settings(int, char **) override; - void coeff(int, char **) override; + public: + PairHbondDreidingLJAngleoffset(class LAMMPS *); + void settings(int, char **) override; + void coeff(int, char **) override; - protected: - double angle_offset_global; - }; + protected: + double angle_offset_global; +}; } // namespace LAMMPS_NS diff --git a/src/MOLECULE/pair_hbond_dreiding_morse.cpp b/src/MOLECULE/pair_hbond_dreiding_morse.cpp index 15e9d49654..34fddc9f89 100644 --- a/src/MOLECULE/pair_hbond_dreiding_morse.cpp +++ b/src/MOLECULE/pair_hbond_dreiding_morse.cpp @@ -303,7 +303,6 @@ void PairHbondDreidingMorse::coeff(int narg, char **arg) (params[nparams].cut_outersq-params[nparams].cut_innersq) * (params[nparams].cut_outersq-params[nparams].cut_innersq) * (params[nparams].cut_outersq-params[nparams].cut_innersq); - // if (angle_offset_flag) params[nparams].angle_offset = angle_offset_one; // flag type2param with either i,j = D,A or j,i = D,A diff --git a/src/OPENMP/pair_hbond_dreiding_morse_angleoffset_omp.h b/src/OPENMP/pair_hbond_dreiding_morse_angleoffset_omp.h index fb95988c6f..5e87908df3 100644 --- a/src/OPENMP/pair_hbond_dreiding_morse_angleoffset_omp.h +++ b/src/OPENMP/pair_hbond_dreiding_morse_angleoffset_omp.h @@ -29,7 +29,9 @@ PairStyle(hbond/dreiding/morse/angleoffset/omp,PairHbondDreidingMorseAngleoffset namespace LAMMPS_NS { -class PairHbondDreidingMorseAngleoffsetOMP : public PairHbondDreidingMorseAngleoffset, public ThrOMP { +class PairHbondDreidingMorseAngleoffsetOMP : + public PairHbondDreidingMorseAngleoffset, + public ThrOMP { public: PairHbondDreidingMorseAngleoffsetOMP(class LAMMPS *); From 66ffb1c39e9a66a62d5f98de08eb5bc14592956b Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Wed, 15 Jan 2025 04:55:36 -0500 Subject: [PATCH 21/29] whitespace --- doc/src/pair_hbond_dreiding.rst | 4 ++-- src/EXTRA-MOLECULE/pair_hbond_dreiding_lj_angleoffset.cpp | 2 +- src/EXTRA-MOLECULE/pair_hbond_dreiding_morse_angleoffset.cpp | 4 ++-- 3 files changed, 5 insertions(+), 5 deletions(-) diff --git a/doc/src/pair_hbond_dreiding.rst b/doc/src/pair_hbond_dreiding.rst index 3b4c6e08a2..e059fdf2ba 100644 --- a/doc/src/pair_hbond_dreiding.rst +++ b/doc/src/pair_hbond_dreiding.rst @@ -40,7 +40,7 @@ Syntax * inner_distance_cutoff = global inner cutoff for Donor-Acceptor interactions (distance units) * outer_distance_cutoff = global cutoff for Donor-Acceptor interactions (distance units) * angle_cutoff = global angle cutoff for Acceptor-Hydrogen-Donor interactions (degrees) -* (with style angleoffset) equilibrium_angle = global equilibrium angle for Acceptor-Hydrogen-Donor interactions (degrees) +* (with style angleoffset) equilibrium_angle = global equilibrium angle for Acceptor-Hydrogen-Donor interactions (degrees) Examples """""""" @@ -207,7 +207,7 @@ follows: * angle cutoff (degrees) For both the *hbond/dreiding/lj/angleoffset* and *hbond/dreiding/morse/angleoffset* styles an additional parameter is added: -* equilibrium angle (degrees) +* equilibrium angle (degrees) For all styles, a single hydrogen atom type K can be specified, or a wild-card asterisk can be used in place of or in conjunction with the K arguments to diff --git a/src/EXTRA-MOLECULE/pair_hbond_dreiding_lj_angleoffset.cpp b/src/EXTRA-MOLECULE/pair_hbond_dreiding_lj_angleoffset.cpp index 403377b3e0..390b72725d 100644 --- a/src/EXTRA-MOLECULE/pair_hbond_dreiding_lj_angleoffset.cpp +++ b/src/EXTRA-MOLECULE/pair_hbond_dreiding_lj_angleoffset.cpp @@ -42,7 +42,7 @@ using namespace MathConst; /* ---------------------------------------------------------------------- */ -PairHbondDreidingLJAngleoffset::PairHbondDreidingLJAngleoffset(LAMMPS *lmp) +PairHbondDreidingLJAngleoffset::PairHbondDreidingLJAngleoffset(LAMMPS *lmp) : PairHbondDreidingLJ(lmp) { angle_offset_flag = 1; diff --git a/src/EXTRA-MOLECULE/pair_hbond_dreiding_morse_angleoffset.cpp b/src/EXTRA-MOLECULE/pair_hbond_dreiding_morse_angleoffset.cpp index 0d2dbf4fcf..37bcdd7ee2 100644 --- a/src/EXTRA-MOLECULE/pair_hbond_dreiding_morse_angleoffset.cpp +++ b/src/EXTRA-MOLECULE/pair_hbond_dreiding_morse_angleoffset.cpp @@ -57,11 +57,11 @@ void PairHbondDreidingMorseAngleoffset::settings(int narg, char **arg) if (narg != 5) error->all(FLERR,"Illegal pair_style command"); // use first four settings args in parent method to update other variables - + PairHbondDreidingMorse::settings(4,arg); angle_offset_global = (180.0 - utils::numeric(FLERR, arg[4], false, lmp)) * MY_PI/180.0; - + } // /* ---------------------------------------------------------------------- From 7fa1bf39f3d2db7488f5a315f770cef01d29cdf3 Mon Sep 17 00:00:00 2001 From: megmcca Date: Wed, 15 Jan 2025 10:24:25 -0700 Subject: [PATCH 22/29] shift settings check from angleoffset to base LJ --- .../pair_hbond_dreiding_lj_angleoffset.cpp | 17 ----------------- .../pair_hbond_dreiding_lj_angleoffset.h | 3 --- .../pair_hbond_dreiding_morse_angleoffset.cpp | 17 ----------------- .../pair_hbond_dreiding_morse_angleoffset.h | 5 +---- src/MOLECULE/pair_hbond_dreiding_lj.cpp | 10 +++++++++- src/MOLECULE/pair_hbond_dreiding_lj.h | 1 + 6 files changed, 11 insertions(+), 42 deletions(-) diff --git a/src/EXTRA-MOLECULE/pair_hbond_dreiding_lj_angleoffset.cpp b/src/EXTRA-MOLECULE/pair_hbond_dreiding_lj_angleoffset.cpp index 390b72725d..02ac009773 100644 --- a/src/EXTRA-MOLECULE/pair_hbond_dreiding_lj_angleoffset.cpp +++ b/src/EXTRA-MOLECULE/pair_hbond_dreiding_lj_angleoffset.cpp @@ -46,23 +46,6 @@ PairHbondDreidingLJAngleoffset::PairHbondDreidingLJAngleoffset(LAMMPS *lmp) : PairHbondDreidingLJ(lmp) { angle_offset_flag = 1; - angle_offset_global = 0.0; -} - -/* ---------------------------------------------------------------------- - global settings -------------------------------------------------------------------------- */ - -void PairHbondDreidingLJAngleoffset::settings(int narg, char **arg) -{ - if (narg != 5) error->all(FLERR,"Illegal pair_style command"); - - // use first four settings args in parent method to update other variables - - PairHbondDreidingLJ::settings(4,arg); - - angle_offset_global = (180.0 - utils::numeric(FLERR, arg[4], false, lmp)) * MY_PI/180.0; - } /* ---------------------------------------------------------------------- diff --git a/src/EXTRA-MOLECULE/pair_hbond_dreiding_lj_angleoffset.h b/src/EXTRA-MOLECULE/pair_hbond_dreiding_lj_angleoffset.h index 8b1a93b1bc..5ae1fcba28 100644 --- a/src/EXTRA-MOLECULE/pair_hbond_dreiding_lj_angleoffset.h +++ b/src/EXTRA-MOLECULE/pair_hbond_dreiding_lj_angleoffset.h @@ -28,11 +28,8 @@ class PairHbondDreidingLJAngleoffset : public PairHbondDreidingLJ { public: PairHbondDreidingLJAngleoffset(class LAMMPS *); - void settings(int, char **) override; void coeff(int, char **) override; - protected: - double angle_offset_global; }; } // namespace LAMMPS_NS diff --git a/src/EXTRA-MOLECULE/pair_hbond_dreiding_morse_angleoffset.cpp b/src/EXTRA-MOLECULE/pair_hbond_dreiding_morse_angleoffset.cpp index 37bcdd7ee2..dedbbb8c9f 100644 --- a/src/EXTRA-MOLECULE/pair_hbond_dreiding_morse_angleoffset.cpp +++ b/src/EXTRA-MOLECULE/pair_hbond_dreiding_morse_angleoffset.cpp @@ -47,23 +47,6 @@ PairHbondDreidingMorseAngleoffset::PairHbondDreidingMorseAngleoffset(LAMMPS *lmp angle_offset_global = 0.0; } - -/* ---------------------------------------------------------------------- - global settings -------------------------------------------------------------------------- */ - -void PairHbondDreidingMorseAngleoffset::settings(int narg, char **arg) -{ - if (narg != 5) error->all(FLERR,"Illegal pair_style command"); - - // use first four settings args in parent method to update other variables - - PairHbondDreidingMorse::settings(4,arg); - - angle_offset_global = (180.0 - utils::numeric(FLERR, arg[4], false, lmp)) * MY_PI/180.0; - -} - // /* ---------------------------------------------------------------------- // set coeffs for one or more type pairs // ------------------------------------------------------------------------- */ diff --git a/src/EXTRA-MOLECULE/pair_hbond_dreiding_morse_angleoffset.h b/src/EXTRA-MOLECULE/pair_hbond_dreiding_morse_angleoffset.h index 82c5b3bd50..5eeb1796a2 100644 --- a/src/EXTRA-MOLECULE/pair_hbond_dreiding_morse_angleoffset.h +++ b/src/EXTRA-MOLECULE/pair_hbond_dreiding_morse_angleoffset.h @@ -28,11 +28,8 @@ class PairHbondDreidingMorseAngleoffset : public PairHbondDreidingMorse { public: PairHbondDreidingMorseAngleoffset(class LAMMPS *); - void settings(int, char **) override; void coeff(int, char **) override; - - protected: - double angle_offset_global; + }; } // namespace LAMMPS_NS diff --git a/src/MOLECULE/pair_hbond_dreiding_lj.cpp b/src/MOLECULE/pair_hbond_dreiding_lj.cpp index 65b333ba02..ebe7fdd123 100644 --- a/src/MOLECULE/pair_hbond_dreiding_lj.cpp +++ b/src/MOLECULE/pair_hbond_dreiding_lj.cpp @@ -57,6 +57,7 @@ PairHbondDreidingLJ::PairHbondDreidingLJ(LAMMPS *lmp) : Pair(lmp) pvector = new double[2]; angle_offset_flag = 0; + angle_offset_global = 0.0; } /* ---------------------------------------------------------------------- */ @@ -309,12 +310,19 @@ void PairHbondDreidingLJ::allocate() void PairHbondDreidingLJ::settings(int narg, char **arg) { - if (narg != 4) error->all(FLERR,"Illegal pair_style command"); + + // narg = 4 for standard form and narg = 5 if angleoffset variant (from EXTRA-MOLECULE) + if (narg != 4 && narg != 5) error->all(FLERR,"Illegal pair_style command"); ap_global = utils::inumeric(FLERR,arg[0],false,lmp); cut_inner_global = utils::numeric(FLERR,arg[1],false,lmp); cut_outer_global = utils::numeric(FLERR,arg[2],false,lmp); cut_angle_global = utils::numeric(FLERR,arg[3],false,lmp) * MY_PI/180.0; + + // update when using angleoffset variant + if (angle_offset_flag) { + angle_offset_global = (180.0 - utils::numeric(FLERR, arg[4], false, lmp)) * MY_PI/180.0; + } } /* ---------------------------------------------------------------------- diff --git a/src/MOLECULE/pair_hbond_dreiding_lj.h b/src/MOLECULE/pair_hbond_dreiding_lj.h index e5cde8da69..131f3946b3 100644 --- a/src/MOLECULE/pair_hbond_dreiding_lj.h +++ b/src/MOLECULE/pair_hbond_dreiding_lj.h @@ -40,6 +40,7 @@ class PairHbondDreidingLJ : public Pair { double cut_inner_global, cut_outer_global, cut_angle_global; int ap_global; int angle_offset_flag; // 1 if angle offset variant used + double angle_offset_global; // updated if angle offset variant used struct Param { double epsilon, sigma; From fecf1c2f692d86b5b620bee9f23af32abee12411 Mon Sep 17 00:00:00 2001 From: megmcca Date: Fri, 17 Jan 2025 15:59:28 -0700 Subject: [PATCH 23/29] update morse coeff method, fix bug in flag --- src/EXTRA-MOLECULE/pair_hbond_dreiding_morse_angleoffset.cpp | 1 - src/MOLECULE/pair_hbond_dreiding_lj.cpp | 2 +- src/MOLECULE/pair_hbond_dreiding_morse.cpp | 4 +++- 3 files changed, 4 insertions(+), 3 deletions(-) diff --git a/src/EXTRA-MOLECULE/pair_hbond_dreiding_morse_angleoffset.cpp b/src/EXTRA-MOLECULE/pair_hbond_dreiding_morse_angleoffset.cpp index dedbbb8c9f..a687d4d7ea 100644 --- a/src/EXTRA-MOLECULE/pair_hbond_dreiding_morse_angleoffset.cpp +++ b/src/EXTRA-MOLECULE/pair_hbond_dreiding_morse_angleoffset.cpp @@ -44,7 +44,6 @@ PairHbondDreidingMorseAngleoffset::PairHbondDreidingMorseAngleoffset(LAMMPS *lmp PairHbondDreidingMorse(lmp) { angle_offset_flag = 1; - angle_offset_global = 0.0; } // /* ---------------------------------------------------------------------- diff --git a/src/MOLECULE/pair_hbond_dreiding_lj.cpp b/src/MOLECULE/pair_hbond_dreiding_lj.cpp index ebe7fdd123..0511ecb021 100644 --- a/src/MOLECULE/pair_hbond_dreiding_lj.cpp +++ b/src/MOLECULE/pair_hbond_dreiding_lj.cpp @@ -311,7 +311,7 @@ void PairHbondDreidingLJ::allocate() void PairHbondDreidingLJ::settings(int narg, char **arg) { - // narg = 4 for standard form and narg = 5 if angleoffset variant (from EXTRA-MOLECULE) + // narg = 4 for standard form, narg = 5 or 6 if angleoffset LJ or Morse variants respectively (from EXTRA-MOLECULE) if (narg != 4 && narg != 5) error->all(FLERR,"Illegal pair_style command"); ap_global = utils::inumeric(FLERR,arg[0],false,lmp); diff --git a/src/MOLECULE/pair_hbond_dreiding_morse.cpp b/src/MOLECULE/pair_hbond_dreiding_morse.cpp index 34fddc9f89..56f9c16442 100644 --- a/src/MOLECULE/pair_hbond_dreiding_morse.cpp +++ b/src/MOLECULE/pair_hbond_dreiding_morse.cpp @@ -246,7 +246,9 @@ void PairHbondDreidingMorse::compute(int eflag, int vflag) void PairHbondDreidingMorse::coeff(int narg, char **arg) { - if (narg < 7 || narg > 11) + int maxarg = 12; + if (angle_offset_flag == 1) maxarg = 12; + if (narg < 7 || narg > maxarg) error->all(FLERR,"Incorrect args for pair coefficients"); if (!allocated) allocate(); From bb83497f6164a8fb7228ead44f1c4dda03d7c1c8 Mon Sep 17 00:00:00 2001 From: megmcca Date: Tue, 28 Jan 2025 12:40:50 -0700 Subject: [PATCH 24/29] add offset code to base omp files --- src/OPENMP/pair_hbond_dreiding_lj_omp.cpp | 7 +++++++ src/OPENMP/pair_hbond_dreiding_morse_omp.cpp | 7 +++++++ 2 files changed, 14 insertions(+) diff --git a/src/OPENMP/pair_hbond_dreiding_lj_omp.cpp b/src/OPENMP/pair_hbond_dreiding_lj_omp.cpp index b0f6dcfb5b..9bc5f5d254 100644 --- a/src/OPENMP/pair_hbond_dreiding_lj_omp.cpp +++ b/src/OPENMP/pair_hbond_dreiding_lj_omp.cpp @@ -222,6 +222,13 @@ void PairHbondDreidingLJOMP::eval(int iifrom, int iito, ThrData * const thr) if (c < -1.0) c = -1.0; ac = acos(c); + if (angle_offset_flag){ + ac = ac + pm.angle_offset; + c = cos(ac); + if (c > 1.0) c = 1.0; + if (c < -1.0) c = -1.0; + } + if (ac > pm.cut_angle && ac < (2.0*MY_PI - pm.cut_angle)) { s = sqrt(1.0 - c*c); if (s < SMALL) s = SMALL; diff --git a/src/OPENMP/pair_hbond_dreiding_morse_omp.cpp b/src/OPENMP/pair_hbond_dreiding_morse_omp.cpp index 0e43e2a037..bb0970dbdc 100644 --- a/src/OPENMP/pair_hbond_dreiding_morse_omp.cpp +++ b/src/OPENMP/pair_hbond_dreiding_morse_omp.cpp @@ -222,6 +222,13 @@ void PairHbondDreidingMorseOMP::eval(int iifrom, int iito, ThrData * const thr) if (c < -1.0) c = -1.0; ac = acos(c); + if (angle_offset_flag){ + ac = ac + pm.angle_offset; + c = cos(ac); + if (c > 1.0) c = 1.0; + if (c < -1.0) c = -1.0; + } + if (ac > pm.cut_angle && ac < (2.0*MY_PI - pm.cut_angle)) { s = sqrt(1.0 - c*c); if (s < SMALL) s = SMALL; From f8ac22ade7ced34117f0755e3b19864d2a9ff5f2 Mon Sep 17 00:00:00 2001 From: megmcca Date: Tue, 28 Jan 2025 12:42:34 -0700 Subject: [PATCH 25/29] add flags and set up inheritance --- ...pair_hbond_dreiding_lj_angleoffset_omp.cpp | 282 +----------------- .../pair_hbond_dreiding_lj_angleoffset_omp.h | 16 +- ...r_hbond_dreiding_morse_angleoffset_omp.cpp | 282 +----------------- ...air_hbond_dreiding_morse_angleoffset_omp.h | 18 +- 4 files changed, 9 insertions(+), 589 deletions(-) diff --git a/src/OPENMP/pair_hbond_dreiding_lj_angleoffset_omp.cpp b/src/OPENMP/pair_hbond_dreiding_lj_angleoffset_omp.cpp index 9d4b06af43..1a5126971a 100644 --- a/src/OPENMP/pair_hbond_dreiding_lj_angleoffset_omp.cpp +++ b/src/OPENMP/pair_hbond_dreiding_lj_angleoffset_omp.cpp @@ -38,285 +38,7 @@ static constexpr double SMALL = 0.001; /* ---------------------------------------------------------------------- */ PairHbondDreidingLJAngleoffsetOMP::PairHbondDreidingLJAngleoffsetOMP(LAMMPS *lmp) : - PairHbondDreidingLJAngleoffset(lmp), ThrOMP(lmp, THR_PAIR) -{ - suffix_flag |= Suffix::OMP; - respa_enable = 0; - hbcount_thr = hbeng_thr = nullptr; + PairHbondDreidingLJOMP(lmp) { + angle_offset_flag = 1; } -/* ---------------------------------------------------------------------- */ - -PairHbondDreidingLJAngleoffsetOMP::~PairHbondDreidingLJAngleoffsetOMP() -{ - if (hbcount_thr) { - delete[] hbcount_thr; - delete[] hbeng_thr; - } -} - -/* ---------------------------------------------------------------------- */ - -void PairHbondDreidingLJAngleoffsetOMP::compute(int eflag, int vflag) -{ - ev_init(eflag,vflag); - - const int nall = atom->nlocal + atom->nghost; - const int nthreads = comm->nthreads; - const int inum = list->inum; - - if (!hbcount_thr) { - hbcount_thr = new double[nthreads]; - hbeng_thr = new double[nthreads]; - } - - for (int i=0; i < nthreads; ++i) { - hbcount_thr[i] = 0.0; - hbeng_thr[i] = 0.0; - } - -#if defined(_OPENMP) -#pragma omp parallel LMP_DEFAULT_NONE LMP_SHARED(eflag,vflag) -#endif - { - int ifrom, ito, tid; - - loop_setup_thr(ifrom, ito, tid, inum, nthreads); - ThrData *thr = fix->get_thr(tid); - thr->timer(Timer::START); - ev_setup_thr(eflag, vflag, nall, eatom, vatom, nullptr, thr); - - if (evflag) { - if (eflag) { - if (force->newton_pair) eval<1,1,1>(ifrom, ito, thr); - else eval<1,1,0>(ifrom, ito, thr); - } else { - if (force->newton_pair) eval<1,0,1>(ifrom, ito, thr); - else eval<1,0,0>(ifrom, ito, thr); - } - } else { - if (force->newton_pair) eval<0,0,1>(ifrom, ito, thr); - else eval<0,0,0>(ifrom, ito, thr); - } - - thr->timer(Timer::PAIR); - reduce_thr(this, eflag, vflag, thr); - } // end of omp parallel region - - // reduce per thread hbond data - if (eflag_global) { - pvector[0] = 0.0; - pvector[1] = 0.0; - for (int i=0; i < nthreads; ++i) { - pvector[0] += hbcount_thr[i]; - pvector[1] += hbeng_thr[i]; - } - } -} - -template -void PairHbondDreidingLJAngleoffsetOMP::eval(int iifrom, int iito, ThrData * const thr) -{ - int i,j,k,m,ii,jj,kk,jnum,knum,itype,jtype,ktype,iatom,imol; - tagint tagprev; - double xtmp,ytmp,ztmp,delx,dely,delz,rsq,rsq1,rsq2,r1,r2; - double factor_hb,force_angle,force_kernel,evdwl,eng_lj; - double c,s,a,b,ac,a11,a12,a22,vx1,vx2,vy1,vy2,vz1,vz2; - double fi[3],fj[3],delr1[3],delr2[3]; - double r2inv,r10inv; - double switch1,switch2; - int *ilist,*jlist,*numneigh,**firstneigh; - const tagint *klist; - - evdwl = 0.0; - - const auto * _noalias const x = (dbl3_t *) atom->x[0]; - auto * _noalias const f = (dbl3_t *) thr->get_f()[0]; - const tagint * _noalias const tag = atom->tag; - const int * _noalias const molindex = atom->molindex; - const int * _noalias const molatom = atom->molatom; - const int * _noalias const type = atom->type; - const double * _noalias const special_lj = force->special_lj; - const int * const * const nspecial = atom->nspecial; - const tagint * const * const special = atom->special; - const int molecular = atom->molecular; - Molecule * const * const onemols = atom->avec->onemols; - double fxtmp,fytmp,fztmp; - - ilist = list->ilist; - numneigh = list->numneigh; - firstneigh = list->firstneigh; - - // ii = loop over donors - // jj = loop over acceptors - // kk = loop over hydrogens bonded to donor - - int hbcount = 0; - double hbeng = 0.0; - - for (ii = iifrom; ii < iito; ++ii) { - i = ilist[ii]; - itype = type[i]; - if (!donor[itype]) continue; - if (molecular == Atom::MOLECULAR) { - klist = special[i]; - knum = nspecial[i][0]; - } else { - if (molindex[i] < 0) continue; - imol = molindex[i]; - iatom = molatom[i]; - klist = onemols[imol]->special[iatom]; - knum = onemols[imol]->nspecial[iatom][0]; - tagprev = tag[i] - iatom - 1; - } - jlist = firstneigh[i]; - jnum = numneigh[i]; - fxtmp=fytmp=fztmp=0.0; - - xtmp = x[i].x; - ytmp = x[i].y; - ztmp = x[i].z; - - for (jj = 0; jj < jnum; jj++) { - j = jlist[jj]; - factor_hb = special_lj[sbmask(j)]; - j &= NEIGHMASK; - - jtype = type[j]; - if (!acceptor[jtype]) continue; - - delx = xtmp - x[j].x; - dely = ytmp - x[j].y; - delz = ztmp - x[j].z; - rsq = delx*delx + dely*dely + delz*delz; - - for (kk = 0; kk < knum; kk++) { - if (molecular == Atom::MOLECULAR) k = atom->map(klist[kk]); - else k = atom->map(klist[kk]+tagprev); - if (k < 0) continue; - ktype = type[k]; - m = type2param[itype][jtype][ktype]; - if (m < 0) continue; - const Param &pm = params[m]; - - if (rsq < pm.cut_outersq) { - delr1[0] = xtmp - x[k].x; - delr1[1] = ytmp - x[k].y; - delr1[2] = ztmp - x[k].z; - domain->minimum_image(delr1); - rsq1 = delr1[0]*delr1[0] + delr1[1]*delr1[1] + delr1[2]*delr1[2]; - r1 = sqrt(rsq1); - - delr2[0] = x[j].x - x[k].x; - delr2[1] = x[j].y - x[k].y; - delr2[2] = x[j].z - x[k].z; - domain->minimum_image(delr2); - rsq2 = delr2[0]*delr2[0] + delr2[1]*delr2[1] + delr2[2]*delr2[2]; - r2 = sqrt(rsq2); - - // angle (cos and sin) - - c = delr1[0]*delr2[0] + delr1[1]*delr2[1] + delr1[2]*delr2[2]; - c /= r1*r2; - if (c > 1.0) c = 1.0; - if (c < -1.0) c = -1.0; - ac = acos(c); - - ac = ac + pm.angle_offset; - c = cos(ac); - if (c > 1.0) c = 1.0; - if (c < -1.0) c = -1.0; - - if (ac > pm.cut_angle && ac < (2.0*MY_PI - pm.cut_angle)) { - s = sqrt(1.0 - c*c); - if (s < SMALL) s = SMALL; - - // LJ-specific kernel - - r2inv = 1.0/rsq; - r10inv = r2inv*r2inv*r2inv*r2inv*r2inv; - force_kernel = r10inv*(pm.lj1*r2inv - pm.lj2)*r2inv * - powint(c,pm.ap); - force_angle = pm.ap * r10inv*(pm.lj3*r2inv - pm.lj4) * - powint(c,pm.ap-1)*s; - - eng_lj = r10inv*(pm.lj3*r2inv - pm.lj4); - if (rsq > pm.cut_innersq) { - switch1 = (pm.cut_outersq-rsq) * (pm.cut_outersq-rsq) * - (pm.cut_outersq + 2.0*rsq - 3.0*pm.cut_innersq) / - pm.denom_vdw; - switch2 = 12.0*rsq * (pm.cut_outersq-rsq) * - (rsq-pm.cut_innersq) / pm.denom_vdw; - force_kernel = force_kernel*switch1 + eng_lj*switch2/rsq; - force_angle *= switch1; - eng_lj *= switch1; - } - - if (EFLAG) { - evdwl = eng_lj * powint(c,pm.ap); - evdwl *= factor_hb; - } - - a = factor_hb*force_angle/s; - b = factor_hb*force_kernel; - - a11 = a*c / rsq1; - a12 = -a / (r1*r2); - a22 = a*c / rsq2; - - vx1 = a11*delr1[0] + a12*delr2[0]; - vx2 = a22*delr2[0] + a12*delr1[0]; - vy1 = a11*delr1[1] + a12*delr2[1]; - vy2 = a22*delr2[1] + a12*delr1[1]; - vz1 = a11*delr1[2] + a12*delr2[2]; - vz2 = a22*delr2[2] + a12*delr1[2]; - - fi[0] = vx1 + b*delx; - fi[1] = vy1 + b*dely; - fi[2] = vz1 + b*delz; - fj[0] = vx2 - b*delx; - fj[1] = vy2 - b*dely; - fj[2] = vz2 - b*delz; - - fxtmp += fi[0]; - fytmp += fi[1]; - fztmp += fi[2]; - - f[j].x += fj[0]; - f[j].y += fj[1]; - f[j].z += fj[2]; - - f[k].x -= vx1 + vx2; - f[k].y -= vy1 + vy2; - f[k].z -= vz1 + vz2; - - // KIJ instead of IJK b/c delr1/delr2 are both with respect to k - - if (EVFLAG) ev_tally3_thr(this,k,i,j,evdwl,0.0,fi,fj,delr1,delr2,thr); - if (EFLAG) { - hbcount++; - hbeng += evdwl; - } - } - } - } - } - f[i].x += fxtmp; - f[i].y += fytmp; - f[i].z += fztmp; - } - const int tid = thr->get_tid(); - hbcount_thr[tid] = static_cast(hbcount); - hbeng_thr[tid] = hbeng; -} - -/* ---------------------------------------------------------------------- */ - -double PairHbondDreidingLJAngleoffsetOMP::memory_usage() -{ - double bytes = memory_usage_thr(); - bytes += (double)comm->nthreads * 2 * sizeof(double); - bytes += PairHbondDreidingLJAngleoffset::memory_usage(); - - return bytes; -} diff --git a/src/OPENMP/pair_hbond_dreiding_lj_angleoffset_omp.h b/src/OPENMP/pair_hbond_dreiding_lj_angleoffset_omp.h index ef2ff9ee2e..ce79a03dee 100644 --- a/src/OPENMP/pair_hbond_dreiding_lj_angleoffset_omp.h +++ b/src/OPENMP/pair_hbond_dreiding_lj_angleoffset_omp.h @@ -24,26 +24,14 @@ PairStyle(hbond/dreiding/lj/angleoffset/omp,PairHbondDreidingLJAngleoffsetOMP); #ifndef LMP_PAIR_HBOND_DREIDING_LJ_ANGLEOFFSET_OMP_H #define LMP_PAIR_HBOND_DREIDING_LJ_ANGLEOFFSET_OMP_H -#include "pair_hbond_dreiding_lj_angleoffset.h" -#include "thr_omp.h" +#include "pair_hbond_dreiding_lj_omp.h" namespace LAMMPS_NS { -class PairHbondDreidingLJAngleoffsetOMP : public PairHbondDreidingLJAngleoffset, public ThrOMP { +class PairHbondDreidingLJAngleoffsetOMP : public PairHbondDreidingLJOMP { public: PairHbondDreidingLJAngleoffsetOMP(class LAMMPS *); - ~PairHbondDreidingLJAngleoffsetOMP() override; - - void compute(int, int) override; - double memory_usage() override; - - protected: - double *hbcount_thr, *hbeng_thr; - - private: - template - void eval(int ifrom, int ito, ThrData *const thr); }; } // namespace LAMMPS_NS diff --git a/src/OPENMP/pair_hbond_dreiding_morse_angleoffset_omp.cpp b/src/OPENMP/pair_hbond_dreiding_morse_angleoffset_omp.cpp index ad5c43f7b7..b0f9cde92f 100644 --- a/src/OPENMP/pair_hbond_dreiding_morse_angleoffset_omp.cpp +++ b/src/OPENMP/pair_hbond_dreiding_morse_angleoffset_omp.cpp @@ -38,284 +38,6 @@ static constexpr double SMALL = 0.001; /* ---------------------------------------------------------------------- */ PairHbondDreidingMorseAngleoffsetOMP::PairHbondDreidingMorseAngleoffsetOMP(LAMMPS *lmp) : - PairHbondDreidingMorseAngleoffset(lmp), ThrOMP(lmp, THR_PAIR) -{ - suffix_flag |= Suffix::OMP; - respa_enable = 0; - hbcount_thr = hbeng_thr = nullptr; -} - -/* ---------------------------------------------------------------------- */ - -PairHbondDreidingMorseAngleoffsetOMP::~PairHbondDreidingMorseAngleoffsetOMP() -{ - if (hbcount_thr) { - delete[] hbcount_thr; - delete[] hbeng_thr; - } -} - -/* ---------------------------------------------------------------------- */ - -void PairHbondDreidingMorseAngleoffsetOMP::compute(int eflag, int vflag) -{ - ev_init(eflag,vflag); - - const int nall = atom->nlocal + atom->nghost; - const int nthreads = comm->nthreads; - const int inum = list->inum; - - if (!hbcount_thr) { - hbcount_thr = new double[nthreads]; - hbeng_thr = new double[nthreads]; - } - - for (int i=0; i < nthreads; ++i) { - hbcount_thr[i] = 0.0; - hbeng_thr[i] = 0.0; - } - -#if defined(_OPENMP) -#pragma omp parallel LMP_DEFAULT_NONE LMP_SHARED(eflag,vflag) -#endif - { - int ifrom, ito, tid; - - loop_setup_thr(ifrom, ito, tid, inum, nthreads); - ThrData *thr = fix->get_thr(tid); - thr->timer(Timer::START); - ev_setup_thr(eflag, vflag, nall, eatom, vatom, nullptr, thr); - - if (evflag) { - if (eflag) { - if (force->newton_pair) eval<1,1,1>(ifrom, ito, thr); - else eval<1,1,0>(ifrom, ito, thr); - } else { - if (force->newton_pair) eval<1,0,1>(ifrom, ito, thr); - else eval<1,0,0>(ifrom, ito, thr); - } - } else { - if (force->newton_pair) eval<0,0,1>(ifrom, ito, thr); - else eval<0,0,0>(ifrom, ito, thr); - } - - thr->timer(Timer::PAIR); - reduce_thr(this, eflag, vflag, thr); - } // end of omp parallel region - - // reduce per thread hbond data - if (eflag_global) { - pvector[0] = 0.0; - pvector[1] = 0.0; - for (int i=0; i < nthreads; ++i) { - pvector[0] += hbcount_thr[i]; - pvector[1] += hbeng_thr[i]; - } - } -} - -template -void PairHbondDreidingMorseAngleoffsetOMP::eval(int iifrom, int iito, ThrData * const thr) -{ - int i,j,k,m,ii,jj,kk,jnum,knum,itype,jtype,ktype,imol,iatom; - tagint tagprev; - double xtmp,ytmp,ztmp,delx,dely,delz,rsq,rsq1,rsq2,r1,r2; - double factor_hb,force_angle,force_kernel,evdwl; - double c,s,a,b,ac,a11,a12,a22,vx1,vx2,vy1,vy2,vz1,vz2; - double fi[3],fj[3],delr1[3],delr2[3]; - double r,dr,dexp,eng_morse,switch1,switch2; - int *ilist,*jlist,*numneigh,**firstneigh; - const tagint *klist; - - evdwl = 0.0; - - const auto * _noalias const x = (dbl3_t *) atom->x[0]; - auto * _noalias const f = (dbl3_t *) thr->get_f()[0]; - const tagint * _noalias const tag = atom->tag; - const int * _noalias const type = atom->type; - const int * _noalias const molindex = atom->molindex; - const int * _noalias const molatom = atom->molatom; - const double * _noalias const special_lj = force->special_lj; - const int * const * const nspecial = atom->nspecial; - const tagint * const * const special = atom->special; - const int molecular = atom->molecular; - Molecule * const * const onemols = atom->avec->onemols; - double fxtmp,fytmp,fztmp; - - ilist = list->ilist; - numneigh = list->numneigh; - firstneigh = list->firstneigh; - - // ii = loop over donors - // jj = loop over acceptors - // kk = loop over hydrogens bonded to donor - - int hbcount = 0; - double hbeng = 0.0; - - for (ii = iifrom; ii < iito; ++ii) { - - i = ilist[ii]; - itype = type[i]; - if (!donor[itype]) continue; - if (molecular == Atom::MOLECULAR) { - klist = special[i]; - knum = nspecial[i][0]; - } else { - if (molindex[i] < 0) continue; - imol = molindex[i]; - iatom = molatom[i]; - klist = onemols[imol]->special[iatom]; - knum = onemols[imol]->nspecial[iatom][0]; - tagprev = tag[i] - iatom - 1; - } - jlist = firstneigh[i]; - jnum = numneigh[i]; - fxtmp=fytmp=fztmp=0.0; - - xtmp = x[i].x; - ytmp = x[i].y; - ztmp = x[i].z; - - for (jj = 0; jj < jnum; jj++) { - j = jlist[jj]; - factor_hb = special_lj[sbmask(j)]; - j &= NEIGHMASK; - - jtype = type[j]; - if (!acceptor[jtype]) continue; - - delx = xtmp - x[j].x; - dely = ytmp - x[j].y; - delz = ztmp - x[j].z; - rsq = delx*delx + dely*dely + delz*delz; - - for (kk = 0; kk < knum; kk++) { - if (molecular == Atom::MOLECULAR) k = atom->map(klist[kk]); - else k = atom->map(klist[kk]+tagprev); - if (k < 0) continue; - ktype = type[k]; - m = type2param[itype][jtype][ktype]; - if (m < 0) continue; - const Param &pm = params[m]; - - if (rsq < pm.cut_outersq) { - delr1[0] = xtmp - x[k].x; - delr1[1] = ytmp - x[k].y; - delr1[2] = ztmp - x[k].z; - domain->minimum_image(delr1); - rsq1 = delr1[0]*delr1[0] + delr1[1]*delr1[1] + delr1[2]*delr1[2]; - r1 = sqrt(rsq1); - - delr2[0] = x[j].x - x[k].x; - delr2[1] = x[j].y - x[k].y; - delr2[2] = x[j].z - x[k].z; - domain->minimum_image(delr2); - rsq2 = delr2[0]*delr2[0] + delr2[1]*delr2[1] + delr2[2]*delr2[2]; - r2 = sqrt(rsq2); - - // angle (cos and sin) - - c = delr1[0]*delr2[0] + delr1[1]*delr2[1] + delr1[2]*delr2[2]; - c /= r1*r2; - if (c > 1.0) c = 1.0; - if (c < -1.0) c = -1.0; - ac = acos(c); - - ac = ac + pm.angle_offset; - c = cos(ac); - if (c > 1.0) c = 1.0; - if (c < -1.0) c = -1.0; - - if (ac > pm.cut_angle && ac < (2.0*MY_PI - pm.cut_angle)) { - s = sqrt(1.0 - c*c); - if (s < SMALL) s = SMALL; - - // Morse-specific kernel - - r = sqrt(rsq); - dr = r - pm.r0; - dexp = exp(-pm.alpha * dr); - eng_morse = pm.d0 * (dexp*dexp - 2.0*dexp); - force_kernel = pm.morse1*(dexp*dexp - dexp)/r * powint(c,pm.ap); - force_angle = pm.ap * eng_morse * powint(c,pm.ap-1)*s; - - if (rsq > pm.cut_innersq) { - switch1 = (pm.cut_outersq-rsq) * (pm.cut_outersq-rsq) * - (pm.cut_outersq + 2.0*rsq - 3.0*pm.cut_innersq) / - pm.denom_vdw; - switch2 = 12.0*rsq * (pm.cut_outersq-rsq) * - (rsq-pm.cut_innersq) / pm.denom_vdw; - force_kernel = force_kernel*switch1 + eng_morse*switch2/rsq; - force_angle *= switch1; - eng_morse *= switch1; - } - - if (EFLAG) { - evdwl = eng_morse * powint(c,pm.ap); - evdwl *= factor_hb; - } - - a = factor_hb*force_angle/s; - b = factor_hb*force_kernel; - - a11 = a*c / rsq1; - a12 = -a / (r1*r2); - a22 = a*c / rsq2; - - vx1 = a11*delr1[0] + a12*delr2[0]; - vx2 = a22*delr2[0] + a12*delr1[0]; - vy1 = a11*delr1[1] + a12*delr2[1]; - vy2 = a22*delr2[1] + a12*delr1[1]; - vz1 = a11*delr1[2] + a12*delr2[2]; - vz2 = a22*delr2[2] + a12*delr1[2]; - - fi[0] = vx1 + b*delx; - fi[1] = vy1 + b*dely; - fi[2] = vz1 + b*delz; - fj[0] = vx2 - b*delx; - fj[1] = vy2 - b*dely; - fj[2] = vz2 - b*delz; - - fxtmp += fi[0]; - fytmp += fi[1]; - fztmp += fi[2]; - - f[j].x += fj[0]; - f[j].y += fj[1]; - f[j].z += fj[2]; - - f[k].x -= vx1 + vx2; - f[k].y -= vy1 + vy2; - f[k].z -= vz1 + vz2; - - // KIJ instead of IJK b/c delr1/delr2 are both with respect to k - - if (EVFLAG) ev_tally3_thr(this,k,i,j,evdwl,0.0,fi,fj,delr1,delr2,thr); - if (EFLAG) { - hbcount++; - hbeng += evdwl; - } - } - } - } - } - f[i].x += fxtmp; - f[i].y += fytmp; - f[i].z += fztmp; - } - const int tid = thr->get_tid(); - hbcount_thr[tid] = static_cast(hbcount); - hbeng_thr[tid] = hbeng; -} - -/* ---------------------------------------------------------------------- */ - -double PairHbondDreidingMorseAngleoffsetOMP::memory_usage() -{ - double bytes = memory_usage_thr(); - bytes += (double)comm->nthreads * 2 * sizeof(double); - bytes += PairHbondDreidingMorseAngleoffset::memory_usage(); - - return bytes; + PairHbondDreidingMorseOMP(lmp) { + angle_offset_flag = 1; } diff --git a/src/OPENMP/pair_hbond_dreiding_morse_angleoffset_omp.h b/src/OPENMP/pair_hbond_dreiding_morse_angleoffset_omp.h index 5e87908df3..2808ed1e68 100644 --- a/src/OPENMP/pair_hbond_dreiding_morse_angleoffset_omp.h +++ b/src/OPENMP/pair_hbond_dreiding_morse_angleoffset_omp.h @@ -24,28 +24,16 @@ PairStyle(hbond/dreiding/morse/angleoffset/omp,PairHbondDreidingMorseAngleoffset #ifndef LMP_PAIR_HBOND_DREIDING_MORSE_ANGLEOFFSET_OMP_H #define LMP_PAIR_HBOND_DREIDING_MORSE_ANGLEOFFSET_OMP_H -#include "pair_hbond_dreiding_morse_angleoffset.h" -#include "thr_omp.h" +#include "pair_hbond_dreiding_morse_omp.h" namespace LAMMPS_NS { class PairHbondDreidingMorseAngleoffsetOMP : - public PairHbondDreidingMorseAngleoffset, - public ThrOMP { + public PairHbondDreidingMorseOMP { public: PairHbondDreidingMorseAngleoffsetOMP(class LAMMPS *); - ~PairHbondDreidingMorseAngleoffsetOMP() override; - - void compute(int, int) override; - double memory_usage() override; - - protected: - double *hbcount_thr, *hbeng_thr; - - private: - template - void eval(int ifrom, int ito, ThrData *const thr); + }; } // namespace LAMMPS_NS From 8b85ee22a31ba9466de9d29adee05415cf663501 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Tue, 28 Jan 2025 21:32:13 -0500 Subject: [PATCH 26/29] use consistent formatting across all hbond/dreiding styles --- .../pair_hbond_dreiding_morse_angleoffset.cpp | 6 +++--- src/MOLECULE/pair_hbond_dreiding_lj.cpp | 11 ++++------- src/MOLECULE/pair_hbond_dreiding_morse.cpp | 14 +++++++------- 3 files changed, 14 insertions(+), 17 deletions(-) diff --git a/src/EXTRA-MOLECULE/pair_hbond_dreiding_morse_angleoffset.cpp b/src/EXTRA-MOLECULE/pair_hbond_dreiding_morse_angleoffset.cpp index a687d4d7ea..21f26ea8d1 100644 --- a/src/EXTRA-MOLECULE/pair_hbond_dreiding_morse_angleoffset.cpp +++ b/src/EXTRA-MOLECULE/pair_hbond_dreiding_morse_angleoffset.cpp @@ -46,9 +46,9 @@ PairHbondDreidingMorseAngleoffset::PairHbondDreidingMorseAngleoffset(LAMMPS *lmp angle_offset_flag = 1; } -// /* ---------------------------------------------------------------------- -// set coeffs for one or more type pairs -// ------------------------------------------------------------------------- */ +/* ---------------------------------------------------------------------- + * set coeffs for one or more type pairs + * ---------------------------------------------------------------------- */ void PairHbondDreidingMorseAngleoffset::coeff(int narg, char **arg) { diff --git a/src/MOLECULE/pair_hbond_dreiding_lj.cpp b/src/MOLECULE/pair_hbond_dreiding_lj.cpp index 0511ecb021..fd0e61edd2 100644 --- a/src/MOLECULE/pair_hbond_dreiding_lj.cpp +++ b/src/MOLECULE/pair_hbond_dreiding_lj.cpp @@ -85,7 +85,7 @@ void PairHbondDreidingLJ::compute(int eflag, int vflag) tagint tagprev; double delx,dely,delz,rsq,rsq1,rsq2,r1,r2; double factor_hb,force_angle,force_kernel,evdwl,eng_lj,ehbond,force_switch; - double c,s,a,b,ac,a11,a12,a22,vx1,vx2,vy1,vy2,vz1,vz2,d; + double c,s,a,b,d,ac,a11,a12,a22,vx1,vx2,vy1,vy2,vz1,vz2; double fi[3],fj[3],delr1[3],delr2[3]; double r2inv,r10inv; double switch1,switch2; @@ -196,15 +196,12 @@ void PairHbondDreidingLJ::compute(int eflag, int vflag) r2inv = 1.0/rsq; r10inv = r2inv*r2inv*r2inv*r2inv*r2inv; - force_kernel = r10inv*(pm.lj1*r2inv - pm.lj2)*r2inv * - powint(c,pm.ap); - force_angle = pm.ap * r10inv*(pm.lj3*r2inv - pm.lj4) * - powint(c,pm.ap-1)*s; + force_kernel = r10inv*(pm.lj1*r2inv - pm.lj2)*r2inv * powint(c,pm.ap); + force_angle = pm.ap * r10inv*(pm.lj3*r2inv - pm.lj4) * powint(c,pm.ap-1)*s; + force_switch = 0.0; eng_lj = r10inv*(pm.lj3*r2inv - pm.lj4); - force_switch=0.0; - if (rsq > pm.cut_innersq) { switch1 = (pm.cut_outersq-rsq) * (pm.cut_outersq-rsq) * (pm.cut_outersq + 2.0*rsq - 3.0*pm.cut_innersq) / diff --git a/src/MOLECULE/pair_hbond_dreiding_morse.cpp b/src/MOLECULE/pair_hbond_dreiding_morse.cpp index 56f9c16442..8506765984 100644 --- a/src/MOLECULE/pair_hbond_dreiding_morse.cpp +++ b/src/MOLECULE/pair_hbond_dreiding_morse.cpp @@ -155,7 +155,6 @@ void PairHbondDreidingMorse::compute(int eflag, int vflag) if (c < -1.0) c = -1.0; } - if (ac > pm.cut_angle && ac < (2.0*MY_PI - pm.cut_angle)) { s = sqrt(1.0 - c*c); if (s < SMALL) s = SMALL; @@ -166,6 +165,7 @@ void PairHbondDreidingMorse::compute(int eflag, int vflag) dr = r - pm.r0; dexp = exp(-pm.alpha * dr); eng_morse = pm.d0 * (dexp*dexp - 2.0*dexp); + force_kernel = pm.morse1*(dexp*dexp - dexp)/r * powint(c,pm.ap); force_angle = pm.ap * eng_morse * powint(c,pm.ap-1)*s; force_switch = 0.0; @@ -204,12 +204,12 @@ void PairHbondDreidingMorse::compute(int eflag, int vflag) vz1 = a11*delr1[2] + a12*delr2[2]; vz2 = a22*delr2[2] + a12*delr1[2]; - fi[0] = vx1 + (b+d)*delx; - fi[1] = vy1 + (b+d)*dely; - fi[2] = vz1 + (b+d)*delz; - fj[0] = vx2 - (b+d)*delx; - fj[1] = vy2 - (b+d)*dely; - fj[2] = vz2 - (b+d)*delz; + fi[0] = vx1 + b*delx + d*delx; + fi[1] = vy1 + b*dely + d*dely; + fi[2] = vz1 + b*delz + d*delz; + fj[0] = vx2 - b*delx - d*delx; + fj[1] = vy2 - b*dely - d*dely; + fj[2] = vz2 - b*delz - d*delz; f[i][0] += fi[0]; f[i][1] += fi[1]; From 759a37cc75aa619f9bfd6f3522256a8078d50db0 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Tue, 28 Jan 2025 21:33:03 -0500 Subject: [PATCH 27/29] update and synchronize with implementation of the non-OPENMP version --- src/OPENMP/pair_hbond_dreiding_lj_omp.cpp | 44 ++++++++++---------- src/OPENMP/pair_hbond_dreiding_morse_omp.cpp | 31 ++++++++------ 2 files changed, 41 insertions(+), 34 deletions(-) diff --git a/src/OPENMP/pair_hbond_dreiding_lj_omp.cpp b/src/OPENMP/pair_hbond_dreiding_lj_omp.cpp index 9bc5f5d254..60ec8938fe 100644 --- a/src/OPENMP/pair_hbond_dreiding_lj_omp.cpp +++ b/src/OPENMP/pair_hbond_dreiding_lj_omp.cpp @@ -120,15 +120,16 @@ void PairHbondDreidingLJOMP::eval(int iifrom, int iito, ThrData * const thr) int i,j,k,m,ii,jj,kk,jnum,knum,itype,jtype,ktype,iatom,imol; tagint tagprev; double xtmp,ytmp,ztmp,delx,dely,delz,rsq,rsq1,rsq2,r1,r2; - double factor_hb,force_angle,force_kernel,evdwl,eng_lj; - double c,s,a,b,ac,a11,a12,a22,vx1,vx2,vy1,vy2,vz1,vz2; + double factor_hb,force_angle,force_kernel,evdwl,eng_lj,ehbond,force_switch; + double c,s,a,b,d,ac,a11,a12,a22,vx1,vx2,vy1,vy2,vz1,vz2; double fi[3],fj[3],delr1[3],delr2[3]; double r2inv,r10inv; double switch1,switch2; int *ilist,*jlist,*numneigh,**firstneigh; const tagint *klist; - evdwl = 0.0; + evdwl = ehbond = 0.0; + int hbcount = 0; const auto * _noalias const x = (dbl3_t *) atom->x[0]; auto * _noalias const f = (dbl3_t *) thr->get_f()[0]; @@ -151,9 +152,6 @@ void PairHbondDreidingLJOMP::eval(int iifrom, int iito, ThrData * const thr) // jj = loop over acceptors // kk = loop over hydrogens bonded to donor - int hbcount = 0; - double hbeng = 0.0; - for (ii = iifrom; ii < iito; ++ii) { i = ilist[ii]; itype = type[i]; @@ -237,30 +235,34 @@ void PairHbondDreidingLJOMP::eval(int iifrom, int iito, ThrData * const thr) r2inv = 1.0/rsq; r10inv = r2inv*r2inv*r2inv*r2inv*r2inv; - force_kernel = r10inv*(pm.lj1*r2inv - pm.lj2)*r2inv * - powint(c,pm.ap); - force_angle = pm.ap * r10inv*(pm.lj3*r2inv - pm.lj4) * - powint(c,pm.ap-1)*s; + force_kernel = r10inv*(pm.lj1*r2inv - pm.lj2)*r2inv * powint(c,pm.ap); + force_angle = pm.ap * r10inv*(pm.lj3*r2inv - pm.lj4) * powint(c,pm.ap-1)*s; + force_switch = 0.0; eng_lj = r10inv*(pm.lj3*r2inv - pm.lj4); + if (rsq > pm.cut_innersq) { switch1 = (pm.cut_outersq-rsq) * (pm.cut_outersq-rsq) * (pm.cut_outersq + 2.0*rsq - 3.0*pm.cut_innersq) / pm.denom_vdw; switch2 = 12.0*rsq * (pm.cut_outersq-rsq) * (rsq-pm.cut_innersq) / pm.denom_vdw; - force_kernel = force_kernel*switch1 + eng_lj*switch2/rsq; - force_angle *= switch1; - eng_lj *= switch1; + + force_kernel *= switch1; + force_angle *= switch1; + force_switch = eng_lj*switch2/rsq; + eng_lj *= switch1; } if (EFLAG) { evdwl = eng_lj * powint(c,pm.ap); evdwl *= factor_hb; + ehbond += evdwl; } a = factor_hb*force_angle/s; b = factor_hb*force_kernel; + d = factor_hb*force_switch; a11 = a*c / rsq1; a12 = -a / (r1*r2); @@ -273,12 +275,12 @@ void PairHbondDreidingLJOMP::eval(int iifrom, int iito, ThrData * const thr) vz1 = a11*delr1[2] + a12*delr2[2]; vz2 = a22*delr2[2] + a12*delr1[2]; - fi[0] = vx1 + b*delx; - fi[1] = vy1 + b*dely; - fi[2] = vz1 + b*delz; - fj[0] = vx2 - b*delx; - fj[1] = vy2 - b*dely; - fj[2] = vz2 - b*delz; + fi[0] = vx1 + b*delx + d*delx; + fi[1] = vy1 + b*dely + d*dely; + fi[2] = vz1 + b*delz + d*delz; + fj[0] = vx2 - b*delx - d*delx; + fj[1] = vy2 - b*dely - d*dely; + fj[2] = vz2 - b*delz - d*delz; fxtmp += fi[0]; fytmp += fi[1]; @@ -297,7 +299,7 @@ void PairHbondDreidingLJOMP::eval(int iifrom, int iito, ThrData * const thr) if (EVFLAG) ev_tally3_thr(this,k,i,j,evdwl,0.0,fi,fj,delr1,delr2,thr); if (EFLAG) { hbcount++; - hbeng += evdwl; + ehbond += evdwl; } } } @@ -309,7 +311,7 @@ void PairHbondDreidingLJOMP::eval(int iifrom, int iito, ThrData * const thr) } const int tid = thr->get_tid(); hbcount_thr[tid] = static_cast(hbcount); - hbeng_thr[tid] = hbeng; + hbeng_thr[tid] = ehbond; } /* ---------------------------------------------------------------------- */ diff --git a/src/OPENMP/pair_hbond_dreiding_morse_omp.cpp b/src/OPENMP/pair_hbond_dreiding_morse_omp.cpp index bb0970dbdc..7772ea69fa 100644 --- a/src/OPENMP/pair_hbond_dreiding_morse_omp.cpp +++ b/src/OPENMP/pair_hbond_dreiding_morse_omp.cpp @@ -120,14 +120,14 @@ void PairHbondDreidingMorseOMP::eval(int iifrom, int iito, ThrData * const thr) int i,j,k,m,ii,jj,kk,jnum,knum,itype,jtype,ktype,imol,iatom; tagint tagprev; double xtmp,ytmp,ztmp,delx,dely,delz,rsq,rsq1,rsq2,r1,r2; - double factor_hb,force_angle,force_kernel,evdwl; - double c,s,a,b,ac,a11,a12,a22,vx1,vx2,vy1,vy2,vz1,vz2; + double factor_hb,force_angle,force_kernel,force_switch,evdwl,ehbond; + double c,s,a,b,d,ac,a11,a12,a22,vx1,vx2,vy1,vy2,vz1,vz2; double fi[3],fj[3],delr1[3],delr2[3]; double r,dr,dexp,eng_morse,switch1,switch2; int *ilist,*jlist,*numneigh,**firstneigh; const tagint *klist; - evdwl = 0.0; + evdwl = ehbond = 0.0; const auto * _noalias const x = (dbl3_t *) atom->x[0]; auto * _noalias const f = (dbl3_t *) thr->get_f()[0]; @@ -151,7 +151,6 @@ void PairHbondDreidingMorseOMP::eval(int iifrom, int iito, ThrData * const thr) // kk = loop over hydrogens bonded to donor int hbcount = 0; - double hbeng = 0.0; for (ii = iifrom; ii < iito; ++ii) { @@ -239,8 +238,10 @@ void PairHbondDreidingMorseOMP::eval(int iifrom, int iito, ThrData * const thr) dr = r - pm.r0; dexp = exp(-pm.alpha * dr); eng_morse = pm.d0 * (dexp*dexp - 2.0*dexp); + force_kernel = pm.morse1*(dexp*dexp - dexp)/r * powint(c,pm.ap); force_angle = pm.ap * eng_morse * powint(c,pm.ap-1)*s; + force_switch = 0.0; if (rsq > pm.cut_innersq) { switch1 = (pm.cut_outersq-rsq) * (pm.cut_outersq-rsq) * @@ -248,18 +249,22 @@ void PairHbondDreidingMorseOMP::eval(int iifrom, int iito, ThrData * const thr) pm.denom_vdw; switch2 = 12.0*rsq * (pm.cut_outersq-rsq) * (rsq-pm.cut_innersq) / pm.denom_vdw; - force_kernel = force_kernel*switch1 + eng_morse*switch2/rsq; + + force_kernel *= switch1; force_angle *= switch1; + force_switch = eng_morse*switch2/rsq; eng_morse *= switch1; } if (EFLAG) { evdwl = eng_morse * powint(c,pm.ap); evdwl *= factor_hb; + ehbond += evdwl; } a = factor_hb*force_angle/s; b = factor_hb*force_kernel; + d = factor_hb*force_switch; a11 = a*c / rsq1; a12 = -a / (r1*r2); @@ -272,12 +277,12 @@ void PairHbondDreidingMorseOMP::eval(int iifrom, int iito, ThrData * const thr) vz1 = a11*delr1[2] + a12*delr2[2]; vz2 = a22*delr2[2] + a12*delr1[2]; - fi[0] = vx1 + b*delx; - fi[1] = vy1 + b*dely; - fi[2] = vz1 + b*delz; - fj[0] = vx2 - b*delx; - fj[1] = vy2 - b*dely; - fj[2] = vz2 - b*delz; + fi[0] = vx1 + b*delx + d*delx; + fi[1] = vy1 + b*dely + d*dely; + fi[2] = vz1 + b*delz + d*delz; + fj[0] = vx2 - b*delx - d*delx; + fj[1] = vy2 - b*dely - d*dely; + fj[2] = vz2 - b*delz - d*delz; fxtmp += fi[0]; fytmp += fi[1]; @@ -296,7 +301,7 @@ void PairHbondDreidingMorseOMP::eval(int iifrom, int iito, ThrData * const thr) if (EVFLAG) ev_tally3_thr(this,k,i,j,evdwl,0.0,fi,fj,delr1,delr2,thr); if (EFLAG) { hbcount++; - hbeng += evdwl; + ehbond += evdwl; } } } @@ -308,7 +313,7 @@ void PairHbondDreidingMorseOMP::eval(int iifrom, int iito, ThrData * const thr) } const int tid = thr->get_tid(); hbcount_thr[tid] = static_cast(hbcount); - hbeng_thr[tid] = hbeng; + hbeng_thr[tid] = ehbond; } /* ---------------------------------------------------------------------- */ From 201d1a59b5cc7550c9ff3b768473a78f66c1fb1f Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Tue, 28 Jan 2025 21:33:44 -0500 Subject: [PATCH 28/29] the /angleoffset versions have their own different parameter file and reader --- ...pair_hbond_dreiding_lj_angleoffset_omp.cpp | 85 +++++++++++++++++- .../pair_hbond_dreiding_lj_angleoffset_omp.h | 1 + ...r_hbond_dreiding_morse_angleoffset_omp.cpp | 86 ++++++++++++++++++- ...air_hbond_dreiding_morse_angleoffset_omp.h | 2 +- 4 files changed, 171 insertions(+), 3 deletions(-) diff --git a/src/OPENMP/pair_hbond_dreiding_lj_angleoffset_omp.cpp b/src/OPENMP/pair_hbond_dreiding_lj_angleoffset_omp.cpp index 1a5126971a..11c09ed549 100644 --- a/src/OPENMP/pair_hbond_dreiding_lj_angleoffset_omp.cpp +++ b/src/OPENMP/pair_hbond_dreiding_lj_angleoffset_omp.cpp @@ -22,6 +22,7 @@ #include "force.h" #include "math_const.h" #include "math_special.h" +#include "memory.h" #include "molecule.h" #include "neigh_list.h" #include "suffix.h" @@ -33,7 +34,7 @@ using namespace LAMMPS_NS; using namespace MathConst; using namespace MathSpecial; -static constexpr double SMALL = 0.001; +static constexpr int CHUNK = 8; /* ---------------------------------------------------------------------- */ @@ -42,3 +43,85 @@ PairHbondDreidingLJAngleoffsetOMP::PairHbondDreidingLJAngleoffsetOMP(LAMMPS *lmp angle_offset_flag = 1; } +/* ---------------------------------------------------------------------- + set coeffs for one or more type pairs +------------------------------------------------------------------------- */ + +void PairHbondDreidingLJAngleoffsetOMP::coeff(int narg, char **arg) +{ + auto mylmp = PairHbondDreidingLJ::lmp; + if (narg < 6 || narg > 11) + error->all(FLERR,"Incorrect args for pair coefficients"); + if (!allocated) allocate(); + + int ilo,ihi,jlo,jhi,klo,khi; + utils::bounds(FLERR, arg[0], 1, atom->ntypes, ilo, ihi, error); + utils::bounds(FLERR, arg[1], 1, atom->ntypes, jlo, jhi, error); + utils::bounds_typelabel(FLERR, arg[2], 1, atom->ntypes, klo, khi, mylmp, Atom::ATOM); + + int donor_flag; + if (strcmp(arg[3],"i") == 0) donor_flag = 0; + else if (strcmp(arg[3],"j") == 0) donor_flag = 1; + else error->all(FLERR,"Incorrect args for pair coefficients"); + + double epsilon_one = utils::numeric(FLERR, arg[4], false, mylmp); + double sigma_one = utils::numeric(FLERR, arg[5], false, mylmp); + + int ap_one = ap_global; + if (narg > 6) ap_one = utils::inumeric(FLERR, arg[6], false, mylmp); + double cut_inner_one = cut_inner_global; + double cut_outer_one = cut_outer_global; + if (narg > 8) { + cut_inner_one = utils::numeric(FLERR, arg[7], false, mylmp); + cut_outer_one = utils::numeric(FLERR, arg[8], false, mylmp); + } + if (cut_inner_one>cut_outer_one) + error->all(FLERR,"Pair inner cutoff >= Pair outer cutoff"); + double cut_angle_one = cut_angle_global; + if (narg > 9) cut_angle_one = utils::numeric(FLERR, arg[9], false, mylmp) * MY_PI/180.0; + double angle_offset_one = angle_offset_global; + if (narg == 11) angle_offset_one = (180.0 - utils::numeric(FLERR, arg[10], false, mylmp)) * MY_PI/180.0; + if (angle_offset_one < 0.0 || angle_offset_one > 90.0 * MY_PI/180.0) + error->all(FLERR,"Illegal angle offset"); + + // grow params array if necessary + + if (nparams == maxparam) { + maxparam += CHUNK; + params = (Param *) memory->srealloc(params, maxparam*sizeof(Param), + "pair:params"); + + // make certain all addional allocated storage is initialized + // to avoid false positives when checking with valgrind + + memset(params + nparams, 0, CHUNK*sizeof(Param)); + } + + params[nparams].epsilon = epsilon_one; + params[nparams].sigma = sigma_one; + params[nparams].ap = ap_one; + params[nparams].cut_inner = cut_inner_one; + params[nparams].cut_outer = cut_outer_one; + params[nparams].cut_innersq = cut_inner_one*cut_inner_one; + params[nparams].cut_outersq = cut_outer_one*cut_outer_one; + params[nparams].cut_angle = cut_angle_one; + params[nparams].angle_offset = angle_offset_one; + params[nparams].denom_vdw = + (params[nparams].cut_outersq-params[nparams].cut_innersq) * + (params[nparams].cut_outersq-params[nparams].cut_innersq) * + (params[nparams].cut_outersq-params[nparams].cut_innersq); + + // flag type2param with either i,j = D,A or j,i = D,A + + int count = 0; + for (int i = ilo; i <= ihi; i++) + for (int j = MAX(jlo,i); j <= jhi; j++) + for (int k = klo; k <= khi; k++) { + if (donor_flag == 0) type2param[i][j][k] = nparams; + else type2param[j][i][k] = nparams; + count++; + } + nparams++; + + if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients"); +} diff --git a/src/OPENMP/pair_hbond_dreiding_lj_angleoffset_omp.h b/src/OPENMP/pair_hbond_dreiding_lj_angleoffset_omp.h index ce79a03dee..03d3392e4d 100644 --- a/src/OPENMP/pair_hbond_dreiding_lj_angleoffset_omp.h +++ b/src/OPENMP/pair_hbond_dreiding_lj_angleoffset_omp.h @@ -32,6 +32,7 @@ class PairHbondDreidingLJAngleoffsetOMP : public PairHbondDreidingLJOMP { public: PairHbondDreidingLJAngleoffsetOMP(class LAMMPS *); + void coeff(int, char **) override; }; } // namespace LAMMPS_NS diff --git a/src/OPENMP/pair_hbond_dreiding_morse_angleoffset_omp.cpp b/src/OPENMP/pair_hbond_dreiding_morse_angleoffset_omp.cpp index b0f9cde92f..e7c75f29e4 100644 --- a/src/OPENMP/pair_hbond_dreiding_morse_angleoffset_omp.cpp +++ b/src/OPENMP/pair_hbond_dreiding_morse_angleoffset_omp.cpp @@ -22,6 +22,7 @@ #include "force.h" #include "math_const.h" #include "math_special.h" +#include "memory.h" #include "molecule.h" #include "neigh_list.h" #include "suffix.h" @@ -33,7 +34,7 @@ using namespace LAMMPS_NS; using namespace MathConst; using namespace MathSpecial; -static constexpr double SMALL = 0.001; +static constexpr int CHUNK = 8; /* ---------------------------------------------------------------------- */ @@ -41,3 +42,86 @@ PairHbondDreidingMorseAngleoffsetOMP::PairHbondDreidingMorseAngleoffsetOMP(LAMMP PairHbondDreidingMorseOMP(lmp) { angle_offset_flag = 1; } + +/* ---------------------------------------------------------------------- + * set coeffs for one or more type pairs + * ---------------------------------------------------------------------- */ + +void PairHbondDreidingMorseAngleoffsetOMP::coeff(int narg, char **arg) +{ + auto mylmp = PairHbondDreidingMorse::lmp; + if (narg < 7 || narg > 12) + error->all(FLERR,"Incorrect args for pair coefficients"); + if (!allocated) allocate(); + + int ilo,ihi,jlo,jhi,klo,khi; + utils::bounds(FLERR, arg[0], 1, atom->ntypes, ilo, ihi, error); + utils::bounds(FLERR, arg[1], 1, atom->ntypes, jlo, jhi, error); + utils::bounds_typelabel(FLERR, arg[2], 1, atom->ntypes, klo, khi, mylmp, Atom::ATOM); + + int donor_flag; + if (strcmp(arg[3],"i") == 0) donor_flag = 0; + else if (strcmp(arg[3],"j") == 0) donor_flag = 1; + else error->all(FLERR,"Incorrect args for pair coefficients"); + + double d0_one = utils::numeric(FLERR, arg[4], false, mylmp); + double alpha_one = utils::numeric(FLERR, arg[5], false, mylmp); + double r0_one = utils::numeric(FLERR, arg[6], false, mylmp); + + int ap_one = ap_global; + if (narg > 7) ap_one = utils::inumeric(FLERR, arg[7], false, mylmp); + double cut_inner_one = cut_inner_global; + double cut_outer_one = cut_outer_global; + if (narg > 9) { + cut_inner_one = utils::numeric(FLERR, arg[8], false, mylmp); + cut_outer_one = utils::numeric(FLERR, arg[9], false, mylmp); + } + if (cut_inner_one>cut_outer_one) + error->all(FLERR,"Pair inner cutoff >= Pair outer cutoff"); + double cut_angle_one = cut_angle_global; + if (narg > 10) cut_angle_one = utils::numeric(FLERR, arg[10], false, mylmp) * MY_PI/180.0; + double angle_offset_one = angle_offset_global; + if (narg == 12) angle_offset_one = (180.0 - utils::numeric(FLERR, arg[11], false, mylmp)) * MY_PI/180.0; + if (angle_offset_one < 0.0 || angle_offset_one > 90.0 * MY_PI/180.0) + error->all(FLERR,"Illegal angle offset {}", angle_offset_one); + + // grow params array if necessary + + if (nparams == maxparam) { + maxparam += CHUNK; + params = (Param *) memory->srealloc(params, maxparam*sizeof(Param),"pair:params"); + + // make certain all addional allocated storage is initialized + // to avoid false positives when checking with valgrind + + memset(params + nparams, 0, CHUNK*sizeof(Param)); + } + + params[nparams].d0 = d0_one; + params[nparams].alpha = alpha_one; + params[nparams].r0 = r0_one; + params[nparams].ap = ap_one; + params[nparams].cut_inner = cut_inner_one; + params[nparams].cut_outer = cut_outer_one; + params[nparams].cut_innersq = cut_inner_one*cut_inner_one; + params[nparams].cut_outersq = cut_outer_one*cut_outer_one; + params[nparams].cut_angle = cut_angle_one; + params[nparams].angle_offset = angle_offset_one; + params[nparams].denom_vdw = (params[nparams].cut_outersq-params[nparams].cut_innersq) * + (params[nparams].cut_outersq-params[nparams].cut_innersq) * + (params[nparams].cut_outersq-params[nparams].cut_innersq); + + // flag type2param with either i,j = D,A or j,i = D,A + + int count = 0; + for (int i = ilo; i <= ihi; i++) + for (int j = MAX(jlo,i); j <= jhi; j++) + for (int k = klo; k <= khi; k++) { + if (donor_flag == 0) type2param[i][j][k] = nparams; + else type2param[j][i][k] = nparams; + count++; + } + nparams++; + + if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients"); +} diff --git a/src/OPENMP/pair_hbond_dreiding_morse_angleoffset_omp.h b/src/OPENMP/pair_hbond_dreiding_morse_angleoffset_omp.h index 2808ed1e68..94b6131f98 100644 --- a/src/OPENMP/pair_hbond_dreiding_morse_angleoffset_omp.h +++ b/src/OPENMP/pair_hbond_dreiding_morse_angleoffset_omp.h @@ -33,7 +33,7 @@ class PairHbondDreidingMorseAngleoffsetOMP : public: PairHbondDreidingMorseAngleoffsetOMP(class LAMMPS *); - + void coeff(int, char **) override; }; } // namespace LAMMPS_NS From 229916e11fdb855c8940fd79d5943fad249dc3d8 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Tue, 28 Jan 2025 21:35:11 -0500 Subject: [PATCH 29/29] whitespace --- src/EXTRA-MOLECULE/pair_hbond_dreiding_morse_angleoffset.h | 2 +- src/OPENMP/pair_hbond_dreiding_morse_angleoffset_omp.h | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/src/EXTRA-MOLECULE/pair_hbond_dreiding_morse_angleoffset.h b/src/EXTRA-MOLECULE/pair_hbond_dreiding_morse_angleoffset.h index 5eeb1796a2..20d2a8698d 100644 --- a/src/EXTRA-MOLECULE/pair_hbond_dreiding_morse_angleoffset.h +++ b/src/EXTRA-MOLECULE/pair_hbond_dreiding_morse_angleoffset.h @@ -29,7 +29,7 @@ class PairHbondDreidingMorseAngleoffset : public PairHbondDreidingMorse { public: PairHbondDreidingMorseAngleoffset(class LAMMPS *); void coeff(int, char **) override; - + }; } // namespace LAMMPS_NS diff --git a/src/OPENMP/pair_hbond_dreiding_morse_angleoffset_omp.h b/src/OPENMP/pair_hbond_dreiding_morse_angleoffset_omp.h index 94b6131f98..40a797ac33 100644 --- a/src/OPENMP/pair_hbond_dreiding_morse_angleoffset_omp.h +++ b/src/OPENMP/pair_hbond_dreiding_morse_angleoffset_omp.h @@ -33,7 +33,7 @@ class PairHbondDreidingMorseAngleoffsetOMP : public: PairHbondDreidingMorseAngleoffsetOMP(class LAMMPS *); - void coeff(int, char **) override; + void coeff(int, char **) override; }; } // namespace LAMMPS_NS