diff --git a/src/pair.h b/src/pair.h index 5d70d8efbb..c00f442679 100644 --- a/src/pair.h +++ b/src/pair.h @@ -170,10 +170,9 @@ class Pair : protected Pointers { } virtual void born_matrix(int /*i*/, int /*j*/, int /*itype*/, int /*jtype*/, double /*rsq*/, - double /*factor_coul*/, double /*factor_lj*/, double& du, double& du2) + double /*factor_coul*/, double /*factor_lj*/, double &du, double &du2) { - du = 0.0; - du2 = 0.0; + du = du2 = 0.0; } virtual void settings(int, char **) = 0; diff --git a/src/pair_hybrid.cpp b/src/pair_hybrid.cpp index 836624cd4b..64a59db4ad 100644 --- a/src/pair_hybrid.cpp +++ b/src/pair_hybrid.cpp @@ -874,6 +874,40 @@ double PairHybrid::single(int i, int j, int itype, int jtype, return esum; } +/* ---------------------------------------------------------------------- + call sub-style to compute born matrix interaction + error if sub-style does not support born_matrix call + since overlay could have multiple sub-styles, sum results explicitly +------------------------------------------------------------------------- */ + +void PairHybrid::born_matrix(int i, int j, int itype, int jtype, double rsq, + double factor_coul, double factor_lj, + double &dupair, double &du2pair) +{ + if (nmap[itype][jtype] == 0) + error->one(FLERR,"Invoked pair born_matrix on pair style none"); + + double du, du2; + dupair = du2pair = 0.0; + + for (int m = 0; m < nmap[itype][jtype]; m++) { + if (rsq < styles[map[itype][jtype][m]]->cutsq[itype][jtype]) { + if (styles[map[itype][jtype][m]]->born_matrix_enable == 0) + error->one(FLERR,"Pair hybrid sub-style does not support born_matrix call"); + + if ((special_lj[map[itype][jtype][m]] != nullptr) || + (special_coul[map[itype][jtype][m]] != nullptr)) + error->one(FLERR,"Pair hybrid born_matrix calls do not support" + " per sub-style special bond values"); + + du = du2 = 0.0; + styles[map[itype][jtype][m]]->born_matrix(i,j,itype,jtype,rsq,factor_coul,factor_lj,du,du2); + dupair += du; + du2pair += du2; + } + } +} + /* ---------------------------------------------------------------------- copy Pair::svector data ------------------------------------------------------------------------- */ diff --git a/src/pair_hybrid.h b/src/pair_hybrid.h index ee56783700..675bd16043 100644 --- a/src/pair_hybrid.h +++ b/src/pair_hybrid.h @@ -49,6 +49,8 @@ class PairHybrid : public Pair { void write_restart(FILE *) override; void read_restart(FILE *) override; double single(int, int, int, int, double, double, double, double &) override; + void born_matrix(int, int, int, int, double, double, double, double &, double &) override; + void modify_params(int narg, char **arg) override; double memory_usage() override; diff --git a/src/pair_hybrid_scaled.cpp b/src/pair_hybrid_scaled.cpp index 24158f46a0..a2a7184ddd 100644 --- a/src/pair_hybrid_scaled.cpp +++ b/src/pair_hybrid_scaled.cpp @@ -413,7 +413,7 @@ double PairHybridScaled::single(int i, int j, int itype, int jtype, double rsq, (special_coul[map[itype][jtype][m]] != nullptr)) error->one(FLERR, "Pair hybrid single() does not support per sub-style special_bond"); - scale = scaleval[map[itype][jtype][m]]; + double scale = scaleval[map[itype][jtype][m]]; esum += scale * pstyle->single(i, j, itype, jtype, rsq, factor_coul, factor_lj, fone); fforce += scale * fone; } @@ -423,6 +423,57 @@ double PairHybridScaled::single(int i, int j, int itype, int jtype, double rsq, return esum; } +/* ---------------------------------------------------------------------- + call sub-style to compute born matrix interaction + error if sub-style does not support born_matrix call + since overlay could have multiple sub-styles, sum results explicitly +------------------------------------------------------------------------- */ + +void PairHybridScaled::born_matrix(int i, int j, int itype, int jtype, double rsq, + double factor_coul, double factor_lj, + double &dupair, double &du2pair) +{ + if (nmap[itype][jtype] == 0) error->one(FLERR, "Invoked pair born_matrix on pair style none"); + + // update scale values from variables where needed + + const int nvars = scalevars.size(); + if (nvars > 0) { + double *vals = new double[nvars]; + for (int k = 0; k < nvars; ++k) { + int m = input->variable->find(scalevars[k].c_str()); + if (m < 0) + error->all(FLERR, "Variable '{}' not found when updating scale factors", scalevars[k]); + vals[k] = input->variable->compute_equal(m); + } + for (int k = 0; k < nstyles; ++k) { + if (scaleidx[k] >= 0) scaleval[k] = vals[scaleidx[k]]; + } + delete[] vals; + } + + double du, du2, scale; + dupair = du2pair = scale = 0.0; + + for (int m = 0; m < nmap[itype][jtype]; m++) { + auto pstyle = styles[map[itype][jtype][m]]; + if (rsq < pstyle->cutsq[itype][jtype]) { + if (pstyle->born_matrix_enable == 0) + error->one(FLERR, "Pair hybrid sub-style does not support born_matrix call"); + + if ((special_lj[map[itype][jtype][m]] != nullptr) || + (special_coul[map[itype][jtype][m]] != nullptr)) + error->one(FLERR, "Pair hybrid born_matrix() does not support per sub-style special_bond"); + + du = du2 = 0.0; + scale = scaleval[map[itype][jtype][m]]; + pstyle->born_matrix(i, j, itype, jtype, rsq, factor_coul, factor_lj, du, du2); + dupair += scale*du; + du2pair += scale*du2; + } + } +} + /* ---------------------------------------------------------------------- set coeffs for one or more type pairs ------------------------------------------------------------------------- */ diff --git a/src/pair_hybrid_scaled.h b/src/pair_hybrid_scaled.h index 9bb8901846..a7789c7164 100644 --- a/src/pair_hybrid_scaled.h +++ b/src/pair_hybrid_scaled.h @@ -38,6 +38,7 @@ class PairHybridScaled : public PairHybrid { void write_restart(FILE *) override; void read_restart(FILE *) override; double single(int, int, int, int, double, double, double, double &) override; + void born_matrix(int, int, int, int, double, double, double, double &, double &) override; void init_svector() override; void copy_svector(int, int) override;