From 201d1a59b5cc7550c9ff3b768473a78f66c1fb1f Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Tue, 28 Jan 2025 21:33:44 -0500 Subject: [PATCH] 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