From d2b5f557379fdd8f7abb4ed6886aeb4028968084 Mon Sep 17 00:00:00 2001 From: EiPi Fun Date: Sun, 1 Sep 2024 16:07:52 +0800 Subject: [PATCH] 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